3.1 The CATHY model

The CATHY model (Camporese et al., 2010) is a physics-based hydrological model capable of simulating integrated subsurface, overland and channel water flow. The model combines a Richards equation solver for the three-dimensional flow in variably saturated porous media with a surface water flow module for the solution of the one-dimensional diffusion wave approximation of the de Saint-Venant equation. In this study, however, there is no surface runoff. Therefore, only subsurface flow is considered, according to the following form of the Richards equation:

In Eq. (1), $S_{\mathrm{w}}=\mathit{\theta }/\mathit{\varphi }$ is water saturation, θ and ϕ being the volumetric soil water content and porosity (∕), respectively, S_s is the specific storage coefficient (L⁻¹), ψ is the pressure head (L), t is time (T), ∇ is the gradient operator, K_s is the saturated hydraulic conductivity tensor (L∕T), K_r is the relative hydraulic conductivity function (∕), ${\mathit{\eta }}_z=(\mathrm{0},\mathrm{0},\mathrm{1}{)}^T$ is the vertical direction vector, z is the vertical coordinate directed upward (L), and q_s represents distributed source or sink terms (L³∕L³T).

The unsaturated hydraulic properties are taken into account by means of the van Genuchten functions (e.g., Wösten and van Genuchten, 1988) S_w(ψ) and K_r(ψ):

where $S_{\mathrm{wr}}={\mathit{\theta }}_{\mathrm{r}}/\mathit{\varphi }$ is the residual water saturation, with θ_r the residual water content, α is an empirical constant (L⁻¹) related to the inverse of the air entry suction, and the dimensionless shape parameters n and m are linked by the expression $m=\mathrm{1}-\mathrm{1}/n$. The model solves Eq. (1) by means of Galerkin finite elements with tetrahedral elements and linear basis functions in space and weighted finite differences for integration in time (Camporese et al., 2010).

It is worth noting that the Richards equation is strongly nonlinear, due to the retention curves (Eqs. 2 and 3). Such nonlinearities are enhanced in this study by the fact that the hillslope is characterized by dominant unsaturated conditions (Fig. 2), as opposed to many other previous applications of data assimilation with physics-based models. This makes our hillslope experiment particularly challenging and distinctive from both modeling and data assimilation perspectives.

3.2 The ensemble Kalman filter

The ensemble Kalman filter (Evensen, 2003, 2009a, b) is a sequential data assimilation scheme, in which states (and parameters) are sequentially updated based on a Monte Carlo approximation of the covariance matrices needed in the standard Kalman filter. The process is Markovian of the first order and the implementation of the EnKF does not require the linearization of the model, making it particularly suitable for handling nonlinear problems.

In this paper the EnKF is implemented according to the numerical formulation proposed by Sakov et al. (2010). Let X be an ensemble matrix of M rows and N columns, where N is the number of realizations and M is the state dimension, i.e., the number of nodes in the finite element grid, augmented by the number of parameters that are subject to update. The main idea behind this type of implementation is that the matrix X can be defined as the sum, x+A, of the ensemble average, x,

and the matrix of ensemble anomalies, A,

where 1 and I are a vector with all elements equal to one and the identity matrix, respectively. As usual, superscript T denotes matrix transposition.

Whenever observed data are available, the EnKF can compute the updated matrix X^u as the sum of the updated ensemble mean, x^u, and the updated anomalies, A^u:

In the following, the lack of a superscript u denotes that the matrix or vector is computed at the forecast stage, i.e., at the previous model time step.

The updated mean can be calculated as

where β is a diagonal matrix of dampening factors (Hendricks Franssen and Kinzelbach, 2008), whose elements vary from 0 to 1, and s is the scaled innovation vector

which depends on the measurement error covariance matrix, R, the difference between the measurements, D, and the ensemble mean of the simulated observations, Hx. Dampening factors help to prevent the occurrence of unstable updates and filter divergence (e.g., Evensen, 2009b) and were chosen as an alternative to covariance inflation due to ease of implementation.

The matrix G is defined as

where S is the matrix of scaled ensemble innovation anomalies

HA being the simulated measurement anomalies

and M being defined as

The updated anomalies, A^u, are computed as

