Using data assimilation to optimize pedotransfer functions using field-scale in situ soil moisture observations

Soil moisture predictions from land surface models are important in hydrological, ecological, and meteorological applications. In recent years, the availability of wide-area soil moisture measurements has increased, but few studies have combined model-based soil moisture predictions with in situ observations beyond the point scale. Here we show that we can markedly improve soil moisture estimates from the Joint UK Land Environment Simulator (JULES) land surface model using field-scale observations and data assimilation techniques. Rather than directly updating soil moisture estimates towards observed values, we optimize constants in the underlying pedotransfer functions, which relate soil texture to JULES soil physics parameters. In this way, we generate a single set of newly calibrated pedotransfer functions based on observations from a number of UK sites with different soil textures. We demonstrate that calibrating a pedotransfer function in this way improves the soil moisture predictions of a land surface model at 16 UK sites, leading to the potential for better flood, drought, and climate projections.

1 Introduction

Soil moisture is an important physical variable, significant in agriculture (Pinnington et al., 2018), flood events (Koster et al., 2010; Berghuijs et al., 2019), and processes related to weather and climate (Seneviratne et al., 2010). Land surface models such as the Joint UK Land Environment Simulator (JULES) can be used to make predictions of soil moisture and generally rely on empirical pedotransfer functions (PTFs) to relate readily available or easy-to-measure soil characteristics such as soil texture to the soil hydraulics parameters required by the model (e.g. Van Looy et al., 2017)

There are a number of different types of pedotransfer function, as noted in Van Looy et al. (2017) and Hodnett and Tomasella (2002), with different inputs and outputs depending partly on the representation of soil physics processes of the chosen land surface model. In “class” approaches, soil types are clustered into groups, and hydraulic model parameters are then obtained from a lookup table (Wösten et al., 1999); this results in discrete soil hydraulics parameter sets. Alternatively, continuous pedotransfer functions take soil characteristic information from each sample of interest and apply the function to produce continuous soil hydraulics parameter sets (e.g. Cosby et al., 1984; Hodnett and Tomasella, 2002; Schaap et al., 2001).

To date, pedotransfer functions have been derived by fitting to results from field or laboratory experiments on point- or small-scale soil samples (centimetre to metre), despite the fact that land surface models are generally applied at larger (field to kilometre) scales. The recent development of novel in situ techniques for measuring soil moisture over field rather than point scale presents an opportunity to test whether land surface models, in conjunction with commonly used pedotransfer functions, are able to reproduce field-scale soil moisture observations.

In this paper, we have compared JULES soil moisture predictions with soil moisture observations from the COSMOS-UK dataset (Stanley et al., 2021); these observations are measured by cosmic ray neutron sensor (CRNS) instruments over a footprint of up to 120 000 m². We have then used the LaVEnDAR four-dimensional ensemble variational data assimilation framework (Pinnington et al., 2020) to combine COSMOS-UK soil moisture observations at 16 sites with equivalent JULES soil moisture estimates. We have thereby optimized constants in the Cosby pedotransfer function (Cosby et al., 1984). This results in a newly calibrated set of pedotransfer functions based on field-scale soil moisture observations across 16 sites with a range of soil types. This approach allows us to test whether we can improve the performance of the model by optimizing the pedotransfer functions for larger scales using field-scale soil moisture observations. Our approach also allows for comparison of the soil hydraulics parameters generated using field-scale (∼ hundreds of metres) soil moisture measurements with those generated by the original pedotransfer functions, which are based on small-scale (∼ centimetres) measurements. We chose to optimize the pedotransfer functions rather than directly optimizing soil physics parameters since this preserves the physical relationships between soil physics parameters that the pedotransfer functions describe. This approach also has the advantage that we can assimilate observations from all sites simultaneously to produce one set of pedotransfer functions applicable at all 16 study sites. The same pedotransfer function could then potentially be applied anywhere in the UK that soil texture information is available.

We use CRNS soil moisture measurements in this study. Larger scale soil moisture measurements are also increasingly available from satellite products, and these have been used to good effect in data assimilation frameworks with land surface models (e.g. Pinnington et al., 2018; Liu et al., 2011; De Lannoy and Reichle, 2016; Yang et al., 2016). The advantage of the CRNS measurements used here is that they provide a more direct soil moisture measurement than those from satellites. CRNS soil measurements are also representative of depths of approximately 10 to 30 cm, compared to the top 5 to 10 cm for satellite retrievals.

