Estimating geostatistical parameters and spatially-variable hydraulic conductivity within a catchment system using an ensemble smoother

Groundwater flow models are important tools in assessing baseline conditions and investigating management alternatives in groundwater systems. The usefulness of these models, however, is often hindered by insufficient knowledge regarding the magnitude and spatial distribution of the spatially-distributed parameters, such as hydraulic conductivity (K), that govern the response of these models. Proposed parameter estimation methods frequently are demonstrated using simplified aquifer representations, when in reality the groundwater regime in a given watershed is influenced by strongly-coupled surface-subsurface processes. Furthermore, parameter estimation methodologies that rely on a geostatistical structure of K often assume the parameter values of the geostatistical model as known or estimate these values from limited data. In this study, we investigate the use of a data assimilation algorithm, the Ensemble Smoother, to provide enhanced estimates ofK within a catchment system using the fullycoupled, surface-subsurface flow model CATHY. Both water table elevation and streamflow data are assimilated to condition the spatial distribution of K. An iterative procedure using the ES update routine, in which geostatistical parameter values defining the true spatial structure of K are identified, is also presented. In this procedure, parameter values are inferred from the updated ensemble of K fields and used in the subsequent iteration to generate the K nsemble, with the process proceeding until parameter values are converged upon. The parameter estimation scheme is demonstrated via a synthetic three-dimensional tilted v-shaped catchment system incorporating stream flow and variably-saturated subsurface flow, with spatio-temporal variability in forcing terms. Results indicate that the method is successful in providing improved estimates of theK field, and that the iterative scheme can be used to identify the geostatistical parameter values of the aquifer system. In general, water table data have a much greater ability than streamflow data to condition K. Future research includes applying the methodology to an actual regional study site.

In this study, we investigate the use of a data assimilation algorithm, the Ensemble Smoother, to provide enhanced estimates of K within a catchment system using the fullycoupled, surface-subsurface flow model CATHY.Both water table elevation and streamflow data are assimilated to condition the spatial distribution of K.An iterative procedure using the ES update routine, in which geostatistical parameter values defining the true spatial structure of K are identified, is also presented.In this procedure, parameter values are inferred from the updated ensemble of K fields and used in the subsequent iteration to generate the K ensemble, with the process proceeding until parameter values are converged upon.The parameter estimation scheme is demonstrated via a synthetic three-dimensional tilted v-shaped catchment system incorporating stream flow and variably-saturated subsurface flow, with spatio-temporal variability in forcing terms.Results indicate that the method is successful in providing improved estimates of the K field, and that the iterative scheme can be used to identify the geostatistical parameter values of the aquifer system.In general, water table data have a much greater ability than streamflow data to condition K. Future research includes applying the methodology to an actual regional study site.

Inverse modeling in groundwater applications
Hydrologic models are important tools in assessing baseline conditions and investigating best-management practices in groundwater and catchment-scale systems.Before reliable hydrologic assessments can be made, however, parameter values that drive the response of the model must be appropriately chosen for a specific aquifer or catchment.Direct measurements of hydrologic parameters, however, are scarce and fraught with uncertainty, and typically only apply locally due to the spatial variability of parameter values.
To address this problem of parameter uncertainty, hydrologic models can be used in applications "opposite" or "inverse" to their original use, i.e., parameter values are treated as system unknowns and are determined by extracting Published by Copernicus Publications on behalf of the European Geosciences Union.R. T. Bailey and D. Ba ù: Estimating geostatistical parameters information from observations of system-response variables (Kitanidis and Vomvoris, 1983).The general approach consists of determining the set of parameter values that yields adequate matches between model results and observations from the true hydrologic system.The treatment of parameter values as unknowns that need to be identified constitutes the inverse problem of groundwater modeling (Kitanidis and Vomvoris, 1983), and in most cases must be incorporated in the modeling process (Carrera et al., 2005).
In recent decades numerous methodologies have been proposed and applied to the inverse modeling problem in groundwater modeling, with the general aim to estimate the spatial distribution of hydraulic conductivity (K) or transmissivity (T ) in an aquifer system.An excellent review of early inverse methods is provided by Carrera and Neuman (1986).A review of more recently-proposed methods is given by Carrera et al. (2005).Broadly, parameter estimation is accomplished either through (i) optimization procedures, in which an objective function is defined (typically minimizing the error between model results and measurements) and minimized in a least-squares approach, and (ii) statistical conditioning, in which covariance between the parameters and system-response variables is utilized to condition the parameter values using measurement information.It should be noted that conditioning methods also incorporate a sense of optimization, although the optimization occurs in the derivation of the conditioning algorithm, e.g., through minimizing the trace of the a posteriori error estimate covariance matrix (e.g., Kalman, 1960).
Proposed methodologies are demonstrated typically using simplified hydrologic systems.For applications to groundwater systems, the majority of methodologies are demonstrated using two-dimensional (2-D) confined groundwater flow models (e.g., Gailey et al., 1991;Hantush and Mariño, 1997;Hendricks Franssen et al., 1999;Gómez-Hernández et al., 2003;Drécourt et al., 2006;Hendricks Franssen and Kinzelbach, 2008;Fu and Gómez-Hernández, 2009;Bailey and Baù, 2010).Several studies have employed threedimensional steady-state flow models (Chen and Zhang, 2006;Liu et al., 2008), and several have estimated hydraulic parameters in variably-saturated flow conditions (Yeh and Zhang, 1996;Zhang and Yeh, 1997;Li and Yeh, 1999), although for the latter applications were limited to small 2-D vertical-plane systems.In general, however, critical components of hydrology in watershed systems, e.g., infiltration and percolation in variably-saturated porous media, ponding and overland flow, and stream channel flow have been neglected.Catchment models such as CATHY (CATchment HYdrology), based on the 3-D Richards equation for variably-saturated porous media and a diffusion wave approximation for overland and channel flow, have been used in data assimilation studies (Camporese et al., 2009(Camporese et al., , 2010)), but not yet in parameter estimation.Estimation of parameters in land-surface models has been performed (e.g., Boulet et al., 2002;Xie and Zhang, 2010), although the models treat groundwater flow using simplified approaches.
In recognition that improved parameter estimation occurs when system-response data from more than one governing equation is used (Gailey et al., 1991), with the implication that each data type contains unique information regarding the parameter, numerous studies have employed two or more sets of dissimilar data to condition the parameter values.Such data sets typically include hydraulic head data as well as another data type such as solute concentration data (Gailey et al., 1991;Li and Yeh, 1999;Hendricks Franssen et al., 2003;Gómez-Hernández et al., 2003;Liu et al., 2008), groundwater temperature (Woodbury and Smith, 1988), groundwater travel time (Fu and Gómez-Hernández, 2009), groundwater discharge to surface water (Bailey and Baù, 2010), and tracer breakthrough data at observation wells (Wen et al., 2002).Streamflow data, which carries information regarding the spatial structure of aquifer K due to groundwater-surface water interactions, has been used in data assimilation to improve model performance (Schreider et al., 2001;Aubert et al., 2003;Clark et al., 2008;Camporese et al., 2009Camporese et al., , 2010)), although as yet has not been used to condition K.