When updating the states only, the elements of X are the pressure heads at each node of the finite element grid, while the state augmentation technique is used when also updating the parameters. In this latter case, the desired parameters (e.g., hydraulic conductivity, parameters of the retention curves), transformed as described in Sect. 4.2, are added to X and updated based on their correlation with the system states (Erdal et al., 2015).

Table 1Perturbation parameters for the generation of the ensemble initial conditions, hydraulic conductivities and atmospheric forcing.

^a Rainfall (positive) and evaporation (negative) rates. ^b σ_VG as described in Sect. 4.2.

Download Print Version | Download XLSX

4.1 CATHY setup

The artificial hillslope is discretized with a surface triangular grid resulting from the subdivision of square cells of 10 cm side. The triangular grid is then replicated vertically for a total of 25 layers to generate the three-dimensional tetrahedral mesh (Fig. 3). Fifteen layers are used to represent the top sand, while 10 layers discretize the clay. No flow boundary conditions are assumed at each boundary, except for the subsurface outflow section, where seepage face boundary conditions are used, and the surface, where time-variable rainfall or evaporation rates are imposed. Finally, the soil hydraulic parameters are assigned as reported in Table 1. The values of saturated hydraulic conductivity and van Genuchten retention parameters are perturbed to generate the ensemble of realizations as described in the following section, whereas soil porosities and specific storages are considered deterministic, the former being well characterized with laboratory tests and the latter having little impact on the CATHY model response in this case.

Figure 3Three-dimensional finite element grid of the hillslope.

Download

4.2 EnKF setup

In order to generate the ensemble of realizations needed for the application of the EnKF, we perturb the atmospheric forcing (i.e., rainfall and evaporation rates), soil properties and initial conditions. Table 1 reports a summary of the perturbed variables, along with their nominal mean values as well as the nature and statistics of the perturbations. The uncertainties on model parameters and boundary conditions have been assigned on the basis of previous modeling experiences and preliminary characterization of the soils in the hillslope (Camporese et al., 2009a, b; Lora et al., 2016b).

The ensemble of time-variable atmospheric forcing rates was generated with a sequence of multiplicative perturbations, q_k, correlated in time as in Evensen (2003):

where the subscript k is the time index, w_k is a sequence of white noise drawn from the standard normal distribution, and the coefficient γ is computed as

Δt being the assimilation interval and τ the specified time decorrelation length, here set equal to 108 000 s, i.e., 30 h.

The initial conditions consist of a uniform value of pressure head, ψ₀, whose nominal ensemble mean is −0.67 m, based on the average provided by the tensiometer measurements. We opted for a uniform value of initial pressure head because, at the beginning of the experiment, the entire hillslope was in a highly unsaturated condition, as derived from the tensiometer data. This means that, even without an initially hydrostatic (i.e., equilibrium) profile, water is basically prevented from flowing in or out of the hillslope due to the very small values of relative hydraulic conductivity. This is one of the peculiarities of our test case compared to similar studies. Preliminary analyses showed that the model is not very sensitive to ψ₀, because in any case the initial hydraulic conductivity is so small that the hillslope responds only when the first rainfall event occurs and starts wetting the soil. Therefore, in this study, no warm-up or spin-up was necessary. The ensemble of ψ₀ values is generated by additive perturbations normally distributed with mean equal to 0 and standard deviation equal to 0.2 m.

Perturbed soil parameters, for both sand and clay, include the saturated hydraulic conductivity as well as the parameters of the van Genuchten retention curves. Table 1 reports the nominal mean values of K_s, based on soil samples analyzed in the laboratory, as well as the parameters used for generating their ensembles. The saturated hydraulic conductivities of sand and clay are perturbed independently of each other with multiplicative perturbations sampled from a lognormal distribution.

The parameters of the van Genuchten retention curves α, n, and θ_r are perturbed taking into account their mutual correlation according to Carsel and Parrish (1988), who described their statistics and transformed them into normally distributed variables via the Johnson system (Johnson, 1970). Three main transformation functions are available, i.e, the lognormal (LN), log-ratio (SB), and hyperbolic (SU):

