Articles | Volume 22, issue 9
https://doi.org/10.5194/hess-22-4921-2018
https://doi.org/10.5194/hess-22-4921-2018
Research article
 | 
21 Sep 2018
Research article |  | 21 Sep 2018

Inflation method for ensemble Kalman filter in soil hydrology

Hannes H. Bauser, Daniel Berg, Ole Klein, and Kurt Roth
Abstract

The ensemble Kalman filter (EnKF) is a popular data assimilation method in soil hydrology. In this context, it is used to estimate states and parameters simultaneously. Due to unrepresented model errors and a limited ensemble size, state and parameter uncertainties can become too small during assimilation. Inflation methods are capable of increasing state uncertainties, but typically struggle with soil hydrologic applications. We propose a multiplicative inflation method specifically designed for the needs in soil hydrology. It employs a Kalman filter within the EnKF to estimate inflation factors based on the difference between measurements and mean forecast state within the EnKF. We demonstrate its capabilities on a small soil hydrologic test case. The method is capable of adjusting inflation factors to spatiotemporally varying model errors. It successfully transfers the inflation to parameters in the augmented state, which leads to an improved estimation.

Dates
1 Introduction

Data assimilation combines information from models and measurements into an optimal estimate of a geophysical field of interest (Reichle2008). It has applications in all branches of the geosciences, with weather forecasting as the driving force behind many recent advances (van Leeuwen et al.2015). The advantage of data assimilation methods (in contrast to, e.g., inverse modeling) is the possibility of considering model errors, which are characteristic of geophysical systems.

The ensemble Kalman filter (EnKF) (Evensen1994; Burgers et al.1998) is a popular data assimilation method due to its simple conceptional formulation and ease of implementation (Evensen2003). It is an extension of the Kalman filter (Kalman1960) for nonlinear models.

In hydrology, the EnKF was used for soil moisture estimation from satellite data (e.g., Reichle et al.2002; Crow and Van Loon2006) or from local measurements (e.g., De Lannoy et al.2007, 2009; Camporese et al.2009). However, the largest uncertainties in hydrology are associated with soil hydraulic material properties. They can neither be measured directly nor can they be transferred from the lab to the field, and are typically parameterized. Thus, including material properties in the estimation can be crucial in hydrology. Liu and Gupta (2007) called for an integrated assimilation framework including not only states, but also parameters and even model structure.

The joint estimation of states and parameters in data assimilation might be one possibility to reduce the influence of model errors on parameter estimation (Liu et al.2012). Such a joint estimation in the EnKF with an augmented state was already demonstrated by Anderson (2001) for an atmospheric model. In hydrology Vrugt et al. (2005) combined an EnKF and the shuffled complex evolution Metropolis algorithm, while Moradkhani et al. (2005) used a dual EnKF approach to estimate states and parameters for a rainfall–runoff model. The joint assimilation of states and parameters in an augmented state was successfully performed for example in groundwater research (e.g., Chen and Zhang2006; Hendricks Franssen and Kinzelbach2008; Kurtz et al.2012, 2014; Erdal and Cirpka2016), but also in soil hydrology for land surface models (e.g., Bateni and Entekhabi2012; Han et al.2014; Zhang et al.2017) and on smaller scales based on the Richards equation (e.g., Li and Ren2011; Wu and Margulis2011, 2013; Song et al.2014; Erdal et al.2014, 2015; Shi et al.2015; Bauser et al.2016; Brandhorst et al.2017; Botto et al.2018).

Due to unrepresented model errors and due to a limited ensemble size, the EnKF underestimates model errors, which can lead to filter inbreeding. Systematic model errors are common for example in land surface models (Vereecken et al.2015). Additionally, in soil hydrology spatially and temporally varying model errors occur due to un- or ill-represented processes like preferential flow or hysteresis. Underestimated errors cause an insufficient ensemble spread in the augmented state. This is especially severe for parameters, which are typically not changed through a forward propagation and consequently cannot increase their uncertainty again. Due to the convergent dynamics in soil hydrology, the uncertainty in the state depends strongly on the parameter spread and becomes too small as well.

Covariance inflation can counteract filter inbreeding. Different methods have been proposed: (i) additive inflation, which adds a model error after the forward propagation. This method is especially useful if prior knowledge about the model error exists. In atmospheric sciences additive inflation has been successfully applied by, e.g., using reanalysis of historical weather prediction errors (Whitaker et al.2008). (ii) Relaxation methods, which relax the analysis back to a prior perturbation or spread, have been proposed with tuning parameters (Zhang et al.2004; Whitaker and Hamill2012) or based on deviations to measurements (Ying and Zhang2015). (iii) Multiplicative covariance inflation, which inflates the complete state with a scalar factor, where the inflation factor is either chosen manually (Anderson and Anderson1999) or is estimated based on deviations from measurements (e.g., Wang and Bishop2003; Anderson2007; Li et al.2009). This method has been further extended to inflate each state component individually (Anderson2009).

All these inflation methods are developed in an atmospheric sciences context. Their transfer to soil hydrology is limited, due to the spatiotemporally varying model errors and the typically employed augmented state. For groundwater research, Kurtz et al. (2012) reported improved results by employing the inflation method by Anderson (2007), and Kurtz et al. (2014) used the constant inflation by Anderson and Anderson (1999). In soil hydrology, however, adjusted methods have been used: for example, Han et al. (2014) and Zhang et al. (2017) apply a special case of the inflation method by Whitaker and Hamill (2012) and keep the parameter spread constant to ensure a sufficient ensemble spread. Bauser et al. (2016) used the method by Anderson (2009), but adjusted the inflation of parameters.

Alternatively, no inflation method is reported (e.g., Li and Ren2011; Shi et al.2015), but instead a damping factor (Hendricks Franssen and Kinzelbach2008), which can alleviate the issue, is employed. This is done by, e.g., Wu and Margulis (2011), Song et al. (2014), Erdal et al. (2014), Brandhorst et al. (2017), and Botto et al. (2018), where Erdal et al. (2014) and Brandhorst et al. (2017) combined this method with additive inflation.

