EnKF with closed-eye period - towards a consistent aggregation of information in soil hydrology

Abstract

The representation of soil water movement exposes uncertainties in all model components. We assess the key uncertainties for the specific hydraulic situation of a 1-D soil profile with TDR (time domain reflectometry)-measured water contents. The uncertainties addressed are initial condition, soil hydraulic parameters, small-scale heterogeneity, upper boundary condition, and the local equilibrium assumption by the Richards equation. We employ an ensemble Kalman filter (EnKF) with an augmented state to represent and estimate all key uncertainties, except for the intermittent violation of the local equilibrium assumption. For the latter, we introduce a closed-eye EnKF to bridge the gap. Due to an iterative approach, the EnKF was capable of estimating soil parameters, Miller scaling factors and upper boundary condition based on TDR measurements during a single rain event. The introduced closed-eye period ensured constant parameters, suggesting that they resemble the believed true material properties. This closed-eye period improves predictions during periods when the local equilibrium assumption is met, but requires a description of the dynamics during local non-equilibrium phases to be able to predict them. Such a description remains an open challenge. Finally, for the given representation our results show the necessity of including small-scale heterogeneity. A simplified representation with Miller scaling already yielded a satisfactory description.

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/ doi:10.5194/hess-20-4999-2016 Author(s) 2016. CC Attribution 3.0 License.

Hannes H. Bauser1,2, Stefan Jaumann1,2, Daniel Berg1,2, and Kurt Roth1,3

1Institute of Environmental Physics (IUP), Heidelberg University, Heidelberg, Germany

2HGS MathComp, Heidelberg University, Heidelberg, Germany

3Interdisciplinary Center for Scientic Computing (IWR), Heidelberg University, Heidelberg, Germany

Correspondence to: Hannes H. Bauser ([email protected])

Received: 10 June 2016 Published in Hydrol. Earth Syst. Sci. Discuss.: 20 June 2016 Revised: 27 September 2016 Accepted: 7 November 2016 Published: 19 December 2016

Abstract. The representation of soil water movement exposes uncertainties in all model components. We assess the key uncertainties for the specic hydraulic situation of a 1-D soil prole with TDR (time domain reectometry)-measured water contents. The uncertainties addressed are initial condition, soil hydraulic parameters, small-scale heterogeneity, upper boundary condition, and the local equilibrium assumption by the Richards equation. We employ an ensemble Kalman lter (EnKF) with an augmented state to represent and estimate all key uncertainties, except for the intermittent violation of the local equilibrium assumption. For the latter, we introduce a closed-eye EnKF to bridge the gap. Due to an iterative approach, the EnKF was capable of estimating soil parameters, Miller scaling factors and upper boundary condition based on TDR measurements during a single rain event. The introduced closed-eye period ensured constant parameters, suggesting that they resemble the believed true material properties. This closed-eye period improves predictions during periods when the local equilibrium assumption is met, but requires a description of the dynamics during local non-equilibrium phases to be able to predict them. Such a description remains an open challenge. Finally, for the given representation our results show the necessity of including small-scale heterogeneity. A simplied representation with Miller scaling already yielded a satisfactory description.

1 Introduction

The description of soil water ow in the vadose zone with a mathematical model requires knowledge about material

EnKF with closed-eye period towards a consistent aggregation of information in soil hydrology

properties (typically characterized by soil hydraulic parameters), initial conditions, and boundary conditions. The material properties are especially difcult to determine, since they can neither be measured directly nor transferred directly from the laboratory to the eld.

Soil hydraulic parameters have been estimated inversely based on measurements of the temporal development of the hydraulic state with reviews, e.g., by Hopmans et al. (2002) and Vrugt et al. (2008). Inversions are based on the assumption of a perfect model and incorporate all errors and uncertainties into the estimated parameters, leading to a suboptimal estimation.

In contrast, data assimilation methods are capable of combining information from measurements and models into an optimal estimate of the geophysical eld of interest, but depend on the correct description of corresponding uncertainties (Reichle, 2008). Data assimilation methods include variational methods, particle lters, and ensemble Kalman lters (EnKF).

The EnKF, introduced by Evensen (1994), became popular in hydrology and was used to estimate, for example, water content states from satellite data (e.g., Reichle et al., 2002) or local water content measurements (e.g., De Lannoy et al., 2009).

Liu and Gupta (2007) pointed out that in hydrologic data assimilation the focus on state estimation is not sufcient and called for the incorporation of the errors in model structures and parameters.

Moradkhani et al. (2005) used a dual EnKF approach, whereas Vrugt et al. (2005) combined an EnKF and the shufed complex evolution metropolis algorithm in an iterative

Published by Copernicus Publications on behalf of the European Geosciences Union.

5000 H. H. Bauser et al.: EnKF with closed-eye period

way, to estimate parameters and states in a rainfallrunoff model. Chen and Zhang (2006) used an EnKF with augmented state to simultaneously estimate pressure and conductivity eld for saturated ow.

Based on the Richards equation, Li and Ren (2011) rst employed the simultaneous estimation of states and soil hydraulic parameters with the EnKF, which has been studied in more detail for purely synthetic cases (e.g., Wu and Margulis, 2011; Song et al., 2014; Erdal et al., 2015; Shi et al., 2015; Man et al., 2016) or with a focus on synthetic cases (e.g., Wu and Margulis, 2013; Erdal et al., 2014).

Synthetic cases offer the advantage of a direct control and knowledge about all uncertainties, which is a challenge in real-world cases. The characterization of uncertainties is critical for the success of a data assimilation scheme (Liu et al., 2012). Assuming the wrong uncertainty sources can lead to poorer performance (e.g., Crow and Van Loon, 2006; Zhang et al., 2015). If the correct uncertainty is unknown it can be better to overestimate uncertainties than to disregard or underestimate them (e.g., Jafarpour and Tarrahi, 2011; Zhang et al., 2015).

In this study, we focus on a real-world case to address the challenge of consistent aggregation of the information. We exercise this aggregation with the EnKF on a small example: a 1-D soil prole equipped with time domain reectometry (TDR) probes measuring water content, during a time period of less than 2 months. We assess all uncertainties in the representation of this particular situation qualitatively and design a three-stage approach to reduce the largest uncertainties or to consider them appropriately: rst, we improve the prior knowledge (initial condition and small-scale heterogeneity) to facilitate the subsequent estimation. Second, we perform an assimilation with a standard EnKF approach with an augmented state to directly reduce uncertainties in states, soil hydraulic parameters, small-scale heterogeneity and upper boundary condition. We introduce iterations over the complete EnKF scheme to cope with the short data set, and to determine times when the underlying local equilibrium assumption by the Richards equation is violated. We dene this specic time as the closed-eye period, because in the third stage, we only estimate states, but keep parameters constant during this time and thus prevent the incorporation of the uncertainties in the dynamics into the parameters. The estimation of soil hydraulic parameters is only performed before and after this closed-eye period.

The remainder of the paper is organized as follows: Sect. 2 describes our understanding of (i) a representation of a physical system in general, (ii) the hydrological system we want to represent, (iii) the specic representation and its uncertainties in detail, (iv) the concept of knowledge fusion and(v) the EnKF. Section 3 presents the three stages: (i) improving the prior, (ii) standard EnKF and (iii) closed-eye EnKF and compares the prediction capabilities of standard and closed-eye EnKF. In Sect. 4 we present our conclusions and their implications.

2 Background

2.1 Representation

In this work, we call the mathematical description of a physical system a representation, which comprises in the most general sense the following four components: dynamics, forcing, subscale physics and states.

Dynamics: the dynamics species how the state of the system is propagated in space and time at the scale of interest;

Forcing: the forcing or coupling of the system to the superscale physics in space (initial state) and time (boundary condition);

Subscale physics: the subscale physics not explicitly described by the dynamics, but with a parameterization comprising one or more parameters and their spatial and temporal distribution;

State: distribution of the variable propagated by the dynamics, at a specic time. Notice that this excludes parameters.

According to this denition a representation is specic for each system. Hence, we rst introduce the system, and the details of the representation are presented afterwards.

2.2 Test site

We aim to represent the water movement in a soil prole at the Grenzhof test site close to Heidelberg, Germany.Since 2003, experiments have been conducted at the test site.In 2004, a weather station was built, which measures precipitation and further atmospheric data (wind, temperature, incoming and outgoing long- and shortwave radiation, relative humidity and air temperature) in 10 min intervals. A detailed description of the test site can be found in Wollschlger et al. (2009).