where V denotes the parameter before transformation, bounded within the range [A B], and Y denotes the transformed parameter with normal distribution. In this work, the prior statistics of the van Genuchten parameters are taken from the soil types “sandy loam” and “silt loam” in Carsel and Parrish (1988), assumed as valid representations of our sand and clay, respectively. In particular, the SB transformation has been used for all the parameters, except for n of the sandy loam and α of the silt loam, for which the LN transformation has been applied. The ensemble of transformed parameters is generated by

where y is the vector containing transformed α, n, and θ_r, u contains the transformed variable means, T is an upper diagonal matrix, obtained from the Choleski factorization of the covariance matrix of the three parameters (see Table 2), and z is a vector of normal deviates with mean equal to 0 and standard deviation σ_VG. The transformed variables can be included in the matrix X when updating the parameters, and can then be back-transformed in order to obtain the updated values of the soil retention parameters.

The EnKF algorithm implemented here is actually an ensemble transform Kalman filter (Bishop et al., 2001) that does not require the perturbation of observations. On the other hand, the measurement error covariance matrix, R, must be assumed to be known a priori. In this work, R was estimated directly from the measurements, taking advantage of the high time resolution of the collected data. Pressure head and water content data were collected every 2 s and averaged every 10 min over a 40 s window to obtain the observations to assimilate. Over the same time window, the data were linearly detrended and the residuals were used to calculate the correlation coefficients between all pairs of observations. The final covariance matrices were then assembled by multiplying the correlation coefficients by the relevant standard deviations, assumed as 0.05 m and 0.025 for pressure head and water content, respectively. The subsurface outflow measurements are assumed to be independent of the pressure head and water content data, with a standard deviation equal to 8 % of the measured discharge. All measurement error standard deviations have been estimated based on the accuracy of the sensors and plausible positioning errors.

When assimilating multiple variables, proper normalization of the measurement error covariance matrices, anomalies of the simulated data, and innovation vectors were performed, using values of 0.6 m, 0.58, and $\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{5}}$ m³ s⁻¹ for pressure head, water content and subsurface outflow, respectively. The normalization ensures that in multivariate assimilation scenarios the covariance matrices in the Kalman gain are not ill-conditioned (Camporese et al., 2009b; Evensen, 2003).

Table 2Factored covariance matrices used for the perturbation of van Genuchten parameters (Carsel and Parrish, 1988).

Download Print Version | Download XLSX

4.3 Data assimilation scenarios

A total of 17 data assimilation scenarios have been simulated, whereby the assimilation interval, the assimilated variables, the updated variables, and the uncertainty on the van Genuchten parameters were varied. Table 3 reports a summary of the main characteristics for each scenario. Assimilated variables may include water content only or with subsurface outflow, pressure head only or with subsurface outflow, and all three variables together. With regard to the updated variables, three cases have been analyzed: update of the state variables only; update of the state variables and saturated hydraulic conductivities for both sand and clay; update of the state variables, hydraulic conductivities and van Genuchten parameters, for both sand and clay. In all the data assimilation scenarios, observations were assimilated only during the first 5 days of simulation, leaving the final 7 days as a validation period, during which the ensemble was left to evolve freely. For comparison, two open loop simulations, i.e., without data assimilation, have also been carried out (Table 3).

Table 3Overview of the open loop and data assimilation scenarios.

¹ OL1 and OL2 indicate open loop scenarios, i.e., simulations without data
assimilation. ² WC, PH, and Q denote water content, pressure head, and subsurface
outflow, respectively. ³ St. var., K_s, and V.G. indicate state variables (in terms of
pressure head), saturated hydraulic conductivity, and van Genuchten parameters,
respectively.

Download Print Version | Download XLSX

The performance of the simulations has been evaluated by means of the root mean square error (RMSE), computed for the different variables, i.e., pressure head, water content, and subsurface outflow. The root mean square error is calculated as

where N_o is the number of observations available (six, four, and one for water content, pressure head, and subsurface outflow, respectively), S_i,j refers to the simulated results of the ith realization of the ensemble at the location of the jth observation and O_j is the corresponding experimental value. To obtain a meaningful comparison between the errors of different variables, we also compute the time-averaged normalized root mean square error, NRMSE,

