A Bayesian Consistent Dual Ensemble Kalman Filter for State-Parameter Estimation in Subsurface Hydrology

Ensemble Kalman filtering (EnKF) is an efficient approach to addressing uncertainties in subsurface groundwater models. The EnKF sequentially integrates field data into simulation models to obtain a better characterization of the model's state and parameters. These are generally estimated following joint and dual filtering strategies, in which, at each assimilation cycle, a forecast step by the model is followed by an update step with incoming observations. The Joint-EnKF directly updates the augmented state-parameter vector while the Dual-EnKF employs two separate filters, first estimating the parameters and then estimating the state based on the updated parameters. In this paper, we reverse the order of the forecast-update steps following the one-step-ahead (OSA) smoothing formulation of the Bayesian filtering problem, based on which we propose a new dual EnKF scheme, the Dual-EnKF$_{\rm OSA}$. Compared to the Dual-EnKF, this introduces a new update step to the state in a fully consistent Bayesian framework, which is shown to enhance the performance of the dual filtering approach without any significant increase in the computational cost. Numerical experiments are conducted with a two-dimensional synthetic groundwater aquifer model to assess the performance and robustness of the proposed Dual-EnKF$_{\rm OSA}$, and to evaluate its results against those of the Joint- and Dual-EnKFs. The proposed scheme is able to successfully recover both the hydraulic head and the aquifer conductivity, further providing reliable estimates of their uncertainties. Compared with the standard Joint- and Dual-EnKFs, the proposed scheme is found more robust to different assimilation settings, such as the spatial and temporal distribution of the observations, and the level of noise in the data. Based on our experimental setups, it yields up to 25% more accurate state and parameters estimates.