In 2009, a soil prole was equipped with 11 TDR probes measuring water content hourly. The soil prole itself, described explicitly by Schenk (2011), is depicted in Fig. 1. It consists of four different soil layers up to the depth of 1.4 m.The top three layers are equipped with three TDR probes each, and the bottom layer with two TDR probes. The prole is kept free of vegetation at the surface.

The complete time period considered comprises 60 days from 1 October 2011 (day 1) until 29 November 2011 (day 60). The boundary condition, along with the water content measured by the topmost TDR, is shown in Fig. 2.

2.3 Representation of soil water ow

For this specic soil prole we can now formulate the representation consisting of the four components described in

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/

H. H. Bauser et al.: EnKF with closed-eye period 5001

Figure 1. Soil prole at the Grenzhof test site close to Heidelberg, Germany (Schenk, 2011). It consists of 4 different soil layers and is equipped with 12 temperature probes and 11 TDR probes. The surface is kept free of vegetation.

10 A B

0.5

C

D

Time [day]

Evaporation / precipitation [mm day]

-1

Water content in 8.5 cm depth

0

0.4

-10

0.3

-20

0.2

-30

0.1

0

-40 0 10 20 30 40 50 60

Figure 2. Boundary condition and topmost TDR from 1 October (day 1) until 29 November (day 60) at the Grenzhof test site. We distinguish 4 sections: A (day 12), B (day 317), C (day 1822) and D (day 2360). The evaporation is calculated using the reference FAO PenmanMonteith equation. Before the rst rain the evaporation is set to 0, due to a previous dry period of over a month. The precipitation is displayed with a 6 h resolution.

Sect. 2.1. For each we give the a priori description, as well as the corresponding assumptions and sources of uncertainty.

Dynamics: for a continuous porous medium, the Richards equation (Richards, 1931) describes the change of volumetric water content ()

@

@t r [K( )[rhm( ) 1]] = 0, (1)

with the isotropic conductivity K (L T1) and matric head hm (L). The Richards equation itself is process-based and is expected to describe the actual physics adequately as long as its underlying assumptions are fullled: it assumes

the soil water movement to be a single-phase process. Hence, the Richards equation is expected to fail in very dry or wet conditions. Both are not observed at the Grenzhof test site during the time period considered. The corresponding uncertainties are small. The Richards equation further assumes local equilibrium. This is assumed to be violated during strong inltration events, e.g., during heavy rainfall leading to preferential ow. We expect local non-equilibrium at the Grenzhof test site, but it is a priori not clear how strong the forcing has to be or how to describe these cases. This creates uncertainty in the dynamics during rain events.

We additionally assume that horizontal ow is negligible and that we can describe the system one-dimensionally without additional sources or sinks. The at terrain and horizontal layering at the Grenzhof test site supports this assumption. Still, heterogeneity, which is not represented, might introduce two-dimensional ow.

Material properties: the subscale physics not explicitly represented in Eq. (1) are the dependencies of K and hm on . Several parameterizations are available.

We follow Wollschlger et al. (2009) and choose the Mualemvan Genuchten parameterization (Mualem, 1976; Van Genuchten, 1980):

K([Theta1]) = K0[Theta1] [bracketleftbigg]

1

h1 [Theta1]n/[n1][bracketrightBig]11/n[bracketrightbigg]2, (2)

h[Theta1]n/[n1] 1[bracketrightBig]1/n

, (3)

with the saturation [Theta1] ()

[Theta1] :=

hm([Theta1]) =

1

r

s r

. (4)

The parameterization reduces the description of the subscale physics to a set of six parameters: s (), r (), (L1),

www.hydrol-earth-syst-sci.net/20/4999/2016/ Hydrol. Earth Syst. Sci., 20, 49995014, 2016

5002 H. H. Bauser et al.: EnKF with closed-eye period

n (), K0 (L T1) and (). We do not expect signicant uncertainties due to assuming this particular parameterization, because we only observe a rather small water content range.

The parameters s and r are the saturated and residual water content. The parameters and n can be associated with an inverse air entry value and the width of the pore size distribution respectively, although there is no direct physical relation.K0 is the saturated conductivity and is a fudge factor.

The spatial distribution of these material properties has to be dened on the scale of the dynamics. The large-scale structure for this soil prole is made up of four soil layers, where we assign an individual parameter set to each layer.

Although these parameters cannot be measured directly and are often assigned with the largest uncertainties in soil hydrology, we can in this case use the soil parameters determined by Wollschlger et al. (2009) through an inversion for a soil prole close by. They are expected to differ from the soil hydraulic parameters at the represented soil prole due to spatial heterogeneity.

We do not include a description of hysteresis, leading to additional uncertainties connected with the parameterization.

We assume the parameters to be constant in time, which appears justied, since we only investigate a short time period of less than 2 months.

The assumption of homogeneous soil layers is presumably wrong. A possibility to describe the texture inside soil layers in a simplied way is Miller scaling (Miller and Miller, 1956). They assume geometrical similarity and scale the material properties at each location (grid cell) according to the following:

K( ) = K ( ) 2, (5) hm( ) = h m( )

have to assume an uncertainty of about 20 % for the precipitation as well.

For the lower boundary condition there are no measurements at the site. However, constant water content measurements in the lowest soil layer during the observed time period indicate that the dynamics is decoupled from the groundwater table. Therefore, we follow Wollschlger et al. (2009) and add a fth layer with soil parameters for sand from Carsel and Parrish (1988) up to a depth of 4 m in our representation.The sand layer decouples the dynamics from the groundwater table that we keep constant at the depth of 4 m. This leads to uncertainties, but they are expected to be negligible.

The forcing in space is the initial condition. The initial condition is difcult to estimate. A simple approach is to use the measured water contents and interpolate them linearly with constant extrapolation to layer boundaries. This leads to large uncertainties in between.

States: with the three components dynamics, material properties and forcing we can predict the temporal development of the states, in this case the water content. These states will have large uncertainties due to the uncertainties in the other components.

An additional source of uncertainty is the simulation of the dynamics with a numerical solver, here with MuPhi (Ippisch et al., 2006). The assumption of local equilibrium in the Richards equation has to be granted on the grid size of the numerical model as well. We choose a resolution of 1 cm. For the represented soil and boundary conditions this resolution is sufcient to resolve the occurring fronts smoothly, so that we can neglect effects caused by the numerical solver.

Water content determined with TDR probes will be used to improve the representation. The water content values are calculated from the temperature-corrected electric permittivity with the complex refraction index model (CRIM) following Roth et al. (1990). The porosities required there were obtained from volumetric soil samples.

We assume an unbiased uncertainty of 0.01 for the water contents. This is the largest noise observed with the TDR probes.

The assumptions made for the representation of the soil water movement at the Grenzhof test site are summarized in Table 1. In the second column we give our qualitative assessment of the uncertainties: small no need to represent the uncertainty explicitly; intermediate it might be necessary to represent these uncertainties, and we decide not to but rather keep them in mind; large uncertainties must be represented with the goal of reducing them. Our ndings depend on this uncertainty assessment and might be affected if it does not hold.

2.4 Knowledge fusion

We dene knowledge fusion as the consistent aggregation of all information pertinent to some observed reality. In the presented situation at the Grenzhof test site this would require

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/

1 , (6)

where denotes the reference material properties of the corresponding layer and () is the Miller scaling parameter at this location. This means that we have to dene a Miller scaling factor for each grid cell. The assumption of geometric similarity does not necessarily hold and leads to some uncertainties. The scaling factor themselves are a priori unknown. We assume no heterogeneity (scaling factors of 1), which is most likely wrong and introduces large uncertainties.

Forcing: the forcing in time is the upper and lower boundary condition. For the upper boundary condition, we can measure rain and use the meteorological data from the on-site weather station to calculate the evaporative water ux using the reference FAO PenmanMonteith equation following Allen et al. (1998). Here in particular the estimation of the evaporation is highly uncertain. The daily mean value can differ 1020 % (Foken, 2006). In dry conditions the evaporation is limited by the soil, which leads to additional uncertainties.

At the Grenzhof test site the rain gauge had not been calibrated during the period of consideration. Because of that we