An alternative approach to assimilate CRNS soil moisture measurements into land surface models is taken in Brunetti et al. (2019), Han et al. (2015), and Mwangi et al. (2020). These studies use neutron counts from CRNS instruments as observations, combined with the COSMIC method presented in Shuttleworth et al. (2013) to map modelled soil moisture estimates into equivalent neutron counts. In this study we instead directly compare modelled and CRNS-derived soil moisture.

The rest of the paper is organized as follows: in Sect. 2, we outline the JULES land surface model and the COSMOS-UK data used in this study; we also describe the data assimilation experiment we have performed and introduce the metric we deployed to measure how well the model fits the observations. In Sect. 3, we present results, showing that we can use COSMOS-UK observations from 2017 to improve the fit between the JULES model output and observations over 2 years at all the sites we included. We discuss our results in the context of changes in the JULES soil physics parameters in Sect. 4. In Sect. 5, we conclude that it is possible to optimize pedotransfer functions with field-scale soil moisture measurements and that this markedly improves the fit of JULES soil moisture estimates to COSMOS-UK observations.

2 Methods

2.1 JULES land surface model

JULES uses the Darcy–Richards equation to model soil hydraulic processes (Best et al., 2011), so that the downward water flux, W, between adjacent soil layers is given by

$$\begin{array}{}\text{(1)}& W=K\left({\displaystyle \frac{\partial \mathrm{\Psi }}{\partial z}}+\mathrm{1}\right),\end{array}$$

where Ψ is the soil matric suction, K is the soil hydraulic conductivity, and z is the distance from the soil surface in the vertical direction.

JULES provides two options for representing the relation between soil water content, θ, matric suction, and hydraulic conductivity; in this paper we use the Brooks and Corey soil physics option (Best et al., 2011; Brooks and Corey, 1964), where we assume

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\text{s}}}}={\left({\displaystyle \frac{\mathrm{\Psi }}{{\mathrm{\Psi }}_{\text{s}}}}\right)}^{-\frac{\mathrm{1}}{b}}\end{array}$$

and

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{K}{K_{\text{s}}}}={\left({\displaystyle \frac{\mathit{\theta }}{{\mathit{\theta }}_{\text{s}}}}\right)}^{\mathrm{2}b+\mathrm{3}}.\end{array}$$

In Eqs. (2) and (3), θ_s, K_s, and Ψ_s are values of soil moisture, hydraulic conductivity, and soil matric suction at saturation; b is a soil-dependent constant with a value usually determined through a pedotransfer function. The soil physics parameters used in the implementation of Brooks and Corey soil physics in JULES are briefly described in Table 1; more details are available in Best et al. (2011) or in the JULES user guide (2020).

Table 1 Soil physics parameters.

The values of the eight soil physics parameters outlined in Table 1 are generally calculated via a set of pedotransfer functions. Here we use the Cosby pedotransfer functions, which have the following mathematical formulation (Cosby et al., 1984; Marthews et al., 2014):

$$\begin{array}{}\text{(4)}& {\displaystyle }b={\mathit{\kappa }}_{\mathrm{1}}+{\mathit{\kappa }}_{\mathrm{2}}{f}_{\text{clay}}-{\mathit{\kappa }}_{\mathrm{3}}{f}_{\text{sand}}\text{(5)}& {\displaystyle }{\mathit{\theta }}_{\text{s}}={\mathit{\kappa }}_{\mathrm{4}}-{\mathit{\kappa }}_{\mathrm{5}}{f}_{\text{clay}}-{\mathit{\kappa }}_{\mathrm{6}}{f}_{\text{sand}}\text{(6)}& {\displaystyle }{\mathrm{\Psi }}_{\text{s}}=\mathrm{0.01}\times {\mathrm{10}}^{{\mathit{\kappa }}_{\mathrm{7}}-{\mathit{\kappa }}_{\mathrm{8}}{f}_{\text{clay}}-{\mathit{\kappa }}_{\mathrm{9}}{f}_{\text{sand}}}\end{array}$$