I. INTRODUCTION
In modern hydrology research, uncertainty quantification studies have focused on field applications, including surface and subsurface water flow, contaminant transport, and reservoir engineering. The motivations behind these studies were driven by the uncertain and stochastic nature of hydrological systems. For instance, surface rainfall-runoff models that account for soil moisture and streamflows are subject to many uncertain parameters, such as free and tension water storage content, water depletion rates, and melting threshold temperatures [1]. Groundwater flow models, on the other hand, depend on our knowledge of spatially variable aquifer properties, such as porosity and hydraulic conductivity, which are often poorly known [2], [3]. In addition, contaminant transport models that investigate the migration of pollutants in subsurface aquifers are quite sensitive to reaction parameters like sorption rates, radioactive decay, and biodegradation [4], [5]. To this end, it is important to study the variability of hydrological parameters and reduce their associated uncertainties in order to obtain reliable simulations. To achieve this goal, hydrologists have used various inverse and Monte Carlo-based statistical techniques with the standard procedure of pinpointing parameter values that, when integrated with the simulation models, allow some system-response variables (e.g., hydraulic head, solute concentration) to fit given observations [6]- [11].
Recently, sequential data assimilation techniques, such as the particle filter (PF), has been proposed to handle any type of statistical distribution, Gaussian or not, to properly deal with strongly nonlinear systems [12]. The PF may require, however, a prohibitive number of particles (and thus model runs) to accurately sample the distribution of the state and parameters, making this scheme computationally intensive for large-scale hydrological applications [13]- [15]. This problem has been partially addressed by the popular ensemble Kalman filter (EnKF), which further provides robustness, efficiency and nonintrusive formulation [16]- [19] to tackle the state-parameters estimation problem.
The EnKF is a filtering technique that is relatively simple to implement, even with complex nonlinear models, requiring only an observation operator that maps the state variables from the model space into the observation space. Compared with traditional inverse and direct optimization techniques, which are generally based on least-squares-like formulations, the EnKF has the advantage of being able to account for model errors that are not only present in the uncertain parameters but also in the model structure and inputs, such as external forcings [3]. In addition, and because of its sequential formulation, it does not require storing all past information about the states and parameters, leading to consequent savings in computational cost [5], [20].
The EnKF is widely used in surface and subsurface hydrological studies to tackle state-parameters estimation problems. Two approaches are usually considered based on the joint and the dual stateparameter estimation strategies. The standard joint approach concurrently estimates the state and the parameters by augmenting the state variables with the unknown parameters, that do not vary in time. The parameters could also be set to follow an artificial evolution (random walk) before they get updated with incoming observations [21]. One of the early applications of the Joint-EnKF to subsurface groundwater flow models was carried in [2]. In their study, a conceptual subsurface flow system was considered and ensemble filtering was performed to estimate the transient pressure field alongside the hydraulic conductivity. In a reservoir engineering application, [22] considered a two-dimensional North Sea field model and considered the joint estimation problem of the dynamic pressure and saturation on top of the static permeability field. Groundwater contamination problems were also tackled using the Joint-EnKF (e.g., [4], [23]), in which the hydraulic head, contaminant concentration and spatially variable permeability and porosity parameters were estimated for coupled groundwater flow and contaminant transport systems.
In contrast with the joint approach, the dual approach runs two separate interactive EnKFs; one for the parameters and the other for the state. This separation is motivated by the fact that it is the marginal posterior distributions, and not the joint posterior distribution of the state and parameters vector, that are the relevant quantities to be estimated [24]. The Dual-EnKF has been applied to streamflow forecasting problems using rainfall-runoff models (e.g., [1], [24]). Recently, the authors in [25] successfully implemented the Dual-EnKF with a subsurface compositional flow model using chemical composite data and compared it to the standard Joint-EnKF. The authors concluded that the dual estimation scheme provides more accurate state and parameters estimates than the joint scheme when implemented with large enough ensembles. In terms of complexity, however, the dual scheme is computationally more intensive, as it requires integrating the filter ensemble twice with the numerical model at every assimilation cycle.
In this study, we introduce a new Dual-EnKF algorithm following the one-step-ahead (OSA) smoothing formulation of the Bayesian filtering problem, in which the measurement-update step precedes the forecast (or time-update) step. Unlike the standard Dual-EnKF, this OSA smoothing-based filter, called Dual-EnKF OSA , is consistent with the Bayesian formulation of the state-parameter estimation problem and exploits the observations in both the state smoothing and forecast steps. It further provides a theoretically sound framework for dual filtering and is shown to be beneficial in terms of estimation accuracy, as compared to the standard Joint-and Dual-EnKFs, without requiring any significant additional computational cost. This is confirmed via synthetic numerical experiments using a groundwater flow model and estimating the hydraulic head and the conductivity parameter field.
The rest of the paper is organized as follows. Section II reviews the standard Joint-and Dual-EnKF strategies. The Dual-EnKF OSA is derived in Section III and its relation with the Joint-and Dual-EnKFs is discussed. Section IV presents a conceptual groundwater flow model and outlines the experimental setup. Numerical results are presented and discussed in Section V. Conclusions are offered in Section VI, followed by an Appendix.