In this paper, we propose a novel multiplicative inflation method, specifically designed for the needs in the soil hydrology community. The inflation method can vary rapidly in space and time to cope with the typically varying model errors and it is capable of a transfer of the inflation in the state to the parameters in the augmented state. The remainder of this paper is organized as follows: Sect. 2 describes (i) the EnKF, (ii) our proposed inflation method and (iii) a soil hydrologic test case. Section 3 shows the results of our method applied to the test case, followed by discussion and conclusion in Sects. 4 and 5.

2 Method

2.1 Ensemble Kalman filter

The EnKF (Evensen1994; Burgers et al.1998) is the Monte Carlo extension of the Kalman filter (Kalman1960) for nonlinear models and assumes unbiased Gaussian error distributions to combine model and measurement information. The filter is a sequential method and alternates between a forecast step and an analysis step. The forecast propagates a state including its uncertainty forward in time. The analysis combines uncertain model information with uncertain measurements at this time into an optimal estimate of the state. These two steps are now explained in more detail.

The forecast propagates an ensemble of states φn forward from time k−1 to time k with a model M,

(1) φ k f , n = M φ k - 1 a , n + β n ,

where the superscripts “f” and “a” denote forecast and analysis respectively, while n denotes the ensemble members with n=1, …, N. The uncertainty in the state is directly represented through the ensemble φk-1a,n and then propagated nonlinearly with the model. Unrepresented model errors can be added through the unbiased Gaussian process noise β. This is also called additive inflation. However, the details of the model error are typically unknown and thus not represented adequately. The propagated uncertainties are directly represented through the new forecast ensemble φkf,n.

The state can be extended by, e.g., model parameters ϕ to an augmented state u=[φ, ϕ]. This requires a forecast for each augmented state component. Parameters are typically assumed to be constant in time:

(2) ϕ k f , n = ϕ k - 1 a , n .

The forecast of the state φkf,n now also depends on the corresponding parameter set ϕk-1a,n. This way, uncertainties in the parameters are propagated as well and can be reduced jointly in the analysis.

Assuming unbiased Gaussian distributions, the ensemble of augmented states is characterized through the forecast error covariance matrix Pf,

(3) P k f = u k f , n - u k f u k f , n - u k f T ,

where the overline denotes the expectation value and ukf is the ensemble mean.

The analysis combines model and measurement information based on the Gaussian error assumption. The measurement error covariance matrix R of the measurements d is defined analogously as

(4) R k = ϵ k n ϵ k n T ,

where ϵ is the measurement error. The measurements are linked to the state through the linear measurement operator H, which maps from the state space to the measurement space:

(5) d k = H k u true , k + ϵ k .

The Kalman gain K weighs the forecast error covariance matrix with the measurement error covariance matrix and maps from the measurement space back to the state space, based on the covariances in the forecast error covariance matrix:

(6) K k = P k f H k T H k P k f H k T + R k - 1 .

Based on the measurements, the Kalman gain updates the forecast ensemble to the analysis ensemble:

(7) u k a , n = u k f , n + K k d k + ϵ k n - H k u k f , n .

This update to the ensemble uka,n minimizes the analysis error covariance Pka, which fulfills

(8) P k a = I - K k H k P k f ,

for infinite ensemble sizes.

Through spurious correlations and non-Gaussian distributions, Pka will become too small, which can lead to filter inbreeding and ultimately filter divergence (e.g., Hamill et al.2001). This is intensified if the model error required in Eq. (1) is unknown.

A common way to alleviate this issue in hydrology is the use of a damping factor γ[0,1] (Hendricks Franssen and Kinzelbach2008), which is multiplied to the correction vector in Eq. (7) and consequently lessens the uncertainty reduction. The damping factor can be extended to a vector γ (and an entrywise multiplication) to treat augmented state components differently (Wu and Margulis2011). Typically, parameters are multiplied by a smaller factor than the state. However, the damping factor can only alleviate and not completely prevent the inbreeding problem.

2.2 Multiplicative inflation for soil hydrology

Multiplicative inflation is another heuristic way to avoid filter inbreeding. Anderson and Anderson (1999) proposed to increase the distance of each ensemble member to the ensemble mean by multiplying this distance by λ for the inflation factor λ≥1:

(9) u inf f , n = λ u f , n - u f + u f .

This inflation factor is applied to the complete augmented state and has to be adjusted to the specific problem. By construction, it does not alter the mean value: uinff=uf. A temporally varying inflation factor can be estimated by comparing uncertainties with the distance of measurement and forecast (e.g., Wang and Bishop2003; Anderson2007; Li et al.2009). A spatiotemporally adaptive inflation has been achieved by estimating a vector λ for the complete augmented state (Anderson2009). The author uses the correlation between measurements and augmented state dimensions and asks this question: How much inflation is required in each dimension to explain the observed differences to the measurements? Anderson (2009) showed that this works excellently for the actual state. However, we experienced possible over-inflation in parameters (which do not have any dynamics to compensate for this), which can lead to filter collapses.

We propose a more conservative inflation method and ask this question: How much of the required change in the inflation are we allowed to transfer to the state dimensions based on the correlation information? This can be achieved by applying a Kalman filter for the inflation within the EnKF.

In this Kalman filter, the inflation vector is treated as the state variable. As for parameters, we choose a constant model for the forecast in time:

(10) λ k f = λ k - 1 a .

For convenience we will drop the time subscript k in the following. Furthermore, we will use the same symbols as for the EnKF, but denote them with the subscript λ. We approximate the forecast error covariance matrix for lambda, Pλf, based on the covariance matrix of the augmented state in the EnKF, Pf, as the normalized absolute correlation matrix of the augmented state ensemble. The matrix component ij is determined as

(11) P λ f i j = σ λ 2 P f i j P f i i P f j j - 1 2 ,