Kalman Filter methods
In Kalman Filtering methods, a priori information, i.e., model parameters and associated model results, are merged with observation data from the true system to produce an a posteriori system estimate honoring the true system data at observation points, while still incorporating physically-based information from the numerical model.The resulting algorithm is used to merge model and measurement data whenever measurement data become available during the course of the model simulation.
In contrast to a filter, which assimilates data sequentially as they become available, a smoother incorporates all past model and measurement information in a single assimilation step.The EnKF, EnKS, and ES all use an ensemble of realizations to represent numerically the measurement error statistics (Evensen, 2003), and are designed for large, nonlinear systems.The EKF, EnKF, EnKS, and ES have all been used in hydrologic modeling applications in both systemresponse updating (e.g., Schreider et al., 2001;Aubert et al., 2003;Dunne and Entakhabi, 2005;Clark et al., 2008;Durand et al., 2008;Camporese et al., 2009Camporese et al., , 2010) ) and system parameter conditioning (Hantush and Mariño, 1997;Boulet et al., 2002;Chen and Zhang, 2006;Hendricks Franssen and Kinzelbach, 2008;Liu et al., 2008;Bailey and Baù, 2010;Xie and Zhang, 2010).Application of the EnKF and ES to highly nonlinear hydrologic systems such as a land surface model (Dunne and Entakhabi, 2005) and a coupled surface and variably-saturated subsurface flow model (Camporese et al., 2009) has proven successful.

Geostatistics in parameter estimation
Many parameter estimation studies employ geostatistical models (GMs) to define the a priori estimate of the spatial distribution of log-K or log-T (e.g., Kitanidis and Vomvoris, 1983;Hantush and Mariño, 1997;Chen and Zhang, 2006;Hendricks Franssen and Kinzelbach, 2008), under the assumption that aquifer K in regional systems can generally be described using such models (Kitanidis and Vomvoris, 1983;Hoeksema and Kitanidis, 1985;Carrera et al., 2005).The values of the parameter (e.g., log-K mean, log-K variance, correlation length) that characterize these GMs often have a strong influence on the response of a groundwater model and parameter estimation results (Jafarpour and Tarrahi, 2011), and yet in practice are estimated from limited geologic information and hence are not known with a high degree of certainty (Gautier and Noetinger, 2004;Jafarpour and Tarrahi, 2011).
As a consequence, several methodologies have aimed at estimating the values of GM parameters, with the general approach of (i) performing "structural analysis", in which the form of the GM is selected, followed by (ii) an estimation of the values of the parameters defining the GM using observation data from the aquifer system.For example, Kitanidis and Vomvoris (1983) and Hoeksema and Kitanidis (1984) used maximum likelihood estimation to estimate values for a two-parameter GM using measurements of log-T and hydraulic head in 1-D and 2-D steady-state flow systems, respectively, in their approach to estimating the spatial distribution of log-T .A more recent review of the technique is given in Kitanidis (1996).More recent studies in the field of petroleum-reservoir engineering (e.g., Yortsos and Al-Afaleg, 1997;Gautier and Noetinger, 2004) have used well test data to estimate parameter values of the permeability variogram.For example, Gautier and Noetinger (2004) expanded on the work of Kitanidis and Vomvoris (1983) to develop a methodology for transient flow.

Objectives of this study
The objectives of this study are three-fold.The first objective is to apply the Kalman Filter parameter estimation methodology within a fully-coupled surface and variably-saturated subsurface flow model to provide more realistic simulation of water table elevation, as well as allow for streamflow to be simulated.To accomplish this, the CATHY model is used in a tilted v-catchment setting, similar in design to the vcatchment used by Camporese et al. (2009), with uncertain initial conditions (i.e., water table elevation) and uncertain patterns of applied water at the ground surface in space and time in a 365-day simulation.An ES is used to assimilate water table elevation data from a reference system to provide an updated estimate of the spatial distribution of log-K.Using uncertain initial conditions and forcing terms provides a stiff test for estimating K (Hendricks Franssen and Kinzelbach, 2008) since values of water table elevation and streamflow are not influenced solely by K.The second objective is to exploit the functionality of CATHY to explore the possibility of using streamflow measurements, solely and jointly with water table elevation data, to condition K.
The third objective is to use the ES in an iterative scheme to identify the parameters of a geostatistical model through assimilation of water table elevation data, and hence provide a new methodology for estimating the value of these parameters.In this study, the ability of the scheme to assess the log-mean and log-variance of a geostatistical model is investigated.Uncertainty in correlation scales is not addressed in this study, but is left to future work.Assessment of the true correlation scale for a given aquifer will likely require the direct assimilation of K measurements, whereas in this study only the model response variables are assimilated.
For the first and second objectives, the influence of the number of measurements and the uncertainty assigned the measurement data on the ability of the ES to provide accurate updates is investigated.Overall, with uncertainty in initial conditions, forcing terms, and geostatistical model parameters, the complexity of real-world systems is approached, providing a key liaison between theory and real-world application.