A. Problem formulation
Consider a discrete-time state-parameter dynamical system: in which x n ∈ R Nx and y n ∈ R Ny denote the system state and the observation at time t n , of dimensions N x and N y respectively, and θ ∈ R Nθ is the parameter vector of dimension N θ . M n is a nonlinear operator integrating the system state from time t n to t n+1 , and the observational operator at time t n , H n , is assumed to be linear for simplicity; the proposed scheme is still valid in the nonlinear case.
The model process noise, η = {η n } n∈N , and the observation process noise, ε = {ε n } n∈N , are assumed to be independent, jointly independent and independent of x 0 and θ, which, in turn, are assumed to be independent. Let also η n and ε n be Gaussian with zero means and covariances Q n and R n , respectively.
Throughout the paper, y 0:n def = {y 0 , y 1 , · · · , y n }, and p(x n ) and p(x n |y 0:l ) stand for the prior probability density function (pdf) of x n and the posterior pdf of x n given y 0:l , respectively. All other used pdfs are defined in a similar way.
We focus on the state-parameter filtering problem, say, the estimation at each time, t n , of the state, x n , as well as the parameters vector, θ, from the history of the observations, y 0:n . The standard solution of this problem is the a posteriori mean (AM): which minimizes the a posteriori mean square error. In practice, analytical computation of (2) and (3) is not feasible, mainly due to the nonlinear character of the system. The Joint-and Dual-EnKFs have been introduced as efficient schemes to compute approximations of (2) and (3). Before reviewing these algorithms, we recall the following classical results of stochastic sampling that we will extensively use in the derivation of the filtering algorithms.
Property 1 (Hierarchical sampling [26]). Assuming that one can sample from p(x 1 ) and p(x 2 |x 1 ), then a sample, x * 2 , from p(x 2 ) can be drawn as follows: Property 2 (Conditional sampling [27]). Consider a Gaussian pdf, p(x, y), with P xy and P y denoting the cross-covariance of x and y and the covariance of y, respectively. Then a sample, x * , from p(x|y), can be drawn as follows: B. The Joint-and Dual-EnKFs 1) The Joint-EnKF: The key idea behind the Joint-EnKF is to transform the state-parameter system (1) into a classical state-space system based on the augmented state, z n = x T n θ T T , on which the classical EnKF can be directly applied. The new augmented state-space system can be written as: representing p(z n−1 |y 0:n−1 ), the EnKF uses the augmented state model (1 st Eq. of (4)) to compute the , approximating p(z n |y 0:n−1 ). The observation model (2 nd Eq. of (4)) is then used to obtain the analysis ensemble, {x , at time t n . Let, for an ensemble ,r denotes its empirical mean and S r a matrix with N e -columns whose m th column is defined as r (m) −r . The Joint-EnKF steps can be summarized as follows: • Forecast step: The parameters vector members, θ An approximation of the forecast state,x n|n−1 , which is, by definition, the mean of p(x n |y 0:n−1 ) (similar to (2)), is given by the empirical mean of the forecast ensemble,x f n . The associated forecast error covariance is estimated as P x f n = (N e − 1) −1 S x f n S T x a,(m) The analysis estimates, (2) and (3), and their error covariances, can thus be approximated by the analysis ensemble means,x a n andθ |n , and P x a n = (N e − 1) −1 S x a n S T x a n and P θ|n = (N e − 1) −1 S θ|n S T θ|n , respectively.
2) The Dual-EnKF: In contrast with the Joint-EnKF, the Dual-EnKF is designed following a conditional estimation strategy, operating as a succession of two EnKF-like filters. First, a (parameter) filter , is applied based on the following two steps. , again in two steps that can be summarized as follows.
x a,(m) Note that we ignore here the term for the process noise, η n , which is often the case in geophysics applications. We further consider the same observation forecast members y f,(m) n in both the Joint-EnKF (10) and Dual-EnKF (11) (which means that we consider the same measurement noise members, ε Such a separation of the update steps was shown to provide more consistent estimates of the parameters, especially for strongly heterogeneous subsurface formations as suggested by [28] in their confirming-step EnKF algorithm. The dual-update framework was further shown to provide better performances than the Joint-EnKF at the cost of increased computational burden (see for instance, [1], [24], [25]).
3) Probabilistic formulation: Following a probabilistic formulation, the augmented state system (4) can be viewed as a continuous state Hidden Markov Chain (HMC) with transition density, and likelihood, where N v (m, C) represents a Gaussian pdf of argument v and parameters (m, C).
Eq. (14) refers to a Markovian step (or time-update step) and uses the transition pdf, p (x n |x n−1 , θ), of the Markov chain, {z n } n , to compute the forecast pdf of z n from the previous analysis pdf. Eq. (16) refers to a Bayesian step (or measurement-update step) since it uses the Bayes' rule to update the forecast pdf of z n using the current observation y n . Now, establishing the link between the Joint-EnKF and Eqs. then provide the analysis ensemble of the state (7) and the parameters (8).
Regarding the Dual-EnKF, the forecast ensemble of the state and observations in the parameter filter can be obtained following the same process as in the Joint-EnKF. This is followed by the computation of the analysis ensemble of the parameters using Property 2 and p (θ|y 0:n ) = p (θ, y n |y 0:n−1 ) p (y n |y 0:n−1 ) .
However, in the state filter, the ensemble, { x f,(m) n } Ne m=1 , obtained via Eq. (9) in the forecast step does not represent the forecast pdf, p(x n |y 0:n−1 ), since Eq. (9) involves θ (m) |n rather than θ (m) |n−1 . Accordingly, the Dual-EnKF is basically a heuristic algorithm in spite of its proven performance compared to the Joint-EnKF.
In the next section, we follow the same probabilistic framework to derive a new Dual-EnKF algorithm, in which the state is updated by the observations at both the measurement-update step and the time-update step, without violating the Bayesian filtering formulation of the state-parameter estimation problem. Our goal is to build a fully Bayesian consistent dual-like EnKF that retains the separate formulation of the state and parameters update steps.
The algorithm of this section is inspired by the fact that the classical (time-update, measurementupdate) path (14)- (16), that involves the forecast pdf p(z n |y 0:n−1 ) when moving from the analysis pdf p(z n−1 |y 0:n−1 ) to the analysis pdf at the next time p(z n |y 0:n ), is not the only possible one. Here, we reverse the order of the time-update and measurement-update steps to derive a new dual algorithm involving the one-step-ahead smoothing pdf, p(z n−1 |y 0:n ), between two successive analysis pdfs, p(z n−1 |y 0:n−1 ) and p(z n |y 0:n ). This is expected to improve the filter estimates since the state is constrained by more observations than the standard joint and dual schemes. The reader may consult, for instance, [29], in which the one-step-ahead smoothing-based filtering problem has been studied and based on which a class of KF-and PF-based algorithms have been derived. The more recent work in [30] also proposes a number of algorithms estimating both the system state and the model noise based on the same strategy.
A. The one-step-ahead smoothing-based dual filtering algorithm The analysis pdf, p(x n , θ|y 0:n ), can be computed from p(x n−1 , θ|y 0:n−1 ) in two steps: • Smoothing step: The one-step-ahead smoothing pdf, p(x n−1 , θ|y 0:n ), is first computed as, with, Eq. (19) is derived from the fact that in the state-parameter model (1), the observation noise, ε n , and the model noise, η n−1 , are independent of (x n−1 , θ) and past observations y 0:n−1 .
The smoothing step (18) is indeed a measurement-update step since conditionally on y 0:n−1 , Eq.
• Forecast step: The smoothing pdf at t n−1 is then used to compute the current analysis pdf, p(x n , θ|y 0:n ), with, which, in turn, arises from the fact that ε n and η n−1 are independent of (x n−1 , θ) and y 0:n−1 (see smoothing step above). We note here that only the (marginal) analysis pdf of x n , p(x n |y 0:n ), is of interest since the analysis pdf of θ has already been computed in the smoothing step.
From (21), p(x n |x n−1 , θ, y 0:n ) = p(x n |x n−1 , θ, y n ). Thereby, there is a similarity between Eq. (20) and the forecast step (14), in the sense that (20) can be seen as a forecast step once the current observation y n is known, i.e., (20) coincides with " (14) given the observation y n ". Accordingly, and without abuse of language, we refer to (20)-(21) as the forecast step.