where σλ denotes the uncertainty of the inflation factors. It is a tuning parameter that is kept constant over time and is assigned to all state dimensions. It influences how fast the inflation factors are adjusted. This follows the idea by Anderson (2007, 2009) to avoid a closure problem, where the inflation estimation would require its own inflation. Instead, the uncertainty is kept constant. Furthermore, only the absolute value of the correlation is considered, since the inflation is based on differences between measurement and model, but ignores their direction. Note, that this presumes that the correlations of the model state can be transferred to the inflation. In the presence of unknown model errors this assumption may or may not be correct. However, the estimation at measurement locations will remain meaningful in any case.

For the analysis, the distance dλ between mean forecast and measurement is used as measurement for λ:

(12) d λ = d - H u inf f .

The measurement error covariance matrix Rλ of dλ can be calculated based on the error covariance matrices of d and Huinff,

(13) R λ i j = ( R ) i j + HP inf f H T i j ,

where the inflated forecast error covariance matrix Pinff can be calculated from the inflation vector and the forecast error covariance matrix by combining Eq. (9) (with vector λ and entrywise multiplication) and Eq. (3): Pinff=Pf[λfλfT]. The entrywise product is denoted by  and the entrywise square root of λ by λ.

The expected distance between measurement and mean forecast based on the current inflation is

(14) h λ λ f i = R λ i i 1 2 ,

which combines the uncertainties of d and Huinff. To be able to determine the Kalman gain, we first calculate the Jacobian matrix Hλ of partial derivatives of hλ with respect to λ:

Hλij=λfjhλλfi,=2λfj12hλλfi-1(15)m(H)ij(H)imPfjmλfm12.

With this approximated measurement operator Hλ, the Kalman gain Kλ and the analysis state λa are obtained as

(16)Kλ=PλfHλTHλPλfHλT+Rλ-1,(17)λa=λf+Kλdλ-hλλf.

Note, that the matrices Pλf and Rλ can possibly become indefinite, due to the absolute value in Eqs. (11) and (13). Consequently, the inverse in Eq. (16) could become unfeasible. However, we never encountered such a case. In a situation with uncorrelated measurements, the issue can be resolved by reducing σλ just for that single time step.

With this Kalman filter, the inflation vector is updated at each time step based on the difference of the mean forecast to the measurements. Following Anderson (2007), we additionally prohibit a deflation by constraining the inflation values to (λ)i≥1.

2.3 Model

We test the proposed inflation method on a small hydrologic test case. We constructed it specifically to require a strong inflation. This makes it possible to explore features of the inflation in detail on a rather short timescale. Due to a small ensemble size, the results vary depending on the seed of the random numbers. This however, is related to different performance of the EnKF itself. In simulations (results are not shown), we found that the behavior of the inflation remains consistent. We have also tested the inflation method with real-world data by reanalyzing the application by Bauser et al. (2016) with the main result shown in Appendix B.

The Richards equation describes the change of volumetric soil water content θ (–) in a continuous porous medium,

(18) θ t - K ( θ ) h m ( θ ) - 1 = 0 ,

where K (L T−1) is the isotropic conductivity and hm (L) is the matric head. Both are related to the water content. This relation is typically described through parameterized material properties. We choose the Mualem–van Genuchten parameterization (Mualem1976; van Genuchten1980),

(19)K(Θ)=K0Θτ1-1-Θn/[n-1]1-1/n2,(20)hm(Θ)=1αΘ-n/[n-1]-11/n,

with the saturation Θ (–),

(21) Θ := θ - θ r θ s - θ r .

The parameterization is described by a set of six parameters: θs (–), θr (–), α (L−1), n (–), K0 (L T−1), and τ (–).

We additionally consider small-scale heterogeneity through Miller scaling. It assumes geometrical similarity. With this the microscopic geometry of the pore space at a macroscopic position is parameterized by a single length scale ξ and the macroscopic heterogeneity field can be generated with a single scalar field of this length scale. Miller and Miller (1956) showed that the hydraulic functions scale with this parameter according to

(22)K(θ)=K*(θ)ξ2,(23)hm(θ)=hm*(θ)1ξ,

where the functions K(θ) and hm(θ) are defined at a reference point * with Miller scaling parameter ξ=1 and from there are projected to all locations.

For the test, we choose a one-dimensional case with a depth of 50 cm for a time of 6 days. We set a groundwater table as the lower boundary condition throughout the whole time and start from equilibrium conditions. The upper boundary condition is no flux, except for a rain event with 2.0×10-7 (m s−1) during the fourth day. As observations we choose two water content measurements at depths of 9.5 and 19.5 cm as they would be available from time domain reflectometry (TDR). We set the measurement uncertainty to a standard deviation of 0.007 (e.g., Jaumann and Roth2017).

As material we choose sandy loam from Carsel and Parrish (1988): θs=0.41, θr=0.065, α=-7.5 m−1, n=1.89, K0=1.23×10-5 m s−1, and τ=0.5. For the Miller scaling we choose ξ1=0.32 at the upper measurement position and ξ2=3.2 at the lower measurement position. We reduce the description of the heterogeneity to these two parameters. The full function of the scaling factors is calculated by linearly interpolating between the measurement positions and constantly extrapolating to the boundaries.

The forward simulations are performed using MuPhi (Ippisch et al.2006) with a spatial resolution of 1 cm. This corresponds to a state with 50 dimensions.

To test the inflation method, we perform a perfect model experiment. With the EnKF we estimate the water content state and four parameters (ξ1, ξ2, K0, and τ) through the augmented state u=[θ, log10(ξ1), log10(ξ2), log10(K0),τ]. We choose to include the logarithm of ξ1, ξ2, and K0, because we expect a more linear relation to the water content state than for the actual parameters. For the water content state, we use the correct initial condition as the mean with an uncertainty of 0.005. The uncertainty is spatially correlated using the fifth-order piecewise rational function by Gaspari and Cohn (1999) with the length scale c=5 cm. As an initial guess for the parameters, we start with unknown heterogeneity log10(ξ1)=log10(ξ2)=0.0±0.25, corresponding to 2 SD (standard deviations) away from the true values of log10(ξ1)=-0.5 and log10(ξ2)=0.5. For the saturated hydraulic conductivity, we choose a too small value of log10(K0)=-5.5±0.5, K0 in (m s−1), which is about 1 SD away from the true value of log10(K0)=-4.9. For the tortuosity τ=0.5±0.5 we start from the true value.