Methodology
In this section, the theoretical development of the ES update algorithm is presented within the context of estimating the spatial distribution of K using observed water table elevation (W T ) and streamflow (Q) data from a reference catchment system.The general forecast and update steps of the Kalman Filter are first discussed, followed by a modification of these steps for the ES scheme.Using an ensemble of n MC system realizations to establish the uncertainty in the system, the state of the system is estimated using the model forecast step: where f indicates forecast, X f t contains the ensemble of realizations of the forecasted estimate of the system at time t, t represents the solution to the numerical model, and P, X 0 , q, and b represent the system parameters, initial conditions, forcing terms, and boundary conditions, respectively.The numerical model employed in this study is the CATHY model, and is used to generate values of W T and Q as well as establish relationships between the system parameter (i.e., K) and the system response variables (i.e., W T and Q).
CATHY simulates subsurface, overland, and channel flow by coupling the 3-D Richards equation for variably saturated porous media with a 1-D diffusion wave approximation of the de Saint Venant equation for surface flow (Bixio et al., 2000;Camporese et al., 2010).The groundwater flow equation is given by Camporese et al. (2010): where S w = θ/θ s , with θ and θ s as volumetric water content [-] and saturated water content (porosity) [-], respectively, S S is the specific storage coefficient K s is the saturated hydraulic conductivity tensor [LT −1 ] with K s treated as a scalar field when conditions of isotropy are hypothesized, K r (ψ) is the relative hydraulic conductivity function [-], η z = (0, 0, 1), z is the vertical coordinate directed upward [L], and q ss represents distributed source or sink terms Using a 1-D coordinate system s [L] to describe the channel network, the surface water flow equation is given by (Camporese et al., 2010): where and q s is the inflow or outflow rate from the subsurface to the surface In CATHY, Eq. ( 2a) is solved using Galerkin finite elements (FE), whereas Eq. ( 2b) is solved using an explicit time discretization based on the Muskingum-Cunge routing scheme (Orlandini and Rosso, 1996).In this study, the K r (ψ) and S w (ψ) relationships are specified using the formulation of van Genuchten and Nielsen (1985), although other capillary curves are available in CATHY (see Camporese et al., 2010).The channel network is identified using the terrain topography from a digital elevation model (DEM) and the hydraulic geometry concept used by Orlandini and Rosso (1996).The DEM cells are then triangulated to generate a 2-D triangular FE mesh, which is replicated vertically to construct a 3-D tetrahedral FE mesh for the subsurface system.Interaction between surface water and groundwater modeled in CATHY is described by Putti and Paniconi (2004).
As will be discussed further in Sect.2.2, the resulting system state X (Eq. 1) is comprised of: (1) a W T value for each node of the triangular FE mesh, calculated using the vertical profile of ψ for each of these nodes; (2) a K value for each DEM cell in the horizontal direction; (3) a Q value for each DEM grid cell along a stream channel.If e and n denote the number of DEM cells and FE nodes in the horizontal direction, respectively, and g the number of DEM cells along the stream channel, then the dimension d of X is equal to [n + e + g].The forcing terms q in Eq. ( 1) are represented by q ss in Eq. ( 2a), and in this study correspond to rates of applied water at the ground surface, with uncertainty established by sampling values from a prescribed frequency distribution.Uncertainty in X 0 is also included, as discussed in Sect.3.1.
The spatially-variable values of Y K = log K are generated using SKSIM (Baù and Mayer, 2008), a sequential Gaussian simulation algorithm, where the spatial distribution and correlation is established by a normal distribution wherein the geostatistical model is a 2-D exponential covariance model in the logarithmic domain: where µ Y K and σ Y K , and σ 2 Y K are the mean, standard deviation, and variance of the logarithmic distribution of the parameters, respectively, d i s are the components of the distance vector d, and λ i s are the spatial correlation scales in the coordinate directions.