$$\begin{array}{}\text{(7)}& {\displaystyle }{K}_{\text{s}}={\mathrm{10}}^{-{\mathit{\kappa }}_{\mathrm{10}}-{\mathit{\kappa }}_{\mathrm{11}}{f}_{\text{clay}}+{\mathit{\kappa }}_{\mathrm{12}}{f}_{\text{sand}}}\times {\displaystyle \frac{\mathrm{25.4}}{\mathrm{3600}}}\text{(8)}& {\displaystyle }{\mathit{\theta }}_{\text{crit}}={\mathit{\theta }}_{\text{s}}{\left({\displaystyle \frac{{\mathrm{\Psi }}_{\text{s}}}{\mathrm{3.364}}}\right)}^{\mathrm{1}/b},\text{(9)}& {\displaystyle }{\mathit{\theta }}_{\text{wilt}}={\mathit{\theta }}_{\text{s}}{\left({\displaystyle \frac{{\mathrm{\Psi }}_{\text{s}}}{\mathrm{152.9}}}\right)}^{\mathrm{1}/b}\end{array}$$

$$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}& h_{\text{cap}}=(\mathrm{1}-{\mathit{\theta }}_{\text{s}})(\mathrm{2.373}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}{f}_{\text{clay}}+\mathrm{2.133}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}{f}_{\text{silt}}\\ & +\mathrm{2.133}\mathrm{\times }{\mathrm{10}}^{\mathrm{6}}{f}_{\text{sand}})\end{array}\text{(11)}& {\displaystyle }\begin{array}{rl}& h_{\text{con}}={\mathrm{0.025}}^{{\mathit{\theta }}_{\text{s}}}\times {\mathrm{1.16}}^{f_{\text{clay}}(\mathrm{1}-{\mathit{\theta }}_{\text{s}})}\times {\mathrm{1.57}}^{f_{\text{sand}}(\mathrm{1}-{\mathit{\theta }}_{\text{s}})}\\ & \times {\mathrm{1.57}}^{f_{\text{silt}}(\mathrm{1}-{\mathit{\theta }}_{\text{s}})},\end{array}\end{array}$$

where f_clay, f_sand, and f_silt are fractions of clay, sand, and silt in the soil, by weight. Equations (8) and (9) are rearrangements of Eq. 2 at fixed values of matric suction corresponding to the wilting and critical points. Equation (10) is a linear combination of the assumed heat capacities of sand, silt, and clay, weighted by their relative fractions, and Eq. 11 is as given in Dharssi et al. (2009). The values of the constants κ₁ to κ₁₂ usually used in Eqs. (4) to (11) are those given in Cosby et al. (1984); we present them in Table 2. These values are empirically determined from 1448 small soil samples (centimetre dimensions) taken from 23 states in the United States (for further details of the soil samples and sampling methods, see Rawls, 1976, and Holtan, 1968). The values of the constants given here match those in Marthews et al. (2014) (with soil fraction multipliers adjusted for fraction, rather than percentage, of soil by weight).

Cosby et al. (1984)

Table 2 Values of the constants commonly used in the Cosby pedotransfer functions.

JULES requires meteorological driving data to produce soil moisture estimates. The required input variables are air pressure, air temperature, humidity, downward fluxes of shortwave and longwave radiation, precipitation, and wind speed. In this paper we have used half-hourly meteorological observations measured at COSMOS-UK sites as driving data; in this way we can use JULES to give soil moisture predictions at any COSMOS-UK sites with sufficiently complete meteorological data.

JULES provides estimates of soil moisture at various depths; in the standard configuration used here these correspond to four layers, with depths [0,10 cm], [10 to 35 cm], [35 to 100 cm], and [100 to 300 cm]. The JULES layers are often referred to by their thicknesses, which are 10, 25, 65, and 200 cm respectively. Here, we refer to the soil moisture estimates for the four layers as SM₁₀, SM₂₅, SM₆₅, and SM₂₀₀.

2.2 COSMOS-UK soil moisture data

The COSMOS-UK project comprises a network of soil moisture monitoring stations across the United Kingdom, providing long-term soil moisture measurements at around 50 sites. Data for 2013 to 2019 are available in the EIDC archive (Stanley et al., 2021). Soil moisture observations are made using an innovative cosmic ray neutron sensor (CRNS) instrument at each site; these provide a measurement of soil moisture over an area of up to 120 000 m² (∼ 30 acres) (Antoniou et al., 2019; Evans et al., 2016). The CRNS at each site counts fast neutrons within the sensor's footprint. These counts are corrected for local meteorological conditions using in situ measurements and also background neutron intensity using data from a neutron monitoring station (Evans et al., 2016). The corrected counts are then calibrated for site-specific soil properties determined from destructive soil sampling conducted after site installation. Soil samples were collected from each site following Köhli et al. (2015) and were returned to UKCEH for laboratory analysis. The results were used to determine reference soil moisture, lattice and bound water, bulk density, and organic matter for the day of sampling and are subsequently used to derive soil water content from the corrected CRNS counts. The majority of sites explored in this study are grasslands, and it is therefore expected that CRNS soil moisture results are not significantly affected by seasonal changes in biomass (Baatz et al., 2014).