Through the unrepresented heterogeneity, we can mimic a model error, leading to a bias towards smaller values for the estimation of K0 during times without dynamics, which may necessitate inflation. The parameter τ is expected to have a smaller influence, since the uncertainty is chosen small and it is already at the true value. This way it can act as an indicator parameter for the inflation as it does not require inflation.

The EnKF is set up with a total of 25 ensemble members and a damping vector of γ=[1.0, 0.3, 0.3, 0.3, 0.3], which we also apply to the inflation. The damping factor of 0.3 is applied to the parameters to alleviate issues of nonlinear relations between observations and parameters. For the uncertainty of the inflation factors we choose σλ=1.0.

3 Results

We estimate the water content state together with the four parameters ξ1, ξ2, K0 and τ with the EnKF as described in Sect. 2.3. The development of the water content at the two measurement locations at a depth of 9.5 and 19.5 cm is shown together with the inflation factor at these locations in Fig. 1. The inflation factor is applied to the forecast ensemble before the analysis. The standard deviation of the inflated ensemble should describe the distance of the estimated mean to the synthetic truth. Note, that the inflation factor is not based on this distance and relies on the noisy measurements. Therefore, it is only an indicator.

https://www.hydrol-earth-syst-sci.net/22/4921/2018/hess-22-4921-2018-f01

Figure 1Water content estimation at the two measurement locations. The standard deviation of the inflated ensemble should be able to explain the differences between the inflated mean and the synthetic truth. The inflation factor is increased when the ensemble uncertainty is too small.

Download

During the first 3 days without any dynamics, the uncertainty for the upper measurement is slightly underestimated, while the uncertainty in the lower measurement is slightly overestimated. This leads to an inflation factor of basically 1 for the lower measurement (factors smaller than 1 are not allowed), while the inflation factor for the upper measurement is larger. However, due to correlations between the measurement locations a stronger inflation to fully explain the difference to the truth is prevented.

The deviation from the synthetic truth is induced through the initial guess of no heterogeneity and can also be seen in the systematic deviation of the inflated mean (which is equal to the forecast mean) from the analysis mean. When the infiltration front reaches the measurements, the deviations from the truth, underestimation of the uncertainty, and inflation factors increase rapidly. All of them are more pronounced for the upper measurement location. After the main peak, the differences and also the inflation factors decrease rapidly again.

The inflation factor for the state is shown in Fig. 2. It shows the strong increase in the inflation factor during the infiltration and its fast decrease afterwards. The inflation is strongest at the measurement location at a depth of 9.5 cm. The inflation factor is transferred to the other state locations through the correlations, which decrease with distance. Directly below the measurement locations the inflation factors are increased less than above. This is due to the chosen interpolation of the Miller scaling factors. Through the interpolation between the measurement locations and extrapolation to the boundaries, the dynamics changes at the measurement locations. During the infiltration the dynamics is mainly influenced by the water content above and the correlations with these locations are stronger.

The development of the Miller scaling factors ξ1 and ξ2 at the two measurement positions (9.5 and 19.5 cm depth) is shown in Fig. 3a and b together with the estimated inflation factor for these parameters. Both initial conditions assume no heterogeneity and start at log10(ξ1)=log10(ξ2)=0.0±0.25, corresponding to 2 SD away from the true value. At the upper location the true value of log10(ξ1)=-0.5 corresponds to a finer material. Consequently, the water content drops, as seen in Fig. 1, leading to a strong correlation with the scaling factor, and log10(ξ1) is adjusted rapidly to lower values. Accordingly, the inflation factor is increased quickly in the beginning and then reduced back to 1 when the estimation of log10(ξ1) reaches and eventually underestimates the true value. The underestimation of the scaling factor corresponds to a too fine material, which leads to slower changes in the water content state and therefore smaller correlations. The scaling factor is corrected during the rain event on the fourth day, which also leads to an inflation.

https://www.hydrol-earth-syst-sci.net/22/4921/2018/hess-22-4921-2018-f02

Figure 2Inflation factor for the water content state. The inflation is strongest at the upper measurement location during the infiltration, when the uncertainty is underestimated the most. The inflation factor is transferred to the other measurement locations through the correlations in the Kalman gain. The used interpolation of Miller scaling factors impacts these correlations and leads to the smaller inflation directly below the measurement locations.

Download

The initial guess for the scaling factor for the depth of 19.5 cm underestimates the scaling factor, which corresponds to a too fine material. Again, the correlations are small. The value increases slowly during the dry period in the beginning, but is inflated and adjusted strongly during the rain event.

The saturated hydraulic conductivity K0 (Fig. 3c) was chosen to start a little more than 1 SD below the true value. Due to the unrepresented heterogeneity in the beginning, the value decreases even further. The inflation remains very small due to correlations with both measurement locations. However, as soon as the infiltration event reaches the first measurement location, the value is corrected towards the true value. At the same time, the inflation factor is increased due to the too small uncertainty. After the rain event the inflation factor drops rapidly back to 1. The hydraulic conductivity remains below the true value. Another rain event would be required to improve the estimation further.

The tortuosity τ (Fig. 3d) also influences the hydraulic conductivity function, but has in this case a much smaller impact and consequently smaller correlations with the measurements than K0. We use it as an indicator parameter and start at the true value. During the infiltration event the value is changed due to its correlation. The corresponding inflation factor is increased as well, but remains small enough and drops back to 1 quickly enough to not cause any over-inflation of the parameter.

https://www.hydrol-earth-syst-sci.net/22/4921/2018/hess-22-4921-2018-f03

Figure 3Development of Miller scaling factors ξ1 and ξ2, saturated hydraulic conductivity K0, and tortuosity τ together with their corresponding inflation factors during estimation with the EnKF.

Download