B. Ensemble Formulation
Since it is not possible to derive the analytical solution of (18)-(21) because of the nonlinear character of the model, M(.), we propose an EnKF-like formulation, assuming that p(y n , z n−1 |y 0:n−1 ) is Gaussian for all n. This assumption implies that p(z n−1 |y 0:n−1 ), p(z n−1 |y 0:n ) and p(y n |y 0:n−1 ) are Gaussian. , as with η , as The (cross-) covariances in equations (24) and (25) are defined and practically evaluated from the ensembles as in the EnKF (see Subsection II-B1): 2 Furthermore, one can verify that Eq. (21) leads to a Gaussian pdf: with K n = Q n−1 H T n H n Q n−1 H T n + R n −1 and Q n−1 = Q n−1 − K n H n Q n−1 . However, when the state dimension, N x , is very large, the computational cost of K n and Q n−1 (which may be an off-diagonal matrix even when Q n−1 is diagonal) may become prohibitive. One way to avoid this problem is to directly sample from p(x n |x n−1 , θ, y n ) without explicit computation of this pdf in (29).
Using Properties 1 and 2, an explicit form of such samples can be obtained as (see Appendix), where the (cross)-covariances, P ξn, yn and P yn , are evaluated from { ξ x a,(m) In contrast with the Dual-EnKF, which uses θ n−1 after an update with the current observation, y n , following (24). Therefore, when including the Kalman-like correction term as well, the observation, y n , is used three times in the the Dual-EnKF OSA in a fully consistent Bayesian formulation, but only twice in the Dual-EnKF. This means that the Dual-EnKF OSA exploits the observations more efficiently than the Dual-EnKF, which should provide more information for improved state and parameters estimates.
Despite the smoothing framework of the Dual-EnKF OSA , this algorithm obviously addresses the state forecast problem as well. As discussed in the smoothing step above, the (one-step-ahead) forecast members are inherently computed. The j-step-ahead forecast member, denoted by x (m) n+j|n for j ≥ 2, can be computed following a recursive procedure, where for = 2, 3, · · · , j, one has  i.e., the number of state variables is much larger than the number of observations. This is generally the case for subsurface flow applications due to budget constraints given the consequent costs needed for preparing, drilling, and maintaining subsurface wells.

A. Transient groundwater flow problem
We adopt in this study the subsurface flow problem of [31].  x is the linear KF), C θ : parameter model cost (usually free ≡ identity), C y : observation operator cost (= N y N x in the linear KF), S x : storage volume for one state vector, S θ : storage volume for one parameter vector.