H. H. Bauser et al.: EnKF with closed-eye period 5003

Table 1. Prior uncertainties of the representation of soil water movement at the Grenzhof test site during October and November 2011. We assess the uncertainty caused by assumptions in a qualitative way: small no need to represent the uncertainty explicitly; intermediate it might be necessary to represent these uncertainties, and we decide not to but rather keep them in mind; large uncertainties must be represented with the goal of reducing them.

Assumption Uncertainty

Dynamics

Single phase description (Richards equation) smallLocal equilibrium (Richards equation) smalllarge 1-D, no sources and sinks intermediate

Material properties

Mualemvan Genuchten parameterization for material properties small Mualemvan Genuchten parameter values from Wollschlger et al. (2009) large Mualemvan Genuchten parameters constant in time smallNo hysteresis intermediate Large-scale architecture described with layers small Small-scale architecture described with Miller scaling smallintermediate Miller scaling factors are 1 large

Forcing

Upper boundary condition from rain gauge and FAO PenmanMonteith equation large Lower boundary condition constant groundwater table decoupled with soil layer small Linear interpolation for initial condition large

Operational aspects

Numerical representation small Water content measurements without biases and uncertainty of 0.01 intermediate

the quantitatively correct description of all uncertainties in the representation and a subsequent optimal reduction of all these uncertainties based on all additionally available information. For the measurement part these are primarily the water content data from the 11 TDR probes, but any other information, even expert knowledge, should be incorporated as well. So far, however, this goal is only partly feasible.

With an inversion, all information and uncertainties would be included in the parameters, which is not a consistent aggregation of the information, since the structural errors in the other components are not represented.

Data assimilation methods are capable of representing all uncertainties and have been expanded to estimate not only states but also parameters. For an aggregation of all information we have to reduce uncertainties in all representation components. In this study we aim to describe and reduce the uncertainties classied as large in Table 1 and decided to neglect the ones classied as intermediate or small with the EnKF.

2.5 Ensemble Kalman lter

The EnKF is a data assimilation method that uses a Monte Carlo approach for an optimal state estimation, based on the assumption of unbiased Gaussian error distributions. The

EnKF incorporates the measurements sequentially by alternating between a forecast step (superscript f), which propagates the state in time, and an analysis step (superscript a), which incorporates the information from the measurements at this time to improve the state.

These two steps are now explained in more detail specifically for the given representation. A general description can be found, for example, in Burgers et al. (1998).

Forecast: rst an ensemble of water content states i is

propagated forward in time from k 1 to k: f,ik = M [parenleftBig]

a,ik1,ik,ik,[notdef]ik[parenrightBig] +

ik, (7)

where M denotes the propagation with the Richards equation and i, i and [notdef]i are soil hydraulic parameters, Miller scaling parameters, and the upper boundary condition of ensemble member i, respectively. In this way, all the uncertainties from , , and [notdef] are explicitly represented and propagated forward in time with the non-linear dynamics leading to the forecast ensemble f, which represents these uncertainties.Uncertainties not explicitly represented are incorporated into the process noise , which is difcult to describe, however.We set the process noise to 0. In this way, we represent all the uncertainties classied as large in Table 1 directly, except for uncertainties caused by the local equilibrium assumption, which is expected to only hold outside of strong rain events.

www.hydrol-earth-syst-sci.net/20/4999/2016/ Hydrol. Earth Syst. Sci., 20, 49995014, 2016

5004 H. H. Bauser et al.: EnKF with closed-eye period

We decide to not represent intermediate and small uncertainties.

The uncertainty of the forecast ensemble is characterized based on the Gaussian assumption with the state error covariance matrix P:

Pfk =

1N 1

As parameters to be incorporated, we choose , n and K0 for each of the four layers as well as the saturated water content s for layers 3 and 4. For the upper two layers, the measured water contents were far away from s and r estimated by Wollschlger et al. (2009). Thus, they are irrelevant for the estimation, as changes in s and r can be compensated by the other parameters. For the bottom two layers measurements exceeded s. Thus, it has to be added in the augmented state.

We do not include the fudge factor , since the observed water content ranges are rather small. For other situations however, incorporation of might improve the results.

Since the EnKF assumes Gaussian distributions and linear correlations, we do not directly use K0 but the logarithm (to which the water content response behaves more linearly).

The dimensions of the Miller scaling factors are reduced as well. Only the scaling factors at measurement locations are added. The whole Miller scaling eld is determined by linear interpolation between measurement locations inside a layer and a constant extrapolation towards layer boundaries.As the measurements yield only a little information about the small-scale architecture away from their location, we expect that this assumption has only a slight inuence on the results. Again, the logarithm of the parameter is used in the augmented state.

The upper boundary condition is already a scalar and is added as the ux at the surface to the augmented state.

In this way, soil parameters, Miller scaling factors and the upper boundary condition are updated in addition to the states in the analysis.

The expansion to an augmented state changes the propagation in time. Each component needs an individual forward propagation. We assume the soil hydraulic parameters and Miller scaling factors to be constant in time. This is not possible for the upper boundary condition where the forward equation is unknown from a soils perspective. However, measurements are available to estimate the evaporation and precipitation. Hence, we assume the forward model constant until a new estimation is available. Then, we switch to the estimated boundary condition. To base the improvement of the upper boundary condition on several measurements, we reduce the temporal resolution of precipitation measurements to change daily and at transitions between precipitation and evaporation. This means that the upper boundary condition is treated like the parameters, except that the value can change in the forward propagation. The original temporal resolution of the precipitation data (10 min) is not required due to the dissipative nature of the Richards equation, which smooths the inltration front before it reaches the rst TDR sensor at a depth of 8.5 cm. The estimation of the averaged boundary condition will ensure that there is no global bias on the parameter estimation during the rain event, but could lead to small short-term parameter drifts within a rain event. The expanded forward propagation of the augmented state is as follows:

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/

N

Xi=1

f,ik fk[parenrightBig][parenleftBig]

f,ik fk[parenrightBig]T

, (8)

with the number of ensemble members N and the ensemble mean .

Analysis: to incorporate measurements, we have to quantify their uncertainties by dening the measurement error covariance matrix R:

Rk =

1N 1

N

Xi=1

[epsilon1]ik

[parenrightBig][parenleftBig]

T

, (9)

with the measurement errors [epsilon1]. We assume the measurement errors to be uncorrelated and drawn from a normal distribution with the assumed TDR standard deviation of 0.01.

The Kalman gain K weights the forecast state error covariance and the measurement error covariance

Kk = PfkHTk [bracketleftBig]

HkPfkHTk + Rk[bracketrightBig]1

[epsilon1]ik

, (10)

with the measurement operator H mapping from the state space to the measurement space, and the Kalman gain K mapping from the measurement space back to the state space weighted with the uncertainties of measurements and states. Based on the Kalman gain the weighted difference of forecast and measurements is added to the forecast states to obtain the analysis states

a,ik = f,ik + Kk [bracketleftBig]

zk Hk f,ik + [epsilon1]ik[bracketrightBig]

, (11)

where z is the measurements. The measurement errors, [epsilon1], have to be added to the difference so that the ensemble describes the analysis error covariance correctly (Burgers et al., 1998). The descriptive explanation is that the measurements, z, already include measurement errors, while the measurements of the forecast state, H f, do not. Hence, the error has to be added to compare equivalents.

By alternating between forecast and analysis, the information of all measurements is incorporated to achieve an improved estimation of the states at each time step.

By replacing the water content state with an augmented state , not only the water content states, but all components of the augmented state can be estimated with the EnKF. We choose to include all the directly represented uncertain components in the augmented state = [ , , , [notdef]] consisting of

water content, soil hydraulic parameters, Miller scaling parameters and the upper boundary condition.

Adding components to the augmented state increases its dimension, which in turn requires a larger number of ensemble members, leading to a higher computational effort. To minimize this effect, we keep the added components as small as possible.

H. H. Bauser et al.: EnKF with closed-eye period 5005

fk = M ak1,ak1,ak1,[notdef]ak1 [parenrightbig]

, (12)

fk = ak1, (13) fk = ak1, (14)

[notdef]fk =

Spurious correlations and non-Gaussian distributions can lead to lter inbreeding and ultimately lter divergence (e.g., Houtekamer and Mitchell, 1998; Hamill et al., 2001). To prevent this, the ensemble can be inated by a multiplicative factor (Anderson and Anderson, 1999):

