Effects of Soil Hydraulic Properties on Soil Moisture Estimation

Accurate soil moisture (SM) data are required to understand land surface processes, which include water, heat, and momentum exchanges between ground and atmosphere, as well as water and heat transports (Chen et al., 2018). Soil moisture is required to calculate the proportion of surface runoff that results from precipitation in hydrology (Zohaib et al., 2017; Kolassa et al., 2018; Deng et al., 2020) and is critical to effective water management in areas of agricultural irrigation (Fu et al., 2014; Huang et al., 2016). However, accurate large scale soil moisture data can be difficult to obtain. Some soil moisture observation sites have been established, such as the Walnut Gulch Experimental Watershed (WGEW) in America (Keefer et al., 2008; Yu et al., 2014) and the Heihe River Basin Field Experiment (HEIFE) in China (Chen et al., 2018), but they are spatially extremely dispersed.

Land surface models (Dickinson et al., 1993; Famiglietti and Wood, 1994; Sellers et al., 1996; Oleson et al., 2004; Fu et al., 2018; Lawrence et al., 2018) can be used to predict soil moisture at different spatio–temporal scales. However, the accuracy of simulated soil moisture is affected by the uncertainty of forcing data, model parameters, and model structure among other influences (Yu et al., 2001; Vereecken et al., 2008; Xu et al., 2018; Yang et al., 2021). In the uncertainties, the soil hydraulic properties are critical to land surface processes (Han et al., 2014; Assouline and Selker, 2017; Brandhorst et al., 2017; Chen et al., 2018; Jia et al., 2019), and they have been investigated in some studies. Based on the revised parameterization schemes in community land model (CLM), field sampling, and laboratory tests, Chen et al. (2018) found out that saturated soil hydraulic conductivity decreased and was positively correlated with organic matter content when soil depth increased in the source region of the Yellow River (SRYR). Jia et al. (2019) found that the calculation of saturated water conductivity in CLM4.5 is inaccurate by testing the soil samples and develops a new format of saturated hydraulic conductivity for soil moisture estimation according to the content of sandy soil at the Maqu area in SRYR. An optimization algorithm [e.g., genetic algorithm, GA (Yu et al., 2014)] can be used to optimize values of soil hydraulic properties in order to increase the accuracy of predictive soil moisture.

In-situ observations provide accurate soil moisture data (Keefer et al., 2008; Zheng et al., 2017; Su et al., 2020), but with spatial gaps. Remote sensing (Zhao et al., 2016; Reichle et al., 2017; Tangdamrongsub et al., 2020) provides large-scale soil moisture data with temporal gaps. These data are not satisfactory for practical application (Huang et al., 2016), but data assimilation technique can make model predictions more accurate by merging multi-source observed data into a land surface model (Han and Li, 2008; Yu et al., 2012; Dong et al., 2016; Fu et al., 2019, 2020). For example, the particle filter (PF) (Weerts and El Serafy, 2006; Dong et al., 2015; Yan et al., 2015; Ju et al., 2019), the Kalman filter (KF) and various KF-based methods (Lyu et al., 2011; Huang et al., 2016; Gevaert et al., 2018; Naz et al., 2020; Seo et al., 2021) have been used in soil moisture research for several decades.

Fu et al. (2022) pointed out that the ranges of calculated values of soil hydraulic properties are small according to the numerical solving method in CLM, which implies a weak influence of calculated soil hydraulic properties on soil moisture simulation. We therefore investigate whether the soil hydraulic properties calculated by CLM can be replaced with optimized values based on soil type or field observation to make the simulation of soil moisture more accurate. In view of the distinctiveness of the soil hydraulic properties in the land surface processes, they are firstly optimized by using GA. Then, we conduct soil moisture assimilation experiments using optimized soil hydraulic properties and the unscented weighted ensemble Kalman filter (UWEnKF) at the ELBARA field site in SRYR, in the northeast of the Tibetan Plateau, China.