where N_T is the number of time steps, k is the time index, and NF is a normalization factor equal to 0.58, 0.60 m, and $\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{5}}$ m³ s⁻¹ (i.e., 2.5 L min⁻¹) for the water content, pressure head, and subsurface discharge, respectively. The NRMSE is computed separately for each variable and for the assimilation (N_T=120) and validation (N_T=167) periods, but also as a global index of performance averaged over all the variables and the two periods.

5.1 Preliminary simulations

A preliminary sensitivity analysis over a number of EnKF parameters has been performed, in order to select a final and satisfactory setup for the subsequent data assimilation scenarios. First, simulations with N equal to 32, 128 and 256 have been performed, and it has been found that an ensemble size of 128 ensures a good tradeoff between performance and computational effort. Some preliminary sensitivity analyses of the value of σ_VG have also been performed, resulting in values of 0.1 and 0.25 for the scenarios with and without update of van Genuchten parameters, respectively (Table 3).

Then, several dampening factor values (β in Eqs. 7 and 13) were tested, including combinations of different values for the update of system state and parameters. Based on this analysis, a value of 1 has been chosen for the update of the system state, whereas a dampening factor equal to 0.5 has been selected for the update of the soil hydraulic parameters (K_s, α, n, and θ_r, for both sand and clay). This choice of dampening factors is consistent with previous studies (Brandhorst et al., 2017) and prevents abrupt changes in the retention curve parameters that could lead to difficulties in model convergence and hence loss of realizations.

According to these preliminary analyses, all scenarios reported in Table 3 have been simulated with an ensemble size of 128 and dampening factors of 1 and 0.5 for system state and parameters, respectively.

5.2 Overall EnKF performance

Table 4 and Fig. 4 summarize the performance of the EnKF in all the data assimilation scenarios, expressed in terms of NRMSE for the three measured variables (water content, WC, pressure head, PH, and subsurface outflow, Q), and averaged separately over the assimilation and validation windows. A comparison between scenarios 1–5 and 8–12 in Table 4, characterized by the same assimilated and updated variables but different assimilation intervals, shows that assimilating more frequently does not always result in significant improvements of model predictions. The variable that benefits the most from more frequent updates is subsurface outflow, especially in the scenarios where Q is assimilated (e.g., compare scenarios 3–5 with 10–12).

Figure 4Normalized root mean square errors for water content, pressure head, and subsurface outflow of the data assimilation scenarios versus corresponding values of the open loop. Symbols in magenta represent values calculated over the assimilation period, while symbols in blue represent values calculated over the validation period.

Download

Figure 4 highlights that most of the data assimilation scenarios result in an improvement of model predictions for pressure head and subsurface outflow, compared to the open loop simulations (data pairs below the 45 ^∘ reference line). However, in some scenarios, the filter performance in predicting the water content is actually worse than in the open loop. A close inspection of the values in Table 4 indicates that such scenarios are those where pressure head is assimilated, alone or in conjunction with water content and subsurface outflow. This is likely due to a combination of two factors: (i) only four out of six tensiometers are available for assimilation, compared to the six WCR probes available for water content, and (ii) the pressure head measurements are characterized by relatively poor quality. This can be appreciated from the pressure head data shown in Fig. 2, where diurnal disturbances caused by temperature fluctuations (Warrick et al., 1998) are apparent.

Table 4Normalized root mean square errors (NRMSEs) for the 17 data assimilation scenarios under analysis and two open loop (OL) simulations. The table reports the NRMSE for three variables, water content, pressure head and outflow discharge, for both the assimilation and the validation periods. The last three columns report the mean values calculated over the three variables (WC, PH and Q) and the global average between assimilation and validation. Bold values indicate exceedance of the corresponding open loop errors.

Download Print Version | Download XLSX

Figure 5Ratios between RMSEs in scenarios S6, S8, and S13 (a, c, e), S7, S10, and S15 (b, d, f), and the corresponding open loop values, for water content, WC (a, b), pressure head, PH (c, d), and subsurface outflow, Q (e, f). The magenta bars refer to the data assimilation period, and the blue ones to the validation period.

Download

5.3 Parameter estimation capabilities