Algorithm
Time-update Measurement-update Storage and 15 m are placed on the west and east ends of the aquifer, respectively, with an average saturated thickness, b, of 25 m. The north and south boundaries are assumed to be Neumann with no flow conditions ( Figure 1). The mesh is discretized using a cell-centered finite difference scheme with 10 m × 20 m rectangles, resulting in 2500 elements. The following 2D saturated groundwater flow system is solved: We consider a dynamically complex experimental setup that is similar to a real-world application and is based on various time-dependent external forcings. The recharge is assumed spatially heterogenous and sampled using the GCOSIM3D toolbox [32] with statistical parameters shown in Table II Prior to assimilation, a reference run is first conducted for each experimental setup using the prescribed parameters above, and is considered as the truth. We simulate the groundwater flow system over a yearand-a-half period using the classical fourth-order Runge Kutta method with a time step of 12 hours. The initial hydraulic head configuration is obtained after a 2-years model spin-up starting from a uniform indicate pumping or groundwater that is being removed from the aquifer. Right-Panel: Reference heterogenous spatial recharge values obtained using the sequential Gaussian simulation code [32] with parameters given in Table   II. The black squares represent the pumping wells whereas the black circles denote the position of 3 monitoring wells. The black diamond is a control well.

B. Assimilation Experiments
To imitate a realistic setting, we impose various perturbations on the reference model and set our goal to estimate the water head and the hydraulic conductivity fields using an imperfect forecast model and perturbed data extracted from the reference (true) run. This experimental framework is known as "twin-experiments". In the forecast model, we perturb both transient pumping and spatial recharge rates.
The perturbed recharge field, as compared to the reference recharge in Figure 2, is sampled with different variogram parameters as shown in Table II. Pumping rates from PW1, PW2 and PW3 are perturbed by adding a Gaussian noise with mean zero and standard deviation equal to 20% of the reference transient rates. The flow field simulated by the forecast (perturbed) model after 18 months is shown in the right panel of Figure 3. Compared to the reference field, there are clear spatial differences in the hydraulic head, especially around the first and second pumping wells. 14    To initialize the filters, we follow [25] and perform a 5-years simulation run using the perturbed forecast model starting from the mean hydraulic head of the reference run solution. Then, we randomly select a set of N e hydraulic head snapshots to form the initial state ensemble. By doing so, the dynamic head changes that may occur in the aquifer are well represented by the initial ensemble. The corresponding parameters' realizations are sampled with the geostatistical software, GCOSIM3D, using the same variogram parameters of the reference conductivity field but conditioned on two hard measurements as indicated by black crosses in Figure 1. The two data points capture some parts of the high conductivity regions in the domain, and thus one should expect a poor representation of the low conductivity areas in the initial log(K) ensemble. This is a challenging case for the filters especially when a sparse observational network is considered. To ensure consistency between the hydraulic heads and the conductivities at the beginning of the assimilation, we conduct a spin-up of the whole state-parameters ensemble for a 6-months period using perturbed recharge time-series for each ensemble member.
The filter estimates resulting from the different filters are evaluated based on their average absolute forecast errors (AAE) and their average ensemble spread (AESP): where x t i is the reference "true" value of the variable at cell i, x f,e j,i is the forecast ensemble value of the variable, andx f,e i is the forecast ensemble mean at location i. AAE measures the estimate-truth misfit and AESP measures the ensemble spread, or the confidence in the estimated values [3]. We further assess the accuracy of the estimates by plotting the resulting field and variance maps of both hydraulic head and conductivities.