To emphasize the need for a fast-adapting inflation factor, we reduce the uncertainty of the inflation factors to σλ=0.5 to slow down their adjustment. The results are summarized in Fig. 4. The inflation of the water content state (Fig. 4a) shows that the inflation factor does not reach as high values as before (see Fig. 2). To compensate for this, the inflation acts over a longer period of time. The same effect is also observed in the inflation of the parameters (Fig. 4b and c). This leads to a smaller inflation during the rain event and consequently a too small uncertainty. At later times, when the cause of the error is not active any more, the correlations with measurement locations are reduced, leading to a slower reduction of the inflation in the parameters. In the indicator parameter τ the beginning of an over-inflation can be seen towards later times. This necessitates a more rapid inflation when correlations are used to update inflation information.

The results for the parameters K0 and τ of a run without inflation (and only damping) are shown in Fig. 5. Again, K0 moves further away from the true value due to the unrepresented heterogeneity and comes closer to the true value during the infiltration event. However, since the Miller scaling factor is not inflated in the beginning, it is adjusted slower. Consequently, the K0 is corrected longer in the wrong direction. The uncertainty eventually becomes too small and in the end the mean is more than 5 SD away from the true value, since the uncertainty cannot be increased any more.

https://www.hydrol-earth-syst-sci.net/22/4921/2018/hess-22-4921-2018-f04

Figure 4Development of the inflation factor for the water content state and saturated hydraulic conductivity K0 and tortuosity τ together with their corresponding inflation factors for an estimation with a reduced inflation factor uncertainty of σλ=0.5.

Download

https://www.hydrol-earth-syst-sci.net/22/4921/2018/hess-22-4921-2018-f05

Figure 5Development of saturated hydraulic conductivity K0 and tortuosity τ for an estimation without inflation.

Download

4 Discussion

The proposed inflation method uses a Kalman filter to estimate inflation factors within the EnKF. It is based on the difference between measurements and mean forecast state. It transfers correlations from the forecast of the augmented state to the inflation. Consequently, the performance will be limited if model errors are structurally not represented in the forecast error covariance matrix. The estimation of the inflation factors with a Kalman filter is, like the EnKF itself, based on a linearized analysis. The use of a damping factor can alleviate issues with estimating nonlinear dependent parameters. To keep the inflation consistent with the analysis in the EnKF, we apply the same damping factor for both.

We designed a small synthetic hydrologic test case for the inflation. This test case mimics a model error through initially unrepresented heterogeneity. We designed the test case so that a strong temporally varying inflation is necessary, as it can occur with real data. We choose a short time so that the details of the behavior of the method can be explored. The method showed that it is capable of inflating states and parameters. The inflation is adjusted quickly and differentiates between parameters with strong and not so strong correlations. No over-inflation of weakly correlated parameters occurred. In this specific test case the estimation with inflation is far superior to an estimation without inflation.

The fast adjustment speed of the inflation factor is important because of the fast-changing model errors and correlations with parameters. The adjustment speed is determined by the uncertainty of the inflation factor. This uncertainty is set to a constant value and has to be adjusted. For all our cases a value of σλ=1 was sufficient, but larger values were possible too. The need for such a fast adjustment is shown by estimating the same case with a reduced uncertainty of σλ=0.5, which leads to a slower adaptation of the inflation factor. This leads to smaller inflation factors, which is compensated for by maintaining them for a longer period of time. In this test case this leads to inflation at times after the infiltration front has passed the measurements already and the model error is small again. This can cause over-inflation of weakly correlated parameters. Too large uncertainties of the inflation (in our test case σλ=4), where the uncertainty is larger than the typical range for the values of lambda, can also lead to over-inflation of weakly correlated parameters. Reasons for this can be the linearizations in the analysis and the calculation of the Jacobian (Eq. 15). This limits the adjustment speed of the inflation.

Fast-dropping correlations between measurements and parameters are a limit for the method. An example could be a parameter only acting on an infiltration boundary condition. After the infiltration is over, correlations with this parameter would drop to zero and the inflation factor for this parameter will not be changed any more. If the inflation factor is not equal to 1 at this time, the parameter spread will keep increasing. In such a case, when there is no correlation, the parameter should be excluded from the estimation and consequently also from the inflation.

The method is in principle capable of compensating unrepresented model errors. However, it relies on correlations calculated from the forecast ensemble of the augmented state. If parameters have correlations with measurement locations with underestimated forecast uncertainties, the inflation will keep increasing the parameter spread until the forecast uncertainties are increased sufficiently. Therefore the correlations have to contain useful information. This means that inflating the parameters based on their correlations with measurement locations has to increase the forecast spread at these measurement locations. If the parameters have an insufficient influence on the state uncertainty, an over-inflation of the parameters can occur. An example are measurements with underestimated measurement uncertainties and short times between measurements compared to the timescale of the dynamics. Then the parameters are not able to increase the state uncertainty in the short forecast time between measurements and the forecast dynamics is not able represent the measurement noise. If such errors occur intermittently, e.g., the closed-eye period as proposed by Bauser et al. (2016) could be used to bridge these times. A rather heuristic solution could be a decay of the inflation factor towards values of 1, as already proposed by Anderson (2009).

5 Conclusions

In this work we propose a novel spatiotemporally adaptive inflation method, specifically designed for soil hydrology, which nevertheless is expected to work in similar systems as well. The inflation method is based on a Kalman filter acting within the EnKF. The method is capable of rapid adjustments of inflation factors, treating each augmented state dimension individually. This rapid adjustment is required due to temporally varying model errors, as they can appear through violation of the local equilibrium assumption of the Richards equation, hysteresis, or unrepresented heterogeneity.

We demonstrate the use of our inflation method in combination with a damping factor on a small hydrologic example. We choose heterogeneity as a possible model error, but allow the heterogeneity to be estimated along with the soil hydrologic parameters K0 and τ of the Mualem–van Genuchten parameterization. Our proposed inflation method proved to be stable in combination with parameter estimation. The performance of the estimation improved and parameter uncertainty remained consistent. The method requires that the correlations from in the forecast ensemble contain useful information for the inflation. However, we demonstrate that it even works for only weakly correlated parameters. We expect the inflation method to generally improve data assimilation with the EnKF and to thus lead to better state and parameter estimations in soil hydrology.