The experiment is conducted at the ELBARA field site (a microwave observation system with the third generation European Space Agency (ESA) L-Band Radiometer), which is part of the Maqu soil moisture and soil temperature (SMST) monitoring network in SRYR, China (Fig. 1). The data at ELBARA are observed from 2016 to August 2019 (Su et al., 2020). The time step is 30 min for observed leaf area index (LAI) and meteorological data (e.g., precipitation, wind speed, air pressure, and air temperature) and is 15 min for observed soil moisture. Considering the need for completeness and continuity in the soil moisture observations, the observed data from 31 July 2017 (day 212) to 31 October 2017 (day 304) are used in this study. The main land cover is alpine meadow, and soil texture includes loam, sandy loam, and silt loam at different soil depths in the Maqu SMST monitoring network (Zheng et al., 2017; Zhao et al., 2018). Table 1 shows the observational ranges of basic soil properties at the Maqu monitoring network.

Figure 1.  ELBARA experimental site at the Maqu monitoring network in SRYR (Fu et al., 2022).

The following one-dimensional (1-D) Richards equation is often used to describe the vertical movement of soil water in the unsaturated zone (Richards, 1931; Chirico et al., 2014) in some land surface models (e.g., CLM):

where $ \theta $, $ K $, and $ \psi $ are volumetric soil moisture content (m³ m⁻³), soil hydraulic conductivity (mm s⁻¹), and soil matric potential (mm), respectively; $ t $, $ z $, and $ e $ are time (s), soil depth (mm), and soil moisture sink term (mm mm⁻¹ s⁻¹), respectively. The following relationship between $ \theta $, $ K $, and $ \psi $ is used to solve the above equation (Clapp and Hornberger, 1978):

where $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B are saturated volumetric soil moisture content (i.e., porosity; m³ m⁻³), saturated soil hydraulic conductivity (mm s⁻¹), saturated soil matric potential (mm), and exponent of the soil porosity distribution, respectively. They are calculated in CLM as follows:

where the length scale $ {{z}}_{\mathrm{*}} $ is 0.5 m, $ P_{{\rm s},{i}} $ and $ P_{{\rm c},{i}} $ are the sand and clay proportions in soil, and $ z_{{\rm h},{i}} $ is the interface between two adjacent soil layers i.

The FAO-56 Penman–Monteith equation (Allen, 2000) is used to calculate evapotranspiration; the schematic diagram in Fig. 2 shows the numerical scheme, which is identical to that in CLM.

Figure 2.  Schematic diagram for numerical scheme (Lawrence et al., 2018; Fu et al., 2022).

In this study, we analyze the effects of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B on soil moisture estimation, so we only use the 1-D Richards equation in the experiment. In the following sections, the model is the 1-D Richards equation.

A novel data assimilation method, the unscented weighted ensemble Kalman filter (UWEnKF), was introduced by Fu et al. (2020). The well performance of UWEnKF has been demonstrated through soil moisture assimilation (Fu et al., 2020, 2022). In UWEnKF, the ensemble members $ {{X}}_{{i}}\; ({i}=0,\cdots ,2{n}{N}) $ are symmetric about the expectation value $ \bar{X} $ and are defined as follows:

where ${q}_{{j}{j}} \;({j}{j}=1,\cdots ,{N})$ conforms to N(0, $ {P}_{x} $) and is randomly generated errors, $ n $ is the state variable dimension, N is the generated error number, and $ \lambda $ is an adjustable scaling parameter (Fu et al., 2020). UWEnKF includes two steps: forecasting and updating.

(1) Forecasting. The state variable with ensemble members ${\bar X}_{{t}|{t}-1, { }{i}}$ (soil moisture $ \theta $ in this study) at time $ t $ can be predicted based on $ {{X}}_{{t}-1,{i}} $ at time $ {t}-1 $ using the following equation:

where F is the state function, here is the 1-D Richards equation [Eq. (1)]; $ {{v}}_{{t}} $ is the state noise.

(2) Updating. Each predicted state value is updated by the following equation:

where $ {{Y}}_{{t},{i}} $ is the generated observation ensemble. Kalman gain $ {{K}}_{{t}} $ is calculated as follows:

where H is the observation function, which links states and observations; $ {{u}}_{{t}} $ is the measurement noise; I is the identity matrix, $ {W}=[{{W}}_{{t},0},\cdots ,{{W}}_{{t},2{n}{N}}] $, and $ {I}.\times {{W}}^{{T}} $ is a diagonal matrix; ${{m}}_{{i}}={\rm fix}\left[\left(2{n}{N}+1\right){{\omega }}_{{t},{i}}\right]$; Q and R are the model and the observation error covariance matrix, respectively.

The updated state variable at time t is:

The genetic algorithm (GA) is a very general robust optimization algorithm (Arifovic, 1994; Bies et al., 2006; Yu et al., 2014) based on the Darwinian evolutionary theory. GA has no limitations on derivation and function continuity; it is inherent parallel and can be globally optimized. GA can probabilistically determine and guide the optimized search space automatically, and adaptively adjusts the search without specific rules. The GA process is as follows:

(1) Initialization: generate M individuals randomly to be the initial population P(0) at time t = 0, and set the maximum genetic times (number of generations).

(2) Individual evaluation: calculate the fitness of each individual in the population P(t) at each generation.

(3) Reproduction: select each individual based on fitness and reproduce the individual to the next generation of the population; the aim of this step is to propagate the optimized individuals to the next generation.

(4) Cross over: different individuals may undergo cross over and the cross over individuals are added into the next generation.

(5) Mutation: select a variant individual and mutate it, and then propagate the generated individual into the next generation.

The population P(t + 1) can be obtained after reproduction, cross over and mutation operations have been performed on P(t).

(6) Termination: if t = T, the individual with maximum fitness obtained in the evolution process is output as the optimal solution.

More details of GA are given at http://www.obitko.com/tutorials/genetic-algorithms/index.php, and the flowchart is given in Yu et al. (2014).

The 1-D Richards equation is run from 31 July (day 212) to 13 October 2017 (day 304). We firstly analyze the effects of soil hydraulic properties or parameters ($ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B) on soil moisture estimation. They are set to values with a step equal to 1% of their interval length (maximum value minus minimum value) for each independent model run. Secondly, the soil hydraulic properties are optimized by using GA from 31 July to 31 August 2017. The assimilation experiments are then conducted by using the optimized soil hydraulic properties for 61 day beginning on 1 September 2017 (day 244).

Soil moisture observations at four soil depths (5, 20, 50, and 80 cm) in ELBARA are used, and the soil layers are redistributed. Based on the work of Fu et al. (2022), we used basic soil properties (clay, silt, sand) at depth 40 cm in place of those at depth 50 cm in this study. Table 1 shows that $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is observed only at soil depths 10, 20, 40, and 80 cm, whereas Section 3.1 shows that $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ should be measured at the interface of two adjacent soil layers in this study. Thus, we assume that the value of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ at the interface between soil depths 5 and 20 cm uses the observed value at 10 cm; the value of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ at the interface between 20- and 50-cm depths uses the observed value at 20 cm; the value of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ at the interface between 50- and 80-cm depths uses the observed value at 40 cm; and the value of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ at the interface between 80 cm and low boundary uses the observed value at 80 cm.

Note that soil moisture simulation or simulated soil moisture refers to the 1-D Richards equation [i.e., Eq. (1)], and the soil moisture assimilation or assimilated soil moisture refers to the data assimilation system [i.e., Eqs. (8)–(23)] in this study.

The root-mean-square error (RMSE) and skill score (SS) (Murphy, 1996) are used to evaluate the effects of soil hydraulic properties or parameters on soil moisture estimation and determine how well UWEnKF increases simulation accuracy:

where $ {{X}}_{\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{i},{t}} $ and $ {{X}}_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u},{t}} $ are the soil moisture assimilations and simulations at time t, respectively, $ {{X}}_{\mathrm{o}\mathrm{b}\mathrm{s},t} $ is soil moisture observations at time t, and N is the total number of time steps.

If SS is positive and SS ≤ 1, data assimilation (DA) performs well; if SS = 1, DA performs perfectly; if SS = 0, DA is ineffective.

In analyzing the effects of soil hydraulic properties or parameters on soil moisture estimation, we first use Eqs. (4)–(7) to calculate the ranges of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B based on the observed values of basic soil properties shown in Table 1, and their calculated values are shown in Table 2.

Soil moisture (SM) simulations obtained by the 1-D Richards equation [Eq. (1)] using the observed or calculated values of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (Tables 1–2) are shown in Fig. 3. In this figure, the other soil hydraulic properties $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B are calculated by using Eqs. (4) and (6)–(7) based on the mean values of observed basic soil properties shown in Table 1. Figure 3 shows that SM simulations differed as $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ varied. The purple shaded SM simulations are small because the range of calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is small. SM simulations using observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ are mostly greater than those using calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for the upper two soil layers. In the third and bottom layers, SM simulations with calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (purple area) are sometimes within the gray shaded area, and sometimes not. This suggests that using observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for the bottom two soil layers may result in better SM simulation than using calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, but not for upper two soil layers. RMSE between the simulated and observed SM with varying $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is shown in Fig. 4. The first column shows RMSE for observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and the second column shows RMSE for calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. It can be seen that RMSE increases as observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ increased for three soil layers (first, second, and bottom layers), but not for the third soil layer. This indicates that better SM simulations can be obtained with smaller observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for the three soil layers and a reasonable observed value of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for the other (third) soil layer. There is no clear change in RMSE when using calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ because the change in calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is small.

Figure 3.  Variations of soil moisture simulations and observations in different depths. The gray shaded area (Simulation range) shows the range of an open loop with observed values of $ {\mathrm{\theta }}_{\mathrm{sat}} $ varied; the purple shaded area (Simulation-cal range) shows the range of an open loop with calculated values of ${\rm{\theta }}_{\rm sat}$ varied using Eq. (5); the purple area appears to be a line, because the range of calculated $ {\mathrm{\theta }}_{\mathrm{sat}} $ is small, which leads to a small difference in soil moisture simulations for different calculated $ {\mathrm{\theta }}_{\mathrm{sat}} $; ${\rm{\theta }}_{{\rm sat}{{\text -}{{\rm{min}}}}}$ and ${\mathrm{\theta }}_{{\mathrm{sat}}{\text-}{\rm max} }$ are the soil moisture simulations for the smallest and largest observed values of $ {\mathrm{\theta }}_{\mathrm{sat}} $.

Figure 4.  RMSE between the soil moisture simulation and observation with different values of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. The left column shows RMSE with observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ varied; the right column shows RMSE with calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in Eq. (5), and RMSE seems a constant because of the small change in calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $.

Figure 5 shows soil moisture (SM) simulations using the observed or calculated values of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (Tables 1–2), and Fig. 6 shows RMSE as observed and calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ varied. The values of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B are calculated by using the equations given in Section 3. Figure 5 shows that the purple areas are close to or included in the gray areas, and that there is little change (purple area) in SM simulations using calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. The SM simulations using the smallest and largest observed values of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in Fig. 5 suggest that the change in RMSE is not linear as observed $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ increased, which can be seen in the first column of Fig. 6. RMSE in the second column of Fig. 6 decreases as calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ increased in the upper three soil layers, and RMSE initially decreases and then becomes stable in the bottom soil layer. The RMSE results imply that a reasonable value of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ among the observed intervals can produce better SM estimation.