i

= p [parenleftBig]

i

[notdef]meas if a measurement is available, [notdef]ak1 else.

[parenrightBig] +

(15)

The analysis information for the additional augmented state components lies in the covariance between states at measurement locations and the components. The covariance linearizes this relation that can be highly non-linear.

To reduce the impact of the linearization of nonlinear dependencies, we employ damping factors (Hendricks Franssen and Kinzelbach, 2008), different for each augmented state component (Wu and Margulis, 2011) d = [d , d, d , du]. The correction vector in Eq. (11) is mul

tiplied by the corresponding factors. Hendricks Franssen and Kinzelbach (2008) investigated damping factors between 0.1 and 1.0, while Erdal et al. (2014) employed a damping factor of 0.3 for the parameters. The smaller the value, the better a non-linearity can be handled, but the slower the estimated value approaches the optimal value. We choose the values d = [1.0, 0.1, 0.1, 0.1]. This means that, based on the water

content measurements, there is no damping for the update of the water content state, but a strong damping for soil hydraulic parameters, Miller scaling factors and upper boundary condition.

Due to a limited ensemble size, the EnKF will show spurious covariances between uncorrelated state components. These can be reduced by increasing the ensemble size, or more computationally efcient by introducing a localization (Houtekamer and Mitchell, 2001), where external information about uncorrelated states is introduced. This can be achieved by multiplying the error covariance matrix by the correlation matrix (with entries from 0 to 1) with an entry-wise product. This alters the calculation of the Kalman gain to

Kt = Pft[parenrightBig]

HTt

hHt

i. (17)

We choose the ination as proposed by Anderson (2009), where the ination factor is adapted temporally and spatially for the states. We expand this for the other augmented state components by simply summing over the ination factors at the measurement locations weighted with the correlations.

We iterate the whole EnKF scheme and start the next iteration with the nal estimation of soil parameters, Miller scaling factors and upper boundary condition of the previous iteration. These iterations differ from the typically applied iterative EnKFs like the restart EnKF or conrming EnKF (e.g., Song et al., 2014) and rather resemble the iterations in the Kalman ensemble generator by Nowak (2009), who used a modied EnKF to estimate parameters only. The full iterations applied are required to estimate constant augmented state components even with a small damping factor and a rather short time period of data. This is especially important for the upper boundary condition, where the time periods of constant values are short.

The operational assumptions of the method lead to a suboptimal estimation of the state in each time step. Furthermore, due to the non-linear dynamics, the assumption of Gaussian distributions does not hold. It is not clear how this affects the EnKF performance in detail. The Gaussian assumption leads to a linearized state update in the analysis step. This induces erroneous updates of those state components with dominant non-linear relation between states and measurements. These errors are alleviated by employing the damping factor, which reduces the update but as a consequence also reduces the incorporation of measurement information. Furthermore, we use the largest observed measurement noise to characterize the TDR uncertainties. This leads to possibly too-large measurement uncertainties, which have a similar effect as an additional damping factor.

We do not expect a strong inuence from this sub-optimal state estimation on the mean value of the results for the soil hydraulic parameters, Miller scaling factors and upper boundary condition. However, the nal value will be approached more slowly than in an ideal estimation. This effect is handled by the iterative approach.

On the downside, incorporating the same measurement information several times will lead, together with the other limitations, to incorrect quantitative uncertainties. However, the spatially and temporally adaptive covariance ination by Anderson (2009) ensures that measurements and states agree with each other within their uncertainties. This is only a partial solution, however, as the information is only propa-

www.hydrol-earth-syst-sci.net/20/4999/2016/ Hydrol. Earth Syst. Sci., 20, 49995014, 2016

i

Pft[parenrightBig]

HTt + Rt[bracketrightBig]1, (16)

where the entrywise product is indicated by .

For the water content state, the covariances are reduced with increasing physical distance. We use the fth-order piecewise rational function dened by Gaspari and Cohn (1999), which is similar to a Gaussian but with local support, which introduces a cutoff length. We chose a very large empirical value of 153 cm (which corresponds to 2c in the fth-order piecewise function).

For the soil parameters, we can localize even more strongly by only allowing covariances (entries of 1) between parameters and measurement locations in the respective soil layer and the rst measurement locations in the neighboring layers. For the Miller scaling factors, only covariances to the corresponding measurement locations are used. All other entries are set to 0.

5006 H. H. Bauser et al.: EnKF with closed-eye period

Stage 1: Improving the prior

EnKF with 40 ensemble members estimation of states and Miller scaling factors

Improved Miller scaling factors (mean and standard deviation)

50 iterations

Improved Miller scaling factors (mean values)

EnKF with 100 ensemble members estimation of states only

Improved initial conditions for time period C

Time

period A

Time

period B

Stage 2: Standard EnKF

EnKF with 100 ensemble members estimation of states, soil parameters, Miller scaling factors and upper

Improved soil parameters, Miller scaling factors and

upper boundary condition (mean and standard deviation)

10 iterations

Time

period C

Improved soil parameters and Miller scaling factors (mean and standard deviation)

time interval for closed-eye period

Stage 3: Closed-eye EnKF

Closed-eye EnKF with 100 ensemble members

estimation of states, soil parameters, Miller scaling factors and upper

Improved soil parameters, Miller scaling factors and upper boundary condition

(mean and standard deviation)

10 iterations

Time

period C

Final estimation of believed true soil parameters,

Miller scaling factors and upper boundary condition

Figure 3. Block diagram illustrating the three-stage approach. The different time periods are shown in Fig. 2.

gated through the correlations to the measurement locations.

Therefore, we will not interpret the uncertainties of the ensemble in a quantitative way. Nevertheless, we do expect them to hold qualitative information that can be interpreted.

3 Results

The complete time period considered comprises 60 days and ranges from 1 October (day 1) until 29 November (day 60). The boundary condition, along with the water content measured by the topmost TDR, is shown in Fig. 2. The whole time period is separated into four sections: A (day 12), B (317), C (1822) and D (2360).

Since soil parameters can only be estimated within the ob-

served water content range and will not be valid outside of this range, a rather large rain event is desirable. On the other hand we do not represent uncertainties associated with the assumption of local equilibrium by the Richards equation, which is violated during strong rain events. Time period C (1822 October 2011) combines both: a rather large total rain amount (18.2 mm) and a small maximal intensity (0.7 mm in 10 min).

We designed a three-stage approach (Fig. 3) to improve the representation by incorporating the information from the water content measurements of the 11 TDR probes: (1) improving the prior information during time periods A and B,(2) a standard EnKF approach, and (3) a closed-eye EnKF approach, both during time period C with the considered rain event. The three stages are discussed in the following sections, with results analyzed through forecasts in time periods B, C and D.

3.1 Stage 1: improving the prior

Highly uncertain properties can exacerbate the performance of the EnKF. For example, Huber et al. (2011) showed that the state estimation with an EnKF is superior in case of properly precalibrated model parameters. In our case, there is no prior knowledge about small-scale heterogeneity, which can lead to difculties in the estimation process: in a uniform layer, the water content is typically higher in deeper locations. Heterogeneity can reverse this and can lead to higher water contents above lower water contents inside a layer. Until the inltration front of the rain event reaches these positions, the EnKF cannot distinguish between heterogeneity and an inltration front propagating at very low speed. This

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/

H. H. Bauser et al.: EnKF with closed-eye period 5007

can easily lead to too-small saturated conductivities, especially in lower soil layers that are not reached by the inltration front. Even upper layers can require a tuned relation of the prior uncertainties of Miller scaling factors and saturated conductivity.

To avoid too-small saturated conductivities, we improve the prior of the heterogeneity using the measurements from time interval A. Prior to A, there had been 1 month without rainfall leading to little dynamics during A. We start with the soil parameters estimated by Wollschlger et al. (2009) and linearly interpolated initial condition (with constant extrapolation to layer boundaries). Parameters and boundary condition are kept constant, but states and Miller scaling factors are updated with an EnKF of 40 ensemble members, where the initial uncertainty of the natural logarithm of the scaling factors was chosen to be 0.4.