To assess the capabilities and benefits of parameter estimation with the EnKF, it is useful to compare scenarios with the same assimilated variables but different updated variables. Figure 5a, c, and e show the ratios between RMSE in data assimilation scenarios S6, S8, and S13 and the corresponding open loop values for water content, pressure head, and subsurface outflow, respectively. In these three scenarios, water content alone is assimilated, but the updated variables are system state only in S6, system state and saturated hydraulic conductivity in S8, and system state, the K_s and van Genuchten parameters α, n, and θ_r in S13. Figure 5a highlights that progressively updating more parameters brings significant improvements in water content prediction over the validation period, with reductions of the NRMSE with respect to the open loop of almost 20 % when updating K_s only and 60 % when also updating the retention curve parameters. Moreover, updating parameters improves significantly pressure head predictions in validation and subsurface outflow in both assimilation and validation, as shown in Fig. 5c and e.

Figure 6Plots of the pressure head in P2 (a, b, c), water content in W6 (d, e, f), and outflow discharge (g, h, i) for scenarios S6 (a, d, g), S8 (b, e, h), and S13 (c, f, i). Solid green lines represent experimental data, while solid black and magenta lines indicate the ensemble mean of the open loop and data assimilation scenarios, respectively. The simulated 90 % confidence bands are also reported (dashed black and magenta lines). The locations of the tensiometer P2 and water content probe W6 are shown in Fig. 1.

Download

The effect of parameter updating on model predictions for scenarios S6, S8, and S13 can be visualized in Fig. 6, which shows ensemble means and 90 % confidence bands of simulated water content, pressure head, and subsurface outflow in comparison with the experimental data and the corresponding open loop values. In scenario S6, without parameter update, model predictions during the validation tend to converge again to the open loop simulations in terms of both mean and uncertainty (Fig. 6a, d, g). Updating the parameters (scenarios S8 and S13) results in decreased uncertainty during validation (Fig. 6b, c, e, f, h, i), due to the reduced variability of saturated hydraulic conductivity and van Genuchten parameters. Note also that the update of K_s is particularly beneficial to pressure head (panel a versus b and c, where it can be seen that, in the validation phase, the ensemble mean with data assimilation departs the open loop and gets closer to the measurements) and subsurface outflow (panel g versus h and i), whereas updating α, n, and θ_r improved significantly the water content (panel e versus f).

Figure 7Relative frequency distributions of the saturated hydraulic conductivity, K_s, of sand and clay in scenario S10. Graphs in magenta denote the prior K_s distributions, while the blue ones indicate the posterior distributions, i.e., at the end of the assimilation period.

Download

Figure 8Relative frequency distributions of the saturated hydraulic conductivity, K_s, and van Genucthen parameters of sand and clay in scenario S15. Graphs in magenta denote prior distributions, while the blue ones indicate posterior distributions, i.e., at the end of the assimilation period.

Download

5.4 Tradeoffs in multi-source data assimilation

We now focus our attention on the scenarios where multi-source data are assimilated. The right panels in Fig. 5 show the ratios between RMSE in data assimilation scenarios S7, S10, and S15 and the corresponding open loop values for water content, pressure head, and subsurface outflow. Whereas in scenarios S6, S8, and S13 (see panels (a), (c) and (e) in Fig. 5) water content alone was assimilated, in scenarios S7, S10, and S15 (see panels (b), (d) and (f) in Fig. 5) water content and subsurface outflow were jointly assimilated. The updated variables are system state only in S7, system state and saturated hydraulic conductivity in S10, and system state plus K_s and van Genuchten parameters α, n, and θ_r in S15. As previously noted for the scenarios with assimilation of water content alone, the effect of parameter updating is significant mainly in the validation phase. However, one can now observe an increase in validation RMSE of water content when also updating K_s, with a value that exceeds the one of the corresponding open loop. At the same time, there is a decrease in pressure head RMSE of more than 40 % with respect to the open loop, indicating that the update of K_s is beneficial to pressure head, but not to water content. Including the update of van Genuchten parameters, in scenario S15, has an opposite effect: model predictions of water content in validation improve dramatically (RMSE of almost 60 % less than in the open loop), whereas pressure head predictions worsen slightly and align with RMSE values of scenario S13. This indicates that the update of the van Genuchten parameters is more important for predictions of water content than pressure head, which represents an interesting example of the kinds of tradeoffs associated with multi-source data assimilation and parameter updating in integrated hydrological models.