Update of system state
Correction to the a priori estimate X f t is accomplished by assimilating observed system-response data from the true system, thereby merging the model-calculated and observed values.This correction depends on the uncertainty attached to both the a priori estimate (i.e., the model results) and the true values (i.e., the observations from the true system), with uncertainty in the model forecast provided by the spread in the ensemble values and uncertainty in the observed values specified according to data-sampling methods.The correction made to the model-calculated values by the observed values is dictated by the ratio of these uncertainties.If there is less uncertainty attached to the observed data, which is typically the case, then the model-calculated value at the observation location will be corrected to approach the observed value.Furthermore, model results can also receive correction from observed data if the model value is correlated with the model value at the observation location.In this way information from the true state at observation points can be "spread" to regions between observation locations, and hence throughout the model domain.
This correction procedure is carried out through the following Kalman Filter update equation, with the forecasted ensemble X f t corrected, or updated, at a time t using m observed data stored in a vector M t [m]: where X u t [dxn MC ] is the updated ensemble with u denoting update; D t [mxn MC ] holds the ensemble of perturbed values of the measurement data, with the ensemble of values for each measurement value calculated by adding a Guassian perturbation (stored in the matrix E [mxn MC ]) to each of the m observations stored in M t ; H [mxd] contains binary constants (0 or 1) resulting in the matrix product HX f t that holds model results at measurement locations, and κ t [dxm] is the so-called "Kalman Gain" matrix.In this study, observation data are sampled from a known reference state to enable assessment of the ES scheme.
In Eq. ( 4), the difference, or residual, between the model values and observed values is represented by D t − HX f t , with the weighting of the correction and the spatial spread of the information dictated by κ t , which holds the ratio of uncertainties as well as the covariance between model values at each model node (Bailey and Baù, 2011).The form of κ t is: where C f [dxd] and R [mxm] are the forecast error and observation error covariance matrices, respectively, and are defined as: where each column of X[n x n MC ] holds the average value of the ensemble at each location in the domain.
In a straightforward application of data assimilation to a catchment system, Eq. ( 4) would correspond to merging observed values of W T (or Q) with the model-calculated W T field (or Q along a stream channel) in order to provide a W T field that honors the observed W T data.However, doing so only corrects the system response of the model -the structural difference between the a priori model state and the true state that yields differences in the system response will persist indefinitely.To temper these structural differences, it is essential to correct the parameters that drive the system response.This can be accomplished by utilizing the relationships between K and the resulting values of W T and Q as established through CATHY.For example, the crosscovariance submatrix between W T and Y K is defined as: By including the expression for Eq. ( 7) into Eq.( 5), the values of Y K in X f t can be corrected by observation data in D t through the spatial correlation between Y K , W T , and Q.

Forecast-Update scheme for the ensemble smoother
Whereas Eqs. ( 1) and ( 4) are run in a sequential manner in the KF and EnKF schemes, with correction via Eq.( 4) occurring whenever observation data are sampled from the true system, the ES algorithm includes all previous model state and observation data up to the final data sampling time t nF , at which time the ES update routine is run to provide updated system states at all previous collection times.At time t nF , the forecast matrix X f t F and the observation matrix D t F hold the model state ensembles and the perturbed observation data from all data sampling times (t 1 ,t 2 ,. . ., t nF ): where nF is the number of times at which measurements are collected.Within the ES scheme, the forecast covariance matrix C f t nF contains both spatial covariance terms and temporal covariance terms between cell values from different collection times (Evensen, 2007).The measurement error covariance matrix R t nF also is established using the perturbations for each of the measurement values for each of the nF collection times.By inserting X f t nF and D t nF into Eq.( 4) and C f t nF and R t nF into Eq.( 5), the updated system state matrix X u t nF contains the updated model state for each assimilation time.
The Keppenne (2000) algorithm, which provides an efficient numerical strategy for updating the system state for the EnKF scheme and designed for high-resolution real-world climate numerical models, was modified to include model states and observation data from each assimilation time (Bailey and Baù, 2010) and used to compute Eq. ( 4) within the ES framework.The techniques employed by Keppenne (2000) and hence inherent in the update algorithm used in this study do not require the direct assemblage of C f , hence saving on computer memory and preventing numerical issues.

Evaluating the updated system state
The ability of the ES algorithm to bring the forecasted ensemble into conformity with the true system state is quantified through two location-specific parameters EE (ensemble error) and EP (ensemble precision) and two global parameters AE (absolute error ) and AEP (average ensemble precision) (Hendricks Franssen and Kinzelbach, 2008;Bailey and Baù, 2011): where X i is the ensemble mean of the ith location (node in a FE discretization of the domain or cell of the surface DEM grid),X i,true is the reference "true" value of the ith location, and X i,j is the variable value of the i location of the j th ensemble realization.Equation (9a) and (c) provide a measure of the deviation between the model state and the reference state, and Eq. ( 9b) and (d) provide a measure of the spread of the values around the ensemble mean of the model state.The performance of the update routine is measured by calculating the difference between performance parameters of the forecasted and updated ensembles.As a second type of performance measure, the ensemble of CATHY simulations can be rerun using the updated ensemble of Y K fields to determine if the updated Y K ensemble produces simulated   results that match the observed data, i.e., to see if the model has been "calibrated" adequately.This latter method will be demonstrated in Sect.3.3.2.

Iterative method to estimate geostatistical parameters
In the forecast step of Sect.2.1, it is assumed that the parameters defining the GM of Eq. ( 3) are known, i.e., that the parameter values used to generate the ensemble of Y K fields for the forecast simulations are the same as the reference system from which the observation data is sampled.In recognition   that this is generally not the case, the ES scheme decribed in Sect.2.1 through Sect.2.3 is employed in an iterative scheme to discover the geostatistical parameter values of the true system, as shown in Fig. 1.Beginning with a set of estimated GM parameter values, an ensemble of Y K fields is generated using SKSIM and the corresponding ensemble of CATHY flow simulations is run.Upon assimilating observation data from the reference system and conditioning the Y K ensemble, the GM parameter values of the updated Y K ensemble are inferred from the updated ensemble and used to produce the forecast Y K ensemble for the subsequent iteration.This process proceeds until GM parameter values are converged upon.At each iteration the model-calculated values of W T can also be compared to the observed W T data from the true aquifer system, to verify that the estimated GM parameter values from the previous iteration yield a spatial structure of Y K that furnishes the system response of the true system.It should be noted that estimation of the correlation scale λ was not pursued extensively in this study, as it became evident during initial uses of the iterative approach that λ could not be estimated using only a model-response variable such as water table elevation.As discussed in Sect.4, the direct assimilation of Y K values are likely required to provide information regarding λ, and will be pursued in future work.
3 Parameter Estimation of Spatially-Variable K

