the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 08 Jan 2025
Improving JULES Soil Moisture Estimates through 4D-En-Var Hybrid Assimilation of COSMOS-UK Soil Moisture Observations
Abstract. Accurate soil moisture estimates are essential for effective management and operational planning in various applications, including flood and drought response. However, soil moisture values derived from land surface models often exhibit significant deviations from in-situ observations. Data assimilation combines model information with observations to enhance prediction accuracy. Previous studies have typically focused on either estimating the initial soil moisture state or optimizing Pedotransfer Function (PTF) constants, which link soil texture to the hydraulic properties of the land surface models. In contrast, in this study, we employ a novel approach by performing joint state-parameter assimilation for the JULES model. We optimized both the PTF constants and the initial soil moisture conditions simultaneously. Using Four-Dimensional Ensemble Variational hybrid data assimilation, we ingested field-scale soil moisture observations from the Cosmic-ray Soil Moisture Monitoring Network across 16 diverse sites in the UK. The results demonstrate that joint state-parameter assimilation significantly enhances the accuracy of soil moisture estimates, improving the average Kling Gupta Efficiency values from 0.33 to 0.72 across different soil characteristics. These findings indicate that our proposed joint state-parameter assimilation framework holds great potential for enhanced predictive accuracy in land surface models.
- Preprint
(7539 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2024-3980', Anonymous Referee #1, 04 Mar 2025
General comments
The topic of the manuscript is important and of actual relevance. It deals with the assimilation of soil moisture values derived from cosmic-ray neutron sensing into a land surface model, and addresses improvements resulting from adapting the initial soil moisture state or adapting parameters within a set of pedotransfer functions or both. The manuscript is well-structured and well-written but in parts repetitive and vague in its description of its substance and potential advances. It lacks a reasonable consideration of its methodology in a way that gives full and clear picture to what has already been developed and published, and what is a new development in this study, accompanied by shortcomings in presenting and discussing results of the current study in respect to previous similar or possibly almost identical results. Though the topic itself offers potential to deeper insight and methodological advances, this potential has not truly been exploited. Some more details in the following.
Specific comments
- A big part of the methodology and the monitoring data set have been presented already in Cooper et al., 2021a (https://doi.org/10.5194/hess-25-2445-2021). The second scenario of the manuscript seems to be very similar to this work published by the manuscript’s second author as main author. It is using the same land surface model JULES, the same pedotransfer functions including the same starting parameter values, the same 16 cosmic-ray neutron locations out of the about 50 sites of the COSMOS-UK network, the same year for the data assimilation and the same following year for the forecast, and the same ensemble size. And the 4D-En-Var assimilation method outlined in the manuscript is depicted already in Pinnington et al. (2020), also including a combination of parameter vector and state vector, and seems to be part of LAVENDAR already together with JULES.
- While the study presented in the manuscript does go somewhat further by now also applying the data assimilation method to include the initial soil moisture state, it does not describe the methods and results of the manuscript on basis of the existing work but as an unclear mix with vague formulations of origin. These publications are cited but often rather as context. A clear distinction between existing work and own work is necessary. Maybe there are reasons to use exactly the same setting as before, though using another data set and approach could contribute more novelty, but this needs to be discussed and rectified. Also results should be discussed on basis of existing work, not in a diffuse way. To give one simple example, in the conclusions, second paragraph, the manuscript presents that KGE has improved from 0.33 to 0.66 for the parameter-only assimilation, which is identical to the statement in the conclusion of Cooper et al. (2021a) that ‘we see an improvement in the average KGE metric from 0.33 (range 0.10 to 0.69) before data assimilation to an average of 0.66 after data assimilation’. It should be referred to the previous result and made clear that this finding is identical to this previous result and why or to be discussed how it nevertheless may be different and why.
- Features of cosmic-ray neutron sensing are partly wrong and even contradictory within the manuscript. For example, in line 30 the horizontal footprint size is specified as ‘approximately 25 to 30 hectares’, in line 112 then ‘an area up to 120,000 m2’, which is only half of the former. Also, the depth specifics are not explained adequately and the fourth, deepest layer of the JULES model actually is not linked to the soil moisture observation by cosmic-ray neutron sensing at all, and the third layer likely only sometimes. And, the equations used for a weighted average of model layer soil moisture values to compare to the observed soil moisture (Cooper et al., 2021a, Pinnington et al., 2021) is a mere average accounting for the different layer thicknesses but not accounting for the strongly decaying weight with depth of cosmic-ray neutron sensing. Furthermore, the error estimate for the observations does not account for the relation between hourly values and a daily value for this cumulative measurement nor the Poisson distribution of its uncertainty instead of a Gaussian.
- In respect to the definition of the Cosby’s pedotransfer functions, there are also shortcomings. The manuscript refers to Cosby et al. (1984), but not everything presented is given there and seems neither developed within the manuscript’s study. Cosby et al. (1984) has reported linear relations between grain fractions and four hydraulic variables, but not the full mathematical equations as presented in 2.3. Therefore, a part of the earlier development seems to be missing. Marthews et al. (2014) could be cited directly in this respect, as one component. But further considerations would be helpful. And some discussion, why this set of pedotransfer functions? Only because they have been used in the similar preceding study (Cooper et al., 2021a)? And why not start with the parameters adjusted there already? How does it compare to other pedotransfer functions as used in other studies, etc.
Technical corrections, just some examples
- The title is full of abbreviations and unclear
- The introduction to monitoring of soil moisture starts with a general list of remote sensing sensors (and references) and rather outdated observation networks reported in 2006 and 2007. This part could be more to the point and up to date.
- In line 56 a reference is needed, as such and also for the claim to be more accurate.
- Line 175Â Â It is a bit uncalled-for to first give p a time index and then declare that it is constant in time.
- Line 245 to 253Â Â Replace Ne by the value chosen here (50), as mentioned anyway several times and as the other particular parameters are also specified as values.
- The references contain a large number of malformed doi links.
References here are named as cited in the manuscript.
Citation: https://doi.org/10.5194/egusphere-2024-3980-RC1 -
RC2: 'Comment on egusphere-2024-3980', Anonymous Referee #2, 22 Apr 2025
The paper presents an application of a data assimilation technique, (Ensemble-4DVar or 4D-En-Var) to the estimation of soil moisture at 16 sites by the JULES model. Importantly, I have concerns about the suitability of this technique in this application, and it is presented without any context - would other techniques be more appropriate? what might be the pros and cons? The 4D-En-Var technique is somewhat involved and has developed its own terminology (not the authors fault) but this is not explained in terms understandable for the general readership of HESS. There are many gaps, inconsistencies and, I think, errors in their description of the model and the algorithm, to the extent that this is on the borderline of automatic rejection. Similar work has been published previously, and the novelty here is not very clear. The context of how this relates to other work is missing e.g. have other models been developed using the COSMOS-UK data? However, correcting the inconsistencies and errors is a errors is a fairly simple matter. Explaining and justifying what they have done is a bigger task and requires substantial additional work and thought.
General points:
1. Since the initial states and parameters are both time-invariant, it would be simplest to consider both as "parameters". This would avoid the confusion about updating states (commonly called "nudging") which is usually associated with the term "data assimilation" but is not happening here. The application here is essentially parameter calibration (which is fine, only the description is confusing). For example "During state assimilation, the initial conditions are updated for each site" (line 286) invites confusion. States are not being assimilated or updated in the usual sense used in the context of data assimilation applied to NWP.Â
2. In what sense is this "4D"? There is only one time window, so there is nothing dynamic going on, carrying information across multiple time windows. The observational data form a single time series, but time is not treated as a special dimension with autocorrelation - they appear to considered as independent data. There is no horizontal spatial dimension to the data - the 16 sites are considered separately. Only a single vertical dimension is represented, so is this not "1DVar"?
3. Related to the above, variational DA is usually applied to initial value problems, where the model is highly sensitive to its starting state, as in numerical weather prediction (NWP). This is not the case here: the model is clearly not sensitive to its initial values (Fig. 3): the lines converge after 2 months and are indistinguishable for the 22 months thereafter. The difference between SC2 & 3 is almost literally zero in terms of absolute accuracy of prediction (Table 1). (I suspect KGE is inflating the differences by expressing relative differences.)
4. Given that the application here is essentially parameter calibration, this begs the question of whether "4D-En-Var" is the best method, perhaps over-kill or maybe inappropriate. Would standard optimisation or ML estimates be any better or worse, or would statistical time series analysis be more appropriate? In the light of this, I think it needs some justification as to why this method was chosen, and perhaps some comparison with more obvious choices.
5. The model is referred to as "the JULES community land surface model" but the focus here is purely a calculation of hydraulic conductivity, which is implicitly used to predict soil moisture via Darcy's law (though this is not actually presented). It is not clear if the other components of the JULES model are relevant. It provides a prediction of evapotranspiration, though perhaps there are local site estimates which might be more accurate, using lysimeters or observed LAI & Penman-Monteith etc. To what extent might the model-observation discrepancy be due to bias in the modelled evapotranspiration?
6. Using the HWSD instead of the locally measured soil texture seems an odd choice. A simpler interpretation of the results might be that the Cosby model parameters are fine, but the soil texture estimates from HWSD are wrong (as they will be for any point location), and the calibration results are simply the model parameters counter-balancing incorrect inputs. A comparison using the locally measured soil texture would discriminate between these options.
7. Essentially the algorithm calibrates 12 parameters, which relate soil texture to hydraulic properties, and so we are indirectly calibrating V_sat, b, sathh and satcon. How do the default & calibrated estimates of V_sat (i.e. porosity, easily calculated from bulk density) relate to the locally measured values at the COSMOS sites? This would be a more informative analysis, telling us whether the DA has produced a more accurate estimate of a mechanistically important and measurable parameter. As it stands, the paper is basically reporting that optimising the parameters to minimise RMSE improves RMSE (this is essentially what the cost function does, albeit weighted and penalised towards the defaults). That in itself is rather circular, so requires some analysis of *how* the parameters are changed, whether this is meaningful or arbitrary. Â
8. Are data for b, sathh and satcon available from the COSMOS sites? Soil water retention curves are fairly standard measurements for soil physics/hydrology. Â As above, the comparison between the measured values and those produced after model calibration would be of interest. Indeed, more interest, since this might reveal something of the underlying mechanisms.
9. To the non-JULES user, the distinction between the Cosby model and the JULES model is very arbitrary. For practical purposes, you are calibrating 12 model parameters. Referring to these as a separate "pedotransfer function" and PTF constants does not help explain conceptually what is going on here.
Specific points:l108. 50 sites? Why use only 16? Computational reasons? How were the 16 sites chosen?
l123. Eqn 1: this is just mass balance, not the Richards' equation. Perhaps leave out the R_bl runoff term since you assume it to be zero, so as not to mislead. Where does precipitation fit in? You show the outputs but not the input.
l124, 128, Eqns 1-3. "soil moisture content" can be defined in many ways. I presume we are talking here about volumetric content in m^3 / m^3, since that is what CNRS is sensitive to. Please define as such, with units.
l131+ throughout. Assuming the above is correct "soil moisture content" is denoted with three different but duplicate symbols: theta, V and SM. Why? This is just horrendously confusing. Just use one - theta is fairly standard. In Eqns 2. "soil moisture content" is defined as theta but its value in saturated conditions is defined with a different symbol, V. This obfuscates the point that the LHS is a fraction-of-the-whole term, i.e. soil moisture content as a fraction of its maximum possible value.
l131+ throughout, Eqns 1-3. Similarly, matric suction is defined as psi, but its value in saturated conditions is defined with as "sathh", producing the same obfuscation as above. (And multi-letter symbols are generally frowned up.) Units?
l131+ throughout, Eqns 1-3. Similarly, hydraulic conductivity is defined as K, but its value in saturated conditions is defined with as "satcon", producing the same obfuscation as above.
l146, Eqn 4. K has already been defined as hydraulic conductivity, and now it is given a different meaning! Use different symbols for different quantities, (but the same symbol when referring to the same quantity).
l150. V_sat is defined as the maximum "water-holding capacity of the soil". Given the different wording, It's not clear if this actually refers to volumetric soil moisture content or not.
l158. Same points apply to V_crit and V_wilt.
l158. Are V_crit and V_wilt relevant? They are only mentioned once hereafter, and no field measurements are presented for comparison.
Eqn 8. Is the denominator 3.364 related to -33 kPa? Make the link clear if so; explain if not.
Eqn 9. Is the denominator 152.9 related to -1500 kPa? Make the link clear if so; explain if not.
l164. Is the soil heat storage relevant here? It is not mentioned again I think.
Section 2.3/4. We never return to Eqn 1 to explain how the derived parameters are used to calculate the water flux between vertical layers. How do the equations and parameters presented actually result in a prediction of soil moisture in different layers. Where does precipitation fit in?l172+, Section 2.4. Can you describe in words for the non-expert what the algorithm is doing. The description here is couched in the terminology of 4D-Var, which is not self-explanatory for the general readership of HESS. e.g. I think Eqn 12 just says that the predicted states at time t are modelled as a function of their states at t-1 and some parameters p, but that is not obvious to the casual reader because of the language used.
l178+. It is clearly stated that the parameters do not vary in time. Why then are they shown throughout with the subscript t, which implies precisely that they *do* vary in time? Very confusing.
l181. The "observation operator" h_t() requires explanation. How does it "map the state vector to the observation space"? What does this mean? Subscript t implies that the function h varies in time. Explain how this is so.
l88+. "encapsulates" and "such as" implies there are other variables as well. I don't think that is so, but needs clarifying. The control vector includes parameters, not just "variables" as stated here.
l191. You don't "aim to estimate the initial soil moisture states" - these are not the quantities of interest. These are the initial states that are allowed to vary / be updated.
l192. By definition, the initial values do not vary in time. Why then are they shown throughout with the subscript t, which implies precisely that they do vary in time?Â
l93. Hazarding a guess at what is actually going on, contrary to what is written, I think x_a is the "control" vectors of parameters which are varied in the algorithm, whilst x_t seems to be the states that are compared with observations, and are thereby different things given the same symbol "x". I think x_t = [z_t]^T in all scenarios, i.e. it is always the current soil moisture states that are compared with observations. It is the make-up of x_a and x_b that differ among scenarios, whether they consist of initial soil moisture states, soil parameters, or both.Â
l200. Don't add synonyms for no reason - use either "background", "before" or "prior"; and  "analysis", "after" or "posterior". I'd avoid the last to avoid confusion with Bayesian estimates (this seems to be a penalised maximum likelihood estimate).
l201-217. Explain what the "background error covariance matrix" is in plain English. In this context, it seems to quantify the uncertainty in the default parameters. How then can it be estimated from "a previous forecast" (line 208) since we can never forecast parameters or initial state values? The explanation in terms of the perturbation matrix, x'_b does not make sense to me (perhaps my maths limitations). "Ensemble members" are mentioned without any explanation of what this is at this point.l203, Eqn 14. Can you give a plain English explanation of the cost function? Looks like a weighted average of (i) deviation of parameters from their default values and (ii) deviation of predictions from observed values, with the weighting coming from the error matrix which indicates how much uncertainty we should attach to each.
l224. "The posterior estimate" is easily confused with "posterior" in the Bayesian sense, but I don't think that is what is meant here. Possibly this is just the terminology used in 4DVar, but I think this is essentially a penalised maximum likelihood estimate). Better to refer to it as something like "optimised" or clarify how this sits with the Bayesian sense of estimating the posterior distribution if it is actually meant in that sense.
l233. Are the 12 parameters independent? Probably not. Is this covariance represented in the B matrix? Â
l234. This seems to be a non sequitur as a justification of guessing +/- 10%.
l236. HWSD is coarse and rather unreliable for estimates at a specific point. Â Why not use the locally measured soil texture values, or better, a comparison of sensitivity among different estimates of these inputs?
l245. So are the "initial values" considered to be 1 Jan 2016 or 1 Jan 2017? The latter would make the spin-up irrelevant. I'm confused.
l260. sigma = 50% in the observations, so these are highly uncertain, and not a strong constraint on the model. This is worth highlighting.
l269. Which "modelled soil moisture values" are used? There are four layers but only the top two correspond to the CNRS measurement, depending on soil wetness. Implied that some weighting is done but not explained.
l271, Eqn 21. How does this differ from h_t (x_t ) − y_t in Eqn 14? Or is just a different symbol for soil moisture?
l285. effect on "JULES performance". Semantics perhaps, but I think here we are looking at the effect of the performance of the calibration algorithm.
Fig. 3. The model is clearly not sensitive to its initial values: the lines converge after 2 months and are indistinguishable thereafter. This is very clear from Table 1 SC2 vs SC3, and this result needs to be emphasised. How is this reconciled with the distributions in Fig 3 and the RMSE values in Table 1 (Prior vs SC1)? At both sites, layer 3 has the opposite trend to the rest. Is there some kind of numerical artefact going on here? KGE gives the wrong impression, being based on relative differences when it is the absolute difference that matters.Â
Citation: https://doi.org/10.5194/egusphere-2024-3980-RC2
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 182 | 43 | 8 | 233 | 8 | 15 |
- HTML: 182
- PDF: 43
- XML: 8
- Total: 233
- BibTeX: 8
- EndNote: 15
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1