A. Sensitivity to the ensemble size
We first study the sensitivity of the three algorithms to the ensemble size, N e . In realistic groundwater applications, we would be restricted to small ensembles due to computational limitations. Obtaining accurate state and parameter estimates with small ensembles is thus desirable. We carry the experiments using three ensemble sizes, N e = 50, 100 and 300, and we fix the frequency of the observations to half a day, the number of wells to nine (Figure 3, left observation network) and the measurement error to 0.50 m. We plot the resulting AAE time series of the state and parameters in Figure 4. As shown, the performance of the Joint-EnKF, Dual-EnKF and Dual-EnKF OSA improves as the ensemble size increases, reaching a mean AAE of 0.161, 0.160, and 0.156 m-water for N e = 300, respectively. The Joint-EnKF and the Dual-EnKF exhibit similar behaviors, with a slight advantage for the Dual-EnKF. As demonstrated in [25], the Dual-EnKF is generally expected to produce more accurate results only when large enough ensembles are used. We have tested the Joint-and the Dual-EnKFs using 1000 members and found that the Dual EnKF is around 9% more accurate in term of AAE. The proposed Dual-EnKF OSA is the most accurate in all tested scenarios. On average, with changing ensemble size, the Dual-EnKF OSA leads to about 7% improvement compared with the joint and dual schemes. In terms of the conductivity estimates, the proposed scheme produces more accurate estimates for all three ensemble sizes. At the early assimilation stage, the three schemes seem to provide similar results, but this eventually changes after 6 months and the Dual-EnKF OSA clearly outperforms the standard schemes as the complexity of the model increases because of the perturbed recharge and pumping rates. We furthermore examined the uncertainties of these forecast estimates by computing the spread of both the hydraulic head and conductivity ensembles. To do this, we evaluated the mean AESP of both variables and tabulated the results for the three ensemble sizes in Table III. For all schemes and as expected, the spread seems to increase as the ensemble size increases. Compared to the joint and the dual schemes, the Dual-EnKF OSA retains the smallest mean AESP for all cases, suggesting more confidence in the head and conductivity estimates. In terms of computational cost, we note that our assimilation results were obtained using a 2.30 GHz MAC workstation and 4 cores for parallel looping while integrating the ensemble members. The Joint-EnKF is the least intensive requiring 70.61 sec to perform a year-and-a-half assimilation run using 50 members. The Dual-EnKF and Dual-EnKF OSA , on the other hand, require 75.37 and 77.04 sec, respectively. The Dual-EnKF is more intensive than the Joint-EnKF because it includes an additional propagation step of the ensemble members as discussed in Section III.C. Likewise, the proposed Dual-EnKF OSA requires both an additional propagation step and an update step of the state members. Its computational complexity is thus greater than the joint scheme and roughly equivalent to that of the Dual-EnKF. Note that in the current setup the cost of integrating the groundwater model is not very significant as compared to the cost of the update step. This might however not be the case in large-scale hydrological applications.