Forecast ensemble of K, W T , and Q
The catchment system used for the numerical experiment in this study is a 4.05 km by 4.05 km tilted v-shaped catchment, as shown in Fig. 2, with a stream flowing north to south along the central depression of the catchment.The DEM of the surface terrain, discretized using 50 m by 50 m grid cells for a total of e = 81 × 81 = 6561 grid cells, is shown in Fig. 3a with a contour plot of the ground surface elevation.Aquifer thickness varies between 7.5 m and 15.5 m, with the thickest portion under the central depression, as shown in Figs. 2 and  3b.The subsurface is discretized by n L = 10 layers of varying thickness, with thicknesses ranging from 0.375 m near the ground surface to 3.0 m near the aquifer base.
The characteristics of the DEM, the 3-D mesh, and the parameters of the model are summarized in Table 1.The number of nodes in the 2-D surface FE mesh is n = 82 × 82 = 6724.The 3-D mesh is obtained from replicating the 2-D FE mesh through the vertical extent of the subsurface and contains 3 × n L × 2e = 393 660 tetrahedral elements and n × (n L + 1) = 73 964 nodes.Lateral spatial distribution of Y K for the forecast ensemble is generated using mean µ Y K = 1.301 log m day −1 (K = 20 m day −1 ), standard deviation σ 2 Y K = (0.250 log m day −1 ) 2 , and correlation length (λ x =λ y ) = 1000 m, resulting in K values ranging from approximately 0.2 m day −1 to 1500 m day −1 .In this study, n MC = 100 realizations are used for the ensemble.Values of Y K are calculated for each DEM cell, resulting in a value of Y K assigned to each 3-D FE under the cell.In other words, the spatial distribution of Y K is the same for each of the 10 layers in the 3-D mesh.Vertical hydraulic conductivity, K v , is set equal to one-third the value of horizontal K.
The simulation period is one year, from January to December.No-flow boundaries are assigned for every edge of the aquifer (see Fig. 2).Forcing terms q ss consist of uniformly-distributed net infiltration from precipitation during the months of January through March and November through December, and spatially-varying rates of net infiltra-tion from applied water (e.g., irrigation water) in addition to net infiltration from precipitation during the months of April through October.For the latter, rates are applied according to the cultivation pattern depicted in Fig. 4a, with cells shown in blue receiving monthly values randomly sampled from an exponential distribution (rate parameter γ = 0.75).The pattern of cultivation shown in Fig. 4a is the same for each realization, but monthly rates vary across the realizations.The resulting rates are not intended to represent particular irrigation practices, but rather to provide a spatio-temporal variation in the forcing terms of the catchment system.As an example, the rates of net infiltration from combined precipitation and applied water for the month of July for one realization are shown in Fig. 4b, with values ranging from 0.000355 m day −1 to 0.006 m day −1 (represented by white and black, respectively).The depth of monthly net infiltration is presented in Fig. 2.
Initial conditions for each simulation are achieved as follows.First, a 10 000-day spin-up simulation with a uniform net infiltration rate of 0.0012 m day −1 and isotropic, homogeneous aquifer with K = 30 m day −1 is run in order to achieve steady-state conditions in the catchment as determined by water table elevation and streamflow rate.Second, each realization of the ensemble is run for 365 days using a different anisotropic Y K field and infiltration pattern to eliminate the bias due to the initial conditions.The results of this 365-day simulation are then used as initial conditions for the final 365-day simulation period for each realization, with time steps ranging between approximately 0.10 to 1.0 day.Stream inflow (Fig. 2) is set to 0.
An additional CATHY simulation representing the true catchment system is run, with the true Y K field and resulting true W T field (at 365 days) depicted in Fig. 5.The streamflow rate at the outlet cell of the catchment is shown in Fig. 5b, indicating the increased discharge during the months of April through October due to increased rates of net infiltration from precipitation as well as applied water.In comparing the Y K forecast ensemble to the true Y K field using Eq. ( 9c) and (d), the AE and AEP values are 0.619 and 0.388, respectively.
The GM parameter values used to generate the Y K field are the same as used for the forecast ensemble.This assumption will be relaxed in Sect.3.3.2,when the iterative approach presented in Sect.2.5 is used to estimate the true GM parameter values.The ability of the ES scheme to recover the spatial distribution of the true Y K field by assimilating W T and Q measurements into the forecast ensemble results is explored in Sect.3.2.
W T data and Q data are collected from observation wells (Fig. 5b) and stream gaging points (Fig. 2), respectively.CPU (Central Processing Unit) time to run a single realization on an Intel® Core™2 Duo CPU @ 3.00 GHz desktop computer range from approximately 20 min to 180 min, depending on the spatial distribution of Y K .