Data availability

The synthetic data are available upon request from the corresponding author.

Appendix A: Jacobian in the inflation method

We briefly show the derivation of the Jacobian matrix Hλ for the inflation (Eq. 15). Again, the entrywise product is denoted by  and the entrywise square root of λ by λ:

Hλij=λfjhλ(λf)i=λfjRii+HPfλfλfTHTii12,withPfsymmetric=λfjRii+mk(H)im(H)ikPfkmλfm12λfk1212=2hλλfi-1mk(H)im(H)ikPfkmλfjλfm12λfk12=2hλλfi-1mk(H)im(H)ikPfkm12δmjλfk12λfm12+δkjλfm12λfk12=2λfj12hλλfi-1m(H)ij(H)im(A1)Pfjmλfm12.
Appendix B: Real-world application

We also applied the inflation method to reanalyze the case presented earlier by Bauser et al. (2016), where measurements from 11 TDR probes were assimilated with an EnKF. There, the inflation method confirmed the behavior observed in the small synthetic case presented in this paper. For the details of the real-world case as well as the concept of the closed-eye period please refer to Bauser et al. (2016) or Bauser (2018, chap. 5), the latter of which also includes the reanalysis of the case.

In this paper, we only show the inflation related to the closed-eye period (Fig. B1), which presents the major challenge to the inflation in that particular application. During this time, preferential flow occurs and the underlying local equilibrium assumption of the Richards equation is violated. With a standard approach, parameters become biased to compensate these errors. To avoid this, Bauser et al. (2016) introduced the closed-eye period, which pauses the parameter estimation and only guides the water content states through times, when assumptions are violated. Compared to the standard approach, this leads to a reduced bias in the parameters, but effectively increases the model errors during the closed-eye period. A strong inflation is required to compensate this error. The inflation method used in Bauser et al. (2016) was just able to accomplish this and the authors were concerned that a too slow adjustment speed of the inflation limits the applicability of the closed-eye period for cases with larger model errors.

Figure B1 confirms the fast adjustment speed of the new inflation method proposed in this paper for the real-world application. The strong required inflation stays well within the closed-eye period. This enables the EnKF to pick up the parameter estimation after the period from a water content state consistent with the TDR measurements and facilitates the use of the closed-eye period.

https://www.hydrol-earth-syst-sci.net/22/4921/2018/hess-22-4921-2018-f06

Figure B1Inflation factor for the water content state in a real-world application for a standard and closed-eye EnKF together with atmospheric forcing and the corresponding response of measured water contents. The closed-eye period starts when the infiltration front reaches the topmost TDR. During this time the local equilibrium assumption of the Richards equation is violated and a strong inflation is required. The new inflation method allows a fast adjustment of the inflation factors, which enables the EnKF to effectively guide the water content states with the TDR measurements through the closed-eye period. Modified from Bauser (2018).

Download

Author contributions

HHB designed, implemented, performed and analyzed the presented study. DB provided computational software. All authors participated in continuous discussions. HHB prepared the manuscript with contributions from all authors.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank editor Insa Neuweiler and two anonymous reviewers for their comments, which helped to improve this paper.

This research is funded by the Deutsche Forschungsgemeinschaft (DFG) through project RO 1080/12-1 and the Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg (Az 33-7533.-30-20/6/2). HGS MathComp provided travel expenses for Hannes H. Bauser and Daniel Berg.

Edited by: Insa Neuweiler
Reviewed by: two anonymous referees

References

Anderson, J. L.: An ensemble adjustment Kalman filter for data assimilation, Mon. Weather Rev., 129, 2884–2903, https://doi.org/10.1175/1520-0493(2001)129<2884:AEAKFF>2.0.CO;2, 2001. a

Anderson, J. L.: An adaptive covariance inflation error correction algorithm for ensemble filters, Tellus A, 59, 210–224, https://doi.org/10.1111/j.1600-0870.2006.00216.x, 2007. a, b, c, d, e

Anderson, J. L.: Spatially and temporally varying adaptive covariance inflation for ensemble filters, Tellus A, 61, 72–83, https://doi.org/10.1111/j.1600-0870.2008.00361.x, 2009. a, b, c, d, e, f

Anderson, J. L. and Anderson, S. L.: A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts, Mon. Weather Rev., 127, 2741–2758, https://doi.org/10.1175/1520-0493(1999)127<2741:AMCIOT>2.0.CO;2, 1999. a, b, c

Bateni, S. M. and Entekhabi, D.: Surface heat flux estimation with the ensemble Kalman smoother: Joint estimation of state and parameters, Water Resour. Res., 48, w08521, https://doi.org/10.1029/2011WR011542, 2012. a

Bauser, H. H.: Knowledge Fusion in Soil Hydrology, PhD thesis, Ruperto-Carola University Heidelberg, Heidelberg, Germany, https://doi.org/10.11588/heidok.00024713, 2018. a, b

Bauser, H. H., Jaumann, S., Berg, D., and Roth, K.: EnKF with closed-eye period – towards a consistent aggregation of information in soil hydrology, Hydrol. Earth Syst. Sci., 20, 4999–5014, https://doi.org/10.5194/hess-20-4999-2016, 2016. a, b, c, d, e, f, g, h

Botto, A., Belluco, E., and Camporese, M.: Multi-source data assimilation for physically based hydrological modeling of an experimental hillslope, Hydrol. Earth Syst. Sci., 22, 4251–4266, https://doi.org/10.5194/hess-22-4251-2018, 2018. a, b

Brandhorst, N., Erdal, D., and Neuweiler, I.: Soil moisture prediction with the ensemble Kalman filter: Handling uncertainty of soil hydraulic parameters, Adv. Water Resour., 110, 360–370, https://doi.org/10.1016/j.advwatres.2017.10.022, 2017. a, b, c