As time interval A is only a 2-day period, the lter will not be able to reach constant Miller scaling values during this short time. Therefore, we iterate 50 times over this period, which then does lead to constant scaling factors. The resulting scaling eld of the improved prior (Fig. 4) shows that there is heterogeneity within each layer that requires a representation. Layer 4 in particular exhibits a strong heterogeneity, which is already observable in Fig. 1.

Additionally, a good initial state can improve the estimation. Therefore, we guide the state with an EnKF (only state estimation) with 100 ensemble members through time period B, and achieve a better representation for the initial state for time period C than interpolating between the measurements there.

3.2 Stage 2: standard EnKF

During time period C, we improve the representation with the following uncertain components: improved initial condition, soil hydraulic parameters estimated by Wollschlger et al. (2009), improved Miller scaling factors, and upper boundary condition from precipitation measurements and estimated with the reference FAO PenmanMonteith equation.We improve soil parameters, Miller scaling factors and the upper boundary condition along with the states. The estimated mean values and standard deviations of soil parameters, Miller scaling factors and upper boundary condition are used to generate the ensemble for the next iterations. The initial condition is kept the same. We iterate 10 times over the time period C. This estimation is performed with an ensemble size of 100 members and will be referred to as the standard EnKF. This approach is comparable to the SODA approach by Vrugt et al. (2005) but is computationally more efcient.

The initial ensembles of initial condition and Miller scaling are determined as described in stage 1. For the Miller scaling factors the uncertainty is increased again compared to the estimations in the prior, since they were estimated under the assumption of xed soil parameters. Now, they have to be

0 Layer

1

2

3

4

PriorImproved priorStandard EnKF Closed-eye EnKF Closed-eye EnKF ensemble

Depth [cm]

140 0 1 2 3 4 5Miller scaling f actor

Figure 4. Mean values of the Miller scaling eld. Soil layers are indicated by different gray scales. The heterogeneity is a priori unknown and the according prior is set to 1. The natural logarithm of the scaling factors is estimated at the measurement locations and interpolated linearly between (with constant extrapolation to layer boundaries). Already the improved prior can describe the main features. The further estimations with the standard and closed-eye EnKF lead to further small changes mainly in the rst layer. For the closed-eye EnKF the actual ensemble indicating the uncertainties is additionally shown.

able to adapt to changing parameters. The uncertainty of the natural logarithm of the scaling factors is chosen to be 0.1.For the boundary condition shown in Fig. 2, an uncertainty of 20 % is assumed. The initial soil hydraulic parameters and their uncertainties are summarized in Table 2.

In order to check the improved results, we do not show the states from the last iteration but actually run another iteration, however, this time, without incorporating the measurement information. This corresponds to an ensemble of forward runs during time period C. It is a much more strict test for the quality and objectivity of the assimilation, because the states now cannot be adjusted. In fact, allowing incorporation of measurements leads to a much better agreement. The results for all four layers are shown in Fig. 5. The graphs show that the measurements can be represented well by the ensemble, but deviations occur, mainly during the inltration.

There is heterogeneity inside the layers. For example, in the rst layer, the middle TDR shows the highest water content. This effect is strongest in layer 4. Here, the water content of the TDR about 25 cm above the next has a water content almost 0.1 larger. The results calculated with soil parameters by Wollschlger et al. (2009) with the assumption of homogeneous soil layers are not able to describe the measured water content.

The results for the estimated Miller scaling factors are shown in Fig. 4. The Miller scaling factors were adjusted along with the hydraulic parameters during this standard EnKF approach. This is necessary, because changing parameters require adjusted Miller scaling factors as well. The largest changes occur in layer 1, where information is introduced by the water content changes during the inltration.

www.hydrol-earth-syst-sci.net/20/4999/2016/ Hydrol. Earth Syst. Sci., 20, 49995014, 2016

5008 H. H. Bauser et al.: EnKF with closed-eye period

Table 2. Soil hydraulic parameters as estimated by Wollschlger et al. (2009) and our initial uncertainties for the standard EnKF.

Layer log10(K0) n s r

[m1] [m s

1

]

1 4.8 2.0 4.6 1.0 1.61 0.20 0.39 0.020 0.5

2 11.0 2.0 5.3 1.0 1.80 0.20 0.40 0.070 0.5

3 7.4 2.0 6.0 1.0 1.43 0.20 0.43 0.02 0.090 0.5

4 2.9 2.0 5.4 1.0 1.36 0.20 0.43 0.02 0.100 0.5

5 14.5 4.1 2.68 0.43 0.045 0.5

Parameters without uncertainty are not included in the augmented state and are not estimated. The saturated water content s of the layers 3 and 4 were already estimated during the prior estimation in time period A, because measurement values exceeded the value of 0.36 determined by Wollschlger et al. (2009).

(a)

Layer 1

(b)

Layer 2

0.3 C

0.3 C

TDR (42.5 cm) TDR (52 cm)

TDR (62 cm)

0.25

0.25

Vol. water content

0.2

Vol. water content

0.2

0.15

0.15

TDR (8.5 cm) TDR (19 cm) TDR (28.5 cm)

0.1

0.1

17 18 19 20 21 22

17 18 19 20 21 22

Time [day]

Time [day]

(c)

Layer 3

(d)

Layer 4

0.4 C

0.5 C

0.45

0.35

Vol. water content

0.4

0.3

Vol. water content

0.35

0.3

0.25

TDR (76 cm) TDR (87 cm) TDR (98 cm)

0.25

TDR (111 cm) TDR (136 cm)

0.2

0.2

17 18 19 20 21 22

17 18 19 20 21 22

Time [day]

Time [day]

Figure 5. Forecast at the measurement positions of the four different layers during time period C. The solid lines show the mean value of a total of 100 different forecasts. The pale colors show the results from 25 of these, with soil parameters, Miller scaling parameters and boundary condition sampled from the distributions estimated with the standard EnKF and the initial condition, which was actually used for the estimation itself. The dashed lines show the measurements, the dotted lines results from the original soil parameters by Wollschlger et al. (2009). The measurements were recorded hourly and occurring data gaps were interpolated linearly.

The estimated soil parameters, including their uncertainties, are summarized in Table 3. Although the uncertainties are quantitatively incorrect, they still indicate how much information about each parameter is in the measurements. This can be seen especially well for K0. The estimation of these parameters requires dynamics. Therefore, the smallest uncertainty is observed in layer 1 and increases with the layers. In layer 4 the uncertainty of log10(K0) was reduced least, but still from 1.0 to 0.4. This shows the difculties of the EnKF.

As we did not observe a change in water content, it is unlikely that any ux was introduced during the rain event. Nevertheless, the uncertainty was reduced, which can lead to difcul-ties when there actually is information. The development of the parameters K0, and n of the rst layer are shown in

Fig. 6. Towards later iterations, all three show the same fea-

ture: The value increases during the rain event and decreases before and after, but does not change signicantly between initial and nal values. As we added the boundary condition into our estimation, we can exclude effects due to a bias in the precipitation to cause this parameter shift. The reduced temporal resolution in the boundary condition could only cause short-term parameter changes within the rain event.

We interpret the parameter shift during the rain event as follows: the assumption of local equilibrium during the rain event is wrong, which leads to preferential ow. The inltration is thus too fast for the actual parameters (Fig. 5). The preferential ow leads to an apparent increase in the saturated hydraulic conductivity during the rain event, which is reduced again afterwards. This also shows that our initial as-

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/

H. H. Bauser et al.: EnKF with closed-eye period 5009

Table 3. Results for the soil hydraulic parameters estimated with the standard EnKF.

Layer log10(K0) n s r

[m1] [m s

1

1 4.7 0.3 5.23 0.04 1.43 0.01 0.39 0.02 0.5

2 10.1 0.7 4.08 0.25 1.66 0.02 0.40 0.07 0.5

3 8.2 1.0 6.04 0.31 1.33 0.03 0.42 0.01 0.09 0.5

4 2.3 0.3 5.71 0.40 1.31 0.01 0.43 0.01 0.10 0.5

Parameters without uncertainty are not included in the augmented state and are consequently not estimated.

]

Figure 6. Evolution of the parameters K0, and n of the rst layer during the estimation with standard EnKF (purple) and the closed-eye EnKF (green) for all iterations (a, c, e) and the last iteration (b, d, f).

sumption of zero process noise is wrong. The uncertainty in

troduced by the dynamics is not represented.

3.3 Stage 3: closed-eye EnKF