Update of K ensemble
Observation data from the true catchment system are collected tri-monthly, resulting in four assimilation times during the year.As such, forecast ensemble model results are also saved every three months for use in the ES algorithm.The first set of update scenarios consists of conditioning the forecast Y K ensemble using W T data from the 24 observation wells shown in Fig. 5b, with variations on (i) the number of assimilation times and (ii) the error, defined using coefficient of variation (CV) of W T data, assigned to the observed W T data.For these scenarios, the CPU run-time of the ES update routine is approximately 30 s.
For the scenario where four assimilation times are used and the observation data are assumed to be error-free, the ensemble mean at each computational point for the updated Y K ensemble is shown in Fig. 6a.The AE value of the updated Y K ensemble, measuring the absolute error from the true Y K field, is 0.380, a reduction of 38.6 % from the forecast value of 0.619.In comparison with the true Y K field shown in Fig. 5a, the mean of the updated Y K ensemble in Fig. 6a captures the overall spatial pattern of the true field, with high values of Y K in the northwest and southeast sections of the aquifer and low values in the north-central and southwest portions.The updated AEP value, measuring the spread of the updated Y K ensemble, is 0.203, a reduction of 47.7 % from the forecast value of 0.388.The spatial distribution of EP for the updated Y K ensemble, as calculated by Eq. (9b), is shown in Fig. 7. Notice that the spread of the ensemble values is lowest at and around the locations of the 24 observation wells.
When a CV value of 0.10 is assigned to the W T observation data, the AE and AEP values of the updated Y K ensemble are 0.405 and 0.246, respectively, reductions of 34.5 % and 36.6 % from the forecast values.When a CV value of 0.30 is assigned, the AE and AEP values are 0.441 and 0.284, reductions of 28.8 % and 26.8 % from the forecast values.The ensemble of perturbed values for a data value collected from the observation well located at location (X = 2500 m, Y = 3500 m) for these two cases is shown in Fig. 8.For the  latter case, the ensemble mean of the updated Y K ensemble is shown in Fig. 6b.In comparison to the case of CV = 0.00 (Fig. 6a), the spatial distribution of Y K does not resemble as well the true Y K field shown in Fig. 5a.
When observation data from only one assimilation time (time = 365 days) is assimilated, the AE value of the updated Y K ensemble is 0.408; when observation data from two assimilation times are used (time = 181 days and 365 days), the AE value is 0.385.In comparison to the AE value of 0.380, when four assimilation times are used, the improvement of the updated Y K ensemble with respect to the true Y K field is not considerable.The usefulness of additional assimilation times, however, is seen in the context of observation data error.Figure 9 shows the increase of AE (i.e., the increase in deviation from the true Y K field) with increasing values of CV of the W T observation data, for the cases when one, two, and four assimilation times are used.Notice that the increase of AE is lessened when observation data from multiple times are assimilated, with the best results occurring when 4 assimilation times are used.The use of additional assimilation times yielded no improvement.For the case of CV = 0.50, the AE value using one, two, and four assimilation times is 0.594, 0.563, and 0.471, respectively.
The second set of update scenarios consists of assimilating observed values of Q from the reference system.The ensemble of values of Q for a given surface grid cell are found to be close to lognormally-distributed (Clark et al., 2008), with the coefficient of determination r 2 = 0.702 and the Kolmogorov-Smirnov statistics KS = 0.270 for a lognormal fit.This requires both the Q observation data and Q forecast ensemble to be log-transformed before use in the ES update routine.
When Y Q = log Q data from 4 gaging stations along stream channel from the four measurement times are assimilated into the forecast ensembles of Y Q and Y K , the resulting AE and AEP values of the updated Y K ensemble are 0.512 and 0.305, providing reductions of 17.3 % and 21.3 % from the forecast values.The ensemble mean for the updated Y K ensemble for this scenario is shown in Fig. 10.When compared with the updated Y K ensemble in Fig. 6a, it is clear that observation data from the 24 wells provide an updated Y K ensemble that fits more closely with the true Y K field than using Y Q data.Still, the Y K ensemble mean in Fig. 10 captures the principal features of the spatial pattern of the true field.However, when error is assigned to the observed Y Q data, the ability of the Y Q data to condition the Y K ensemble is reduced dramatically.Table 2 shows the values of AE and AEP of the updated Y K ensemble for scenarios in which the CV of the Y Q data ranges from 0.00 to 1.00.When CV = 0.1,  the reduction in the AE value from the forecast ensemble is 7.3 %; for CV = 0.30, the reduction is only 2.0 %.
Table 3 presents results of scenarios wherein W T and Y Q data are jointly assimilated, with the number of observation wells ranging from 2 to 24.In each case, the observation wells are positioned in a grid network.In order to assess the influence of assimilating Y Q data, four scenarios (1-4) are run with only W T observation data, with the same four scenarios rerun (scenarios 5-8) with the inclusion of assimilating Y Q data from 4 gaging stations.Results indicate that the inclusion of Y Q data only provides enhanced conditioning of the Y K ensemble when the number of observation wells used is less than 8.For example, when 4 observation wells are used, the AE value for the scenarios with and without Y Q data is 0.455 and 0.434, respectively; when two wells are used, the AE value with and without Y Q data is 0.545 and 0.482, respectively.Hence, when sparse W T data are available, the additional Y Q data are able to provide information regarding the spatial distribution of Y K and hence partially maintain the value of AE.

Log-Mean of true Y K field is uncertain
In this section, the ability of the ES update routine to condition the ensemble of Y K using W T observation data sampled from a catchment system where µ Y K is different from the one used in generating the forecast ensemble is explored.In the first scenario, µ Y K of the true Y K field is 1.801 log m day −1 (K = 63.2 m day −1 , one-half order of magnitude higher than the forecast value of 20 m day −1 ); for the second scenario, µ Y K of the true Y K field is 0.801 log m day −1 (K = 6.3 m day −1 , one-half order of magnitude lower than the forecast value of 20 m day −1 ).The value of σ 2 Y K in the true systems and the forecast ensemble is the same [(0.250log m day −1 ) 2 ].In Sect.3.3.2,a more severe test is used to demonstrate the iterative approach described in Sect.2.5.
The true Y K field for the first and second scenarios are shown in Fig. 11a and b, respectively.For the two scenarios, the AE of the forecast Y K ensemble is 0.690 and 0.692, respectively.By assimilating W T data from 24 observation wells for four assimilation times, the AE of the updated Y K ensemble for the two scenario is 0.365 and 0.395, respectively, resulting in reductions of 47.1 % and 43.0 % from the forecast AE values.The spatial distribution of the updated Y K ensemble mean for the two scenarios is shown in Fig. 11c and d.For both scenarios, the updated ensembles capture the general magnitude and spatial distribution of the true Y K field.
For the first scenario, the initially too-low (compared to the true system) values of Y K are conditioned to the higher values present in the true system; for the second scenario, the initially too-high values are conditioned to the lower values.The same effect is observed through a comparison between the reference, forecast, and updated values of Y K across a west-east transect located at Y = 2000 m, as shown in Fig. 12.In both scenarios, the forecast Y K values along the transect are approximately equal to the forecast µ Y K value of 1.301, whereas the updated Y K values have been conditioned to resemble the true profiles.