B. Sensitivity to the frequency of observations
In the second set of experiments, we change only the temporal frequency of observations, i.e., the times at which head observations are assimilated. We implement the three filters with 100 members and use data from nine observation wells perturbed with 0.10 m noise. Figure 5 plots the mean AAE of the hydraulic conductivity estimated using the three filters for the six different observation frequency scenarios. The Dual-and Joint-EnKFs lead to comparable performances, but the latter performs slightly better when data are assimilated more frequently, i.e., every five and three days. The performance of the proposed Dual-EnKF OSA , as seen from the plot, is rather effective and results in more accurate estimates than those obtained by the other two filters. The largest improvements resulting from this scheme are obtained when assimilating data every 1, 3, and 5 days. The improvements over the joint and the dual schemes decrease as the frequency of observations in time decreases. The reason for this is related to the nature of the Dual-EnKF OSA algorithm, which adds a one-step-aheadsmoothing to the analyzed head ensemble members before updating the forecast parameters and the state samples. Therefore, the more data available in time, the greater the number of smoothing steps applied, and hence the better the characterization of the state and parameters. To illustrate, the smoothing step of the state ensemble enhances its statistics and eventually provides more consistent state-parameters crosscorrelations to better predict the data at the next update step. When assimilating hydraulic head data on a daily basis, the proposed Dual-EnKF OSA leads to about 24% more accurate conductivity estimates than the Joint and Dual-EnKFs.

Dual-EnKF
Dual-EnKF-OSA We also compared the hydraulic head estimates when changing the temporal frequency of observations.
Similar to the parameters, the improvements of the Dual-EnKF OSA algorithm over the other schemes become significant when more data are assimilated over time. Overall, the benefits of the proposed scheme are more pronounced for the estimation of the parameters because the conductivity values at all aquifer cells are indirectly updated using hydraulic head data, requiring more observations for efficient estimation.
One effective way to evaluate the estimates of the state is to examine the evolution of the reference heads and the forecast ensemble members at various aquifer locations. For this, we plot in Figure 6 the true and the estimated time-series change in hydraulic head at the assigned monitoring wells as they result from the Joint-EnKF, Dual-EnKF and the Dual-EnKF OSA . We use 100 ensemble members and assume the 9 data points are available every five days. At MW1, the performance of the three filters is quite similar and they all successfully reduce the uncertainties and recover the true evolution of the hydraulic head at that location. We note that between the 5 th and the 9 th month, the Dual-EnKF seems to To further assess the performance of the filters we analyze the spatial patterns of the estimated fields. To do so, we plot and interpret the ensemble mean of the conductivity as it results from the three filters using

D. Sensitivity to measurement errors
In the last set of sensitivity experiments, we fix the number of wells to nine, the observation frequency to 5 days, and we use different measurement errors to perturb the observations. We plot the results of nine different observational error scenarios in Figure 9 and compare the conductivity estimates obtained using the Joint-EnKF, Dual-EnKF and the Dual-EnKF OSA . In general, the performance of the filters appears to degrade as the observations are perturbed with larger degree of noise. All three filters exhibit similar With 0.10 m measurement error, the estimate of the Dual-EnKF OSA is also approximately 12% better.
Finally, we investigated the time-evolution of the ensemble variance of the conductivity estimates as they result from the Dual-EnKF and the Dual-EnKF OSA with 0.10 m-water measurement noise. Spatially, the ensemble variance maps provide insight about the uncertainty reduction due to data assimilation. The initial map (left panel, Figure 10) exhibits zero variance at the sampled two locations and increasing variance away from these locations. The ensemble spread of conductivity field from the two filters (right panels, Figure 10) after 6 and 18 months is quite small and comparable. The Dual-EnKF OSA , however, tends to maintain a larger variance at the edges than the Dual-EnKF, which in turn increases the impact of the observations.
E. Further assessments of the Dual-EnKF OSA scheme To further assess the system performance in terms of parameters retrieval, we have integrated the model in prediction mode (without assimilation) for an additional period of 18 months starting from the end of the assimilation period. We plot in Figure 11, using the final estimates of the conductivity as they