We enhance the EnKF with a closed-eye period during times when the local equilibrium assumption of the Richards equation does not hold. During this time, soil hydraulic parameters and Miller scaling factors are not adjusted but kept constant. In contrast, the water content state is continuously updated. In this way, the state is guided on the basis of measurements through times with uncertain dynamics without incorporating the dynamics uncertainties into the parameter estimation. The estimation of soil hydraulic parameters is only performed before and after this closed-eye period.

We use an identied non-constant parameter to dene the closed-eye period. We choose K0, because it shows the changes most prominently and we can directly relate it to ob-served physics. The closed-eye period is dened as the time period between the minimum and the maximum value of K0 during the last iteration with the standard EnKF (Fig. 7).

To perform the closed-eye EnKF estimation, we add another 10 iterations to the previous standard EnKF iterations, but with the closed-eye period for the rst layer. We continue with the previously estimated mean values, but increase the uncertainty for the rst layer again (to half the initial uncertainty) to allow that the new values can possibly deviate from the standard EnKF results.

www.hydrol-earth-syst-sci.net/20/4999/2016/ Hydrol. Earth Syst. Sci., 20, 49995014, 2016

5010 H. H. Bauser et al.: EnKF with closed-eye period

Time [day]

10 C

0.6

0.3

C

Evaporation / precipitation [mm day]

-1

0

0.5

Water content in 8.5 cm depth

0.25

Vol. water content

-10

0.4

-20

0.3

0.2

-30

0.2

0.15

-40

0.1

TDR (8.5 cm) TDR (19 cm) TDR (28.5 cm)

0

-50 17 18 19 20 21 22

0.1

Figure 7. Boundary condition along with the topmost TDR water content measurements during time period C. The closed-eye time period was determined based on the changes in parameter K0 (Fig. 6). During this time the dynamics cannot describe the water content measurements any more, most likely because of the violation of the local equilibrium assumption. The boundary condition is displayed in the reduced resolution assumed during the estimation.

The additional iterations result in new soil parameters and Miller scaling factors. The scaling factors changed only marginally compared to their regular estimates (Fig. 4). Changes in the parameters of the rst layer are larger (Fig. 6). We recognize two things: rst, the values differ from each other. For n and the changes are rather small, but for K0 the difference is signicant. The closed-eye estimation yields approximately half the hydraulic conductivity as before. Second, the parameters now are constant over time and do not show uctuations any more.

This indicates that we could extract times when the Richards equation is actually valid and we are able to determine soil hydraulic parameters that resemble the believed true material properties. This means, we believe, that if the parameters are constant in time, the parameters can represent reality in the observed water content range during times when the underlying assumptions hold. There is also an apparent downside, though. We cannot use the measurements during the closed-eye period to estimate the parameters. This leads to a smaller observed water content range, limiting the parameter estimation possibilities. We call this apparent, because the corresponding interval does not contain valid information about the modeled system in the rst place.

As a related aspect, the estimation of the boundary condition is not meaningful during the closed-eye period any more. Consequently only an estimation of the evaporation before and after the rain event is possible. Due to the short time interval and small evaporation ux, in our case the effect is negligible.

The forecast during time period C (analogous to the standard EnKF) is shown in Fig. 8. The performance is acceptable, but worse than the standard EnKF (Fig. 5). This is exactly what is expected. In the standard EnKF we incorporated errors in the dynamics into the parameters, e.g., by estimating a too-large hydraulic conductivity, which can represent

17 18 19 20 21 22 Time [day]

Figure 8. Forecast at the measurement positions of the rst layer during time period C with the results from the closed-eye EnKF.The solid lines show the mean value of a total of 100 different forecasts. The pale colors show the results from 25 of these, with soil parameters, Miller scaling parameters and boundary condition sampled from the distributions estimated with the closed-eye EnKF and the initial condition, which was actually used for the estimation itself. The dashed line shows the measurements.

preferential ow up to a point. By introducing the closed-eye period, we do not incorporate these errors into the parameters any more. The parameters are thus believed to be closer to the believed true material properties. Obviously, however, these parameters fail to describe rain events that lead to hydraulic non-equilibrium.

Guiding the states through the closed-eye period is a challenge. A representation of the dynamics uncertainty would be required to estimate optimal states. We did not do this. The adaptive covariance ination reduces the issue by increasing the ensemble spread when measurements cannot be explained by combined measurement and ensemble uncertainty (Fig. 9). Still, it cannot solve the problem completely, since it was not designed to follow rapid changes of the ination factor and it is merely based on the correlations in the ensemble.These correlations are not necessarily valid any more, when a major, spatially structured uncertainty is not represented.In our case this did not lead to complications, because the deviations are still small and could be adjusted towards the end of the closed-eye period. In cases with stronger rain (and stronger non-equilibrium) and larger deviations of representation and measurements, we expect that a description of the dynamics uncertainties becomes mandatory.

3.4 Prediction capabilities

We investigate the predictive capabilities of the standard EnKF and closed-eye EnKF during the combined time period of B and C as well as time period D. As the closed-eye EnKF added another 10 iterations for the estimation, we compare it to standard EnKF results with a total of 20 iterations as well.

During time period B the total amount of rain was 23.2 mm with a maximal ux of 2.5 mm in 10 min (compared to a total of 18.2 mm with a maximal ux of 0.7 mm in 10 min for

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/

H. H. Bauser et al.: EnKF with closed-eye period 5011

17 18 19 20 21 22

Time [day]

(a)

Standard EnKF

(b)

Closed-eye EnKF

0 5 10 15 20 25 30 35

1.05

Depth [cm]

Depth [cm]

0 5 10 15 20 25 30 35

1.04

1.03

Inflation factor

1.02

1.01

1

17 18 19 20 21 22

Time [day]

Figure 9. Ination factor in the rst layer during the last iteration of the standard EnKF and closed-eye EnKF. The adaptive covariance ination increases the ination factor around a measurement location when the deviations due to the inltration reach the corresponding TDR. The deviations and consequently the ination are stronger for the closed-eye EnKF, as the preferential ow is not incorporated into the parameters any more. The values at the end of an iteration are used as the initial ination factors for the next iteration.

Time period B and C with standard EnKF results

(a)

0.3 Time period B and C with closed-eye EnKF results

(b)

0.3

0.25

0.25

Vol. water content

0.2

Vol. water content

0.2

0.15

0.15

0.1

0.1

5 10 15 20

5 10 15 20

Time [day]

Time [day]

(c)

(d)

0.3 D

Time period D with standard EnKF Results

0.3 D

Time Period D with closed-eye EnKF Results

0.25

0.25

Vol. water content

Vol. water content

0.2

0.2

0.15

0.15

TDR (8.5 cm) TDR (19 cm) TDR (28.5 cm)

TDR (19 cm) TDR (28.5 cm)

0.1

0.1

30 Time [day]

40 50 60

30 Time [day]

40 50 60

Figure 10. Forecast at the measurement positions of the rst layer during time period B and C (a, b) as well as D (c, d). The solid lines show the mean value of a total of 100 different forecasts. The pale colors show the results from 25 of these, with soil parameters and Miller scaling parameters sampled from the distributions estimated with the standard EnKF (a, c) and closed-eye EnKF (b, d). The dashed line shows the measurements, the dotted line the lowest water content measured during time period C, indicating the lower limit of the water content range, for which the parameters are deemed valid.

time period C). As we already assume local non-equilibrium during the rain event in time period C, we expect an even stronger effect during time period B.

The forecast results for time period B and C for soil parameters and Miller scaling parameters from the standard EnKF and closed-eye EnKF are depicted in Fig. 10a and b (there is no improved estimation of the upper boundary condition during time period B and we did not use the improved boundary condition from time period C). We performed a spin-up with an EnKF with state estimation during time period A to achieve a good initial condition. Neither of them is able to predict the water contents during the rain event in time pe-

riod B in detail. However, the standard EnKF describes the rain event in time period C better. The closed-eye EnKF is not able to predict the water contents during rain events as, for both, local non-equilibrium is expected, which is not incorporated into the parameter estimation. For the standard EnKF only the rain event from time period C is incorporated into the parameters and it predicts this rain event better.However, it fails to predict the stronger rain event in time period B showing that the local non-equilibrium in different rain events affects the dynamics differently. Hence, their effects should not be incorporated into parameters, but rather considered separately. The uncertainties are large for both