Burgers, G., van Leeuwen, P. J., and Evensen, G.: Analysis scheme in the ensemble Kalman filter, Mon. Weather Rev., 126, 1719–1724, https://doi.org/10.1175/1520-0493(1998)126<1719:ASITEK>2.0.CO;2, 1998. a, b

Camporese, M., Paniconi, C., Putti, M., and Salandin, P.: Ensemble Kalman filter data assimilation for a process-based catchment scale model of surface and subsurface flow, Water Resour. Res., 45, W10421, https://doi.org/10.1029/2008WR007031, 2009. a

Carsel, R. F. and Parrish, R. S.: Developing joint probability distributions of soil water retention characteristics, Water Resour. Res., 24, 755–769, https://doi.org/10.1029/WR024i005p00755, 1988. a

Chen, Y. and Zhang, D.: Data assimilation for transient flow in geologic formations via ensemble Kalman filter, Adv. Water Resour., 29, 1107–1122, https://doi.org/10.1016/j.advwatres.2005.09.007, 2006. a

Crow, W. T. and Van Loon, E.: Impact of incorrect model error assumptions on the sequential assimilation of remotely sensed surface soil moisture, J. Hydrometeorol., 7, 421–432, https://doi.org/10.1175/JHM499.1, 2006. a

De Lannoy, G. J. M., Houser, P. R., Pauwels, V. R. N., and Verhoest, N. E. C.: State and bias estimation for soil moisture profiles by an ensemble Kalman filter: Effect of assimilation depth and frequency, Water Resour. Res., 43, w06401, https://doi.org/10.1029/2006WR005100, 2007. a

De Lannoy, G. J. M., Houser, P. R., Verhoest, N. E. C., and Pauwels, V. R. N.: Adaptive soil moisture profile filtering for horizontal information propagation in the independent column-based CLM2.0, J. Hydrometeorol., 10, 766–779, https://doi.org/10.1175/2008JHM1037.1, 2009. a

Erdal, D. and Cirpka, O. A.: Joint inference of groundwater-recharge and hydraulic-conductivity fields from head data using the ensemble Kalman filter, Hydrol. Earth Syst. Sci., 20, 555–569, https://doi.org/10.5194/hess-20-555-2016, 2016. a

Erdal, D., Neuweiler, I., and Wollschläger, U.: Using a bias aware EnKF to account for unresolved structure in an unsaturated zone model, Water Resour. Res., 50, 132–147, https://doi.org/10.1002/2012WR013443, 2014. a, b, c

Erdal, D., Rahman, M., and Neuweiler, I.: The importance of state transformations when using the ensemble Kalman filter for unsaturated flow modeling: Dealing with strong nonlinearities, Adv. Water Resour., 86, 354–365, https://doi.org/10.1016/j.advwatres.2015.09.008, 2015. a

Evensen, G.: Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics, J. Geophys. Res.-Oceans, 99, 10143–10162, https://doi.org/10.1029/94JC00572, 1994. a, b

Evensen, G.: The ensemble Kalman filter: theoretical formulation and practical implementation, Ocean Dynam., 53, 343–367, https://doi.org/10.1007/s10236-003-0036-9, 2003. a

Gaspari, G. and Cohn, S. E.: Construction of correlation functions in two and three dimensions, Q. J. Roy. Meteorol. Soc., 125, 723–757, https://doi.org/10.1002/qj.49712555417, 1999. a

Hamill, T. M., Whitaker, J. S., and Snyder, C.: Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter, Mon. Weather Rev., 129, 2776–2790, https://doi.org/10.1175/1520-0493(2001)129<2776:DDFOBE>2.0.CO;2, 2001. a

Han, X., Hendricks Franssen, H.-J., Montzka, C., and Vereecken, H.: Soil moisture and soil properties estimation in the Community Land Model with synthetic brightness temperature observations, Water Resour. Res., 50, 6081–6105, https://doi.org/10.1002/2013WR014586, 2014. a, b

Hendricks Franssen, H. J. and Kinzelbach, W.: Real-time groundwater flow modeling with the ensemble Kalman filter: Joint estimation of states and parameters and the filter inbreeding problem, Water Resour. Res., 44, w09408, https://doi.org/10.1029/2007WR006505, 2008. a, b, c

Ippisch, O., Vogel, H.-J., and Bastian, P.: Validity limits for the van Genuchten–Mualem model and implications for parameter estimation and numerical simulation, Adv. Water Resour., 29, 1780–1789, https://doi.org/10.1016/j.advwatres.2005.12.011, 2006. a

Jaumann, S. and Roth, K.: Effect of unrepresented model errors on estimated soil hydraulic material properties, Hydrol. Earth Syst. Sci., 21, 4301–4322, https://doi.org/10.5194/hess-21-4301-2017, 2017. a

Kalman, R. E.: A new approach to linear filtering and prediction problems, J. Basic Eng., 82, 35–45, 1960. a, b

Kurtz, W., Hendricks Franssen, H.-J., and Vereecken, H.: Identification of time-variant river bed properties with the ensemble Kalman filter, Water Resour. Res., 48, w10534, https://doi.org/10.1029/2011WR011743, 2012. a, b

Kurtz, W., Hendricks Franssen, H.-J., Kaiser, H.-P., and Vereecken, H.: Joint assimilation of piezometric heads and groundwater temperatures for improved modeling of river-aquifer interactions, Water Resour. Res., 50, 1665–1688, https://doi.org/10.1002/2013WR014823, 2014. a, b

Li, C. and Ren, L.: Estimation of unsaturated soil hydraulic parameters using the ensemble Kalman filter, Vadose Zone J., 10, 1205–1227, https://doi.org/10.2136/vzj2010.0159, 2011. a, b

Li, H., Kalnay, E., and Miyoshi, T.: Simultaneous estimation of covariance inflation and observation errors within an ensemble Kalman filter, Q. J. Roy. Meteorol. Soc., 135, 523–533, https://doi.org/10.1002/qj.371, 2009. a, b