Figure 5.  As in Fig. 3, but for $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $.

Figure 6.  RMSE between soil moisture simulation and observation for various ${{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}$. The left column shows RMSE for varied observed $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $; the right column shows RMSE for varied calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in Eq. (4).

Simulated soil moisture using calculated values of $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and B are shown in Fig. 7. Other soil hydraulic properties or parameters are calculated by using the mean values of observed basic soil properties given in Table 1. The figure shows that the gray area is small for both $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and B, which is due to the small ranges of calculated $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and B.

Figure 7.  Variations of soil moisture simulations and observations. The gray shaded area (simulation-cal range) shows the range of an open loop with various calculated values of $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in the left column of Eq. (7) and various calculated values of B in the right column of Eq. (6).

Only calculated values of $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and B are available, and their ranges are small (Table 2), which lead to only small changes in soil moisture simulations (Fig. 7). We therefore do not optimize $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ or B, and instead calculated them using the mean values of observed basic soil properties (Table 1) for the following experiments.

The variable K_sat at each interface of two adjacent soil layers is optimized by GA based on the range of each layer (Table 1) from 31 July (day 212) to 31 August 2017 (day 243). The ranges of K_sat for the top to bottom soil interfaces are respectively 1.14 × 10⁻³ to 8.53 × 10⁻³, 7.44 × 10⁻⁴ to 6.27 × 10⁻³, 1.57 × 10⁻⁴ to 7.33 × 10⁻⁴, and 2.67 × 10⁻⁴ to 2.63 × 10⁻² mm s⁻¹. We run the GA independently several times with 100 genetic times (number of generations), and values of K_sat converge to the same optimized K_sat values. Figure 8 shows one optimized process, and the optimized values of K_sat are respectively 0.00853, 0.00627, 0.000733, and 0.000267 mm s⁻¹ for the four interfaces.

Figure 8.  Optimized K_sat with different ranges depending on the values for four soil interfaces (Table 1): K_sat ∈ (1.14 × 10⁻³, 8.53 × 10⁻³) mm s⁻¹ for interface 1, K_sat ∈ (7.44 × 10⁻⁴, 6.27 × 10⁻³) mm s⁻¹ for interface 2, K_sat ∈ (1.57 × 10⁻⁴, 7.33 × 10⁻⁴) mm s⁻¹ for interface 3, and K_sat ∈ (2.67 × 10⁻⁴, 2.63 × 10⁻²) mm s⁻¹ for interface 4.

Only the optimized values of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ are used in this assimilation experiment. Values of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B are calculated using the equations in Section 3. The model is validated by using the optimized values of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ before the assimilation experiment. Figures 9a1–a4 show the results of using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ from days 244 to 304, and Figs. 9b1−b4 show the results of calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (calculated the mean values of observed basic soil properties in Table 1). It is noted that the simulated results using the optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ well reflected soil moisture changes and are closer to the observed values than those by calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. RMSE between soil moisture simulations and observations in Table 3 also indicates that the 1-D Richards equation performs well using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $.

Figure 9.  Soil moisture assimilation using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (left column; a1–a4) and calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (right column; b1–b4) with an assimilation interval 12 ∆t.

After verification, soil moisture assimilation experiments with optimized and calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ are conducted. Ensemble size is set to 16,001 based on the work of Fu et al. (2020). It is difficult to determine the explicit model and observation error statistics (Han and Li, 2008) for data assimilation in practice, so we assume that those errors are known in advance for perfect in-situ observations and refer to those in Fu et al. (2022).

Assimilated soil moisture (SM) results [obtained by Eqs. (8)–(23)] with a 12-∆t assimilation interval (here ∆t = 1 h is the simulation time step) are shown in Fig. 9. It is noted that the simulated SM values are not as close to observed values as assimilated values whether optimized or calculated (without optimization) $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is used in assimilation experiment. SM assimilation using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ produced values that are little different from observed values, in contrast to assimilation by calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. RMSE (Table 3) shows that UWEnKF performs well whether using calculated or optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. The absolute error (AE) between RMSE for simulations and for assimilations shows that UWEnKF with optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ provides more accurate SM estimation than with calculated $ { {K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. SS for assimilation by optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is significantly larger than for by calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, which suggests that UWEnKF performs better when model performance is better.

Figure 10 shows the soil moisture (SM) assimilations with daily assimilation (24 ∆t assimilation interval). The figure shows that UWEnKF increased the accuracy of SM simulations even when assimilation frequency decreased. We also note that SM assimilations using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ are closer to observations than using calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ because better SM simulations are obtained by using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $.

Figure 10.  Soil moisture assimilation using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (left column; a1–a4) and calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $; (right column; b1–b4) with 24 ∆t assimilation interval.

Figure 5 shows that the soil moisture simulations using calculated K_sat (purple area) are close to or included in those using observed K_sat (gray area), and the first column of Fig. 6 shows that RMSE does not change linearly as observed K_sat increases. We therefore used GA to optimize K_sat from the same range, 1.57 × 10⁻⁴ mm s⁻¹ to 2.63 × 10⁻² mm s⁻¹ for every soil layer with several independent runs.

Figure 11 shows the optimization of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in four independent runs, each with 1000 genetic times (number of generations). The optimized values of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (Fig. 11) for the four soil layers are all respectively 0.00853, 0.00627, 0.000733, and 0.000267 mm s⁻¹ in the four independent runs. The optimized values of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ here are the same as in Fig. 8, but the optimization processes differs. Figure 11 shows that optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ may be obtained within 100, 500, 700, 900, or even more genetic times. There is thus considerable uncertainty around the number of generations required for optimization. This is because the same range (1.57 × 10⁻⁴ to 2.63 × 10⁻² mm s⁻¹) for all soil layers is larger than the individual ranges for each layer shown in Fig. 8, which increases the randomness of initial values used in GA. We therefore need to increase the genetic times (number of generations) to ensure optimized values when the optimization range increases.

Figure 11.  Optimized ${{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}$ for the same range 1.57 × 10⁻⁴ to 2.63 × 10⁻² mm s⁻¹ at the four interfaces between two adjacent soil layers, i.e., ${{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}$ ∈ (1.57 × 10⁻⁴, 2.63 × 10⁻²) mm s⁻¹ for interfaces 1 to 4.

The soil texture at the Maqu SMST monitoring network includes loam, sandy loam, and silt loam at different soil depths. This suggests that $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ may be optimized within a certain range for other sites without the need for observed $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ values but with only knowledge of the soil texture and the range of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (for example, 1.57 × 10⁻⁴ to 2.63 × 10⁻² mm s⁻¹ at Maqu for loam soil texture).

The gray area in Fig. 3 shows that if observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is selected from the individual range of each layer, soil moisture simulations using observed $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (gray area) are less accurate than those using calculated $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ (purple area), especially in the upper two soil layers. Thus, we decide not to optimize $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ from the ranges of individual layers, but from the overall range (0.37 to 0.85 m³ m⁻³) for the entire soil depth. To ensure reasonable optimized values of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for the entire soil depth, we assume that an individual-layer optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ would decrease from top layer to bottom layer (i.e., $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t},\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}1}\geqslant {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t},\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}2}\geqslant {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t},\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}3}\geqslant $${\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t},\mathrm{l}\mathrm{a}\mathrm{y}\mathrm{e}\mathrm{r}4} $). Figure 12 shows optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ from day 212 to day 243 with 1000 genetic times (number of generations), and optimized values of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for each soil layers are, in descending order, 0.4452, 0.4452, 0.4452, and 0.4452 m³ m⁻³.