www.hydrol-earth-syst-sci.net/20/4999/2016/ Hydrol. Earth Syst. Sci., 20, 49995014, 2016

5012 H. H. Bauser et al.: EnKF with closed-eye period

EnKFs. The reason for this is the large uncertainty in the boundary conditions.

We emphasize that the standard EnKF excels if the goal is to guide some partly awed representation with a stream of data. However, if in addition it is the goal to estimate objective parameters for the representation within its range of applicability, then the closed-eye EnKF must be used, at the cost of a lower heuristic performance.

The forecast during the dry time period D is depicted in Fig. 10c and d. Here, we used the nal state of the last iteration in the standard EnKF and closed-eye EnKF estimation as the initial condition. In the observed water content range the closed-eye EnKF performs better than the standard EnKF.This is the result that we expected. We were able to estimate the believed true material properties in the observed water content range and can make predictions in this range as long as the basic assumptions of the Richards equation are met.

Both EnKF can describe the topmost TDR better than the lower TDRs. This is expected as well, since the rst TDR measures the strongest dynamics and hence dominates the parameter estimation. The small deviations in the other TDR sensors show that the assumption of Miller similarity, including the linear interpolation in between, is capable of improving the representation. It cannot represent heterogeneity completely, however.

Below the observed soil water content we do not expect predictive capabilities, since the material properties are heuristic and hence cannot be applied outside of the calibrated range. Nevertheless, we can see that the closed-eye EnKF does describe the topmost TDR acceptably. Again, we see deviations in the two lower TDR sensors, this time even stronger. Again, we attribute this to the limited description of the heterogeneity with Miller scaling.

As a last comment, we did investigate the possibility of improving the results by estimating a multiplicative factor to the uncertain evaporation in analogy to Wollschlger et al. (2009). This factor can improve the estimation of the topmost TDR towards later times, but does not change the qualitative results above.

4 Conclusions

In this study we improved the representation of soil water movement in a 1-D soil prole at the Grenzhof test site close to Heidelberg, Germany with an EnKF based on water content measurements by TDR probes.

We assess key uncertainties of this specic representation.

These are initial conditions, soil hydraulic parameters, Miller scaling factors (describing small-scale heterogeneity), upper boundary condition, and the local equilibrium assumption by the Richards equation. These were accounted for in the process. Other components are deemed to be of lesser importance and were neglected. The most noteworthy neglected

uncertainties are hysteresis, violations of the 1-D assumption and measurement biases, which might affect our results.

We designed a three-stage approach to directly represent and estimate all key uncertainties, except for errors caused by the local hydraulic equilibrium assumption. These intermittent errors are handled by introducing a closed-eye period.Using an iterative approach we can perform the estimation on a single rain event.

In the rst stage, the prior is improved. In our case this is the initial water content distribution and Miller scaling factors. A good prior for the Miller scaling proved particularly benecial, because during phases without dynamics, water content distributions generated by heterogeneity cannot be distinguished from low hydraulic conductivities.

The second stage is the iterative application of a standard EnKF with an augmented state to improve the representation. This state consists of soil water contents, soil hydraulic parameters, Miller scaling factors, and the upper boundary condition. This approach incorporates possible errors in the dynamics into the augmented state components. Since parameters are assumed to be constant in time, we can detect intermittent errors in the dynamics from uctuating parameters.

The third stage is the application of a closed-eye EnKF that only estimates the full augmented states when and where we assume the dynamics to be correctly represented. Outside of this range we close an eye and do not estimate the previously varying parameters. In this way, the state is guided through this ill-represented phase and the full estimation is picked up afterwards again.

Our study showed that assuming diagnostic soil layers to be homogeneous is not sufcient to represent heterogeneity.We were forced to also assume small-scale heterogeneity by observations of consistently higher water contents at probes that are located higher up within a soil layer. This emphasizes the role of heterogeneity, which must be considered more extensively in studies that rely on local measurements like the used TDR measurements.

We employed a simplied representation of heterogeneity by estimating Miller scaling factors at measurement locations and interpolation between them. This captures the main features. Miller scaling is not capable of representing the heterogeneity completely, however. Predictions for the rst layer showed that measurements from the topmost TDR can be predicted better than the measurements of the second and third TDR probes. We attribute this to limitations of the Miller scaling. The parameter estimation is primarily inuenced by the topmost TDR. Describing the other measurements with the same parameter set scaled by the Miller factor cannot fully describe the material properties there.

An iterative standard EnKF was used to successfully estimate soil parameters, Miller scaling parameters and upper boundary condition. It was capable of predicting the rain event that was used for its calibration very well. However, it fails at the prediction of a different rain event. The rea-

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/

H. H. Bauser et al.: EnKF with closed-eye period 5013

son is the violation of the local equilibrium assumption by the Richards equation (Eq. 1) during strong rain events. Non-represented errors are incorporated into the parameters. Since these errors are different for different rain events (with different strengths and different local non-equilibrium), they cannot be predicted. Due to the iterative approach, we could detect the times when the local-equilibrium assumption is violated: the variations of the parameters are larger than the change from initial to nal value during later iterations. The changes in the direction of the parameter update then determine the closed-eye period.

Generally, the closed-eye period can be detected if the operational limits of the model are known. In our case, we base this on the changing parameters. However, a direct detection in the state or forcing could be possible as well.

The closed-eye EnKF omits the incorporation of the model structural errors in the parameters and is a generally applicable concept. In our current study, it yields better predictions during periods when the underlying assumptions are fullled: the drying period after a rain event when there is local equilibrium, showing the strength of the Richards equation there. As a consequence, however, its predictions are worse during rain events, when the local equilibrium assumption is violated.

During the closed-eye phase, a description of these uncertainties caused by the non-equilibrium is desirable to be able to optimally guide the states through this phase. In our case, the method still performed well without such a description, as the errors were small enough to be compensated by the adaptive covariance ination.

Our approach is capable of nding parameters closer to the believed true material properties of the soil than a standard EnKF. Predictions during rain events would require an additional representation of the fast dynamics during the event, though. Still, our approach shows a way to limit the incorporation of errors into parameters and is one step towards the concept of knowledge fusion the consistent aggregation of all information pertinent to some observed reality.

5 Data availability

The underlying measurement data is available at http://ts.iup.uni-heidelberg.de/data/bauser-etal-2016-hess.zip

Web End =http://ts.iup. http://ts.iup.uni-heidelberg.de/data/bauser-etal-2016-hess.zip

Web End =uni-heidelberg.de/data/bauser-etal-2016-hess.zip .

Author contributions. Hannes H. Bauser designed, implemented, performed and analyzed the presented study. Stefan Jaumann provided computational software. Daniel Berg provided discussions on the statistical foundations from a particle lter point of view.Kurt Roth came up with the idea of knowledge fusion and the closed-eye period. All authors participated in continuous discussions. Hannes H. Bauser prepared the manuscript with contributions from all authors.

Acknowledgements. We thank the two anonymous reviewers for their constructive comments, which helped to improve this paper.

This research is funded by Deutsche Forschungsgemeinschaft (DFG) through project RO 1080/12-1.

Edited by: I. NeuweilerReviewed by: two anonymous referees

References

Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: Crop evapotranspiration: Guidelines for computing crop requirements, FAO Irrigation and Drainage Paper No. 56, FAO Food and Agriculture Organization of the United Nations, Rome, Italy, 1998.

Anderson, J. L.: Spatially and temporally varying adaptive covariance ination for ensemble lters, Tellus A, 61, 7283, 2009.Anderson, J. L. and Anderson, S. L.: A Monte Carlo implementation of the nonlinear ltering problem to produce ensemble assimilations and forecasts, Mon. Weather Rev., 127, 27412758, 1999.

Burgers, G., van Leeuwen, P. J., and Evensen, G.: Analysis scheme in the ensemble Kalman lter, Mon. Weather Rev., 126, 1719 1724, 1998.

Carsel, R. F. and Parrish, R. S.: Developing joint probability distributions of soil water retention characteristics, Water Resour.Res., 24, 755769, 1988.

Chen, Y. and Zhang, D.: Data assimilation for transient ow in geologic formations via ensemble Kalman lter, Adv. Water Resour., 29, 11071122, 2006.

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, 421432, 2006.

De Lannoy, G. J., Houser, P. R., Verhoest, N. E., and Pauwels, V. R.:

Adaptive soil moisture prole ltering for horizontal information propagation in the independent column-based CLM2.0, J.Hydrometeorol., 10, 766779, 2009.

Erdal, D., Neuweiler, I., and Wollschlger, U.: Using a bias aware

EnKF to account for unresolved structure in an unsaturated zone model, Water Resour. Res., 50, 132147, 2014.

Erdal, D., Rahman, M., and Neuweiler, I.: The importance of state transformations when using the ensemble Kalman lter for unsaturated ow modeling: Dealing with strong nonlinearities, Adv. Water Resour., 86, 354365, 2015.

Evensen, G.: Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics, J. Geophys. Res.-Oceans, 99, 1014310162, 1994.Foken, T.: Angewandte Meteorologie, Springer, Berlin, Heidelberg,

2006.

Gaspari, G. and Cohn, S. E.: Construction of correlation functions in two and three dimensions, Q. J. Roy. Meteorol. Soc., 125, 723 757, 1999.

Hamill, T. M., Whitaker, J. S., and Snyder, C.: Distance-dependent ltering of background error covariance estimates in an ensemble Kalman lter, Mon. Weather Rev., 129, 27762790, 2001.Hendricks Franssen, H.-J. and Kinzelbach, W.: Real-time ground-water ow modeling with the ensemble Kalman lter: Joint estimation of states and parameters and the l-

www.hydrol-earth-syst-sci.net/20/4999/2016/ Hydrol. Earth Syst. Sci., 20, 49995014, 2016

5014 H. H. Bauser et al.: EnKF with closed-eye period

ter inbreeding problem, Water Resour. Res., 44, W09408, doi:http://dx.doi.org/10.1029/2007WR006505

Web End =10.1029/2007WR006505 http://dx.doi.org/10.1029/2007WR006505

Web End = , 2008.

Hopmans, J. W., imnek, J., Romano, N., and Durner, W.: Simultaneous determination of water transmission and retention properties. Inverse Methods, in: Methods of Soil Analysis. Part 4.Physical Methods, edited by: Dane, J. H. and Topp, G. H., Soil Science Society of America Book Series No. 5., Soil Science Society of America, Madison, WI, 9631008, 2002.Houtekamer, P. L. and Mitchell, H. L.: Data assimilation using an ensemble Kalman lter technique, Mon. Weather Rev., 126, 796 811, 1998.

Houtekamer, P. L. and Mitchell, H. L.: A sequential ensemble Kalman lter for atmospheric data assimilation, Mon. Weather Rev., 129, 123137, 2001.

Huber, E., Hendricks-Franssen, H., Kaiser, H., and Stauffer, F.: The role of prior model calibration on predictions with ensemble Kalman lter, Ground Water, 49, 845858, 2011.

Ippisch, O., Vogel, H.-J., and Bastian, P.: Validity limits for the van GenuchtenMualem model and implications for parameter estimation and numerical simulation, Adv. Water Resour., 29, 17801789, 2006.

Jafarpour, B. and Tarrahi, M.: Assessing the performance of the ensemble Kalman lter for subsurface ow data integration under variogram uncertainty, Water Resour. Res., 47, W05537, doi:http://dx.doi.org/10.1029/2010WR009090

Web End =10.1029/2010WR009090 http://dx.doi.org/10.1029/2010WR009090

Web End = , 2011.

Li, C. and Ren, L.: Estimation of unsaturated soil hydraulic parameters using the ensemble Kalman lter, Vadose Zone J., 10, 1205 1227, 2011.

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

Web End =10.1029/2006WR005756 http://dx.doi.org/10.1029/2006WR005756

Web End = , 2007.

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, 38633887, doi:http://dx.doi.org/10.5194/hess-16-3863-2012

Web End =10.5194/hess-16-3863-2012 http://dx.doi.org/10.5194/hess-16-3863-2012

Web End = , 2012.

Man, J., Li, W., Zeng, L., and Wu, L.: Data assimilation for unsatu-rated ow models with restart adaptive probabilistic collocation based Kalman lter, Adv. Water Resour., 92, 258270, 2016.Miller, E. and Miller, R.: Physical theory for capillary ow phenomena, J. Appl. Phys., 27, 324332, 1956.

Moradkhani, H., Sorooshian, S., Gupta, H. V., and Houser, P. R.:

Dual stateparameter estimation of hydrological models using ensemble Kalman lter, Adv. Water Resour., 28, 135147, 2005.Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513522, 1976.

Nowak, W.: Best unbiased ensemble linearization and the quasi-linear Kalman ensemble generator, Water Resour. Res., 45, W04431, doi:http://dx.doi.org/10.1029/2008WR007328

Web End =10.1029/2008WR007328 http://dx.doi.org/10.1029/2008WR007328

Web End = , 2009.

Reichle, R. H.: Data assimilation methods in the Earth sciences,Adv. Water Resour., 31, 14111418, 2008.

Reichle, R. H., McLaughlin, D. B., and Entekhabi, D.: Hydrologic data assimilation with the ensemble Kalman lter, Mon. Weather Rev., 130, 103114, 2002.

Richards, L. A.: Capillary conduction of liquids in porous mediums,Physics, 1, 318333, 1931.

Roth, K., Schulin, R., Flhler, H., and Attinger, W.: Calibration of time domain reectometry for water content measurement using a composite dielectric approach, Water Resour. Res., 26, 2267 2273, 1990.

Schenk, G. U.: Soil-Atmosphere coupling: Looking from within the soil towards the surface, PhD thesis, Ruperto-Carola University Heidelberg, Heidelberg, Germany, http://www.ub.uni-heidelberg.de/archiv/12840

Web End =http://www.ub. http://www.ub.uni-heidelberg.de/archiv/12840

Web End =uni-heidelberg.de/archiv/12840 (last access: 7 December 2016), 2011.

Shi, L., Song, X., Tong, J., Zhu, Y., and Zhang, Q.: Impacts of different types of measurements on estimating unsaturated ow parameters, J. Hydrol., 524, 549561, 2015.

Song, X., Shi, L., Ye, M., Yang, J., and Navon, I. M.: Numerical comparison of iterative ensemble Kalman lters for unsaturated ow inverse modeling, Vadose Zone J., 13, doi:http://dx.doi.org/10.2136/vzj2013.05.0083

Web End =10.2136/vzj2013.05.0083 http://dx.doi.org/10.2136/vzj2013.05.0083

Web End = , 2014.

Van Genuchten, M. T.: A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am.J., 44, 892898, 1980.

Vrugt, J. A., Diks, C. G., Gupta, H. V., Bouten, W., and Verstraten, J. M.: Improved treatment of uncertainty in hydro-logic modeling: Combining the strengths of global optimization and data assimilation, Water Resour. Res., 41, W01017, doi:http://dx.doi.org/10.1029/2004WR003059

Web End =10.1029/2004WR003059 http://dx.doi.org/10.1029/2004WR003059

Web End = , 2005.

Vrugt, J. A., Stauffer, P. H., Whling, T., Robinson, B. A., and Vesselinov, V. V.: Inverse modeling of subsurface ow and transport properties: A review with new developments, Vadose Zone J., 7, 843864, 2008.

Wollschlger, U., Pfaff, T., and Roth, K.: Field-scale apparent hydraulic parameterisation obtained from TDR time series and inverse modelling, Hydrol. Earth Syst. Sci., 13, 19531966, doi:http://dx.doi.org/10.5194/hess-13-1953-2009

Web End =10.5194/hess-13-1953-2009 http://dx.doi.org/10.5194/hess-13-1953-2009

Web End = , 2009.

Wu, C.-C. and Margulis, S. A.: Feasibility of real-time soil state and ux characterization for wastewater reuse using an embedded sensor network data assimilation approach, J. Hydrol., 399, 313325, 2011.

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

Web End =10.2136/vzj2011.0176 http://dx.doi.org/10.2136/vzj2011.0176

Web End = , 2013. Zhang, D., Madsen, H., Ridler, M. E., Refsgaard, J. C., and Jensen,K. H.: Impact of uncertainty description on assimilating hydraulic head in the MIKE SHE distributed hydrological model, Adv. Water Resour., 86, 400413, 2015.

Hydrol. Earth Syst. Sci., 20, 49995014, 2016 www.hydrol-earth-syst-sci.net/20/4999/2016/