We have used daily averaged soil moisture data from 16 COSMOS-UK sites as observations in this paper. The sites were selected based on completeness of soil moisture and meteorological data over a 3-year period from 2016–2018 and are listed in Table 3, with details of land cover and broad soil descriptions taken from Antoniou et al. (2019). Locations of the sites are shown in Fig. 1. For more details of each of the sites, see Antoniou et al. (2019). The Cosby pedotransfer function was designed to work for mineral soils, and the CRNS calibration is most reliable at sites with minimal vegetation. We therefore consider that the first seven sites listed in Table 3 are those at which the JULES model can be expected to provide a good match to observations via our chosen PTF; soil types and land cover at the remaining sites mean that JULES may not be able to represent the observed soil moisture time series as accurately.

Table 3 COSMOS-UK sites selected for this study. “Heath” indicates some shrubs are present at the site.

Figure 1 Locations of COSMOS-UK sites used in this study.

Both the depth and the footprint over which the CRNS measures soil moisture change with soil moisture (Evans et al., 2016; Köhli et al., 2015; Antoniou et al., 2019), with the footprint and depth of the measurement both becoming smaller as soil moisture increases. The COSMOS-UK dataset (Stanley et al., 2021) includes estimates of the depth over which each daily soil moisture value is valid, known as a D86 value. Measurements of several other environmental variables are made at COSMOS-UK sites, using a suite of instrumentation. These include point soil moisture and temperature measurements at various depths in the soil and meteorological variables. We have used half-hourly in situ meteorological data from the COSMOS-UK dataset as driving data for the JULES model.

2.3 Data assimilation

Data assimilation is a group of methods in which information from models and observations is combined in order to give the best estimate of the state of a physical system and/or model parameter values. In this paper, we have used the four-dimensional ensemble variational data assimilation technique, LaVEnDAR, which is introduced in Pinnington et al. (2020) and is based on Liu et al. (2008). We use LaVEnDAR to optimize 12 constants, κ₁ to κ₁₂, in the Cosby pedotransfer functions (Eqs. 4 to 11) based on estimates of soil moisture from JULES and corresponding field-scale observations of soil moisture from COSMOS-UK. LaVEnDAR optimizes κ₁ to κ₁₂ here by minimizing a cost function with two terms. The first term is a measure of the difference between the observed and modelled soil moisture, and the second term is a measure of the difference between prior and posterior values of κ₁ to κ₁₂.

The values of κ₁ to κ₁₂ are assumed to be constant in time and space; the same values are used across all sites to generate soil JULES moisture estimates via the pedotransfer functions.

2.4 Experimental details