Figure 12.  Optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ with the range 0.37 to 0.85 m³ m⁻³ by GA for four soil layers.

Soil moisture (SM) assimilations using optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ with an assimilation interval 12 ∆t are shown in Fig. 13. Comparison of Fig. 13 and Figs. 9b1–b4 shows that using optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in the assimilation system increased the accuracy of SM estimations, as indicated by RMSE (Tables 3 and 4). Using optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ provides better SM estimations, but optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ is not consistent to the ranges of measured values (Table 1) for the corresponding soil layers; this may reduce the utility of the optimized results in some extent.

Figure 13.  Soil moisture assimilation using optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ with assimilation interval 12 ∆t.

The results of using both above optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in the assimilation system with 12 ∆t assimilation interval are shown in Fig. 14. RMSE and SS are shown in Table 4. Figures 9a1–a4, 13, and 14 show that the soil moisture (SM) estimations are more accurate using both optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ than when using either of them singly. UWEnKF still improves the accuracy of simulated SM according to SS, but not significantly over the results of using only optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ or optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $. This is because the 1-D Richards equation performs better than when using only optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ or optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $.

Figure 14.  Soil moisture assimilation using optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and optimized ${{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}}$ with an assimilation interval 12 ∆t.

The effects of $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B on soil moisture (SM) estimation are investigated by using the 1-D Richards equation with data from ELBARA at the Maqu SMST monitoring network in SRYR, China. Before assimilation, GA is used to optimize the soil hydraulic properties, and then, SM assimilation experiments using optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ and $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ are conducted.