Figure 9Time evolution of the saturated hydraulic conductivity and van Genuchten parameters (mean values, solid line, together with minimum and maximum values, in dashed lines, to indicate the ensemble spread) for the two types of soil, sand and clay, in scenario S17.

Download

Further insights into the differences between scenarios S10 and S15 can be gained from Figs. 7 and 8, which report the prior and posterior (i.e., at the end of the assimilation period) distributions of soil parameters. Fig. 7 reports the results for scenario S10, where K_s only was updated, and shows that the sand K_s is clearly identifiable, whereas the clay K_s is not, with a large residual uncertainty and a mean value that is not consistent with the actual soil type in the hillslope. As no data are available in the clay layer, the large residual variability should be expected, but the bias in the mean value is probably caused by spurious correlations and, perhaps, also by the low sensitivity of the assimilated variables to the clay permeability. This is likely the reason for the poor model prediction of water content.

Figure 8 shows the prior and posterior distributions of K_s, α, n, and θ_r in S15. Again, for the sand, the saturated hydraulic conductivity can be clearly identified, as well as the exponent n of the van Genuchten retention function, while parameters θ_r and α are more difficult to estimate. As for the clay, posterior uncertainty is large for all the parameters and the mean value of K_s shows again a bias with respect to the prior value, although in this case the final value is more consistent with the actual soil type and this could explain why the water content model predictions improve significantly compared to scenario S10. This can also be explained by the fact that all the sensors are located in the sand layer, while no experimental data from the clay layer are assimilated.

Figure 10Plots of the pressure head in P2, (a, b, c), water content in W6, (d, e, f) and outflow discharge (g, h, i) for scenarios S15 (a, d, g), S16 (b, e, h), and S17 (c, f, i). Solid green lines represent experimental data, while solid black and magenta lines indicate the ensemble mean of the open loop and data assimilation scenarios, respectively. The simulated 90 % confidence bands are also reported (dashed black and magenta lines). The locations of the tensiometer P2 and water content probe W6 are shown in Fig. 1.

Download

An additional perspective on parameter estimation is given by Fig. 9, which shows the time evolution of ensemble mean and spread of the saturated hydraulic conductivity and van Genuchten parameters for both sand and clay in scenario S17 (one of the scenarios reported in Fig. 10). The convergence toward stable values in the data assimilation phase is only achieved for the parameters K_s and n of the sand layer, while the dispersion remains generally higher for θ_r (to which the model is typically not very sensitive) and α, as well as for all the parameters of the clay layer.

In summary, the results of parameter updating for the clay show that data would be needed in all the soil layers. However, when dealing with large heterogeneous structures, it is very expensive to have every soil zone properly probed, which is why it is important to assess whether multivariate data assimilation approaches are capable of compensating for the lack of distributed observations with alternative sources of information. Here, an integrated measurement such as the subsurface outflow does not seem to be sufficient to compensate for this lack of representativeness.

Finally, we analyze the tradeoffs in system state predictions associated with multi-source data assimilation for scenarios S15, S16, and S17. Figure 10 shows pressure head in P2, water content in W6, and subsurface outflow as simulated in the data assimilation and open loop scenarios, in terms of ensemble mean and 90 % confidence bands, compared to the measurements. In scenario S15 (Fig. 10a, d, g), where water content and subsurface discharge were assimilated, model results are very good for these variables but not so for pressure head (see also NRMSE values in Table 4). On the other hand, in scenario S16, where pressure head and subsurface outflow were assimilated, pressure head and discharge are well simulated, but not water content (Fig. 10b, e, h, and Table 4). Finally, in scenario S17, where all the available data were assimilated, the model predicts well both pressure head and water content, but at the cost of a slightly degraded prediction of subsurface outflow compared to scenarios S15 and S16 (Fig. 10c, f, i, and Table 4).

Similar issues were reported by Zhang et al. (2016), who found that the joint assimilation of soil moisture and groundwater head does not improve model predictions, unless update localization is used. In their study, this was likely caused by unrealistic cross-variable correlations due to limited ensemble sizes, whereas in our case it might be that the relatively poor quality of pressure head measurements, compared to water content observations, and the lack of observations in the clay layer do not allow us to obtain accurate estimates of the van Genuchten parameters.