In order to use COSMOS-UK data with JULES outputs in the LaVEnDAR scheme, we require both sets of soil moisture values to correspond to the same soil depth. We have therefore devised a weighted depth approach, in which we extract from each JULES prediction an average soil moisture corresponding to the UK-COSMOS observed depth. The observed depth changes with soil moisture and with horizontal distance from the CRNS instrument; here we have used the reported observation depth at 75 m from the CRNS (in the horizontal direction). For each day, we calculate a depth-adjusted JULES soil moisture estimate, SM_depth, depending on the 75 m observation depth value, D86, provided for that day, such that

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\text{SM}}_{\text{depth}}=\\ & \left\{\begin{array}{ll}{\text{SM}}_{\mathrm{10}},& \text{if}\mathrm{D}\mathrm{86}\le \mathrm{10}\text{cm},\\ \frac{\mathrm{10}}{\mathrm{D}\mathrm{86}}{\text{SM}}_{\mathrm{10}}+\frac{(\mathrm{D}\mathrm{86}-\mathrm{10})}{\mathrm{D}\mathrm{86}}{\text{SM}}_{\mathrm{25}}& \text{if}\mathrm{10}\text{cm}<\mathrm{D}\mathrm{86}\le \mathrm{35}\text{cm},\\ \frac{\mathrm{10}}{\mathrm{D}\mathrm{86}}{\text{SM}}_{\mathrm{10}}+\frac{\mathrm{25}}{\mathrm{D}\mathrm{86}}{\text{SM}}_{\mathrm{25}}+\frac{(\mathrm{D}\mathrm{86}-\mathrm{35})}{\mathrm{D}\mathrm{86}}{\text{SM}}_{\mathrm{65}},& \text{if}\mathrm{35}\text{cm}<\mathrm{D}\mathrm{86}\le \mathrm{65}\text{cm},\end{array}\right.\end{array}\end{array}$$

where SM₁₀, SM₂₅, and SM₆₅ are the JULES-predicted soil moisture values from the [0,10 cm], [10 to 35 cm], and [35 to 100 cm] layers respectively, and the D86 value is given in centimetres. In this way, thickness-weighted contributions to the soil moisture are taken from every JULES layer which would be wholly or partly contained within the D86 depth. We have not taken the COSMOS-UK variable footprint into account in this study.

In this paper we have used an ensemble size of 50, as in related experiments in Pinnington et al. (2020) and Liu et al. (2008). In order to implement the LaVEnDAR scheme we conducted the following:

1. We generated a 50-member ensemble of each of the 12 PTF constants κ₁ to κ₁₂. These were obtained by sampling from a Gaussian distribution centred on the value given in Table 2, with standard deviation equal to 10 % of the mean. This standard deviation value was chosen fairly arbitrarily; future work could assess the sensitivity of the results to the values chosen for each PTF constant.

2. We assembled 50 unique sets of 12 constants κ₁ to κ₁₂.

3. We used each unique set of constants in Eqs. (4) to (11) to generate 50 sets of soil physics parameters for each site. Soil texture information for each site was taken from the Harmonized World Soil Database (HWSD) (Fischer et al., 2008).

4. We used the soil parameter sets to run 50 realizations of JULES at each of our selected sites over a 2-year time window to create a prior ensemble of 50 soil moisture time series per site.

5. We used the LaVEnDAR scheme to generate a new, posterior ensemble of values for each of the 12 PTF constants, taking into account COSMOS-UK soil moisture observations from 2017. Here, we assumed uncorrelated observation errors of 50 % of the mean soil moisture value at each site.

6. We used the new posterior ensemble of PTF constants to generate 50 posterior sets of soil physics variables at each site.

7. We ran 50 posterior realizations of JULES at each site to create posterior soil moisture time series.

These steps are also shown in schematic form in Fig. 2.

We assume that the soil texture values from the HWSD are correct; they are not changed during the data assimilation process. We used a global soil texture dataset since site-specific soil texture observations were not available; using a global soil texture dataset also has the advantage that our method then has the potential for extension to other UK areas without local soil measurements. Other open-source global soil texture products are also available, e.g. SoilGrids (Hengl et al., 2017). We acknowledge that there can be discrepancies between the HWSD and local measurements (e.g. Zhao et al., 2018), but our decision to use the HWSD here follows recent successful integration of soil texture data from the HWSD with JULES in studies such as Martínez-de la Torre et al. (2019), Ritchie et al. (2019), and Ehsan Bhuiyan et al. (2019).

Figure 2 Schematic showing data assimilation experimental design; % ssc refers to site-specific fractions of sand, silt, and clay in the soil. In this study, only observations from 2017 (at each site) were used in the assimilation algorithm.

We have assumed a high observation error value in this experiment. The daily soil moisture measurements we use are averaged from hourly soil moisture measurements; analysis of these data shows that the standard deviation of the hourly data around the daily mean is approximately 20 %. We have inflated this here to 50 % observation error; we note that similar observation error covariance inflation techniques have been used in, for example, assimilation of satellite observations in numerical weather prediction (Fowler et al., 2018; Hilton et al., 2009). The reason for inflating the observation error is essentially because we found that smaller observation error values impacted negatively on the posterior soil moisture results. We suggest that inflation of the observation error is necessary here to compensate for otherwise neglected sources of error (e.g. the error in converting neutron counts to soil moisture) and for the assumption of uncorrelated observation error; in fact there will likely be intra-site correlations between observation errors due to site-specific instrument calibration.

We have used COSMOS-UK measurements from 2017 only in our data assimilation experiments but compared the prior and posterior JULES runs from 2017 and 2018 with observations.

Figure 3 Observed and modelled (ensemble mean) soil moisture time series at Bickley Hall (BICKL). The dotted line separates the period over which observations have been used for assimilation (2017) from the period in which no observations have been assimilated (2018).

2.5 Metrics

In order to assess how well our prior and posterior JULES runs match COSMOS-UK observations, we require a metric. Here we have used the Kling–Gupta efficiency metric, as described in Gupta et al. (2009) and Knoben et al. (2019), to compare the goodness of fit between observed and modelled (ensemble mean) soil moisture times series. The Kling–Gupta efficiency (KGE) is given by

$$\begin{array}{}\text{(13)}& \text{KGE}=\mathrm{1}-\sqrt{(\mathrm{1}-r{)}^{\mathrm{2}}+(\mathrm{1}-\mathit{\alpha }{)}^{\mathrm{2}}+(\mathrm{1}-\mathit{\beta }{)}^{\mathrm{2}}},\end{array}$$

where

$$\begin{array}{}\text{(14)}& \mathit{\alpha }={\displaystyle \frac{{\mathit{\sigma }}_{\text{model}}}{{\mathit{\sigma }}_{\text{obs}}}}\end{array}$$

and

$$\begin{array}{}\text{(15)}& \mathit{\beta }={\displaystyle \frac{{\mathit{\mu }}_{\text{model}}}{{\mathit{\mu }}_{\text{obs}}}}.\end{array}$$

In Eqs. (14) and (15), μ_model and μ_obs are mean values of the modelled and measured soil moisture time series respectively; σ_model and σ_obs are the standard deviations in the modelled and observed soil moisture time series. The value of r is the Pearson correlation coefficient between the model and the observation time series data and can vary between −1 (anti-correlation) and 1 (perfect correlation), with a score of 0 indicating no correlation. The value of α reflects how well the spread in the modelled soil moisture values matches that of the observations, with a value of 1 corresponding to perfect matching. Equation (15) shows that the value of β represents bias between the model and observations, with a value of 1 indicating zero bias. Since α and β can be larger or smaller than 1, the value of the KGE can range between 1 (perfect model fit to data) and very large negative values. In Knoben et al. (2019), the authors argue that while in some studies a threshold of KGE ≥ 0 has been used to denote “good” model performance, a lower threshold of KGE ≥ −0.41 is required for the model to perform better than a mean persistence forecast. We used Python 3.7.1 to calculate metrics and prepare plots.

Figure 4 Observed vs. modelled (ensemble mean) soil moisture at Bickley Hall (BICKL) for prior and posterior JULES runs. The diagonal line shows the 1:1 perfect correspondence line. The correlation coefficient at this site changed from 0.93 (prior) to 0.94 (posterior), and the RMSE reduced from 0.13 (prior) to 0.03 (posterior).

3 Results

3.1 Effect of data assimilation on JULES soil moisture predictions

Figures 3 to 6 show measured and modelled soil moisture time series for 2017 and 2018 at two representative COSMOS-UK sites. In all cases,, the modelled soil moisture series is the ensemble mean. These figures show that the JULES runs using posterior PTF constants produce soil moisture estimates which are a better match to the observations than the JULES runs using the prior PTF constants. Figures 3 and 4 show results from Bickley Hall (BICKL), which is a site at which we expect soil moisture to be well represented by JULES via the Cosby PTF (this site has a typical mineral soil). Figures 5 and 6 represent results from a site at which the high organic content of the soil and the presence of trees mean that we do not expect our JULES setup to match the observations so successfully.

Figure 5 Observed and modelled (ensemble mean) soil moisture time series at Gisburn Forest (GISBN). The dotted line separates the period over which observations have been used for assimilation (2017) from the period in which no observations have been assimilated (2018).

Figure 7a shows the KGE values for prior and posterior JULES runs at all 16 sites included in our study. These metrics show how closely the prior and posterior JULES runs match the observations over the period of 2017 and 2018 before and after assimilation of observations from 2017. Figure 7a shows that data assimilation markedly improves the fit to observations at all sites according to the Kling–Gupta metric; all the analysis Kling–Gupta efficiency scores are closer to the ideal value of 1 than the prior values. We note that for all sites, the match between model and measurements is better in 2017 and 2018, even though only observations from 2017 were used in the optimization process. This indicates that the new values for the PTF constants allow JULES to simulate field-scale soil moisture measurements better than the original (prior) PTF constants. Figure 7b shows that the prior and posterior correlation coefficients, r, are very similar at most sites, although there is a slight deterioration of the correlation coefficient at the majority of the sites. Despite this, the reduction in r is very small compared to the overall improvement in the KGE metric at all sites, and the prior and posterior r values are all greater than 0.8 at sites with a typical mineral soil. The r value stays low at Moorhouse (MOORH), perhaps because the soil at this site is too highly organic for the Cosby parameters to really be applicable and for the COSMOS-UK measurements to be reliable. The r value also stays low at Gisburn Forest (GISBN), which is likely due to the fact that there are a large number of trees at this site. The presence of above-ground biomass may make the site-specific calibration less reliable than at other sites (Baatz et al., 2014). The high organic carbon content of the soil at Gisburn Forest likely also contributes to this as our chosen PTF is designed to work best with mineral soils. Interception is another process which potentially complicates the calibration at sites with vegetation, although the authors of Bogena et al. (2013) report that water intercepted by the canopy constitutes a negligible amount of the water detected in the CRNS footprint, even in coniferous forests.

Figure 6 Observed vs. modelled (ensemble mean) soil moisture at Gisburn Forest (GISBN) for prior and posterior JULES runs. The diagonal line shows the 1:1 perfect correspondence line. The correlation coefficient at this site changed from 0.73 (prior) to 0.69 (posterior), and the RMSE reduced from 0.25 (prior) to 0.15 (posterior).

Figure 7c shows that a significant contribution to improved KGE at all sites comes from improvement in the alpha component, which is much closer to the ideal value of 1 for all of the posterior JULES runs than the prior JULES runs. The alpha component represents how well the spread in the model matches the spread in the observations. We saw in time series plots such as Figs. 3 and 5 that the spread in JULES soil moisture was too small at all sites; our results show that the data assimilation has acted to correct this by updating the value of the PTF constants. Figure 7d shows that the beta parameter is closer to the ideal value of 1 after data assimilation than before at all sites except for Cardington; i.e. data assimilation is correcting a bias in the JULES outputs at all but one site. The prior bias at Cardington is in the opposite direction to bias at all of the other sites.

Figure 7 Kling–Gupta efficiency scores for JULES runs using prior and posterior PTF variable values. Horizontal dotted lines show the value of the metric for a perfect match between the model and observations.

Figure 8 Ensemble prior (orange) and posterior (blue) parameter values at each site. These are topsoil results, which we have assumed to correspond to the top two soil layers in JULES (0–35 cm depth from the surface).

3.2 Effect of data assimilation on JULES soil physics parameters

The data assimilation algorithm in this study acts directly on the PTF constants κ₁–κ₁₂ which make up the state vector. The resulting changes to the JULES soil physics parameters through Eqs. (4)–(11) are presented here in Sect. 3.2. Figures 8 and 9 show changes to the eight JULES soil physics parameters used for the topsoil and subsoil layers respectively (Sect. 3.3 shows how the underlying PTF constants are updated).

Figure 9 Ensemble prior (orange) and posterior (blue) parameter values at each site. These are “subsoil” results, which we have assumed to correspond to the deeper two soil layers in JULES (35–300 cm depth from the surface).

Figure 10 Prior and posterior PTF constant value distributions. Orange shows prior and blue posterior. The blue line shows the original value of the constant as in Table 2.

Figures 8 and 9 show that the mean value of K_s (satcon) gets smaller (4 to 5 times smaller) at each site after data assimilation and that the posterior distribution of the K_s (satcon) parameter is narrower than the prior distribution. The results in Figs. 8 and 9 also show that the site-to-site variability of the b parameter reduces following data assimilation; the largest mean prior values of b are reduced, and the distributions with the smallest mean values are shifted to slightly larger values. Figures 8 and 9 show that the mean value θ_s (vsat) has increased at all the sites following data assimilation, and the distribution of θ_s (vsat) at each site has become much narrower. The mean values of the θ_crit (vcrit) and θ_wilt (vwilt) distributions have stayed broadly similar or increased slightly after data assimilation. We also see that at all sites Ψ_s (sathh) becomes very small (∼ 30 times smaller) after data assimilation.

Figures 8 and 9 show that h_cap and h_con change through data assimilation. However, this translates into minimal differences between the prior and post soil temperatures; both prior and post data assimilation temperature estimates are close to the in situ COSMOS-UK measurements (not shown).

3.3 Effect of data assimilation on pedotransfer function constants

In this section we present the changes to the 12 PTF constants κ₁–κ₁₂. These updates are the direct result of applying the data assimilation algorithm.

Figure 10 shows prior (orange) and posterior (blue) distributions of the 12 PTF constants, κ₁ to κ₁₂. These plots demonstrate how the dependence of the soil physics parameters on texture is changed in Eqs. (4) to (11) via data assimilation. The values of κ₁, κ₂, and κ₃ control the magnitude of the soil physics parameter b through Eq. (4). The decreases of κ₂ and κ₃ after data assimilation translate to a decreased dependence of b on clay and sand fractions through Eq. (4). Changes to κ₄, κ₅, and κ₆ contribute to changes to θ_s through Eq. (5). The large increase in κ₄ values allows larger values of θ_s to be realized after data assimilation. The parameter Ψ_s is controlled by κ₇, κ₈, and κ₉. The mean value of κ₇ is greatly reduced following data assimilation, and this leads to the much smaller posterior values of Ψ_s seen in Figs. 8 and 9. The constants κ₁₀, κ₁₁, and κ₁₂ determine the values of K_s through Eq. (7). The shift in the κ₁₀ distribution to larger values leads to the reduction in values of K_s seen in Figs. 8 and 9.

4 Discussion

The results in Sect. 3.1 show that we have been able to successfully update the constants in a Cosby-like pedotransfer function based on field-scale in situ soil moisture measurements. The new set of constants obtained in this way generate soil physics parameters at each studied COSMOS-UK site such that there is a large improvement in the match between modelled and observed field-scale soil moisture at all sites.

Our results suggest that it is primarily a combination of the changes to θ_s, Ψ_s, and K_s distributions which results in a better match to the observations after data assimilation. Calibrating the PTF using field-scale soil moisture observations allows the model to access higher soil moisture values. We suggest that the data assimilation is effectively acting to slow the drainage of water in JULES, especially close to saturation, by increasing θ_s and decreasing K_s.

Representation of soil physics processes in land surface models is fundamentally important in modelling soil moisture, and it is important to note that though the soil physics parameter values calculated here fall within a physically reasonable range, they may not exactly match physically expected values for a number of reasons. Firstly, we have generated new soil physics parameters based on field-scale COSMOS-UK measurements; differences in parameter values from the prior values may therefore reflect the different spatial scales over which they were calculated. Additionally, the COSMOS-UK soil moisture observations likely include contributions from processes which are important to soil moisture but we have not taken account of here with JULES, such as ponding of water on the soil surface, interception of water on vegetation, groundwater processes, and local soil compaction. Therefore, we may be effectively correcting for these processes (and others not included in JULES) through our new soil physics parameters. In this experiment we have mainly used grass sites to minimize the impact of vegetation in the daily averaged moisture measurements (JULES outputs show the amount of water intercepted to be, at most, of the order 100 times smaller than the amount of water in the topsoil layer).

The data assimilation is potentially correcting for deficiencies in supporting datasets (such as soil texture information or driving meteorological data) as well as parameter values or process representation in JULES, and we acknowledge that our optimization of the Cosby PTF here relies on consistent soil texture data from the HWSD. When using a land surface model, there is inevitable uncertainty in the soil texture used to generate soil physics parameters. Textures taken from any global dataset, as in this study, are likely inappropriately coarse. On the other hand, soil texture measurements taken at a point will also be unrepresentative of the scales on which land surface models are run. In order to strengthen our conclusions, we repeated our experiment using the SoilGrids soil texture database (Hengl et al., 2017). This gave similar results to the ones shown; optimizing the Cosby PTF produced a better match to the observations at all sites, and we saw a resultant increase in θ_s (vsat) and reduction in both K_s (satcon) and Ψ_s (sathh).

Despite these potential limitations, the improvements in soil moisture seen here were obtained by assimilating all the soil moisture values across 16 sites simultaneously rather than on a per-site basis. This strengthens our implicit assumption that the same physical processes can be modelled (through JULES and the Cosby pedotransfer function) for a range of different UK sites and soil types. The fact that one newly optimized PTF improves the fit to data across all 16 sites suggests that this is a systematic improvement to the PTF; i.e. we are improving the mapping between soil texture as reported in the HWSD and soil physics parameters relevant to field-scale application of JULES.

5 Conclusions

We have shown that it is possible to use the LaVEnDAR data assimilation framework to improve JULES estimates of soil moisture based on 1 year’s worth of field-scale COSMOS-UK soil moisture measurements across 16 sites. We have demonstrated improved fit to observations over a 2-year period at all 16 sites by adjusting the values of constants in the underlying pedotransfer function. Averaging across all the sites, 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 (range 0.31 to 0.89).

The method we propose here for calibrating a PTF using a data assimilation approach could be used for any different choice of land surface model, soil texture data, and/or PTF; our choice of PTF here was motivated by the fact that it is widely used and has a relatively simple mathematical formulation. Calibrating PTFs for the soils on which they are to be used and at the scales at which they are applied, rather than on small-scale field or lab soil samples, will ultimately improve the performance of land surface models. This will allow for better estimates from flood forecasting models, earth system models, and numerical weather prediction.