Iterative approach to discover geostatistical parameter values
To demonstrate the iterative approach, the true aquifer system has µ Y K and σ 2 Y K values of 0.301 log m day −1 (K = 2.0 m day −1 ) and 0.500 (log m day −1 ) 2 .The resulting true Y K field is shown in Fig. 13.For each iteration, observation data from 24 observation wells and four assimilations are used to condition the Y K ensemble.Beginning with a forecast Y K ensemble generated using µ Y K = 1.301 log m day −1 and σ 2 Y K = 0.250 (log m day −1 ) 2 , eight iterations are performed, with the µ Y K and σ 2 Y K values of the updated Y K ensemble after each iteration shown in Fig. 14.As seen in Fig. 14, the value of µ Y K reaches the parameter value from the true system within three iterations, but eight iterations are required to determine if convergence has been achieved.For σ 2 Y K , the value from the updated Y K decreases during the first several iterations, but eventually converges upon a value (0.580) slightly higher than the true value of 0.500 (log m day −1 ) 2 .If the GM parameter values of the true system were unknown, then it would be assumed that µ Y K is just under 0.300 log m day −1 and σ 2 Y K is equal to 0.580 (log m day −1 ) 2 .Besides the convergence to the true GM parameter values, the approach of the updated Y K ensemble to the spatial distribution of the true Y K field is demonstrated in Figs. 15 and  16.The ensemble mean of the updated Y K ensemble for iterations 1 through 4 is shown in Fig. 15a-d, with the AE value generally decreasing from the forecast value of 1.106 (0.755, 0.565, 0.507, and 0.518, respectively), although a slight increase occurs between iterations 3 and 4. The structure of the Y K spatial distribution, however, progressively approaches the pattern of the true Y K field shown in Fig. 13 with each successive iteration.
Figure 16  Finally, comparisons are made between observed W T data from the true system and model-calculated W T values generated by CATHY using the updated values of µ Y K and σ 2 Y K from the previous iteration.This is especially important since such a comparison, i.e., re-running the numerical model using the estimated parameter values and comparing model results with observed data at observation locations, is generally the only means by which the parameter estimation method can be verified.Figure 17 shows the comparison between observed W T data and model-calculated W T values from the forecast ensemble, the W T ensemble generated using the updated GM parameter values from the first iteration, and the W T ensemble using the updated GM parameters from the second iteration.Comparisons for times = 273 days and 365 days are shown in Fig. 17a and b, respectively.
The match between the forecast values and the true value is much improved upon using the results from the first iteration, and an excellent match occurs using the results from the second iteration.Quantitatively, the sum of squared differences between the model results and true system values is 95.32, 4.29, and 0.59, respectively for time = 273 days, and 99.16, 4.98, and 1.83, respectively for time = 365 days.Notice that the forecasted W T values are lower than the observed W T values from the true aquifer system, since the forecast Y K values are generated using a higher value of µ Y K .However, using the lower updated value of µ Y K from the first and second iterations, the W T values become higher and more in accordance with the observed W T values.

Conclusions
The ES update routine, a derivative of the Kalman Filter approach, has been evaluated for the estimation of spatially-variable K in a catchment system using the fullycoupled, surface-subsurface flow model CATHY.A 4.05 km by 4.05 km tilted v-catchment was used in demonstration, with spatio-temporal variability in forcing terms to provide increased uncertainty in the system and to strive to mimic real-world conditions.
Both W T data and Q data were collected from a reference catchment system and assimilated into the ensemble of model results to condition the spatial distribution of K to approach the reference K field.Assimilating W T from a network of observation wells provided a distinct improvement in the K ensemble in relation to the true K field, with sets of data from multiple collection times tempering the decrease in improvement when error was assigned to the observed W T data.Assimilating Q data only slightly improved the K ensemble in relation to the true K field.Jointly assimilating Q and W T data only improved the estimate of K when data from a small number (2,4) of observation wells were assimilated.This is due to the region of influence of Q, i.e., the regions of the aquifer where the K values directly influence Q and hence can be conditioned by observed values of Q, being small compared to the collective region of influence of W T values at the observation wells.
For cases in which the parameter values defining the geostatistical structure of the aquifer system are uncertain to a small degree (i.e., mean of true K field is one-half order of magnitude different than the assumed mean), the methodology is still able to condition adequately the forecast K ensemble to approach the magnitude and spatial structure of the true aquifer system.For more severe cases, i.e., the true and assumed means are different by an order of magnitude and the true and assumed variance of the K field is different, an iterative process using the ES is used to converge upon the true geostatistical parameter values.Results indicate that the process is successful in approximating the true values.
For the present study uncertainty in the correlation length of the K field is not investigated, and an amendment to the iterative scheme to converge upon unknown correlation length is left to future research, with assimilation of measurements of K likely necessary.Future studies also include an application of the methodology to an actual catchment system in order to estimate the geostatistical parameter values as well as the spatial distribution of K.

Fig. 1 .
Fig. 1.Flow chart for the iterative approach using an Ensemble Smoother to estimate geostatistical values.

Fig. 2 .
Fig. 2. Conceptual model of tilted v basin through the catchment outlet.The monthly depth of -shaped catchment, with groundwater feeding the river flowing out of the net infiltration from precipitation is also shown.