(1) Based on the solving numerical scheme for the 1-D Richards equation in Section 3.1, and the preceding analysis and discussion, $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\psi }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, and B all affect the SM estimation. RMSE for SM estimation changes linearly for some soil layers when $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ varies, but this is not the case when $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ varies.

(2) According to the above numerical experiment, more accurate SM estimations may be obtained by using GA to optimize $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ based on its observation. However, the optimized $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for all layers may not be consistent with the ranges of measured $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for individual layers, which may impact its utility. For $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, the model performance improves significantly using GA for optimization according to its observation, and UWEnKF performs better with optimized $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ than with calculated $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ in SM assimilation.

(3) Based on the analysis in Sections 4.2 and 4.4, when individual ranges of observed $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ are used for each soil layer, $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ can be optimized by using GA with relatively few genetic times (e.g., 100). However, if the range of observed $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ increases, the uncertainty of genetic times required for GA to optimize $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ increases. Thus, if we use only the range of $ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $ for the entire soil depth, we need to substantially increase the genetic times required to ensure optimization, which increases computational cost.

Overall, optimizing soil hydraulic properties or parameters using GA increases the accuracy of model estimations. UWEnKF is effective and practical for increasing the accuracy of SM simulations.

In this study, we consider the effects of soil hydraulic properties on soil moisture estimation when using the solving numerical scheme for 1-D Richards equation, which adopted in CLM, not use CLM to simulate soil moisture directly. However, despite that the 1-D Richards equation is solved with four soil layers in our study, based on the four layers’ observed soil hydraulic properties ($ {{K}}_{\mathrm{s}\mathrm{a}\mathrm{t}} $, $ {\mathrm{\theta }}_{\mathrm{s}\mathrm{a}\mathrm{t}} $), the conclusions maybe have some references for some land surface models, because 1-D Richards equation has been used to simulate soil moisture in some of them. Thus, based on the above analysis, further research and field experiments are required to find suitable ranges of soil hydraulic properties or parameters for other soil textures (e.g., clay, silt clay, and sand clay), and to optimize them to increase model predictive accuracy for different soil textures over large areas. Additionally, the influence of soil moisture properties on soil moisture estimation with different solving numerical schemes for 1-D Richards equation with multi-layers would be investigated in future, and a numerical experiment in other areas with long period observations will be conducted to address the study limitation in single site.