Liu, Y. and Gupta, H. V.: Uncertainty in hydrologic modeling: Toward an integrated data assimilation framework, Water Resour. Res., 43, w07401, https://doi.org/10.1029/2006WR005756, 2007. a

Liu, Y., Weerts, A. H., Clark, M., Hendricks Franssen, H.-J., Kumar, S., Moradkhani, H., Seo, D.-J., Schwanenberg, D., Smith, P., van Dijk, A. I. J. M., van Velzen, N., He, M., Lee, H., Noh, S. J., Rakovec, O., and Restrepo, P.: Advancing data assimilation in operational hydrologic forecasting: progresses, challenges, and emerging opportunities, Hydrol. Earth Syst. Sci., 16, 3863–3887, https://doi.org/10.5194/hess-16-3863-2012, 2012. a

Miller, E. E. and Miller, R. D.: Physical theory for capillary flow phenomena, J. Appl. Phys., 27, 324–332, https://doi.org/10.1063/1.1722370, 1956. a

Moradkhani, H., Sorooshian, S., Gupta, H. V., and Houser, P. R.: Dual state–parameter estimation of hydrological models using ensemble Kalman filter, Adv. Water Resour., 28, 135–147, https://doi.org/10.1016/j.advwatres.2004.09.002, 2005. a

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522, https://doi.org/10.1029/WR012i003p00513, 1976. a

Reichle, R. H.: Data assimilation methods in the Earth sciences, Adv. Water Resour., 31, 1411–1418, https://doi.org/10.1016/j.advwatres.2008.01.001, 2008. a

Reichle, R. H., McLaughlin, D. B., and Entekhabi, D.: Hydrologic data assimilation with the ensemble Kalman filter, Mon. Weather Rev., 130, 103–114, https://doi.org/10.1175/1520-0493(2002)130<0103:HDAWTE>2.0.CO;2, 2002. a

Shi, L., Song, X., Tong, J., Zhu, Y., and Zhang, Q.: Impacts of different types of measurements on estimating unsaturated flow parameters, J. Hydrol., 524, 549–561, https://doi.org/10.1016/j.jhydrol.2015.01.078, 2015. a, b

Song, X., Shi, L., Ye, M., Yang, J., and Navon, I. M.: Numerical comparison of iterative ensemble Kalman filters for unsaturated flow inverse modeling, Vadose Zone J., 13, https://doi.org/10.2136/vzj2013.05.0083, 2014. a, b

van Genuchten, M. T.: A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892–898, https://doi.org/10.2136/sssaj1980.03615995004400050002x, 1980. a

van Leeuwen, P. J., Cheng, Y., and Reich, S.: Nonlinear data assimilation, in: vol. 2, Springer International Publishing, Switzerland, https://doi.org/10.1007/978-3-319-18347-3, 2015. a

Vereecken, H., Huisman, J. A., Hendricks Franssen, H. J., Brüggemann, N., Bogena, H. R., Kollet, S., Javaux, M., van der Kruk, J., and Vanderborght, J.: Soil hydrology: Recent methodological advances, challenges, and perspectives, Water Resour. Res., 51, 2616–2633, https://doi.org/10.1002/2014WR016852, 2015. a

Vrugt, J. A., Diks, C. G. H., Gupta, H. V., Bouten, W., and Verstraten, J. M.: Improved treatment of uncertainty in hydrologic modeling: Combining the strengths of global optimization and data assimilation, Water Resour. Res., 41, w01017, https://doi.org/10.1029/2004WR003059, 2005. a

Wang, X. and Bishop, C. H.: A Comparison of breeding and ensemble transform Kalman filter ensemble forecast schemes, J. Atmos. Sci., 60, 1140–1158, https://doi.org/10.1175/1520-0469(2003)060<1140:ACOBAE>2.0.CO;2, 2003. a, b

Whitaker, J. S. and Hamill, T. M.: Evaluating methods to account for system errors in ensemble data assimilation, Mon. Weather Rev., 140, 3078–3089, https://doi.org/10.1175/MWR-D-11-00276.1, 2012. a, b

Whitaker, J. S., Hamill, T. M., Wei, X., Song, Y., and Toth, Z.: Ensemble data assimilation with the NCEP global forecast system, Mon. Weather Rev., 136, 463–482, https://doi.org/10.1175/2007MWR2018.1, 2008. a

Wu, C.-C. and Margulis, S. A.: Feasibility of real-time soil state and flux characterization for wastewater reuse using an embedded sensor network data assimilation approach, J. Hydrol., 399, 313–325, https://doi.org/10.1016/j.jhydrol.2011.01.011, 2011. a, b, c

Wu, C.-C. and Margulis, S. A.: Real-time soil moisture and salinity profile estimation using assimilation of embedded sensor datastreams, Vadose Zone J., 12, https://doi.org/10.2136/vzj2011.0176, 2013. a

Ying, Y. and Zhang, F.: An adaptive covariance relaxation method for ensemble data assimilation, Q. J. Roy. Meteorol. Soc., 141, 2898–2906, https://doi.org/10.1002/qj.2576, 2015.  a

Zhang, F., Snyder, C., and Sun, J.: Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble Kalman filter, Mon. Weather Rev., 132, 1238–1253, https://doi.org/10.1175/1520-0493(2004)132<1238:IOIEAO>2.0.CO;2, 2004.  a

Zhang, H., Hendricks Franssen, H.-J., Han, X., Vrugt, J. A., and Vereecken, H.: State and parameter estimation of two land surface models using the ensemble Kalman filter and the particle filter, Hydrol. Earth Syst. Sci., 21, 4927–4958, https://doi.org/10.5194/hess-21-4927-2017, 2017. a, b

Download
Short summary
Data assimilation methods like the ensemble Kalman filter (EnKF) can combine models and measurements to estimate states and parameters, but require a proper representation of uncertainties. In soil hydrology, model errors typically vary rapidly in space and time, which is difficult to represent. Inflation methods can account for unrepresented model errors. To improve estimations in soil hydrology, we designed a method that can adjust the inflation of states and parameters to fast varying errors.