Fig. 3 .
Fig. 3. Contour representation of (A) ground surface elevation and (B) aquifer thickness.Both datasets are used to create the threedimensional subsurface finite-element mesh.

Table 1 .
Frequency distribution of the ensemble of perturbed values for the observed W T value of 18.5 m collected on day 91 at location (X = 2500 m, Y = 3500 m) of the "true aquifer", for CV values of 0.10 and 0.30.hydraulic conductivity K s GM: Eq. (3) Mean µ of K fields log 1.30 (m d −1 ) Variance σ 2 of K fields log 0.25 (m d −1 ) 2 Correlation Length λ of K fields 1000 m Specific storage S s 0.01 m −1 Porosity n 0.35 Residual moisture content θ r 0.061 van Genuchten parameters α = 0.43 m −1 , n = 1.70 Simulation details Monte Carlo Ensemble Size 100 Simulation period days 730.0

Fig. 4 .
Fig. 4. (A) Cultivated (blue) and non-cultivated fields (white), with cultivated fields receiving additional applied water during the months of April through October, and (B) Net infiltration in the month of July (of the second year) for the true system, with values ranging from 0.000355 m day −1 to 0.006 m day −1 (represented by white and black, respectively).

Fig. 5 .
Fig. 5. (A) Reference field of true Y K and (B) corresponding W T field at time = 365 days as calculated by CATHY.In (B), red crosses indicate the location of 24 observation wells.The streamflow at the outlet cell during the 365-day simulation is shown in the subpanel.

Fig. 6 .
Fig. 6.Spatial distribution of the updated Y K ensemble using data from 24 observation wells sampled 4 times during the 365-day period, with the CV of observation data set to (A) 0.0 and (B) 0.3.Compare to the true state shown in Figure 5A.

Fig. 7 .Figure 7 . 4 Figure 8 ..
Fig. 7. Spatial distribution of EP for the updated Y K ensemble using data from 24 observation wells sampled 4 times during the 365-day period, with the CV of observation data set to 0.0.

Fig. 8 .
Fig. 8. Frequency distribution of the ensemble of perturbed values for the observed W T value of 18.5 m collected on day 91 at location (X = 2500 m, Y = 3500 m) of the "true aquifer", for CV values of 0.10 and 0.30.

1Figure 9 .
Figure 9.The increase of AE (i.e., the decrease in conditioning to the true 2

Fig. 9 .
Fig. 9.The increase of AE (i.e., the decrease in conditioning to the true Y K field) with increasing values of CV of the W T observation data, for the cases when one, two, and four assimilation times are used.

Fig. 10 .
Fig. 10.Spatial distribution of the mean of the updated Y K ensemble using Y Q data from four gaging stations along the stream channel sampled four times during the 365-day period, with the CV of observation data set to 0.0.Compare to the true state shown in Fig. 5a.

Fig. 11 .
Fig. 11.(A) True Y K field using µ y K = 1.801 log m day −1 (K = 63.2 m day-1), (B) True Y K field using µ y K = 0.801 log m day −1 (K = 6.3 m day-1), (C) mean distribution of the updated Y K ensemble for the first scenario (µ y K = 1.801), and (D) mean distribution of the updated Y K ensemble for the second scenario (µ y K = 0.801).For both scenarios scenario, 24 W T data from four assimilation times are used.

Fig. 12 .
Fig. 12.Comparison of reference and updated values of Y K along the Y = 2000 m transect where the true is higher (1.801 log m day −1 ) and lower (0.801 log m day −1 ) than the of the forecast ensemble (1.301 log m day −1 ).The values from the reference state are shown in red, and the updated values are shown in blue.For both scenarios, 24 W T measurements are assimilated.

Fig. 14 .
Fig. 14.Progression of the estimated GM parameters and , demonstrating the convergence of the parameter values to ∼0.300 log m day −1 and 0.580 (log m day −1 ) 2 , respectively.

Fig. 15 .Figure 16 .
Fig. 15.Ensemble mean of updated Y K ensemble for the (A) 1st iteration, (B) 2nd iteration, (C) 3rd iteration, and (D) 4th iteration, for the case where the true values of and are 0.301 log m day −1 and 0.500(log m day −1 ) 2 , respectively.Compare to the true state in Fig. 13.The forecast AE values was 1.106.

Fig. 16 .
Fig. 16.Value of Y K along the Y = 2000 m transect for the reference state (red), the foreacast ensemble mean (dotted black line) and the update ensemble mean for the 1st (solid black line) and 2nd (solid blue line) iteration, for the case where the true values of and are 0.301 log m day −1 and 0.500 (log m day −1 ) 2 , respectively.

Fig. 17 .
Fig. 17.Comparison between model-calculated W T values and observed W T data from the true system for the forecast W T ensemble, the W T ensemble generated using the updated GM parameter values from the 1st iteration, and the W T ensemble using the updated GM parameters from the 2nd iteration.Comparison are made for (A) time = 273 days and (B) time = 365 days.
shows the true Y K values, forecast Y K values, and the updated values of Y K for the 1st and 3rd iterations at the Y = 2000 west-east transect.Dramatic improvement in the updated Y K values in relation to the true Y K values occurs from the forecast to the first iteration, and from the first iteration to the third iteration.

Table 2 .
The increase of AE (i.e., the decrease in conditioning to the true Y K field) with increasing values of CV of the W T observation data, for the cases when one, two, and four assimilation times are used.

Table 3 .
Spatial distribution of the mean of the updated Y K ensemble using Y Q data from four gaging stations along the stream channel sampled four times during the 365-day period, with the CV of observation data set to 0.0.Compare to the true state shown in Fig.5a.