the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 23 Mar 2026
Regularisation of the 4DEnVar Data Assimilation method for Calibration of Land Surface Models
Abstract. The four-dimensional ensemble variational (4DEnVar) data assimilation method is an attractive choice for land surface model calibration due to its ease of use, speed, and its circumvention of tangent linear model and adjoint calculations. However, in certain circumstances, this method may prove ineffective when implemented with a highly non-linear model, leading to out-of-range or non-physical posterior parameter values. To address this, the 4DEnVar cost function can be adapted through the introduction of a hyperparameter, which inflates the weight of the background term. In this study, we explore the anticipated challenges of applying 4DEnVar with in situ eddy-covariance flux measurements and present some explanations of the expected behaviour. We show that, when using a hyperparameter to regularise the optimisation, 4DEnVar is able to successfully calibrate two complex land surface models, JULES and ORCHIDEE, with comparable accuracy in its ability to produce model runs with an improved match to the observations. To our knowledge, this is the first study of its kind to compare parameter calibration of two different land surface models in the same experimental setting. In addition to the aforementioned benefits of using 4DEnVar, we show that the method also exhibits considerable versatility, not only with regard to the land surface model to which it is applied, but also in terms of parameter set selection and ensemble size.
- Preprint
(818 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-132', Anonymous Referee #1, 29 Apr 2026
-
AC1: 'Reply on RC1', Natalie Douglas, 14 May 2026
We thank the referee for the comments given. Our replies to each comment are below:
- L. 74-76: This sounds like a result or conclusion. In the introduction, it should be rephrased as a research question to be addressed.
We have changed:
‘As we demonstrate in this study, when not set up with sufficient care, the method is capable of retrieving unlikely (out of an expected range) or even impossible (e.g. implausibly negative) parameter values, which can result in failed model runs using posterior parameters. This is due to a control variable transform of the original 4DVar cost function that it implements, affording the method limited restriction in its search for optimal and physically valid parameters (see section 2.2.3). In this article, we examine these behaviours and describe the key changes necessary to enable the assimilation of data that was not previously possible, giving special attention to the nuances of 4DEnVar’.To:
‘When not set up with sufficient care, the method is capable of retrieving unlikely (out of an expected range) or even impossible (e.g. implausibly negative) parameter values, which can result in failed model runs using posterior parameters. This behaviour arises from the control variable transform of the original 4DVar cost function that it implements, which poses limited constraints on the search for physically valid solutions (see section 2.2.3). This raises the question of under what conditions such behaviour occurs, and how it can be mitigated when applying 4DEnVar to complex land surface models. In this study, we therefore investigate these behaviours and explore the key modifications required to enable robust assimilation of real observations, with particular attention to the role of regularisation within the 4DEnVar framework- L. 76-78: Results should be interpreted in the Discussion section, not the Introduction. This sentence could be replaced with a research question that needs answering.
See above.
- L. 126 (Table 1): Replace "N/A" by "-" or by "Unitless". I am surprised that a key parameter such as rooting depth has not been included in the analysis. In what way does each model describe soil moisture stress?
We have replaced ‘N/A’ in Table 1 with ‘-‘. In JULES, rooting depth is not described as a single parameter, rather, root water uptake is governed by a prescribed root distribution profile across soil layers. Soil moisture stress is calculated as a function of soil moisture availability weighted by this root distribution. We chose to keep the parameter set small for JULES, focusing on photosynthetic parameter calibration in JULES, rather than hydrological parameter estimation, and so we did not include these aspects in the calibration. However, for ORCHIDEE, the hum_cst parameter speaks directly to root distribution and is included in the calibration. In ORCHIDEE, soil moisture stress is represented through its effects on stomatal conductance and photosynthesis, mediated by soil hydraulic properties and root water uptake.
We have changed ‘Soil hydraulic and thermal properties are determined from soil moisture and soil texture, with textures classified into 12 categories based on the USDA soil map…’ to ‘Soil hydraulic and thermal properties are determined from soil moisture, soil texture, and the vertical distribution of roots, with textures classified into 12 categories based on the USDA soil map…’ to aid clarity.
- L. 127 (Section 2.2): Does this method account for errors in the atmospheric variables that the models use as input?
Errors in the atmospheric variables are not explicitly represented in the 4DEnVar method, but they are implicitly included in the model-observation error term represented by the R-matrix. We have changed ‘...and a fixed diagonal ${\bf R}$-matrix with standard deviation of 0.35 (gC m$^{-2}$ d$^{-1}$) for both LSMs representing the joint model-observation error.’ to ‘...and a fixed diagonal ${\bf R}$-matrix with standard deviation of 0.35 (gC m$^{-2}$ d$^{-1}$) for both LSMs representing the joint model-observation error (including uncertainties in the atmospheric driving variables).’- L. 130: "observations given the target" is unclear. Do you mean "given the value of the target"?
Yes, we intended to mean “given the value of the target variables”. We have changed ‘given the target’ to ‘given the value of the target variables’.
- L. 260: The reason why this Fluxnet site has been chosen over another should be explained. Are soil moisture deficits observed at this site?
We have changed ‘. Both models are run offline at the Harvard Forest site (FLUXNET2015 Dataset, 2015) using the meteorological forcing data provided by the FluxNet network for this site (Pastorello et al., 2020). For both models, we assimilate GPP observations derived from Net Ecosystem Exchange eddy covariance tower measurements.’ to ‘. Both models are run offline at the Harvard Forest site (FLUXNET2015 Dataset, 2015) using the meteorological forcing data provided by the FluxNet network for this site (Pastorello et al., 2020). The Harvard Forest site was selected because it provides a long, continuous, and well-characterised eddy-covariance and meteorological record and experiences seasonal soil moisture limitations during summer periods, making it suitable for both photosynthetic and hydrological parameter calibration. For both models, we assimilate GPP estimates derived from eddy-covariance Net Ecosystem Exchange (NEE) measurements using standard FLUXNET partitioning methods that separate ecosystem respiration and gross primary productivity.’ in order to address this comment and the next one.
- L. 262: 'GPP observations' are used. It should be made clear that GPP is not directly observed. I suppose the GPP values you use are derived from NEE observations and an ecosystem respiration model?
See previous comment.
- L. 327 (Fig. 5): RMSD values should have units. Year 2010 is interesting because both calibrated models are wrong. Why is there such a large delay in peak GPP for the two models?
We have added units to the RMSD values in Figure 5 which shows the validation period. Both models were calibrated against years where peak GPP occurred in July, whereas the 2010 observations exhibit an earlier peak. This discrepancy likely reflects interannual variability in phenology or environmental conditions that is not fully represented in the calibrated models or forcing data. We have changed ‘Notably, the 2010 observations exhibit a monthly maximum GPP that is offset by one month compared to the simulations from both models. In this case, the models have been calibrated with observations where GPP peaks are happening in July. In 2010, however, whether down to observational error, changes in meteorological forcing data or a genuine early peak in photosynthesis, the models are unable to reach this peak seen in observations when run with the optimised parameter values’ to ‘Notably, the 2010 observations exhibit a monthly maximum GPP that is offset by one month compared to the posterior simulations from both models. In this case, the models have been calibrated against observations where peak GPP occurred in July and the optimised parameter sets therefore favour this seasonal timing. The mismatch in 2010 likely reflects interannual variability in environmental or phenological conditions that is not fully captured by the calibrated models or forcing data, leading both models to delay the onset and peak of photosynthetic activity for this year.
- L. 386 (Fig. A1): Figure A1 is not mentioned in the text. The current caption of Figure A1 includes an interpretation. This interpretation should be incorporated into the main text. The reasons why the L-curve technique cannot be used to select gamma should be clearly explained.
We have now referred to Figure A1 in the text and have removed the interpretation from the Figure A1 caption. We have provided a more coherent explanation in the main text by changing ‘However, we note that h (xb) + Y′bw∗ serves as a linear approximation to h (x∗) and so these norms are calculated without the need for posterior
model runs. Consequently, any γ chosen this way does not translate well to the calculation of genuine posterior residual norms. As an example to illustrate this point, one can simply consider the case where γ = 1 that results in implausible parameter values - the residual norm can be calculated in w-space but it is not possible to calculate the true residual norm as the model can not be run on the corresponding parameters in x-space. A deterioration of the L-curve, when plotting the prior norm against the
residual norm using the linear approximation versus using the posterior model run, can be seen in Appendix A.’ to ‘However, h (xb) + Y′bw∗ only serves as a linear approximation to h (x∗) meaning that these residual norms are evaluated in the transformed control-variable space rather than from genuine posterior model simulations. Consequently, standard regularisation selection techniques such as the L-curve criterion can become misleading in the context of 4DEnVar. In particular, the approximate residual norms may suggest improved agreement with observations even when the corresponding posterior parameter values are physically implausible or produce degraded posterior model runs. Figure A1 demonstrates this behaviour for the JULES experiment: the L-curve constructed using the linear approximation differs substantially from that obtained using posterior model runs.
Citation: https://doi.org/10.5194/egusphere-2026-132-AC1 -
AC2: 'Edit to Reply on AC1', Natalie Douglas, 14 May 2026
'rooting depth is not described as a single parameter, rather, root water uptake is governed..' should say 'rooting depth is described as a single parameter but root water uptake is governed...'
Citation: https://doi.org/10.5194/egusphere-2026-132-AC2
-
AC2: 'Edit to Reply on AC1', Natalie Douglas, 14 May 2026
-
AC1: 'Reply on RC1', Natalie Douglas, 14 May 2026
-
RC2: 'Comment on egusphere-2026-132', Anonymous Referee #2, 10 Sep 2026
Thank you for the opportunity to review this paper. The authors address an important problem in applying 4DEnVar to parameter calibration in complex land surface models, and demonstrate that the proposed regularisation can prevent implausibly large and/or unphysical parameter updates and improve the subsequent model simulations in both JULES and ORCHIDEE.
I have two main questions about the interpretation of the gamma hyperparameter. The mathematics here is close to the edge of my own comfort zone, so I would genuinely welcome clarification if I have misunderstood things.
1. Relationship between gamma and the observation error covariance R
As I understand Eq. 9, increasing gamma increases the relative weight given to staying close to the prior. However, if the entire cost function is divided by gamma, the location of the minimum is unchanged, and the resulting expression appears to be exactly equivalent to standard 4DEnVar with gamma equal to 1 but with R multiplied by gamma.
If I have followed this correctly, the final analysis parameter vector obtained using gamma and R should therefore be identical to that obtained by leaving gamma at 1 and inflating R by the same factor. The authors note this relationship in the Discussion, but I wonder whether its implications deserve somewhat more attention.
For example, with the prescribed error standard deviation of 0.35 gC m-2 d-1, gamma values of 60 and 100 correspond, for determining the analysis mean, to effective standard deviations of about 2.7 and 3.5 gC m-2 d-1, respectively.
These are very large changes in the relative weighting of the prior and observational terms. I therefore wonder whether the experiments as presented can distinguish between gamma acting primarily to restrict the update to a region where the local linear approximation remains valid, and gamma compensating for the weighting implied by the prescribed R. These seem to have somewhat different implications for how a future user should choose R and gamma.
Could the authors discuss this distinction more explicitly? The simple nonlinear example earlier in the manuscript might perhaps provide a useful controlled setting in which to demonstrate specifically how regularisation addresses the nonlinear or linearisation problem when R is known.
2. Interpretation of the posterior uncertainty
My second question concerns Eq. 11 and the narrow posterior distributions shown in Fig. 3.
As a simple check, consider the limiting case in which the observations provide no information about the parameters. In that case, Eq. 11 appears to imply that the posterior perturbations are reduced by a factor of the square root of gamma, and therefore that the posterior covariance is reduced by a factor of gamma.
If this interpretation is correct, then with gamma values of 60 or 100 the posterior variance would be reduced by factors of 60 or 100 even in the absence of any observational information. More generally, it appears that the posterior covariance from Eq. 11 can never be larger than the original prior covariance divided by gamma.
This makes me uncertain about how to interpret the narrow posterior distributions in Fig. 3. A substantial reduction in spread relative to the stated prior appears to be imposed directly by the choice of gamma, irrespective of any additional constraint provided by the GPP observations.
I can see that Eq. 11 would have a consistent probabilistic interpretation if introducing gamma is equivalent to redefining the prior covariance as the original prior covariance divided by gamma. However, I had understood gamma in the manuscript primarily as an algorithmic regularisation parameter used to keep the update within the region where the local linear approximation is useful.
Could the authors therefore clarify what statistical quantity the posterior uncertainties in Fig. 3 represent? In particular, I am not sure that the narrowness of these distributions can straightforwardly be interpreted as evidence that the observations have tightly constrained the parameters unless the probabilistic role of gamma is made clearer.
Overall, I found the paper interesting and the practical improvement from regularisation convincing. My questions are mainly about how that improvement should be interpreted statistically. A clearer discussion of the relationship between gamma and R, and particularly of the meaning of the posterior uncertainties when gamma is greater than 1, would substantially strengthen the manuscript.
Citation: https://doi.org/10.5194/egusphere-2026-132-RC2 -
AC3: 'Reply on RC2', Natalie Douglas, 21 Sep 2026
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2026/egusphere-2026-132/egusphere-2026-132-AC3-supplement.pdf
-
AC3: 'Reply on RC2', Natalie Douglas, 21 Sep 2026
Model code and software
Code for processing JULES model runs Natalie Douglas https://doi.org/10.5281/zenodo.19001565
Code for ORCHIDEE model Simon Beylat https://doi.org/10.14768/c68bc728-da71-4383-84df-dcde31d9c006
Code for processing ORCHIDEE model runs Simon Beylat https://doi.org/10.5281/zenodo.14609416
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 611 | 224 | 65 | 900 | 86 | 84 |
- HTML: 611
- PDF: 224
- XML: 65
- Total: 900
- BibTeX: 86
- EndNote: 84
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Review of
Regularisation of the 4DEnVar Data Assimilation method for Calibration of Land Surface Models.
by Douglas et al.
General comments:
This methodological paper focuses on improving a land surface model calibration data assimilation scheme. The authors attempt to enhance an existing method by incorporating a tunable parameter (gamma), which can be adjusted to optimise the efficiency of model parameter retrieval. The performance of the new method is evaluated using data from a Fluxnet forest site in the context of GPP (photosynthesis) observation assimilation. The method is applied to two different models: JULES and ORCHIDEE. Photosynthesis parameter values are retrieved for both models. In the case of ORCHIDEE, the retrieval also includes phenology and hydrology parameters. This well-written paper could be accepted after minor revisions.
Particular comment:
- L. 74-76: This sounds like a result or conclusion. In the introduction, it should be rephrased as a research question to be addressed.
- L. 76-78: Results should be interpreted in the Discussion section, not the Introduction. This sentence could be replaced with a research question that needs answering.
- L. 126 (Table 1): Replace "N/A" by "-" or by "Unitless". I am surprised that a key parameter such as rooting depth has not been included in the analysis. In what way does each model describe soil moisture stress?
- L. 127 (Section 2.2): Does this method account for errors in the atmospheric variables that the models use as input?
- L. 130: "observations given the target" is unclear. Do you mean "given the value of the target"?
- L. 260: The reason why this Fluxnet site has been chosen over another should be explained. Are soil moisture deficits observed at this site?
- L. 262: 'GPP observations' are used. It should be made clear that GPP is not directly observed. I suppose the GPP values you use are derived from NEE observations and an ecosystem respiration model?
- L. 327 (Fig. 5): RMSD values should have units. Year 2010 is interesting because both calibrated models are wrong. Why is there such a large delay in peak GPP for the two models?
- L. 386 (Fig. A1): Figure A1 is not mentioned in the text. The current caption of Figure A1 includes an interpretation. This interpretation should be incorporated into the main text. The reasons why the L-curve technique cannot be used to select gamma should be clearly explained.