Dual-EnKF
Dual-EnKF-OSA Finally, in order to demonstrate that our results are statistically robust, 10 other test cases with different reference conductivity and heterogeneous recharge maps were investigated. In each of these cases, we sampled the reference fields by varying the variogram parameters, such as variance, x and y ranges, etc.
The pumping rates and the initial head configuration among the cases were also altered. For all 10 test cases, we fixed the ensemble size to 100 and used data from nine observation wells every 3 days. We set the measurement error to 0.10 m. We plot the resulting conductivity estimates (mean AAE) from each case in Figure 12. The estimates of the three filters, as shown on the plot, give a statistical evidence that the proposed scheme always provides more accurate estimates than the Joint-/Dual-EnKF and is more robust to changing dynamics and experimental setups. Similar results were obtained for the hydraulic head estimates. Averaging over all test cases, the proposed scheme provides about 17% more accurate estimates in term of AAE than the standard Joint-and Dual-EnKFs.   TC1  TC2  TC3  TC4  TC5  TC6  TC7  TC8  TC9  TC10 Mean AAE (log-m/s)

Conductivity Estimates for Different Test Cases
Joint-EnKF

Dual-EnKF
Dual-EnKF-OSA is inverted. Compared with the Dual-EnKF, this introduces a smoothing step to the state by future observations, which seems to provide the model, at the time of forecasting, with better and rather physically-consistent state and parameters ensembles.
We tested the proposed Dual-EnKF OSA on a conceptual groundwater flow model in which we estimated the hydraulic head and spatially variable conductivity parameters. We conducted a number of sensitivity experiments to evaluate the accuracy and the robustness of the proposed scheme and to compare its performance against those of the standard Joint and Dual EnKFs. The experimental results suggest that the Dual-EnKF OSA is more robust, successfully estimating the hydraulic head and the conductivity field under different modeling scenarios. Sensitivity analyses demonstrate that with the assimilation of more observations, the Dual-EnKF OSA becomes more effective and significantly outperforms the standard Jointand Dual-EnKF schemes. In addition, when using a sparse observation network in the aquifer domain, the accuracy of the Dual-EnKF OSA estimates is better preserved, unlike the Dual-EnKF, which seems to be more sensitive to the number of hydraulic wells. Moreover, the Dual-EnKF OSA results are shown to be more robust against observation noise. On average, the Dual-EnKF OSA scheme leads to around 10% more accurate state and parameters solutions than those resulting from the standard Joint-and Dual-EnKFs.
The proposed scheme is easy to implement and only requires minimal modifications to a standard EnKF code. It is further computationally feasible, requiring only a marginal increase in the computational cost compared to the Dual-EnKF. This scheme should therefore be beneficial to the hydrology community given its consistency, high accuracy, and robustness to changing modeling conditions. It could serve as an efficient estimation tool for real-world problems, such as groundwater, contaminant transport and reservoir monitoring, in which the available data are often sparse and noisy. Potential future research includes testing the Dual-EnKF OSA with realistic large-scale groundwater, contaminant transport and reservoir monitoring problems. Furthermore, integrating the proposed state-parameter estimation scheme with other iterative and hybrid ensemble approaches may be a promising direction for further improvements.

ACKNOWLEDGMENTS
Research reported in this publication was supported by King Abdullah University of Science and Technology (KAUST). Datasets and codes of the assimilation experiments can be accessed by contacting the second author: mohamad.gharamti@nersc.no.