Extended Kalman Filters, Least Squares Methods and Linear Successive Estimators for Inverse Problems in Geophysical Applications

1. Introduction

Reliable characterization of hydraulic conductivity is a central task in groundwater hydrology and subsurface flow modeling. Hydraulic conductivity controls the spatial distribution of hydraulic head, Darcy flux, solute transport pathways, and travel times, but it is rarely observed directly at the spatial resolution required by predictive groundwater models. Field measurements of conductivity are typically sparse, indirect, and scale-dependent, whereas hydraulic-head observations are more readily available. However, the hydraulic head value depends on conductivity through a nonlinear forward model, known as the richards equation. As a result, the reconstruction of hydraulic conductivity from hydraulic-head data is a classical inverse problem in hydrogeology, combining partial differential equation models, spatial statistics, optimization, and uncertainty quantification [1,2,3,4,5,6,7]. In saturated soil, the hydraulic head is commonly modeled by the steady groundwater flow equation, i.e. the saturated form of Richards’ equation in steady states, in which the hydraulic conductivity field enters as a spatially heterogeneous coefficient. Since conductivity is positive and often varies over several orders of magnitude, it is standard to estimate the logarithmic conductivity, or logarithmic transmissivity, as a random field. In unsaturated soil however, the inverse problem is more complex because the hydraulic conductivity depends on the water content or pressure head through constitutive relationships. The Richards equation couples pressure head, water content, and hydraulic conductivity in a nonlinear and often strongly degenerate manner in this case [8,9,10,11]. These features make the identification of hydraulic conductivity in both saturated and unsaturated soils ill-posed, sensitive to measurement error, and dependent on appropriate prior information. A class of approaches for hydrogeological inverse problems uses geostatistical priors to describe spatial variability and to regularize the estimation problem. Kriging and co-kriging provide best linear unbiased estimators based on covariance and cross-covariance information, allowing direct conductivity measurements and auxiliary hydraulic-head observations to be incorporated in a unified spatial estimator [6,12,13,14,15,16,17,18]. In this setting, the prior mean and covariance of the logarithmic conductivity field encode both geological structure and measurement information. Such priors are particularly useful in inverse problems because they reduce nonuniqueness and provide a statistically meaningful initial estimate for iterative reconstruction methods. Among iterative stochastic inverse methods, the successive linear estimator (SLE) introduced by [19] has been influential for estimating conditional effective transmissivity and hydraulic-head fields [20,21,22,23,24,25,26,27]. The essential idea is to begin with a geostatistical estimate of the log-conductivity field and then update it successively using discrepancies between measured and simulated hydraulic heads. At each iteration, the nonlinear map from conductivity to hydraulic head is linearized, and the update is computed from covariance and cross-covariance relationships between parameter perturbations and head perturbations. The resulting formula has the structure of a gain-times-innovation correction: the current estimate is adjusted by a linear combination of hydraulic-head residuals. This structure suggests a close connection between SLE and recursive filtering methods. Kalman filtering provides a natural mathematical framework for such recursive estimation. For linear dynamical systems with quadratic loss functions, the Kalman filter gives the optimal linear unbiased update and can be derived equivalently from minimum-variance estimation, Bayesian Gaussian conditioning, or recursive least-squares optimization [28,29,30,31,32]. For nonlinear systems, the extended Kalman filter (EKF) applies the same prediction-correction mechanism after local linearization of the nonlinear state and observation maps [33,34]. In static parameter-estimation problems, the parameter can be treated as a state with identity dynamics, and the observation operator is the nonlinear forward map from parameters to measured quantities. This viewpoint has been widely used in subsurface data assimilation and groundwater model calibration, including ensemble and sequential extensions for joint state-parameter estimation [35,36,37,38]. The purpose of this paper is to establish the connection between SLE and EKF explicitly in the setting of hydraulic-conductivity estimation. We formulate the successive linear estimation of the hydraulic conductivity from measured hydraulic-head data as a parameter-estimation problem within the framework of extended Kalman filtering and least-squares optimization. First, we review the Kalman filter and the EKF from the perspective of weighted least-squares minimization. In particular, the Kalman correction can be interpreted as a one-step Newton update applied to a quadratic least-squares objective, while the EKF corresponds to applying this same update after linearizing the nonlinear observation operator. This least-squares viewpoint clarifies the roles of the prior covariance, observation covariance, sensitivity matrix, and innovation term. We then show that, for saturated soil, the SLE update is precisely an EKF-type correction for a static parameter-estimation problem in which the state variable is the logarithmic conductivity field and the observation operator is the groundwater flow solver mapping log-conductivity to hydraulic head. The co-kriging estimate supplies the initial prior state and covariance, sensitivity of the hydraulic head with respect to log-conductivity supplies the EKF Jacobian, and the SLE weights correspond to the Kalman gain. Thus, SLE can be interpreted not only as a geostatistical iterative estimator, but also as a particular implementation of the extended Kalman filter with a geostatistical initialization and a linearized least-squares update. Finally, we extend the same formulation to unsaturated soil to provide a framework for estimating hydraulic conductivity using spatial-temporal data. In this case, the observation operator is defined by the transient Richards equation that depends on the nonlinear relationship among pressure head, water content, and hydraulic conductivity. By parameterizing the unknown conductivity function or its logarithm and applying the EKF to the resulting static or pseudo-dynamic parameter-estimation problem, one obtains an SLE-type update for unsaturated flow. This formulation places saturated and unsaturated conductivity estimation in a common recursive least-squares framework and provides a basis for incorporating measurement noise, regularization, covariance modeling, and nonlinear sensitivity information in a systematic way. The remainder of the paper is organized as follows. § 2 reviews the Kalman filter, its least-squares interpretation, and the EKF for nonlinear systems. § 3 outlines the SLE formulation for saturated soil, beginning with kriging and co-kriging priors and then showing the equivalence between SLE and an EKF update. § 4 extends the formulation to unsaturated soil governed by the time-dependent Richards’ equation. § 5 summarizes the main conclusions and discusses implications for stochastic inverse modeling in groundwater applications.

2. Kalman Filters and Least Squares Methods

We review the standard Kalman filter for a discrete linear system identification problem nd its extension to a nonlinear system in the context of weighted least squares methods.

2.1. Kalman Filters for Discrete Linear Systems

The Kalman filter is an optimal linear state estimation algorithm in the context of unbiased estimators. It fuses a prediction model with observation data to estimate the true state of a dynamic system in the presence of noise and uncertainty [28]. We use the following discrete, linear time invariant (LTI) system, as an example.

$$\begin{array}{c}\hfill \begin{array}{c}\begin{array}{cc}\hfill \mathrm{state}\mathrm{equation}\text{:}{\mathbf{x}}_k& =\mathbf{A}{\mathbf{x}}_{k-1}+\mathbf{B}{\mathbf{u}}_k+{\mathbf{w}}_k,\hfill \\ \hfill \mathrm{observational}\mathrm{equation}\text{:}{\mathbf{z}}_k& =\mathbf{H}{\mathbf{x}}_k+{\mathbf{v}}_k,k=1,\dots ,\hfill \end{array}\hfill \end{array}\end{array}$$

(2.1)

where ${\mathbf{x}}_k\in {\mathbb{R}}^{n\times 1}$ denotes the state vector of the dynamic system at time step k, which is the true n-dimensional state to be estimated. $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ serves as the state transition matrix, illustrating how the state evolves from time step $k-1$ to k in the absence of control inputs, ${\mathbf{u}}_k$, and process noise, ${\mathbf{w}}_k$. $\mathbf{B}\in {\mathbb{R}}^{n\times m}$ is the control input matrix, mapping the m-dimensional control input, ${\mathbf{u}}_k\in {\mathbb{R}}^{m\times 1}$, to its impact on the state transition in ${\mathbf{x}}_{k-1}\to {\mathbf{x}}_k$ in ${\mathbf{R}}^n$. ${\mathbf{z}}_k\in {\mathbb{R}}^{p\times 1}$ represents the observation (measurement) vector at time step k, acquired from sensors or other measurement devices, where p is the number of observations. $\mathbf{H}\in {\mathbb{R}}^{p\times n}$ is the observation matrix, establishing the relationship between the n-dimensional system state ${\mathbf{x}}_k$ and the p-dimensional observation ${\mathbf{z}}_k$. ${\mathbf{u}}_k\in {\mathbb{R}}^{m\times 1}$ is the known control input, ${\mathbf{w}}_k\in {\mathbb{R}}^{n\times 1}\sim \mathcal{N}(0,{\mathbf{Q}}_k)$ is the process noise with a zero mean and covariance ${\mathbf{Q}}_k=\mathbb{E}\left[{\mathbf{w}}_k{\mathbf{w}}_k^T\right]\in {\mathbb{R}}^{n\times n}$, and ${\mathbf{v}}_k\in {\mathbb{R}}^{p\times 1}\sim \mathcal{N}(0,{\mathbf{R}}_k)$ is the measurement noise with a zero mean and covariance ${\mathbf{R}}_k=\mathbb{E}\left[{\mathbf{v}}_k{\mathbf{v}}_k^T\right]\in {\mathbb{R}}^{p\times p}$. Assume that ${\mathbf{w}}_k$ and ${\mathbf{v}}_k$ are independent Gaussian white noises and uncorrelated with each other. We note that $n,m$ are independent positive integers. However, $n,m$ and p are not completely independent. If $p+m<n$, i.e., the number of control variables and the measurements is smaller than the dimension of the state vector to be estimated, there could be some problems in identifying the state properly. In fact, the state estimation problem of system (1) may not have a solution or unique solution of ${\mathbf{x}}_k$, given ${\mathbf{x}}_{k-1},{\mathbf{u}}_k,{\mathbf{z}}_k,{\mathbf{w}}_k,{\mathbf{v}}_k$. We also need to supply an initial state

$$\begin{array}{c}\hfill {\mathbf{x}}_0={\mathbf{x}}_0^{-}+{\mathbf{v}}_0,\end{array}$$

(2.2)

where ${\mathbf{x}}_0^{-}$ is a given initial value and ${\mathbf{v}}_0$ the noise. Let ${\widehat{\mathbf{x}}}_k^{-}$ denote the a priori state estimate at time k, i.e., the prediction before incorporating the observation ${\mathbf{z}}_k$, which we will be specified later. Given ${\widehat{\mathbf{x}}}^{-}$, the posterior estimate at time k is

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{\mathbf{x}}}_k={\widehat{\mathbf{x}}}_k^{-}+{\mathbf{K}}_k({\mathbf{z}}_k-\mathbf{H}{\widehat{\mathbf{x}}}_k^{-}),\hfill \end{array}\end{array}$$

(2.3)

where ${\mathbf{K}}_k$ is the Kalman gain given below. The objective is to find the Kalman gain matrix, ${\mathbf{K}}_k$, such that the estimate ${\widehat{\mathbf{x}}}_k$ approaches the true state ${\mathbf{x}}_k$ as an unbiased estimator, i.e., $\mathbb{E}\left[{\mathbf{x}}_k\right]=\mathbb{E}\left[{\widehat{\mathbf{x}}}_k\right]$. Define the estimation error as

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{e}}_k={\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k\hfill \\ ={\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k^{-}-{\mathbf{K}}_k(\mathbf{H}{\mathbf{x}}_k+{\mathbf{v}}_k-\mathbf{H}{\widehat{\mathbf{x}}}_k^{-})\hfill \\ =(\mathbf{I}-{\mathbf{K}}_k\mathbf{H})({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k^{-})-{\mathbf{K}}_k{\mathbf{v}}_k\hfill \\ =(\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{e}}_k^{-}-{\mathbf{K}}_k{\mathbf{v}}_k,\hfill \end{array}\end{array}$$

(2.4)

where ${\mathbf{e}}_k^{-}={\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k^{-}.$ Then, $\mathbb{E}\left[{\mathbf{e}}_k\right]=0$ and

$$\begin{array}{c}\hfill \begin{array}{c}\mathbb{E}\left[{\mathbf{e}}_k{\mathbf{e}}_k^T\right]=\mathbb{E}\left[((\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{e}}_k^{-}-{\mathbf{K}}_k{\mathbf{v}}_k){((\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{e}}_k^{-}-{\mathbf{K}}_k{\mathbf{v}}_k)}^T\right]\hfill \\ =(\mathbf{I}-{\mathbf{K}}_k\mathbf{H})E\left[{\mathbf{e}}_k^{-}{\mathbf{e}}_k^{-T}\right]{(\mathbf{I}-{\mathbf{K}}_k\mathbf{H})}^T+{\mathbf{K}}_{k}E\left[{\mathbf{v}}_k{\mathbf{v}}_k^T\right]{\mathbf{K}}_k^T\hfill \\ =(\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{P}}_k^{-}{(\mathbf{I}-{\mathbf{K}}_k\mathbf{H})}^T+{\mathbf{K}}_k{\mathbf{R}}_k{\mathbf{K}}_k^T,\hfill \end{array}\end{array}$$

(2.5)

where

$${\mathbf{P}}_k^{-}=\mathbb{E}\left[{\mathbf{e}}_k^{-}{\mathbf{e}}_k^{-T}\right],{\mathbf{e}}_k^{-}={\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k^{-}.$$

The objective is to minimize the trace of the expected error covariance matrix:

$$\underset{{\mathbf{K}}_k}{min}\mathrm{tr}\left(\mathbb{E}\left[{\mathbf{e}}_k{\mathbf{e}}_k^T\right]\right),$$

which leads to the optimal gain ${\mathbf{K}}_k$.

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial {\mathbf{K}}_k}}tr\left(\mathbb{E}\left[{\mathbf{e}}_k{\mathbf{e}}_k^T\right]\right)=2({\mathbf{K}}_k\mathbf{H}{\mathbf{P}}_k^{-}{\mathbf{H}}^T+{\mathbf{K}}_k{\mathbf{R}}_k)-2{\mathbf{P}}_k^T{\mathbf{H}}^T=0.\end{array}$$

(2.6)

This implies

$$\begin{array}{c}\hfill {\mathbf{K}}_k={\mathbf{P}}_k^{-}{\mathbf{H}}^T{(\mathbf{H}{\mathbf{P}}_k^{-}{\mathbf{H}}^T+{\mathbf{R}}_k)}^{-1}.\end{array}$$

(2.7)

This is the Kalman gain. At this value,

$$\begin{array}{c}\hfill \begin{array}{c}\mathbb{E}\left[{\mathbf{e}}_k{\mathbf{e}}_k^T\right]=(\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{P}}_k^{-}{(\mathbf{I}-{\mathbf{K}}_k\mathbf{H})}^T+{\mathbf{K}}_k{\mathbf{R}}_k{\mathbf{K}}_k^T\hfill \\ =(\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{P}}_k^{-},\hfill \end{array}\end{array}$$

(2.8)

${\mathbf{P}}_k^{-}$ is the a priori estimation error covariance matrix. Here, we define

$$\begin{array}{c}\hfill {\mathbf{P}}_k=\mathbb{E}\left[{\mathbf{e}}_k{\mathbf{e}}_k^T\right]=(\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{P}}_k^{-}.\end{array}$$

(2.9)

We want to have the guess close to the state. So, we use

$$\begin{array}{c}\hfill {\widehat{\mathbf{x}}}_k^{-}=\mathbf{A}{\widehat{\mathbf{x}}}_{k-1}+\mathbf{B}{\mathbf{u}}_k.\end{array}$$

(2.10)

$$\begin{array}{c}\hfill {\mathbf{e}}_k^{-}={\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k^{-}=\mathbf{A}{\mathbf{e}}_{k-1}+{\mathbf{w}}_k.\end{array}$$

(2.11)

With this choice, we have a formula for ${\mathbf{P}}_k^{-}$

$$\begin{array}{cc}\hfill {\mathbf{P}}_k^{-}& =\mathbf{A}{\mathbf{P}}_{k-1}{\mathbf{A}}^T+{\mathbf{Q}}_k.\end{array}$$

(2.12)

The Kalman filter can be written into a prediction-correction scheme as follows.

Algorithm 2.1

(Kalman Filter). Given an initial estimate, ${\widehat{\mathbf{x}}}_0$, and initial error covariance, ${\mathbf{P}}_0$, the Kalman filter proceeds as follows for each time step $k=1,2,\dots$

1. Prediction

Given the posterior estimate ${\widehat{\mathbf{x}}}_{k-1}$ and covariance ${\mathbf{P}}_{k-1}$, the prediction step computes:

$$\begin{array}{cc}\hfill \mathrm{A}\mathrm{priori}\mathrm{state}\mathrm{estimate}\text{:}{\widehat{\mathbf{x}}}_k^{-}& =\mathbf{A}{\widehat{\mathbf{x}}}_{k-1}+\mathbf{B}{\mathbf{u}}_k,\hfill \end{array}$$

(2.13)

$$\begin{array}{cc}\hfill \mathrm{A}\mathrm{priori}\mathrm{error}\mathrm{covariance}\text{:}{\mathbf{P}}_k^{-}& =\mathbf{A}{\mathbf{P}}_{k-1}{\mathbf{A}}^T+{\mathbf{Q}}_k.\hfill \end{array}$$

(2.14)

2. Correction (Update)

$$\begin{array}{cc}\hfill \mathrm{Kalman}\mathrm{gain}\text{:}{\mathbf{K}}_k& ={\mathbf{P}}_k^{-}{\mathbf{H}}^T{\left(\mathbf{H}{\mathbf{P}}_k^{-}{\mathbf{H}}^T+{\mathbf{R}}_k\right)}^{-1},\hfill \end{array}$$

(2.15)

$$\begin{array}{cc}\hfill \mathrm{Measurement}\mathrm{residual}\text{:}{\mathbf{y}}_k& ={\mathbf{z}}_k-\mathbf{H}{\widehat{\mathbf{x}}}_k^{-},\hfill \end{array}$$

(2.16)

$$\begin{array}{cc}\hfill \mathrm{A}\mathrm{posteriori}\mathrm{state}\mathrm{estimate}\text{:}{\widehat{\mathbf{x}}}_k& ={\widehat{\mathbf{x}}}_k^{-}+{\mathbf{K}}_k{\mathbf{y}}_k,\hfill \end{array}$$

(2.17)

$$\begin{array}{cc}\hfill \mathrm{A}\mathrm{posteriori}\mathrm{error}\mathrm{covariance}\text{:}{\mathbf{P}}_k& =(\mathbf{I}-{\mathbf{K}}_k\mathbf{H}){\mathbf{P}}_k^{-}.\hfill \end{array}$$

(2.18)

Lemma 1.

Given $f\left(\mathbf{x}\right)=f_1\left(\mathbf{x}\right)+f_2\left(\mathbf{x}\right)$, where $f_1$ and $f_2$ are quadratic functions,

$$\begin{array}{c}\hfill {\mathbf{x}}_m=arg\underset{\mathbf{x}}{min}f\left(\mathbf{x}\right)={\mathbf{x}}^{*}-{\left({\nabla }^{2}f\left({\mathbf{x}}^{*}\right)\right)}^{-1}\nabla f\left({\mathbf{x}}^{*}\right),\end{array}$$

(2.19)

for any ${\mathbf{x}}^{*}$. If we choose ${\mathbf{x}}^{*}$ such that $\nabla f_1\left({\mathbf{x}}^{*}\right)=0$,

$$\begin{array}{c}\hfill {\mathbf{x}}_m=arg\underset{\mathbf{x}}{min}f\left(\mathbf{x}\right)={\mathbf{x}}^{*}-{\left({\nabla }^{2}f\left({\mathbf{x}}^{*}\right)\right)}^{-1}\nabla f_2\left({\mathbf{x}}^{*}\right).\end{array}$$

(2.20)

This means that one can optimize $f_1\left(\mathbf{x}\right)$ firstly. The overall optimization can be obtained using the above one-step Newton’s formula. In general,

$$\begin{array}{c}\hfill F_m\left(\mathbf{x}\right)={\displaystyle \sum _{i=1}^m}{f}_i\left(\mathbf{x}\right),\end{array}$$

(2.21)

where $f_i\left(\mathbf{x}\right)$ are quadratic functions. The optimizer of $F_m\left(\mathbf{x}\right)$ can be obtained sequentially from the optimizer of $F_{m-1}\left(\mathbf{x}\right)$, ${\mathbf{x}}_{m-1}$, via

$$\begin{array}{c}\hfill {\mathbf{x}}_m={\mathbf{x}}_{m-1}-{\left[{\nabla }^2{F}_m\left({\mathbf{x}}_{m-1}\right)\right]}^{-1}\nabla f_m\left({\mathbf{x}}_{m-1}\right).\end{array}$$

(2.22)

With this lemma, we reformulate the Kalman filter problem into a least square minimization problem that can be carried out sequentially. We denote

$$\begin{array}{c}\hfill {\mathbf{Z}}_k={({\mathbf{x}}_0,\dots ,{\mathbf{x}}_k)}^T.\end{array}$$

(2.23)

i.e. ${\mathbf{Z}}_k$ contains all the state variables from the initial time 0 to the current time k. Then, we define the objective function as follows

$$\begin{array}{c}\hfill J\left({\mathbf{Z}}_k\right)={\displaystyle \sum _{i=1}^k}{\displaystyle \frac{1}{2}}\left[\parallel {\mathbf{x}}_i-\mathbf{A}{\mathbf{x}}_{i-1}-\mathbf{B}{\mathbf{u}}_i{\parallel }_{{{\mathbf{Q}}_i}^{-1}}^2+{\parallel {\mathbf{z}}_i-\mathbf{H}{\mathbf{x}}_i\parallel }_{{{\mathbf{R}}_i}^{-1}}^2\right]+{\displaystyle \frac{1}{2}}{\parallel {\mathbf{x}}_0-{\mu }_0\parallel }_{{{\mathbf{Q}}_0}^{-1}}^2.\end{array}$$

(2.24)

${\mathbf{x}}_0$ denotes the initial state and

$$\begin{array}{c}\hfill {\parallel \mathbf{x}\parallel }_{{\mathbf{Q}}^{-1}}=\sqrt{{\mathbf{x}}^T·{\mathbf{Q}}^{-1}·\mathbf{x}},\end{array}$$

(2.25)

with a positive definite matrix $\mathbf{Q}>0$. It is assumed that the true value of ${\mathbf{x}}_0$ follows a distribution around a certain known prior mean ${\mu }_0$, and its uncertainty is characterized by the covariance matrix ${\mathbf{Q}}_0$.

Based on the recursiveness of J, we have

$$\begin{array}{c}\hfill J\left({\mathbf{Z}}_k\right)=J\left({\mathbf{Z}}_{k-1}\right)+{\displaystyle \frac{1}{2}}\left[\parallel {\mathbf{x}}_k-\mathbf{A}{\mathbf{x}}_{k-1}-\mathbf{B}{\mathbf{u}}_k{\parallel }_{{{\mathbf{Q}}_k}^{-1}}^2+{\parallel {\mathbf{z}}_k-\mathbf{H}{\mathbf{x}}_k\parallel }_{{{\mathbf{R}}_k}^{-1}}^2\right].\end{array}$$

(2.26)

Its gradient with respect to ${\mathbf{Z}}_k$ is

$$\begin{array}{c}\hfill \nabla J\left({\mathbf{Z}}_k\right)=\left(\begin{array}{c}\nabla J\left({\mathbf{Z}}_{k-1}\right)+{(0,\dots ,0,\mathbf{A})}^T{\mathbf{Q}}_k^{-1}\left(\mathbf{A}{\mathbf{x}}_{k-1}-{\mathbf{x}}_k+\mathbf{B}{\mathbf{u}}_k\right)\hfill \\ -{\mathbf{Q}}_k^{-1}\left(\mathbf{A}{\mathbf{x}}_{k-1}-{\mathbf{x}}_k+\mathbf{B}{\mathbf{u}}_k\right)+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}(\mathbf{H}{\mathbf{x}}_k-{\mathbf{z}}_k)\hfill \end{array}\right),\end{array}$$

(2.27)

and its Hessian is

$${\nabla }^{2}J\left({\mathbf{Z}}_k\right)=\left(\begin{array}{cc}{\nabla }^{2}J\left({\mathbf{Z}}_{k-1}\right)+{(0,\dots ,0,\mathbf{A})}^T{\mathbf{Q}}_k^{-1}(0,\dots ,0,\mathbf{A})\hfill & \hfill -{(0,\dots ,0,\mathbf{A})}^T{\mathbf{Q}}_k^{-1}\\ -{\mathbf{Q}}_k^{-1}(0,\dots ,0,\mathbf{A})\hfill & \hfill {\mathbf{Q}}_k^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H}\end{array}\right).$$

(2.28)

The least square solution of the quadratic function is given by the one-step Newton’s formula

$$\begin{array}{c}\hfill {\mathbf{Z}}_k={\mathbf{Z}}_k^{*}-{\left[{\nabla }^{2}J\left({\mathbf{Z}}_k^{*}\right)\right]}^{-1}\nabla J\left({\mathbf{Z}}_k^{*}\right),\end{array}$$

(2.29)

where ${\mathbf{Z}}_k^{*}$ is arbitrary.

Specifically, we choose

$$\begin{array}{c}\hfill {\mathbf{Z}}_k^{*}=\left(\begin{array}{c}{\widehat{\mathbf{Z}}}_{k-1}\hfill \\ {\widehat{\mathbf{x}}}_k^{*}\hfill \end{array}\right),\end{array}$$

(2.30)

where ${\widehat{\mathbf{x}}}_k^{*}$, a prior, satisfies

$$\begin{array}{c}\hfill {\widehat{\mathbf{x}}}_k^{*}-\mathbf{A}{\widehat{\mathbf{x}}}_{k-1}-\mathbf{B}{\mathbf{u}}_k=0,\end{array}$$

(2.31)

and

$$\begin{array}{c}\hfill \nabla J\left({\widehat{\mathbf{Z}}}_{k-1}\right)=0.\end{array}$$

(2.32)

With this choice, the gradient reduces to

$$\begin{array}{c}\hfill \nabla J\left({\mathbf{Z}}_k^{*}\right)=\left(\begin{array}{c}0\hfill \\ {\mathbf{H}}^T{\mathbf{R}}_k^{-1}(\mathbf{H}{\widehat{\mathbf{x}}}_k^{*}-{\mathbf{z}}_k)\hfill \end{array}\right),\end{array}$$

(2.33)

and the Hessian is given by

$${\nabla }^{2}J\left({\widehat{\mathbf{Z}}}_k^{*}\right)=\left(\begin{array}{cc}{\nabla }^{2}J\left({\widehat{\mathbf{Z}}}_{k-1}\right)+{(0,\dots ,0,\mathbf{A})}^T{\mathbf{Q}}_k^{-1}(0,\dots ,0,\mathbf{A})\hfill & \hfill -{(0,\dots ,0,\mathbf{A})}^T{\mathbf{Q}}_k^{-1}\\ -{\mathbf{Q}}_k^{-1}(0,\dots ,0,\mathbf{A})\hfill & \hfill {\mathbf{Q}}_k^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H}\end{array}\right).$$

(2.34)

We then apply the one-step Newton update

$${\mathbf{Z}}_k={\mathbf{Z}}_k^{*}-{\left({\nabla }^{2}J\left({\mathbf{Z}}_k^{*}\right)\right)}^{-1}\nabla J\left({\mathbf{Z}}_k^{*}\right)$$

yields the following block system

$$\left[\begin{array}{cc}{\nabla }^{2}J\left({\widehat{\mathbf{Z}}}_{k-1}\right)+{\mathbf{A}}_{\mathrm{off}}^T{\mathbf{Q}}_k^{-1}{\mathbf{A}}_{\mathrm{off}}& -{\mathbf{A}}_{\mathrm{off}}^T{\mathbf{Q}}_k^{-1}\\ -{\mathbf{Q}}_k^{-1}{\mathbf{A}}_{\mathrm{off}}& {\mathbf{Q}}_k^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H}\end{array}\right]\left[\begin{array}{c}\Delta {\mathbf{Z}}_{\mathrm{hist}}\\ \Delta {\mathbf{x}}_k\end{array}\right]=-\left[\begin{array}{c}0\\ {\mathbf{H}}^T{\mathbf{R}}_k^{-1}(\mathbf{H}{\mathbf{x}}_k^{*}-{\mathbf{z}}_k)\end{array}\right],$$

(2.35)

where ${\mathbf{A}}_{\mathrm{off}}=(0,\dots ,0,\mathbf{A})$ denotes the off-diagonal block of the state transition.

From the first block row of this system, we have

$$({\nabla }^{2}J\left({\widehat{\mathbf{Z}}}_{k-1}\right)+{\mathbf{A}}_{\mathrm{off}}^T{\mathbf{Q}}_k^{-1}{\mathbf{A}}_{\mathrm{off}})\Delta {\mathbf{Z}}_{\mathrm{hist}}-{\mathbf{A}}_{\mathrm{off}}^T{\mathbf{Q}}_k^{-1}\Delta {\mathbf{x}}_k=0.$$

(2.36)

Solving for $\Delta {\mathbf{Z}}_{\mathrm{hist}}$ gives

$$\Delta {\mathbf{Z}}_{\mathrm{hist}}={[{\nabla }^{2}J\left({\widehat{\mathbf{Z}}}_{k-1}\right)+{\mathbf{A}}_{\mathrm{off}}^T{Q}_k^{-1}{\mathbf{A}}_{\mathrm{off}}]}^{-1}{\mathbf{A}}_{\mathrm{off}}^T{\mathbf{Q}}_k^{-1}\Delta {\mathbf{x}}_k.$$

(2.37)

From the second block row of (35), we have

$$-{\mathbf{Q}}_k^{-1}{\mathbf{A}}_{\mathrm{off}}\Delta {\mathbf{Z}}_{\mathrm{hist}}+({\mathbf{Q}}_k^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H})\Delta {\mathbf{x}}_k=-{\mathbf{H}}^T{\mathbf{R}}_k^{-1}(\mathbf{H}{\mathbf{x}}_k^{*}-{\mathbf{z}}_k).$$

(2.38)

By eliminating $\Delta {\mathbf{Z}}_{\mathrm{hist}}$ via the Schur complement, the effective Hessian for $\Delta {\mathbf{x}}_k$ becomes

$$\begin{array}{cc}\hfill {\mathbf{S}}_k& =({\mathbf{Q}}_k^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H})-{\mathbf{Q}}_k^{-1}{\mathbf{A}}_{\mathrm{off}}{[{\nabla }^{2}J\left({\widehat{\mathbf{Z}}}_{k-1}\right)+{\mathbf{A}}_{\mathrm{off}}^T{\mathbf{Q}}_k^{-1}{\mathbf{A}}_{\mathrm{off}}]}^{-1}{\mathbf{A}}_{\mathrm{off}}^T{\mathbf{Q}}_k^{-1}\hfill \\ \hfill & ={\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H}+{({\mathbf{Q}}_k+{\mathbf{A}}_{\mathrm{off}}{\left({\nabla }^{2}J\left({\widehat{\mathbf{Z}}}_{k-1}\right)\right)}^{-1}{\mathbf{A}}_{\mathrm{off}}^T)}^{-1}.\hfill \end{array}$$

(2.39)

Note that ${\mathbf{A}}_{\mathrm{off}}{\left({\nabla }^2{J}_{k-1}\right)}^{-1}{\mathbf{A}}_{\mathrm{off}}^T$ corresponds precisely to the predicted prior error covariance of the previous time step. Adopting the standard Kalman filter notation, this term is identified as $\mathbf{A}{\mathbf{P}}_{k-1}{\mathbf{A}}^T$. Consequently, the composite term $({\mathbf{Q}}_k+\mathbf{A}{\mathbf{P}}_{k-1}{\mathbf{A}}^T)$ is denoted as the predicted covariance ${\mathbf{P}}_k^{-}$. Thus, the Schur complement ${\mathbf{S}}_k$ simplifies to

$${\mathbf{S}}_k={\left({\mathbf{P}}_k^{-}\right)}^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H},$$

(2.40)

where ${\mathbf{P}}_k^{-}$ denotes the predicted error covariance. Consequently, the Newton update $\Delta {\mathbf{x}}_k$ is given by

$$\Delta {\mathbf{x}}_k={\mathbf{S}}_k^{-1}{\mathbf{H}}^T{\mathbf{R}}_k^{-1}({\mathbf{z}}_k-\mathbf{H}{\mathbf{x}}_k^{*})={[{\left({\mathbf{P}}_k^{-}\right)}^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H}]}^{-1}{\mathbf{H}}^T{\mathbf{R}}_k^{-1}({\mathbf{z}}_k-\mathbf{H}{\mathbf{x}}_k^{*}).$$

(2.41)

Applying the SMW identity again to the gain term, we obtain the standard Kalman gain form

$${\mathbf{K}}_k={[{\left({\mathbf{P}}_k^{-}\right)}^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H}]}^{-1}{\mathbf{H}}^T{\mathbf{R}}_k^{-1}={\mathbf{P}}_k^{-}{\mathbf{H}}^T{(\mathbf{H}{\mathbf{P}}_k^{-}{\mathbf{H}}^T+{\mathbf{R}}_k)}^{-1}.$$

(2.42)

Therefore, the update equation becomes

$$\Delta {\mathbf{x}}_k={\mathbf{K}}_k({\mathbf{z}}_k-\mathbf{H}{\mathbf{x}}_k^{*}),$$

(2.43)

which is exactly the Kalman filter measurement update. Since the cost function $J\left({\mathbf{Z}}_k\right)$ is quadratic, the Newton step converges to the global minimum in a single iteration, making the derivation exact. The historical states ${\mathbf{Z}}_{\mathrm{hist}}$ do not require updates because the chosen initial point ${\widehat{\mathbf{Z}}}_{k-1}$ already satisfies the optimality condition for the past subsystem i.e., $\nabla J\left({\widehat{\mathbf{Z}}}_{k-1}\right)=0$.

However, substituting this back into the full Newton step, or recognizing that the objective function $J\left(Z_k\right)$ is quadratic, we observe that the update for the historical states $\Delta {\mathbf{Z}}_{\mathrm{hist}}$ does not need to be explicitly calculated for the recursion. The key insight is that the current state update $\Delta {\mathbf{x}}_k$ can be computed independently by considering the Schur complement of the Hessian, leading to the standard Kalman gain formulation.

Specifically, while a non-zero innovation $({\mathbf{z}}_k-\mathbf{H}{\widehat{\mathbf{x}}}_k^{*})$ generally induces a correction term $\Delta {\mathbf{Z}}_{\mathrm{hist}}$—a process characteristic of Kalman smoothers that refines historical trajectories—the standard Kalman filter operates under the Markov assumption. This assumption implies that for real-time recursion, the estimate ${\widehat{\mathbf{Z}}}_{k-1}$ is fixed, effectively imposing the constraint $\Delta {\mathbf{Z}}_{\mathrm{hist}}=0$.

This constraint does not imply that historical estimates are error-free; rather, it signifies that their uncertainty, characterized by the covariance ${\mathbf{P}}_{k-1}$, is propagated to the current step via the system dynamics matrix $\mathbf{A}$ instead of being retroactively optimized. Consequently, the optimization problem is reduced to minimizing a cost function $J_k\left({\mathbf{x}}_k\right)$ solely with respect to the current state ${\mathbf{x}}_k$, treating the predicted state ${\widehat{\mathbf{x}}}_k^{-}=\mathbf{A}{\widehat{\mathbf{x}}}_{k-1}+\mathbf{B}{\mathbf{u}}_k$ as a prior measurement with covariance ${\mathbf{P}}_k^{-}=\mathbf{A}{\mathbf{P}}_{k-1}{\mathbf{A}}^T+{\mathbf{Q}}_k$. The resulting cost function is given by

$$J_k\left({\mathbf{x}}_k\right)={\displaystyle \frac{1}{2}}\parallel {\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k^{-}{\parallel }_{{\left({\mathbf{P}}_k^{-}\right)}^{-1}}^2+{\displaystyle \frac{1}{2}}{\parallel \mathbf{H}{\mathbf{x}}_k-{\mathbf{z}}_k\parallel }_{{\mathbf{R}}_k^{-1}}^2.$$

(2.44)

Note that the weight of the prediction term is determined by the inverse of the predicted covariance ${\left({\mathbf{P}}_k^{-}\right)}^{-1}$, rather than the process noise covariance ${\mathbf{Q}}_k^{-1}$, as the former accounts for the accumulated uncertainty from previous steps. Calculating the gradient and Hessian of this reduced cost function yields:

$$\begin{array}{cc}\hfill \nabla J_k\left({\mathbf{x}}_k\right)& ={\left({\mathbf{P}}_k^{-}\right)}^{-1}({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k^{-})+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}(\mathbf{H}{\mathbf{x}}_k-{\mathbf{z}}_k),\hfill \end{array}$$

(2.45)

$$\begin{array}{cc}\hfill {\nabla }^2{J}_k\left({\mathbf{x}}_k\right)& ={\left({\mathbf{P}}_k^{-}\right)}^{-1}+{\mathbf{H}}^T{\mathbf{R}}_k^{-1}\mathbf{H}={\mathbf{S}}_k.\hfill \end{array}$$

(2.46)

By initializing the guess at the predicted state ${\mathbf{x}}_k^{*}={\widehat{\mathbf{x}}}_k^{-}$, the Newton step becomes

$$\Delta {\mathbf{x}}_k=-{\left[{\nabla }^2{J}_k\left({\mathbf{x}}_k\right)\right]}^{-1}\nabla J_k\left({\mathbf{x}}_k^{*}\right)={\mathbf{S}}_k^{-1}{\mathbf{H}}^T{\mathbf{R}}_k^{-1}({\mathbf{z}}_k-\mathbf{H}{\widehat{\mathbf{x}}}_k^{-}),$$

(2.47)

which is mathematically identical to the standard Kalman filter update derived in (43).

Proposition 1.

The Kalman filer for a linear system, (1), is given by the solution to the following equivalent least squares minimization problems:

1.

$$\begin{array}{c}\hfill \underset{{\mathbf{Z}}_k}{min}J\left({\mathbf{Z}}_k\right),\end{array}$$

(2.48)

where $J\left({\mathbf{Z}}_k\right)$ is given in (24);

2.

$$\begin{array}{c}\hfill \underset{{\mathbf{x}}_k}{min}{J}_k\left({\mathbf{x}}_k\right),\end{array}$$

(2.49)

where $J_k\left({\mathbf{Z}}_k\right)$ is given in (44).

2.2. Extended Kalman Filters for discrete Nonlinear Systems

The Extended Kalman filter (EKF) is an extension of the Kalman Filter (KF) to nonlinear systems. When system state equations or observation equations are nonlinear, the linear Kalman Filter alluded to above is no longer applicable. Thus, one proposed EKF. It approximates the nonlinear system as a linear system through a local linearization approach, thereby retaining the recursive framework of the Kalman Filter and achieving optimal estimation (minimum mean squared error) of the state of nonlinear dynamic systems[33].

We once again use a discrete nonlinear system as an example. In the discrete nonlinear system, the state and measurement equations are given as follows

$$\begin{array}{c}\hfill \begin{array}{c}\begin{array}{cc}\hfill {\mathbf{x}}_k& =f({\mathbf{x}}_{k-1},{\mathbf{u}}_k,{\mathbf{w}}_k),\hfill \\ \hfill {\mathbf{z}}_k& =h({\mathbf{x}}_k,{\mathbf{v}}_k),\hfill \end{array}\hfill \end{array}\end{array}$$

(2.50)

where ${\mathbf{x}}_k\in {\mathbb{R}}^{n\times 1}$ is the true state vector at step k, ${\mathbf{u}}_k\in {\mathbb{R}}^{m\times 1}$ is the known control input vector, ${\mathbf{w}}_k\in {\mathbb{R}}^{q\times 1}$ is the process noise, assumed to be zero-mean Gaussian with covariance ${\mathbf{Q}}_k\in {\mathbb{R}}^{q\times q}$, i.e., ${\mathbf{w}}_k\sim \mathcal{N}(0,{\mathbf{Q}}_k)$, ${\mathbf{z}}_k\in {\mathbb{R}}_k^{p\times 1}$ is the observed measurement vector, ${\mathbf{v}}_k\in {\mathbb{R}}_k^{r\times 1}$: the measurement noise, assumed to be zero-mean Gaussian with covariance ${\mathbf{R}}_k\in {\mathbb{R}}^{r\times r}$, i.e., ${\mathbf{v}}_k\sim \mathcal{N}(0,{\mathbf{R}}_k)$, $f:{\mathbb{R}}^n\times {\mathbb{R}}^m\times {\mathbb{R}}^q\to {\mathbb{R}}^n$ is the nonlinear state transition function, $h:{\mathbb{R}}^n\times {\mathbb{R}}^r\to {\mathbb{R}}^p$ is the nonlinear observation function.

To apply the EKF, the nonlinear system is linearized using first-order Taylor expansion around the current state estimate ${\widehat{\mathbf{x}}}_{k-1}$. Let

$$\begin{array}{c}\hfill {\tilde{\mathbf{x}}}_k=f({\widehat{\mathbf{x}}}_{k-1},{\mathbf{u}}_k,\mathbf{0}),\tilde{\mathbf{z}}=h({\tilde{\mathbf{x}}}_k,\mathbf{0}).\end{array}$$

(2.51)

We note that the nonlinear control terms in the discrete system does not impact the Kalman filter. So, we do not need to linearize the nonlinear functions with respect to control variables.

Then, the linearized system has the following form

$$\begin{array}{cc}\hfill {\mathbf{x}}_k& ={\tilde{\mathbf{x}}}_k+{\mathbf{A}}_k({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_{k-1})+{\mathbf{W}}_k{\mathbf{w}}_k,\hfill \\ \hfill {\mathbf{z}}_k& ={\tilde{\mathbf{z}}}_k+{\mathbf{H}}_k({\mathbf{x}}_k-{\tilde{\mathbf{x}}}_k)+{\mathbf{V}}_k{\mathbf{v}}_k,\hfill \end{array}$$

(2.52)

where the Jacobian of f w.r.t. $\mathbf{x}$ is ${\mathbf{A}}_k={\left.{\displaystyle \frac{\partial f}{\partial \mathbf{x}}}\right|}_{({\widehat{\mathbf{x}}}_{k-1},{\mathbf{u}}_k,0)}\in {\mathbb{R}}^{n\times n}$. the Jacobian of f w.r.t. $\mathbf{w}$ is ${\mathbf{W}}_k={\left.{\displaystyle \frac{\partial f}{\partial \mathbf{w}}}\right|}_{({\widehat{\mathbf{x}}}_{k-1},{\mathbf{u}}_k,0)}\in {\mathbb{R}}^{n\times q}$, the Jacobian of h w.r.t. $\mathbf{x}$ is ${\mathbf{H}}_k={\left.{\displaystyle \frac{\partial h}{\partial \mathbf{x}}}\right|}_{({\tilde{\mathbf{x}}}_k,0)}\in {\mathbb{R}}^{p\times n}$, and the Jacobian of h w.r.t. $\mathbf{v}$ is ${\mathbf{V}}_k={\left.{\displaystyle \frac{\partial h}{\partial \mathbf{v}}}\right|}_{({\tilde{\mathbf{x}}}_k,0)}\in {\mathbb{R}}^{p\times r}$.

The noise terms in the linearized system have distributions ${\mathbf{W}}_k{\mathbf{w}}_k\sim \mathcal{N}(\mathbf{0},{\mathbf{W}}_k{\mathbf{Q}}_k{\mathbf{W}}_k^T)$, ${\mathbf{V}}_k{\mathbf{v}}_k\sim \mathcal{N}(\mathbf{0},{\mathbf{V}}_k{\mathbf{R}}_k{\mathbf{V}}_k^T)$.

The EKF algorithm that proceeds iteratively in two steps is formulated as follows.

Algorithm 2.2.

1.

Prediction Step

Compute the a priori state estimate and error covariance:

$$\begin{array}{cc}\hfill \mathrm{A}\mathrm{priori}\mathrm{state}\mathrm{estimate}\text{:}& {\widehat{\mathbf{x}}}_k^{-}=f({\widehat{\mathbf{x}}}_{k-1},{\mathbf{u}}_k,\mathbf{0}),\hfill \\ \hfill \mathrm{A}\mathrm{priori}\mathrm{error}\mathrm{covariance}\text{:}& {\mathbf{P}}_k^{-}={\mathbf{A}}_k{\mathbf{P}}_{k-1}{\mathbf{A}}_k^T+{\mathbf{W}}_k{\mathbf{Q}}_k{\mathbf{W}}_k^T.\hfill \end{array}$$

(2.53)

where ${\widehat{\mathbf{x}}}_k^{-}\in {\mathbb{R}}^{n\times 1}$ is the predicted state before incorporating new measurement, and ${\mathbf{P}}_k^{-}\in {\mathbb{R}}^{n\times n}$ the predicted error covariance matrix.

2.

Correction Step

Update the estimate using the latest measurement ${\mathbf{z}}_k$:

$$\begin{array}{cc}\hfill \mathrm{Kalman}\mathrm{gain}\text{:}& {\mathbf{K}}_k={\mathbf{P}}_k^{-}{\mathbf{H}}_k^T{({\mathbf{H}}_k{\mathbf{P}}_k^{-}{\mathbf{H}}_k^T+{\mathbf{V}}_k{\mathbf{R}}_k{\mathbf{V}}_k^T)}^{-1},\hfill \\ \hfill \mathrm{A}\mathrm{posteriori}\mathrm{state}\mathrm{estimate}\text{:}& {\widehat{\mathbf{x}}}_k={\widehat{\mathbf{x}}}_k^{-}+{\mathbf{K}}_k\left({\mathbf{z}}_k-h({\widehat{\mathbf{x}}}_k^{-},\mathbf{0})\right),\hfill \\ \hfill \mathrm{A}\mathrm{posteriori}\mathrm{error}\mathrm{covariance}\text{:}& {\mathbf{P}}_k=({\mathbf{I}}_n-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_k^{-}.\hfill \end{array}$$

(2.54)

where ${\mathbf{K}}_k\in {\mathbb{R}}^{n\times p}$ is the Kalman gain matrix, ${\mathbf{I}}_n$ is the identity matrix of size $n\times n$, and $h({\widehat{\mathbf{x}}}_k^{-},\mathbf{0})$ the predicted measurement based on the a priori state estimate.

For this, we cannot establish that this is an unbiased estimator for the state.

3. Successive Linear Estimators for Saturated Soil

We first review the method of successive linear estimators [19] and then formulate it in the context of the extended Kalman filter with nonlinear state and observation equations. We consider the steady state Richards equation in a spatially heterogenous media

$$\begin{array}{c}\hfill \nabla ·\left(K\right(\mathbf{x})\nabla \varphi (\mathbf{x}\left)\right)=S\left(\mathbf{x}\right),\end{array}$$

(3.1)

This equation describes steady-state saturated flow, i.e., groundwater flow under dynamic equilibrium conditions, such as aquifer flow under long-term stable recharge or discharge. The hydraulic conductivity, $K\left(x\right)\ge 0$, depends solely on properties of the medium and is independent of the water content. $\varphi \left(\mathbf{x}\right)$ is the water head and $S\left(\mathbf{x}\right)$ is the source. In anisotropic media, $K\left(\mathbf{x}\right)$ is a second order tensor. Our objective is to approximate $K\left(\mathbf{x}\right)$ using limited measurements of $\varphi \left(\mathbf{x}\right)$ and $K\left(\mathbf{x}\right)$ spatially.

To ensure $K\left(\mathbf{x}\right)$ is nonnegative, we use

$$\begin{array}{c}\hfill lnK\left(\mathbf{x}\right)=F\left(\mathbf{x}\right)+f\left(\mathbf{x}\right),\varphi =H+h(\mathbf{x},t),\end{array}$$

(3.2)

where $F\left(\mathbf{x}\right)=\mathbb{E}[lnK]$ is the expected value of $lnK$ when K is viewed as a random variable, and $H=\mathbb{E}\left[\varphi \right]$ is the expected value of the water head at $(\mathbf{x},t)$. Both F and $\varphi$ are regarded as random variables. The interpretation here is that K is a random variable so that $\varphi$ as a solution that depends on K is random passively. As the result, f and h are random variables with zero means $\mathbb{E}\left[f\right]=\mathbb{E}\left[h\right]=0$.

Given a finite number of measurements of $K\left(\mathbf{x}\right)$ and $\varphi \left(\mathbf{x}\right)$, can one determine $F\left(\mathbf{x}\right)$? Before we formulate the problem rigorously, let’s first review Kriging and co-Kriging method for obtaining linear estimators.

3.1. Kriging

Kriging is an elegant geostatistical method, used for the optimal unbiased estimation of unknown points based on their spatial correlations [39]. Assuming $f\left(\mathbf{x}\right)$ are measured at ${\mathbf{x}}_i,i=1,\dots ,n$, where n is the number of the observation points, one would like to estimate $f\left(\mathbf{z}\right)$ at another point $\mathbf{z}$ via a function as follows

$$\begin{array}{c}\hfill f\left(\mathbf{z}\right)=F(f\left({\mathbf{x}}_1\right),\dots ,f\left({\mathbf{x}}_n\right)),\end{array}$$

(3.3)

where F belongs to a class of functions such as linear, quadratic, or more complicated form. If F is linear, it is called a linear estimator and is also known as Kriging. In this case, we posit the approximation as follows

$$\begin{array}{c}\hfill f_e\left(\mathbf{z}\right)={\displaystyle \sum _{i=1}^n}{\lambda }_{i}f\left({\mathbf{x}}_i\right),\end{array}$$

(3.4)

where ${\lambda }_i$ are parameters to be estimated.

Unbiasedness requires that the expected estimated value be equal to the expected true value

$$\begin{array}{c}\hfill \begin{array}{c}\mathbb{E}\left[f_e\left(\mathbf{z}\right)\right]=\mathbb{E}\left[f\left(\mathbf{z}\right)\right].\hfill \end{array}\end{array}$$

(3.5)

It implies

$${\displaystyle \sum _{i=1}^n}{\lambda }_i=1,\mathrm{if}\mathbb{E}\left[f\left(\mathbf{z}\right)\right]\ne 0.$$

We hope to find parameter $\Lambda =\left({\lambda }_i\right)$ to minimize the variance of the estimation error

$$\begin{array}{c}\hfill \begin{array}{c}\underset{\Lambda }{min}{\sigma }_e^2\left(\mathbf{z}\right)=\underset{\Lambda }{min}\mathbb{E}\left[{\left(f\left(\mathbf{z}\right)-f_e\left(\mathbf{z}\right)\right)}^2\right]\hfill \\ =\underset{\Lambda }{min}[C(\mathbf{z},\mathbf{z})+{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^n}{\lambda }_i{\lambda }_{j}C({\mathbf{x}}_i,{\mathbf{x}}_j)-2{\displaystyle \sum _{i=1}^n}{\lambda }_{i}C(\mathbf{z},{\mathbf{x}}_i)],\hfill \end{array}\end{array}$$

(3.6)

where $C({\mathbf{x}}_i,{\mathbf{x}}_j)=\mathrm{Cov}(f\left({\mathbf{x}}_i\right),f\left({\mathbf{x}}_j\right))$ is the covariance of $f\left({\mathbf{x}}_i\right)$ and $f\left({\mathbf{x}}_j\right)$. Notice that $f_e\left(z\right)$ is a weight sum of the observed data, $f_e\left(\mathbf{z}\right)$ is the measured or estimated mean of $f\left(\mathbf{z}\right)$.

To estimate $\Lambda$, we introduce Lagrange multipliers $\eta$ and construct the Lagrange function as follows

$$\begin{array}{c}\hfill \begin{array}{c}\mathcal{L}(\Lambda ,\eta )=C(\mathbf{z},\mathbf{z})+{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^n}{\lambda }_i{\lambda }_{j}C({\mathbf{x}}_i,{\mathbf{x}}_j)-2{\displaystyle \sum _{i=1}^n}{\lambda }_{i}C(\mathbf{z},{\mathbf{x}}_i)-2\eta \left(1-{\displaystyle \sum _{i=1}^n}{\lambda }_i\right).\hfill \end{array}\end{array}$$

(3.7)

By utilizing the standard procedures for extremum seeking, the extremum conditions are derived as follows

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \sum _{j=1}^n}{\lambda }_{j}C({\mathbf{x}}_k,{\mathbf{x}}_j)=C(\mathbf{z},{\mathbf{x}}_k)-\eta ,\forall k=1,\dots ,n,1-{\displaystyle \sum _{i=1}^n}{\lambda }_i=0.\hfill \end{array}\end{array}$$

(3.8)

Substituting (62) into (61), we obtain the value of the Lagrange function at the extremum:

$$\begin{array}{c}\hfill \begin{array}{c}\mathcal{L}(\Lambda ,\eta )=C(\mathbf{z},\mathbf{z})-{\mathbf{c}}^T\Lambda -\eta ,\hfill \end{array}\end{array}$$

(3.9)

where $\mathbf{c}$ is an $n\times 1$ vector with elements $c_i=C(\mathbf{z},{\mathbf{x}}_i).$ We denote $\mathbf{C}$ as an $n\times n$ covariance matrix with elements $C_{ij}=C({\mathbf{x}}_i,{\mathbf{x}}_j),$ $\Lambda ={[{\lambda }_1,\dots ,{\lambda }_n]}^T$ and $\mathbf{1}={[1,\dots ,1]}^T.$

The complete system is given by

$$\begin{array}{c}\hfill \begin{array}{c}\mathbf{A}{\Lambda }_1={\mathbf{c}}_1,\hfill \end{array}\end{array}$$

(3.10)

where

$$\mathbf{A}=\left[\begin{array}{cc}\mathbf{C}& \mathbf{1}\\ {\mathbf{1}}^T& 0\end{array}\right],{\Lambda }_1=\left[\begin{array}{c}\Lambda \\ \eta \end{array}\right],{\mathbf{c}}_1=\left[\begin{array}{c}\mathbf{c}\\ 1\end{array}\right].$$

The least square solution is

$$\begin{array}{c}\hfill \begin{array}{c}{\Lambda }_1={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T{\mathbf{c}}_1.\hfill \end{array}\end{array}$$

(3.11)

We note that the covariance $\mathbf{c}$ needs to be known or estimated in order to get the solution. I.e., in order to have a linear approximation using the observed data, the covariance between the estimate function value and the observed values must be known a priori.

3.2. Co-Kriging

Co-Kriging constructs a linear unbiased optimal estimator by jointly utilizing the spatial correlation between the target variable (primary variable) and auxiliary variables [15]. For a more accurate estimation, we use both measured values, $f\left({\mathbf{x}}_i\right),i=1,\dots ,n_f$ (primary variable) and $h\left(x_j\right),j=1,\dots ,n_h$ (auxiliary variable) to estimate $f\left(\mathbf{z}\right)$ linearly as follows

$$\begin{array}{c}\hfill f_e\left(\mathbf{z}\right)={\displaystyle \sum _{i=1}^{n_f}}{\lambda }_{i}f\left({\mathbf{x}}_i\right)+{\displaystyle \sum _{j=1}^{n_h}}{\beta }_{j}h\left({\mathbf{x}}_j\right).\end{array}$$

(3.12)

The least squre solution is given by

$$\begin{array}{c}\hfill (\Lambda ,\beta )=arg\underset{\Lambda ,\beta }{min}\mathbb{E}\left[{\left(f\left(\mathbf{z}\right)-f_e\left(\mathbf{z}\right)\right)}^2\right].\end{array}$$

(3.13)

The solution is governed by the following linear equation system

$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\mathbf{C}}_{ff}& {\mathbf{C}}_{fh}\\ {\mathbf{C}}_{hf}& {\mathbf{C}}_{hh}\end{array}\right]\left[\begin{array}{c}\lambda \\ \beta \end{array}\right]=\left[\begin{array}{c}{\mathbf{c}}_f\\ {\mathbf{c}}_h\end{array}\right].\end{array}$$

(3.14)

Here, the covariance matrices are defined by

$$\begin{array}{c}\hfill \begin{array}{c}{C}_{fh}({\mathbf{x}}_i,{\mathbf{x}}_j)=\mathrm{Cov}[f\left({\mathbf{x}}_i\right),h\left({\mathbf{x}}_j\right)].\hfill \end{array}\end{array}$$

(3.15)

The covariance matrices with subscripts $ff,hf,hh$ are defined analogously. The covariance vectors on the right hand side are defined by

$$\begin{array}{c}\hfill \begin{array}{c}{\left({\mathbf{c}}_f\right)}_i=\mathrm{Cov}[f\left(\mathbf{z}\right),f\left({\mathbf{x}}_i\right)],i=1,\dots ,n_f,\hfill \\ {\left({\mathbf{c}}_h\right)}_j=\mathrm{Cov}[f\left(\mathbf{z}\right),h\left({\mathbf{x}}_j\right)],j=1,\dots ,n_h.\hfill \end{array}\end{array}$$

(3.16)

To make it an unbiased estimator, constraints can be employed like in the case of Kriging. Two constraints in Co-Kriging are given by

$$\begin{array}{c}\hfill {\displaystyle \sum _{i=1}^{n_f}}{\lambda }_i=1,{\displaystyle \sum _{j=1}^{n_h}}{\beta }_j=0.\end{array}$$

(3.17)

Once again, the covariances must be known or estimated to get the solution.

The normal equations are

$$\begin{array}{c}\hfill \left[\begin{array}{cccc}{\mathbf{C}}_{ff}& \mathbf{1}& {\mathbf{C}}_{fh}& \mathbf{0}\\ {\mathbf{1}}^T& 0& {\mathbf{0}}^T& 0\\ {\mathbf{C}}_{hf}& \mathbf{0}& {\mathbf{C}}_{hh}& \mathbf{1}\\ {\mathbf{0}}^T& 0& {\mathbf{1}}^T& 0\end{array}\right]\left[\begin{array}{c}\lambda \\ {\eta }_f\\ \beta \\ {\eta }_h\end{array}\right]=\left[\begin{array}{c}{\mathbf{c}}_f\\ 1\\ {\mathbf{c}}_h\\ 0\end{array}\right],\end{array}$$

(3.18)

where ${\eta }_f$ and ${\eta }_h$ are two Lagrange multipliers used for enforcing that the linear estimator is unbiased.

3.3. Successive Linear Estimators

In the context of SLE for hydraulic conductivity, the state equation is given by

$$lnK={\widehat{Y}}_c+f,f\sim \mathcal{N}(0,R_{ff}),R_{ff}=\mathbb{E}\left[f{f}^T\right],$$

(3.19)

where ${\widehat{Y}}_c=\mathbb{E}[lnK]$ and $\mathbb{E}\left[f\right]=0$. The observation equation is given implicitly by

$$\begin{array}{c}\hfill \begin{array}{c}\varphi =\mathcal{F}(lnK),\iff \nabla ·(e^{lnK}\nabla \varphi \left(x\right))-S=0,\hfill \end{array}\end{array}$$

(3.20)

where $\mathcal{F}$ is the mapping that maps $lnK$ to the solution of the Richards equation, $K\left(x\right)$ is the saturated Transmissivity, $\varphi \left(x\right)$ is the total water head, and

$$h=\varphi -\mathbb{E}\left[\varphi \right]\sim \mathcal{N}(0,R_{hh}),R_{hh}=\mathbb{E}\left[h{h}^T\right].$$

The successive linear estimator at a location, $x_o$, for $lnK$ is given iteratively by

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{Y}}_c^{(r+1)}\left(x_o\right)={\widehat{Y}}_c^{\left(r\right)}\left(x_o\right)+{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\omega }_j^{\left(r\right)}\left[{\varphi }^{*}\left(x_j\right)-{\varphi }^{\left(r\right)}\left(x_j\right)\right],\hfill \end{array}\end{array}$$

(3.21)

where ${\widehat{Y}}_c^{\left(r\right)}\left(x_o\right)$ is the rth iteration of the estimated conditional mean of $lnK$, r is the iteration index, ${\omega }_j$ is the model weight parameter, ${\varphi }^{*}$ denotes the observation, and ${\varphi }^{\left(r\right)}$ is the outcome of the observation equation or the solution of the Richards equation (model-computed value) at the r-th iteration.

To ensure the estimator having minimal variance, the goal is to minimize the variance with respect to $\omega$’s

$$\begin{array}{c}\hfill \begin{array}{c}min{\mathbb{E}}_{{\omega }^{\left(r\right)}}\left[{\left(lnK-{\widehat{\mathbf{Y}}}_c^{(r+1)}\right)}^2\right].\hfill \end{array}\end{array}$$

(3.22)

We differentiate the squared error with respect to ${\omega }^{\left(r\right)}$ and set the resultant to zero to yield

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\omega }_j^{\left(r\right)}{\epsilon }_{hh}^{\left(r\right)}(x_j,x_l)={\epsilon }_{yh}^{\left(r\right)}(x_o,x_l),\hfill \end{array}\end{array}$$

(3.23)

where $h={\varphi }^{*}-\varphi ,y=lnK-{\widehat{Y}}_c,{\epsilon }_{hh}=\mathbb{E}\left[hh\right],{\epsilon }_{yh}=\mathbb{E}\left[yh\right].$ In the following, we set ${\widehat{Y}}_c=0$.

The iteration process of SLE is summarized in the following algorithm.

Algorithm 3.1

(SLE).

1. Initialization

The prior estimate is obtained using the classical cokriging technique:

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{Y}}_c^{\left(0\right)}\left({\mathbf{x}}_o\right)={\displaystyle \sum _{i=1}^{n_f}}{\lambda }_{io}{f}_i^{*}\left({\mathbf{x}}_i\right)+{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\mu }_{jo}{h}_j^{*}\left({\mathbf{x}}_j\right).\hfill \end{array}\end{array}$$

(3.24)

where λ and μ are calculated by the following equations:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \sum _{i=1}^{n_f}}{\lambda }_{io}{R}_{ff}({\mathbf{x}}_i,{\mathbf{x}}_l)+{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\mu }_{jo}{R}_{hf}({\mathbf{x}}_j,{\mathbf{x}}_l)=R_{ff}({\mathbf{x}}_o,{\mathbf{x}}_l),l=1,2,\dots ,n_f,\hfill \\ {\displaystyle \sum _{i=1}^{n_f}}{\lambda }_{io}{R}_{fh}({\mathbf{x}}_i,{\mathbf{x}}_l)+{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\mu }_{jo}{R}_{hh}({\mathbf{x}}_j,{\mathbf{x}}_l)=R_{fh}({\mathbf{x}}_o,{\mathbf{x}}_l),l=n_f+1,n_f+2,\dots ,n_f+n_h.\hfill \end{array}\end{array}$$

(3.25)

The prior covariance is:

$$\begin{array}{c}\hfill \begin{array}{c}{\epsilon }_{yy}^{\left(1\right)}({\mathbf{x}}_o,{\mathbf{x}}_k)=R_{ff}({\mathbf{x}}_o,{\mathbf{x}}_k)-{\displaystyle \sum _{i=1}^{n_f}}{\lambda }_{io}{R}_{ff}({\mathbf{x}}_o,{\mathbf{x}}_i)-{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\mu }_{jo}{R}_{fh}({\mathbf{x}}_k,{\mathbf{x}}_j).\hfill \end{array}\end{array}$$

(3.26)

2. Prediction ((r+1)-th step)

At the current (r+1)-th step, the optimal solution ${\widehat{Y}}_c^{\left(r\right)}$ from the previous step (i.e., the r-th step) is treated as the state prior for this step. Meanwhile, the error covariance ${\epsilon }_{yy}^{\left(r\right)}$ from the previous step is taken as the error covariance prior for this step.

3. Correction (r+1)-th step

${\omega }^{\left(r\right)}$, similar to the Kalman gain, is given by

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\omega }_j^{\left(r\right)}{\epsilon }_{hh}^{\left(r\right)}(x_j,x_l)={\epsilon }_{yh}^{\left(r\right)}(x_o,x_l),l=n_f+1,\dots ,n_f+n_h.\hfill \end{array}\end{array}$$

(3.27)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\epsilon }_{hh}^{\left(r\right)}=J^{\left(r\right)}{\epsilon }_{yy}^{\left(r\right)}{J}^{{\left(r\right)}^T},\hfill \\ {\epsilon }_{yh}^{\left(r\right)}=J^{\left(r\right)}{\epsilon }_{yy}^{\left(r\right)}.\hfill \end{array}\end{array}$$

(3.28)

and the Jacobian matrix (Sensitivity matrix) J of the observation equation is given by:

$$\begin{array}{c}\hfill \begin{array}{c}{J}^{\left(r\right)}={\left.{\displaystyle \frac{\partial \mathcal{F}}{\partial lnK}}\right|}_{{\widehat{Y}}_c^{\left(r\right)}}.\hfill \end{array}\end{array}$$

(3.29)

In practical calculations, an additional relaxation term, Θ, is usually added to the diagonal of the matrix in (81) during each iteration:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\omega }_j^{\left(r\right)}{\epsilon }_{hh}^{\left(r\right)}(x_j,x_l)+{\Theta }^{\left(r\right)}{\delta }_{ll}={\epsilon }_{yh}^{\left(r\right)}(x_o,x_l),l=n_f+1,\dots ,n_f+n_h.\hfill \end{array}\end{array}$$

(3.30)

This is actually an $L_2$ regularization (also known as the Tikhonov regularization or the Ridge Regression). From a Bayesian perspective, it is equivalent to the parameter ${\omega }_j^{\left(r\right)}$ following a Gaussian distribution with a mean of 0 and a variance proportional to $1/\Theta$. Functionally, Θ can be interpreted as a “virtual observation noise," which makes the system perceive the observed data as noisier than it actually is, thereby suppressing overfitting and numerical oscillations.

The posterior estimate is then given by

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{Y}}_c^{(r+1)}\left(x_o\right)={\widehat{Y}}_c^{\left(r\right)}\left(x_o\right)+{\displaystyle \sum _{j=n_f+1}^{n_f+n_h}}{\omega }_j^{\left(r\right)}\left[{\varphi }_j^{*}\left(x_j\right)-{\varphi }_j^{\left(r\right)}\left(x_j\right)\right],\hfill \end{array}\end{array}$$

(3.31)

where ${\varphi }_j^{*}\left(x_j\right)$ are given. The error covariance is updated as follows

$$\begin{array}{c}\hfill \begin{array}{c}{\epsilon }_{yy}^{(r+1)}(x_o,x_k)={\epsilon }_{yy}^{\left(r\right)}(x_o,x_k)-{\displaystyle \sum _{l=n_f+1}^{n_f+n_h}}{\omega }_l^{\left(r\right)}{\epsilon }_{yh}^{\left(r\right)}(x_k,x_l).\hfill \end{array}\end{array}$$

(3.32)

3.4. Reformulation of the SLE as an EKF

In this subsection, we show that the Successive Linear Estimator (SLE) is essentially a specific application of the Extended Kalman Filter (EKF) in static parameter estimations using the Richards equation.

In the problem of estimating the hydraulic conductivity coefficient, the logarithmic conductivity coefficient $lnK\left(x\right)$ is defined as the state variable and the hydraulic head $\varphi \left(x\right)$ is defined as the observation variable. A nonlinear model system is established based on the solution solver of the Richards equation with the given hydraulic conductivity, $lnK\left(x\right)$.

This is a steady state parameter estimation problem, so the state transition function is an identity mapping augmented with a Gaussian noise:

$${\mathbf{x}}_k={\mathbf{x}}_{k-1}+{\mathbf{w}}_k,(\mathbf{x}\leftarrow {\widehat{Y}}_c).$$

(3.33)

Within the SLE framework, the process noise ${\mathbf{w}}_k\equiv 0$, i.e., $\mathbf{Q}=\mathbb{E}\left[{\mathbf{w}}_k{\mathbf{w}}_k^T\right]=0$. The condition of $\mathbf{Q}=0$ indicates a state of complete confidence that the state (i.e., $lnK$) will not undergo unpredictable changes during the iteration process. Its uncertainty propagates solely through the initial covariance matrix $P_0^{-}$ and the Jacobian matrix $\mathbf{J}$ presented below. Here, ${\mathbf{x}}_k={\widehat{Y}}_c^{\left(k\right)}\in {\mathbb{R}}^{N\times 1}$ is the conditional mean of $lnK\left(x\right)$ at the k-th iteration (the state to be estimated).

The observation equation with a Gaussian noise is given by

$${\mathbf{z}}_k=\mathcal{F}\left({\mathbf{x}}_k\right)+{\mathbf{v}}_k,(\mathbf{z}\leftarrow \varphi ),$$

(3.34)

where ${\mathbf{z}}_k\in {\mathbb{R}}^{n_h\times 1}$ and $\mathcal{F}$ is the observation operator. I.e., ${\mathbf{z}}_k$ is the solution of the steady state Richards equation $\nabla ·\left(e^{{\mathbf{x}}_k}\nabla {\mathbf{z}}_k\right)=S\left(x\right)$. without loss of generality, we assume $S\left(x\right)=0$ in the following. The measurements are assumed error-free, with the observation noise ${\mathbf{v}}_k\equiv 0$. As a result, the covariance of the observation noise ${\mathbf{R}}_k=\mathbb{E}\left[{\mathbf{v}}_k{\mathbf{v}}_k^{\mathrm{T}}\right]=0$.

One defines the state error as ${\mathbf{e}}_k={\mathbf{x}}_k-{\widehat{\mathbf{x}}}_k=lnK-{\widehat{Y}}_c$. The objective of the error estimation is to minimize the covariance of the state error with respect to parameter ${\mathbf{K}}_k$, i.e.,

$$\begin{array}{c}\hfill \underset{{\mathbf{K}}_k}{min}tr\left(\mathbb{E}\left[{\mathbf{e}}_k{\mathbf{e}}_k^T\right]\right).\end{array}$$

(3.35)

One denote ${\widehat{Y}}_c$ as the unbiased estimator of $lnK$, that is, $\mathbb{E}\left[{\widehat{Y}}_c\right]=lnK$; it follows that $\mathbb{E}\left[{\mathbf{e}}_k\right]=0$.

Initialization

The prior state ${\widehat{\mathbf{x}}}_0^{-}={\widehat{Y}}_c^{\left(0\right)}\in {\mathbb{R}}^{N\times 1}$ is determined by the cokriging method using both the conductivity and hydraulic head observations through the following equation system

$${\widehat{\mathbf{x}}}_0^{-}=\Lambda {\mathbf{f}}^{\prime }+\mathbf{M}{\mathbf{h}}^{\prime }.$$

(3.36)

Here, ${\mathbf{f}}^{\prime }\in {\mathbb{R}}^{n_f\times 1}$ is the de-meaned logarithmic conductivity coefficient vector, with ${\mathbf{f}}_i^{\prime }=lnK\left(x_i\right)-\mathbb{E}[lnK]$ ($i=1,\dots ,n_f$); ${\mathbf{h}}^{\prime }\in {\mathbb{R}}^{n_h\times 1}$ is the de-meaned hydraulic head vector, with ${\mathbf{h}}_j^{\prime }=\varphi \left(\mathbf{x}\right)-\mathbb{E}\left[\varphi \right]$ ($j=1,\dots ,n_h$). $\Lambda \in {\mathbb{R}}^{N\times n_f}$ is the global kriging weight matrix for transmissivity observations, with ${[\Lambda ]}_{o,i}={\lambda }_{io}$ ($o=1,\dots ,N$; $i=1,\dots ,n_f$); $\mathbf{M}\in {\mathbb{R}}^{N\times n_h}$ is the global kriging weight matrix for hydraulic head observations, with ${\left[\mathbf{M}\right]}_{o,j}={\mu }_{jo}$ ($o=1,\dots ,N$; $j=1,\dots ,n_h$). N is the total number of spatial nodes to be estimated, $n_f$ is the number of transmissivity observation points, and $n_h$ is the number of hydraulic head observation points.

$\Lambda$ and M are determined by the following system of linear equations:

$$\left[\begin{array}{cc}{\mathbf{R}}_{ff,obs}& {\mathbf{R}}_{hf,obs}\\ {\mathbf{R}}_{fh,obs}& {\mathbf{R}}_{hh,obs}\end{array}\right]\left[\begin{array}{c}{\lambda }_o\\ {\mu }_o\end{array}\right]=\left[\begin{array}{c}{\mathbf{r}}_{ff,o}\\ {\mathbf{r}}_{fh,o}\end{array}\right],$$

(3.37)

where ${\lambda }_o\in {\mathbb{R}}^{n_f\times 1}$ is the weight sub-vector corresponding to the transmissivity observations for the target node $x_o$ (the o-th row of $\Lambda$), and ${\mu }_o\in {\mathbb{R}}^{n_h\times 1}$ is the weight sub-vector corresponding to the hydraulic head observations for the target node $x_o$ (the o-th row of $\mathbf{M}$). ${\mathbf{R}}_{ff,obs}\in {\mathbb{R}}^{n_f\times n_f}$ is the autocovariance matrix among the transmissivity observation points, with ${\left[{\mathbf{R}}_{ff,obs}\right]}_{i,l}={\mathbf{R}}_{ff}(x_i,x_l)$ ($i,l=1,\dots ,n_f$); ${\mathbf{R}}_{hh,obs}\in {\mathbb{R}}^{n_h\times n_h}$ is the autocovariance matrix among the hydraulic head observation points, with ${\left[{\mathbf{R}}_{hh,obs}\right]}_{j,l}={\mathbf{R}}_{hh}(x_j,x_l)$ ($j,l=1,\dots ,n_h$); ${\mathbf{R}}_{fh,obs}\in {\mathbb{R}}^{n_f\times n_h}$ is the cross-covariance matrix between the transmissivity and hydraulic head observation points, with ${\left[{\mathbf{R}}_{fh,obs}\right]}_{i,j}={\mathbf{R}}_{fh}(x_i,x_j)$ ($i=1,\dots ,n_f$; $j=1,\dots ,n_h$), and ${\mathbf{R}}_{hf,obs}={\mathbf{R}}_{fh,obs}^T$. ${\mathbf{r}}_{ff,o}\in {\mathbb{R}}^{n_f\times 1}$ is the autocovariance vector between the target node $x_o$ and the transmissivity observation points, with ${\left[{\mathbf{r}}_{ff,o}\right]}_l={\mathbf{R}}_{ff}(x_o,x_l)$ ($l=1,\dots ,n_f$); ${\mathbf{r}}_{fh,o}\in {\mathbb{R}}^{n_h\times 1}$ is the cross-covariance vector between the target node $x_o$ and the hydraulic head observation points, with ${\left[{\mathbf{r}}_{fh,o}\right]}_l={\mathbf{R}}_{fh}(x_o,x_l)$ ($l=1,\dots ,n_h$).

The value of the initial prior covariance ${\mathbf{P}}_0^{-}={\epsilon }_{yy}^{\left(0\right)}\in {\mathbb{R}}^{N\times N}$ is the conditional covariance after incorporating some observations ($n_f$ conductivity observations and $n_h$ hydraulic head observations) through one-step classical co-kriging:

$${\mathbf{P}}_0^{-}={\mathbf{R}}_{ff}-\Lambda {\mathbf{R}}_{ff,obs}{\Lambda }^T-\mathbf{M}{\mathbf{R}}_{fh}{\mathbf{M}}^T,$$

(3.38)

where ${\mathbf{R}}_{ff}\in {\mathbb{R}}^{N\times N}$ is the global autocovariance matrix of the logarithmic conductivity coefficient perturbation term f, with ${\left[{\mathbf{R}}_{ff}\right]}_{a,b}={\mathbf{R}}_{ff}(x_a,x_b)$ ($a,b=1,\dots ,N$); ${\mathbf{R}}_{fh}\in {\mathbb{R}}^{N\times n_h}$ is the global cross-covariance matrix between the logarithmic conductivity coefficient perturbation term f and the hydraulic head perturbation term h, with ${\left[{\mathbf{R}}_{fh}\right]}_{a,j}={\mathbf{R}}_{fh}(x_a,x_j)$ ($a=1,\dots ,N$; $j=1,\dots ,n_h$).

SLE Algorithm formulated in EKF framework

1.

Prediction Step Based on the previous state estimate ${\widehat{\mathbf{x}}}_{k-1}$ (posterior estimate), predict the prior state and prior error covariance of the current iteration.

$$\begin{array}{cc}\hfill \mathrm{Prior}\mathrm{state}\mathrm{estimate}\text{:}& {\widehat{\mathbf{x}}}_k^{-}={\widehat{\mathbf{x}}}_{k-1},\hfill \\ \hfill \mathrm{Prior}\mathrm{error}\mathrm{covariance}\text{:}& {\mathbf{P}}_k^{-}={\mathbf{P}}_{k-1}+{\mathbf{Q}}_k\in {\mathbb{R}}^{N\times N}.\hfill \end{array}$$

(3.39)

When ${\mathbf{Q}}_k=0$, ${\mathbf{P}}_k^{-}={\epsilon }_{yy}^{\left(k\right)}.$ Since the state equation has no dynamic evolution and is only a static estimate, the state estimate remains the result of the previous optimization in the absence of new observations, and the state uncertainty remains unchanged in the absence of new information. However, in practical groundwater systems, the hydraulic conductivity may exhibit unmodeled minor dynamic perturbations due to factors such as medium heterogeneity, temperature variations, and pressure fluctuations—these perturbations are referred to as process noise.

2.

Correction Step Combined with the hydraulic head observation data ${\mathbf{z}}_k={\varphi }^{*}\left(\mathbf{x}\right)$ of the current iteration, complete an iterative optimization by linearizing the observation equation, calculating the Kalman gain, and correcting the state estimate and error covariance.

Perform a first-order Taylor expansion of the nonlinear observation operator $\mathcal{F}\left(\mathbf{x}\right)$ at the current prior state ${\widehat{\mathbf{x}}}_k^{-}$ (i.e., the posterior ${\widehat{\mathbf{x}}}_{k-1}$ of the previous step). The linearized observation equation is

$${\mathbf{z}}_k={\tilde{\mathbf{z}}}_k+{\mathbf{J}}_k({\mathbf{x}}_k-{\widehat{\mathbf{x}}}_{k-1})+{\mathbf{v}}_k,$$

(3.40)

where ${\tilde{\mathbf{z}}}_k=\mathcal{F}\left({\widehat{\mathbf{x}}}_{k-1}\right)$ and ${\mathbf{J}}_k={\left.{\displaystyle \frac{\partial \mathcal{F}}{\partial \mathbf{x}}}\right|}_{{\widehat{\mathbf{x}}}_{k-1}}\in {\mathbb{R}}^{n_h\times N}$ is the Jacobian matrix (sensitivity matrix) of the observation equation.

$$\begin{array}{cc}\hfill \mathrm{Kalman}\mathrm{gain}\text{:}& {\mathbf{K}}_k={\mathbf{P}}_k^{-}{\mathbf{J}}_k^T{({\mathbf{J}}_k{\mathbf{P}}_k^{-}{\mathbf{J}}_k^T+{\mathbf{R}}_k)}^{-1},\hfill \\ \hfill \mathrm{Posterior}\mathrm{state}\mathrm{estimate}\text{:}& {\widehat{\mathbf{x}}}_k={\widehat{\mathbf{x}}}_k^{-}+{\mathbf{K}}_k({\mathbf{z}}_k-{\tilde{\mathbf{z}}}_k),\hfill \\ \hfill \mathrm{Posterior}\mathrm{error}\mathrm{covariance}\text{:}& {\mathbf{P}}_k=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{J}}_k){\mathbf{P}}_k^{-}.\hfill \end{array}$$

(3.41)

where ${\mathbf{K}}_k\in {\mathbb{R}}^{N\times n_h},{\left[{\mathbf{K}}_k\right]}_{o,j}={\omega }_{jo}^{\left(k\right)}\in {\mathbb{R}}^{N\times n_h}$, ${\widehat{\mathbf{x}}}_k={\widehat{Y}}_c^{\left(k\right)}$, ${\widehat{\mathbf{x}}}_k^{-}={\widehat{Y}}_c^{(k-1)}$, ${\mathbf{z}}_k={\varphi }^{*}\left(x\right)$, ${\tilde{\mathbf{z}}}_k={\varphi }^{\left(k\right)}\left(x\right)$, ${\mathbf{P}}_k^{-}={\epsilon }_{yy}^{(k-1)}$, and ${\mathbf{P}}_k={\epsilon }_{yy}^{\left(k\right)}$.

Here, ${\epsilon }_{yy}^{(k-1)}$ denotes the error covariance matrix of the log-conductivity estimate at the $(k-1)$-th iteration, defined as

$${\epsilon }_{yy}^{(k-1)}=\mathbb{E}\left[\left(lnK-{\widehat{Y}}_c^{(k-1)}\right){\left(lnK-{\widehat{Y}}_c^{(k-1)}\right)}^T\right],$$

where ${\widehat{Y}}_c^{(k-1)}$ is the conditional mean estimate from the previous iteration. This corresponds to the prior error covariance before the k-th correction step. However, in practical measurements, observation errors are inevitable due to sensor precision limitations and environmental interferences (e.g., water level fluctuations). The larger the value of ${\mathbf{R}}_k$, the less accurate the observation; correspondingly, the Kalman gain ${\mathbf{K}}_k$ becomes smaller, leading to a greater reliance on the a priori estimate.

From a mathematical perspective, both ${\mathbf{R}}_k$ and the regularization parameter, $\Theta$, that can be added as a $L_2$ regularization of the least squares optimization serve to stabilize matrix inversion. The presence of ${\mathbf{R}}_k$ ensures that the matrix ${\mathbf{J}}_k{\mathbf{P}}_k^{-}{\mathbf{J}}_k^T+{\mathbf{R}}_k$ is positive definite, thereby allowing safe inversion and avoiding numerical instability. In other words, when solving a linear system, both $\Theta$ and ${\mathbf{R}}_k$ play the role of regularization terms, which are used to prevent the matrix from being ill-conditioned. From a Bayesian viewpoint, both ${\mathbf{R}}_k$ and $\Theta$ represent prior information. The $\Theta$ term implicitly assumes that the Kalman gain ${\omega }_{jo}$ itself follows a Gaussian distribution centered at 0 with variance $1/\Theta$; this is a type of prior distribution on the parameters. The ${\mathbf{R}}_k$ term represents the assumption that the observation ${\mathbf{z}}_k$ follows a Gaussian distribution centered at the true value with variance ${\mathbf{R}}_k$; this is a type of prior distribution on the data. Although ${\mathbf{R}}_k$ and $\Theta$ have equivalent roles, their origins are not the same. ${\mathbf{R}}_k$ models the uncertainty in the data within the Bayesian framework, whereas $\Theta$ is an artificially introduced regularization term that models the uncertainty in the parameters.

By all means, we claim that the EKF outcome defined here is exactly the iterative solution of the SLE. This is the primary objective of this paper.

4. Parameter Estimation in Unsaturated Soil Using the Extended Kalman Filter

The general form of the Richards equation for flows in unsaturated soil is given by

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\partial \theta \left(\varphi \right(\mathbf{x},t\left)\right)}{\partial t}}+S(x,t)=\nabla ·\left[K(\theta ,\mathbf{x})\nabla \left(\varphi (\mathbf{x},t)+z\right)\right],\hfill \end{array}\end{array}$$

(4.1)

where $\theta (\mathbf{x},t)$ is the volumetric water content, $\varphi (\mathbf{x},t)$ is the hydraulic head, z is the elevation (the third coordinate in $\mathbf{x}$), $K(\theta ,\mathbf{x})$ is the unsaturated hydraulic conductivity, which is strongly coupled to the water content, $\theta$. $\theta \left(\varphi \right)$ is a monotonically increasing function of $\varphi$ and $K(\theta ,\mathbf{x})$ is a monotonically increasing function of $\theta$ as well as a monotonically increasing function of $\varphi$. In unsaturated soil, lower water content reduces pore-water connectivity, resulting in a smaller K; at saturation, K reaches its maximum value $K_s\left(\mathbf{x}\right)$ (the saturated hydraulic conductivity). $\varphi$ as the pressure head is normally negative, z is the elevation head, $\nabla (\varphi +z)$ is the gradient of the total head, which includes both the pressure gradient and the gravitational gradient. The pressure head $\varphi$ and water content $\theta$ are in one-to-one correspondence via the soil–water retention curve (SWRC), expressed as $\theta =\theta \left(\varphi \right).$

Given the Richards equation with conductivity, $K(\theta ,\mathbf{x})>0$, one would like to estimate $K(\theta ,\mathbf{x})$ through a few measurements of $K(\theta ({\mathbf{x}}_i,t_k),{\mathbf{x}}_i),i=1\dots ,n_f$ and $\varphi ({\mathbf{x}}_j,t_k),j=n_f+1,\dots ,n_f+n_h,k=1,\dots ,N_t$ subject to given initial and boundary conditions. Since $K>0$, one normally uses

$$\begin{array}{c}\hfill lnK=F+f\end{array}$$

(4.2)

as the effective variable to ensure its positivity, where $F=\mathbb{E}[lnK]$ is the expectation of $lnK$ and f is the small fluctuation with $\mathbb{E}\left[f\right]=0$ and $\mathbb{E}\left[\right(f,f\left)\right]=Q\left(\mathbf{x}\right)>0$. Analogously, we set

$$\begin{array}{c}\hfill \varphi =H+h,\end{array}$$

(4.3)

where $H=\mathbb{E}\left[\varphi \right]$ is the mean of $\varphi$ and his the perturbation with $\mathbb{E}\left[h\right]=0$ and $\mathbb{E}[h,h]=R(\mathbf{x},t)>0$.

This inverse problem could also be formulated using EKF. The objective here is to estimate $K(\theta ,\mathbf{x})$ using the spatial-temporal observation data The state equation for K and the observation equation for $\varphi$ is defined, respectively, by

$$\begin{array}{c}\hfill \begin{array}{c}K(\theta ,\mathbf{x})=K(\theta ,\mathbf{x}),\mathcal{F}\left(K\right)=\varphi ,\hfill \end{array}\end{array}$$

(4.4)

where $\mathcal{F}$ represents the time-dependent, inverse operator defined by the solution solver of the Richards equation. However, exploit the time-dependence of $\mathcal{F}$ separates this inverse problem from the saturate case alluded to in the previous section.

4.1. Parametric Inversion Aided by the Van Genuchten–Mualem Model

To make the problem tractable, we assume the hydraulic conductivity is admits a separable form that separates the spatial variability of the saturated hydraulic conductivity $K_s\left(\mathbf{x}\right)$ from the universal dependence of the relative permeability on pressure head $K_r\left(\theta \right)$ as follows,

$$\begin{array}{c}\hfill \begin{array}{c}K(\theta ,\mathbf{x})=K_s\left(\mathbf{x}\right)K_r\left(\theta \right),\hfill \end{array}\end{array}$$

(4.5)

where $K_s\left(\mathbf{x}\right)$ is the maximum hydraulic conductivity and $K_r\left(\theta \right)$ is given by one of the empirical constitutive laws such as the van Genuchten-Mualem model. For de-meaned $ln{K}_s\left(\mathbf{x}\right)$, we assume it can be expanded in a set of orthogonal basis functions ${\left\{{\phi }_i\left(x\right)\right\}}_{i=1}^M$. The choice of orthogonal system is not restricted a priori: on regular domains (e.g., rectangular or layered), classical bases such as Fourier series, cosine or sine bases (suited to specific boundary conditions), and orthogonal polynomials (e.g., Legendre, Chebyshev) are commonly employed. For problems where multi-scale or localized spatial features are anticipated, wavelet bases (e.g., Haar, Daubechies) provide alternatives with compact support in both space and scale. Of course, we could also consider the Karhunen–Loève expansion.

After the orthogonal basis is determined, we expand $ln{K}_s\left(\mathbf{x}\right)$ as follows

$$\begin{array}{c}\hfill \begin{array}{c}ln{K}_s\left(\mathbf{x}\right)=F\left(\mathbf{x}\right)+{\displaystyle \sum _{i=1}^M}{c}_i{\phi }_i\left(\mathbf{x}\right),\mathbb{E}\left[c_i\right]=0,\hfill \end{array}\end{array}$$

(4.6)

where $F\left(\mathbf{x}\right)=\mathbb{E}[ln{K}_s\left(\mathbf{x}\right)]$ is the prior mean field, prescribed as a constant or a spatial trend based on available information (e.g., geological interpretation or direct measurements of $ln{K}_s$). $\mathbf{c}={[c_1,c_2,\dots ,c_M]}^T\in {\mathbb{R}}^M$ is the vector of unknown coefficients. The prior covariance of $\mathbf{c}$ is denoted by ${\Sigma }_c\in {\mathbb{R}}^{M\times M}$.

The expansion, (101), serves two purposes. First, it reduces the dimensionality of the estimation problem from the number of grid nodes N (typically ${10}^3$–${10}^5$) to a much smaller number of coefficients M (typically ${10}^1$–${10}^2$). This compression is possible because the spatial correlation of $ln{K}_s$ implied by the prior covariance ${\Sigma }_c$ allows the random field to be accurately represented by the leading M basis functions. Second, the orthogonality of $\left\{{\phi }_i\right\}$ simplifies the subsequent Co-Kriging and EKF computations by providing a diagonal or nearly diagonal prior covariance structure for the coefficients. The truncation level M is chosen such that the retained basis functions capture a prescribed fraction (e.g., 90%–95%) of the total prior variance, i.e.,

$${\displaystyle \frac{{\sum }_{i=1}^M{\lambda }_i}{{\sum }_{i=1}^{\infty }{\lambda }_i}}\approx 0.9\text{–}0.95,$$

where ${\lambda }_i$ are the eigenvalues of the prior covariance operator associated with $\left\{{\phi }_i\right\}$. In practice, M is determined by the decay rate of the eigenvalue spectrum and is typically much smaller than N, making the inversion computationally feasible.

The prior mean field $F\left(\mathbf{x}\right)$ is kept fixed throughout the inversion; it does not appear in the state vector and is not updated by the EKF. At each iteration, the total $ln{K}_s$ field is reconstructed by adding $F\left(\mathbf{x}\right)$ to the fluctuation estimated from the current state ${\mathbf{x}}_k$, and the Richards equation is solved using this total field to compute the predicted hydraulic head. The objective of the inverse problem is essentially to estimate the unknown expansion coefficients $\left\{c_i\right\}$. The relative hydraulic conductivity $K_r\left(\theta \right)$ and the soil-water retention curve $\theta \left(\varphi \right)$ are prescribed by the van Genuchten-Mualem model, detailed below. The complete unsaturated hydraulic conductivity is then given by

$$\begin{array}{c}\hfill \begin{array}{c}K(\theta ,\mathbf{x})=K_s\left(\mathbf{x}\right)K_r\left(\theta \right)=exp\left(F\left(\mathbf{x}\right)+{\displaystyle \sum _{i=1}^M}{c}_i{\phi }_i\left(\mathbf{x}\right)\right)·K_r\left(\theta \right).\hfill \end{array}\end{array}$$

(4.7)

Because $K_r\left(\theta \right)$ is a monotonically increasing function of $\theta$ and hence of $\varphi$, and the prefactor $exp(F\left(\mathbf{x}\right)+\sum c_i{\phi }_i\left(\mathbf{x}\right))$ is strictly positive, the product $K(\theta ,\mathbf{x})$ inherits the monotonicity with respect to $\varphi$ required by physics.

The relative hydraulic conductivity $K_r\left(\theta \right)$ appearing in (102) is defined by the van Genuchten-Mualem model [9,10]. In this model, the effective water saturation $S_e$ is expressed as a function of the pressure head $\varphi$

$$\begin{array}{c}\hfill \begin{array}{c}{S}_e\left(\varphi \right)={\displaystyle \frac{\theta \left(\varphi \right)-{\theta }_r}{{\theta }_s-{\theta }_r}}=\left\{\begin{array}{cc}{\displaystyle \frac{1}{{[1+(\alpha \left|\varphi \right|)}^n{]}^m}},\hfill & \varphi <0,\hfill \\ 1,\hfill & \varphi \ge 0,\hfill \end{array}\right.\hfill \end{array}\end{array}$$

(4.8)

where ${\theta }_s$ and ${\theta }_r$ are the saturated and residual water contents, $\alpha$ and n are empirical shape parameters, and $m=1-1/n$. The relative hydraulic conductivity is then given by the Mualem model [9],

$$\begin{array}{c}\hfill \begin{array}{c}{K}_r\left(S_e\right)=S_e^{\ell }{\left[1-{\left(1-S_e^{1/m}\right)}^m\right]}^2,\hfill \end{array}\end{array}$$

(4.9)

where ℓ is a pore-connectivity parameter. The parameters $\{\alpha ,n,{\theta }_s,{\theta }_r,\ell \}$ are treated as known constants throughout the EKF inversion; they can be determined a priori from laboratory experiments or pedotransfer functions. If one wants, these parameter set could also be estimated using EKF subject to some predetermined constraints on their ranges.

With this parameterization, the state vector of the EKF is defined as

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{x}}_k={[c_1,c_2,\dots ,c_M]}^T\in {\mathbb{R}}^M,\hfill \end{array}\end{array}$$

(4.10)

and the observation operator $\mathcal{F}\left({\mathbf{x}}_k\right)$ maps the coefficients to the hydraulic head by solving the Richards equation (96) with the conductivity field constructed from (102). The inverse problem thus reduces to estimating the finite-dimensional coefficient vector ${\mathbf{x}}_k$ from head observations.

Specifically, let the state vector ${\mathbf{x}}_k\in {\mathbb{R}}^M$ represent the parameters to be estimated. The state evolution equation is given by

$${\mathbf{x}}_k={\mathbf{x}}_{k-1}+{\mathbf{w}}_k,{\mathbf{w}}_k\sim \mathcal{N}(\mathbf{0},{\mathbf{Q}}_k),$$

where ${\mathbf{Q}}_k$ is the process noise covariance matrix.

The observation equation relates the state to the measured quantities through a nonlinear function $\mathcal{F}(·)$

$${\mathbf{z}}_k=\mathcal{F}\left({\mathbf{x}}_k\right)+{\mathbf{v}}_k,{\mathbf{v}}_k\sim \mathcal{N}(\mathbf{0},{\mathbf{R}}_k),$$

where ${\mathbf{z}}_k\in {\mathbb{R}}^N$ is the observation vector (e.g., hydraulic head measurements), and ${\mathbf{R}}_k$ is the observation noise covariance matrix. The nonlinear function $\mathcal{F}\left({\mathbf{x}}_k\right)$ to ${\mathbf{z}}_k$ is defined by the solution solver of the Richards equation. Or, one can learn this mapping using an operator learning method such as the DeepONet approach.

This inverse problem differs from the one discussed in the previous section in the choice of the observation variable ${\mathbf{z}}_k$. Given the time dependent nature of the pressure head values in unsaturated soil, we present two choices. In the first choice, we lump all the available, observed pressure head values at $({\mathbf{x}}_j,t_k),j=n_f+1,\dots ,N,k=0,\dots ,N_t$ to one vector ${\mathbf{z}}_l\in {\mathbf{R}}^{n_h{N}_t}$, where $l=1,2,\dots ,$ is the EKF iteration index. Each time, $\mathcal{F}$ maps ${\mathbf{x}}_l$ into ${\mathbf{z}}_l$ for $l=0,1,2,\dots .$ We hope ${\mathbf{x}}_l$ converges to a limit in the EKF recursion, from which the hydraulic conductivity can be calculated. The second choice is to align the time step index in $t_k$ with the iteration index so that ${\mathbf{z}}_k={\left(h({\mathbf{x}}_{j-n_f},t_k)\right)}^T$. This will require either the measurement to be consecutive without a limit, i.e., $k=1,2,\dots .$ This continuous measurement is the prerequisite for the convergence of ${\mathbf{x}}_k$. Or alternatively, one must introduce a second index for internal iteration using EKF: ${\mathbf{z}}_k^{\left(l\right)},l=1,\dots$. In the following, our proposed methods apply to all the above cases. However, we will use the last alternative choice as an example.

In practical applications to flows in unsaturated soil, the choice of the covariance matrices ${\mathbf{Q}}_k$ and ${\mathbf{R}}_k$ requires careful consideration. The process noise covariance ${\mathbf{Q}}_k$ is typically set to a small diagonal matrix, ${\mathbf{Q}}_k={\sigma }_q^2\mathbf{I}$, where ${\sigma }_q^2$ is chosen to reflect the expected magnitude of unmodeled temporal variations in the log-conductivity field. A common practice is to take ${\sigma }_q^2$ as a fraction of the diagonal entries of the prior covariance ${\mathbf{P}}_0^{-}$, so that the EKF iterations remain numerically stable without artificially inflating the uncertainty. The observation noise covariance ${\mathbf{R}}_k$ should account for both measurement error and model error arising from the simplified representation of the soil-water retention curve, the finite representation of the saturated hydraulic conductivity, and the numerical solver. When hydraulic head observations are collected at multiple time instants, ${\mathbf{R}}_k$ may be taken as a diagonal matrix whose entries are proportional to the variance of the head measurements at each location, or as a constant multiple of the identity matrix if the measurement accuracy is uniform across sensors. In either case, the magnitude of ${\mathbf{R}}_k$ controls the relative weight assigned to the observation term in the cost function (44) and directly influences the convergence rate and the smoothness of the reconstructed conductivity field.

Initialization of Expansion Coefficients Using Co-Kriging

Recall from Section 3.4 that the SLE initial estimate is obtained by Co-Kriging the log-conductivity and hydraulic head measurements. In the present setting, the unknown quantity is the coefficient vector $\mathbf{c}={[c_1,\dots ,c_M]}^T$ of the expansion. Because $\mathbb{E}\left[\mathbf{c}\right]=\mathbf{0}$, the estimator takes the form of a Co-Kriging

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{\mathbf{c}}}_0={\Lambda }^T{\mathbf{f}}^{*}+{\mathbf{M}}^T{\mathbf{h}}^{*},\hfill \end{array}\end{array}$$

(4.11)

where ${\mathbf{f}}^{*}\in {\mathbb{R}}^{n_f}$ contains the de-meaned log-conductivity measurements,

$$f_i^{*}=ln{K}_s^{\mathrm{obs}}\left({\mathbf{x}}_i\right)-F\left({\mathbf{x}}_i\right),i=1,\dots ,n_f,$$

and ${\mathbf{h}}^{*}\in {\mathbb{R}}^{n_h}$ contains the de-meaned hydraulic head measurements,

$$h_j^{*}={\varphi }^{\mathrm{obs}}\left({\mathbf{x}}_j\right)-\mathbb{E}\left[\varphi \left({\mathbf{x}}_j\right)\right],j=1,\dots ,n_h.$$

where $\mathbb{E}\left[\varphi \left({\mathbf{x}}_j\right)\right]$ is the prior mean hydraulic head at ${\mathbf{x}}_j$, obtained by solving the Richards equation (96) with the prior mean conductivity field $K_s\left(\mathbf{x}\right)=exp\left(F\left(\mathbf{x}\right)\right)$ and the prescribed retention curve.

In practice, the total hydraulic conductivity $K(\theta ,x)$ is measured at a set of locations ${\left\{x_i\right\}}_{i=1}^{n_f}$. At these locations, the pressure head $\varphi$ is not necessarily measured simultaneously. To extract the corresponding saturated hydraulic conductivity, the following procedure is used. First, the Richards equation (96) is solved with the prior mean conductivity field $K_s\left(x\right)=exp\left(F\left(x\right)\right)$ and the van Genuchten-Mualem model (103)–(104) to obtain the prior pressure head field ${\varphi }^{\left(0\right)}\left(x\right)$. At each measurement location $x_i$, the relative conductivity is then evaluated as

$$K_r^{\mathrm{VG}}\left({\varphi }^{\left(0\right)}\left(x_i\right)\right)=S_e^{\ell }{\left[1-{\left(1-S_e^{1/m}\right)}^m\right]}^2,$$

where $S_e$ is computed from ${\varphi }^{\left(0\right)}\left(x_i\right)$ via (103). The corresponding saturated hydraulic conductivity is obtained as

$$\begin{array}{c}\hfill \begin{array}{c}{K}_s^{\mathrm{obs}}\left(x_i\right)={\displaystyle \frac{K^{\mathrm{obs}}({\theta }_i,x_i)}{K_r^{\mathrm{VG}}\left({\varphi }^{\left(0\right)}\left(x_i\right)\right)}},ln{K}_s^{\mathrm{obs}}\left(x_i\right)=ln{K}^{\mathrm{obs}}({\theta }_i,x_i)-ln{K}_r^{\mathrm{VG}}\left({\varphi }^{\left(0\right)}\left(x_i\right)\right).\hfill \end{array}\end{array}$$

(4.12)

These $ln{K}_s^{\mathrm{obs}}$ values are then de-meaned by subtracting the prior mean $F\left(x_i\right)$ to obtain the fluctuation measurements $f_i^{*}$ used in the Co-Kriging initialization. This procedure is physically consistent: the pressure head used to compute $K_r$ is obtained from the same forward model that governs the flow, rather than from a purely statistical interpolation. If pressure head measurements are available at the conductivity measurement locations, they can be used directly in place of ${\varphi }^{\left(0\right)}\left(x_i\right)$ to avoid the dependence on the prior mean field.

The weight matrices $\Lambda \in {\mathbb{R}}^{n_f\times M}$ and $\mathbf{M}\in {\mathbb{R}}^{n_h\times M}$ are obtained by minimizing the trace of the estimation error covariance $\mathbb{E}\left[(\mathbf{c}-{\widehat{\mathbf{c}}}_0){(\mathbf{c}-{\widehat{\mathbf{c}}}_0)}^T\right]$. The resulting system is

$$\begin{array}{c}\hfill \begin{array}{c}\left[\begin{array}{cc}{\mathbf{C}}_{ff}^{\mathrm{obs}}& {\mathbf{C}}_{fh}^{\mathrm{obs}}\\ {\mathbf{C}}_{hf}^{\mathrm{obs}}& {\mathbf{C}}_{hh}^{\mathrm{obs}}\end{array}\right]\left[\begin{array}{c}\Lambda \\ \mathbf{M}\end{array}\right]=\left[\begin{array}{c}{\Phi }_f{\Sigma }_c\\ {\mathbf{C}}_{ch}^T\end{array}\right],\hfill \end{array}\end{array}$$

(4.13)

where ${\mathbf{C}}_{ff}^{\mathrm{obs}}\in {\mathbb{R}}^{n_f\times n_f}$ is the auto-covariance matrix of the log-conductivity fluctuations at the measurement locations,

$${\mathbf{C}}_{ff}^{\mathrm{obs}}={\Phi }_f{\Sigma }_c{\Phi }_f^T,$$

with ${\Phi }_f\in {\mathbb{R}}^{n_f\times M}$ defined by ${\left[{\Phi }_f\right]}_{ik}={\phi }_k\left({\mathbf{x}}_i\right)$; ${\mathbf{C}}_{hh}^{\mathrm{obs}}\in {\mathbb{R}}^{n_h\times n_h}$ is the auto-covariance matrix of the hydraulic head perturbations; ${\mathbf{C}}_{fh}^{\mathrm{obs}}\in {\mathbb{R}}^{n_f\times n_h}$ is the cross-covariance matrix between the log-conductivity fluctuations and the head perturbations. ${\mathbf{C}}_{fh}^{\mathrm{obs}}$ can be evaluated via first-order sensitivity as

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{C}}_{fh}^{\mathrm{obs}}\approx {\Phi }_f{\Sigma }_c{\mathbf{S}}^T,\hfill \end{array}\end{array}$$

(4.14)

where $\mathbf{S}\in {\mathbb{R}}^{n_h\times M}$ is the sensitivity matrix with entries $S_{jk}=\partial h_j/\partial c_k$. Note that this expression is consistent with ${\mathbf{C}}_{ch}={\Sigma }_c{\mathbf{S}}^T$ below and the relation ${\mathbf{C}}_{fh}^{\mathrm{obs}}={\Phi }_f{\mathbf{C}}_{ch}$. ${\Sigma }_c\in {\mathbb{R}}^{M\times M}$ is the prior covariance of the expansion coefficients; ${\mathbf{C}}_{ch}\in {\mathbb{R}}^{M\times n_h}$ is the cross-covariance matrix between the coefficients $\mathbf{c}$ and the head perturbations. Its $(i,j)$ entry is

$${\left[{\mathbf{C}}_{ch}\right]}_{ij}=\mathrm{Cov}[c_i,h_j^{*}]\approx {\displaystyle \sum _{k=1}^M}{\left[{\Sigma }_c\right]}_{ik}{\displaystyle \frac{\partial h_j}{\partial c_k}},$$

where $\partial h_j/\partial c_k$ is the sensitivity of the hydraulic head at observation point ${\mathbf{x}}_j$ with respect to the coefficient $c_k$, evaluated by solving the adjoint of the discretized Richards equation.

The system (108) is solved once. The corresponding initial error covariance is

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{P}}_{0,c}={\Sigma }_c-\left[\begin{array}{cc}{\Lambda }^T& {\mathbf{M}}^T\end{array}\right]\left[\begin{array}{cc}{\mathbf{C}}_{ff}^{\mathrm{obs}}& {\mathbf{C}}_{fh}^{\mathrm{obs}}\\ {\mathbf{C}}_{hf}^{\mathrm{obs}}& {\mathbf{C}}_{hh}^{\mathrm{obs}}\end{array}\right]\left[\begin{array}{c}\Lambda \\ \mathbf{M}\end{array}\right]\in {\mathbb{R}}^{M\times M}.\hfill \end{array}\end{array}$$

(4.15)

Expanding (110) gives the four-term expression

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{P}}_{0,c}={\Sigma }_c-{\Lambda }^T{\mathbf{C}}_{ff}^{\mathrm{obs}}\Lambda -{\Lambda }^T{\mathbf{C}}_{fh}^{\mathrm{obs}}\mathbf{M}-{\mathbf{M}}^T{\mathbf{C}}_{hf}^{\mathrm{obs}}\Lambda -{\mathbf{M}}^T{\mathbf{C}}_{hh}^{\mathrm{obs}}\mathbf{M}.\hfill \end{array}\end{array}$$

(4.16)

Equations (106) and (111) provide the initial state and error covariance for the EKF recursion described below. If the covariance of the log-conductivity field on the computational grid is needed (e.g., for visualization or uncertainty quantification), it can be recovered by a post-processing step

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{P}}_{ln{K}_s}=\Phi {\mathbf{P}}_{0,c}{\Phi }^T,\hfill \end{array}\end{array}$$

(4.17)

where $\Phi \in {\mathbb{R}}^{N\times M}$ is the basis matrix evaluated at all grid nodes, with entries ${[\Phi ]}_{ai}={\phi }_i\left({\mathbf{x}}_a\right)$.

It is worth noting that the EKF recursion could in principle be started from an uninformative prior, e.g., ${\widehat{c}}_0=\mathbf{0}$ and $P_{0,c}={\Sigma }_c$, without performing the Co-Kriging initialization. In that case, the initial fluctuation field would be zero everywhere, and the EKF would rely entirely on the first few correction steps to adjust the coefficients toward plausible values. However, the Co-Kriging-based initialization incorporates the available $ln{K}_s$ and hydraulic head measurements into the initial estimate of c and reduces the prior covariance according to the explanatory power of these observations. This yields a starting field that is physically more consistent with the observed data and typically accelerates the convergence of the iterative estimator. The advantage is most pronounced when direct conductivity measurements are sparse, because the head observations, which are more abundant, help constrain the initial coefficients through the cross-covariance structure encoded in $C_{fh}^{\mathrm{obs}}$ and $C_{ch}$.

At each EKF iteration, the Richards equation is solved numerically with the reconstructed $K_s^{\left(k\right)}\left(x\right)K_r\left(\theta (x,t)\right)$ to obtain the pressure head field ${\varphi }^{\left(k\right)}(x,t)$. The van Genuchten-Mualem model provides the constitutive relations $\theta \left(\varphi \right)$ and $K_r\left(\varphi \right)$, which are evaluated at the current iterate of $\varphi$ within the nonlinear solver.

Within the EKF framework, the dependence of the observation operator $\mathcal{F}$ on the state vector ${\mathbf{x}}_k$ is mediated by the parameterization in (101)-(104). At each EKF iteration, the Jacobian matrix

$${\mathbf{J}}_k={\left.{\displaystyle \frac{\partial \mathcal{F}}{\partial \mathbf{x}}}\right|}_{{\widehat{\mathbf{x}}}_k^{-}}\in {\mathbb{R}}^{n_h\times M}$$

is evaluated by applying the chain rule through the Richards equation solver. Specifically, for the i-th component of the state vector,

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\partial \mathcal{F}}{\partial c_i}}={\displaystyle \frac{\partial \mathcal{F}}{\partial K}}·{\displaystyle \frac{\partial K}{\partial (ln{K}_s)}}·{\displaystyle \frac{\partial (ln{K}_s)}{\partial c_i}}={\displaystyle \frac{\partial \mathcal{F}}{\partial K}}·K(\theta ,\mathbf{x})·{\phi }_i\left(\mathbf{x}\right),\hfill \end{array}\end{array}$$

(4.18)

where $\partial \mathcal{F}/\partial K$ is obtained from the adjoint of the discretized Richards equation (96). The derivatives $d\theta /d\varphi$ and $\partial K/\partial \theta$ required by the adjoint solver are supplied analytically from the van Genuchten-Mualem formulas (103)–(104). Specifically, when the Richards equation (96) is discretized, $d\theta /d\varphi$ enters through the time derivative term $\partial \theta /\partial t$, and $\partial K/\partial \theta$ enters through the flux term $\nabla ·\left[K\right(\theta ,\mathbf{x})\nabla (\varphi +z\left)\right]$; both are required by the adjoint solver when computing ${\mathbf{J}}_k$. These analytical derivatives improve the accuracy and efficiency of the Jacobian computation compared with purely numerical differentiation.

Extending the static EKF formulation to the transient unsaturated case in the first choice of the observation variable ${\mathbf{z}}_k$ is straightforward. However, it requires assimilating observations at multiple time instants $t_k$ in the second choice of ${\mathbf{z}}_k$. Here the iteration index k is aligned with the temporal index of observations, and inner/outer iterations are introduced to handle the strong nonlinearity of the Richards equation[40]. Then the EKF prediction-correction procedure with inner iterations for the unsaturated case proceeds as follows. Let the subscript k denote the time step index, and let the superscript $(\ell )$ denote the inner iteration index within the same time step. At each time step k, hydraulic head observations ${\mathbf{z}}_k={\left[{\varphi }^{\mathrm{obs}}(x_j,t_k)\right]}_{j=1}^{n_h}$ collected at time $t_k$ are assimilated. The observation operator ${\mathcal{F}}_k(·)$ solves the transient Richards equation over the interval $[t_{k-1},t_k]$ with the initial condition inherited from the previous assimilation step, and the Jacobian $J_k$ is evaluated at the current state estimate and time $t_k$.

1. Prediction Step For each time step $k=1,2,\dots$

$$\begin{array}{cc}\hfill \mathrm{Prior}\mathrm{state}\mathrm{estimate}\text{:}& {\widehat{\mathbf{x}}}_k^{-}={\widehat{\mathbf{x}}}_{k-1},\hfill \\ \hfill \mathrm{Prior}\mathrm{error}\mathrm{covariance}\text{:}& {\mathbf{P}}_k^{-}={\mathbf{P}}_{k-1}+{\mathbf{Q}}_k.\hfill \end{array}$$

(4.19)

2. Correction Step with inner iterations Set ${\widehat{x}}_k^{\left(0\right)}={\widehat{x}}_k^{-}$. For $\ell =1,2,\dots ,L_{max}$, Compute the Jacobian matrix of the observation function

$${\mathbf{J}}_k^{(\ell )}={\displaystyle \frac{\partial {\mathcal{F}}_k}{\partial \mathbf{x}}}{|}_{{\widehat{\mathbf{x}}}_k^{(\ell -1)}},$$

$$\begin{array}{cc}\hfill \mathrm{Kalman}\mathrm{gain}\text{:}& {\mathbf{K}}_k^{(\ell )}={\mathbf{P}}_k^{-}{\left[{\mathbf{J}}_k^{(\ell )}\right]}^T{\left({\mathbf{J}}_k^{(\ell )}{\mathbf{P}}_k^{-}{\left[{\mathbf{J}}_k^{(\ell )}\right]}^T+{\mathbf{R}}_k\right)}^{-1},\hfill \\ \hfill \mathrm{Posterior}\mathrm{state}\mathrm{estimate}\text{:}& {\widehat{\mathbf{x}}}_k^{(\ell )}={\widehat{\mathbf{x}}}_k^{(\ell -1)}+{\mathbf{K}}_k^{(\ell )}\left({\mathbf{z}}_k-{\mathcal{F}}_k\left({\widehat{\mathbf{x}}}_k^{(\ell -1)}\right)\right),\hfill \\ \hfill \mathrm{Posterior}\mathrm{error}\mathrm{covariance}\text{:}& {\mathbf{P}}_k^{(\ell )}=(\mathbf{I}-{\mathbf{K}}_k^{(\ell )}{\mathbf{J}}_k^{(\ell )}){\mathbf{P}}_k^{-}.\hfill \end{array}$$

(4.20)

At each inner iteration, the current ${\varphi }^{(\ell )}$ is used to evaluate $K_r\left(\varphi \right)$ via (103)–(104), so that $K=K_s·K_r\left(\varphi \right)$ is self-consistently updated within the nonlinear solver.

The inner iteration terminates when

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\parallel {\widehat{x}}_k^{(\ell )}-{\widehat{x}}_k^{(\ell -1)}{\parallel }_2}{\parallel {\widehat{x}}_k^{(\ell -1)}{\parallel }_2}}<{\epsilon }_{\mathrm{inner}},\hfill \end{array}\end{array}$$

(4.21)

or when $\ell =L_{max}$. Denote by L ($L\le L_{max}$) the actual number of inner iterations executed before the stopping criterion (116) is met. The posterior estimate and covariance at time step k are then

$$\begin{array}{c}\hfill \begin{array}{c}{\widehat{\mathbf{x}}}_k={\widehat{\mathbf{x}}}_k^{\left(L\right)},{\mathbf{P}}_k={\mathbf{P}}_k^{\left(L\right)},\hfill \end{array}\end{array}$$

(4.22)

In practice, ${\mathbf{J}}_k^{(\ell )}$ may be recomputed only when the relative change of ${\widehat{\mathbf{x}}}_k^{(\ell )}$ exceeds a threshold, to balance accuracy and computational cost.

Because the observation operator ${\mathcal{F}}_k$ defined by the transient Richards equation is strongly nonlinear, a single sequential pass over $k=1,\dots ,N_t$ may terminate in a locally linearized solution. To mitigate this, one may perform multiple passes over the same spatiotemporal dataset, each pass using the posterior of the previous pass as the prior, until the estimate stabilizes.

3. Stopping criterion for the outer loop. Let ${\widehat{\mathbf{x}}}_{N_t}^{\left(p\right)}$ denote the posterior estimate of the state vector x after the p-th complete pass over all $N_t$ time steps. In each pass, the EKF proceeds sequentially from $k=1$ to $k=N_t$ according to the prediction–correction procedure described above, and ${\widehat{\mathbf{x}}}_{N_t}^{\left(p\right)}$ is the final estimate produced at the end of the p-th pass. The subscript $N_t$ therefore indicates that the estimate is obtained after assimilating the entire set of spatiotemporal observations, not merely the observation at the single time instant $t_{N_t}$.

The outer loop is terminated when the relative change between two consecutive passes falls below a prescribed tolerance,

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\parallel {\widehat{\mathbf{x}}}_{N_t}^{\left(p\right)}-{\widehat{\mathbf{x}}}_{N_t}^{(p-1)}{\parallel }_2}{\parallel {\widehat{\mathbf{x}}}_{N_t}^{(p-1)}{\parallel }_2}}<{\epsilon }_{\mathrm{outer}},\hfill \end{array}\end{array}$$

(4.23)

or when a preset maximum number of passes $P_{max}$ is reached. The first pass ($p=1$) starts from the Co-Kriging initialization ${\widehat{\mathbf{x}}}_0={\widehat{\mathbf{c}}}_0$ and ${\mathbf{P}}_0={\mathbf{P}}_{0,c}$, and proceeds through $k=1,\dots ,N_t$ to produce ${\widehat{\mathbf{x}}}_{N_t}^{\left(1\right)}$. A typical choice for the inner tolerance is ${\epsilon }_{\mathrm{inner}}={10}^{-3}$ to ${10}^{-4}$, with a maximum of $L_{max}=5$–10 inner iterations per time step. The outer-loop tolerance is set to the same range, ${\epsilon }_{\mathrm{outer}}={10}^{-3}$ to ${10}^{-4}$. These values are standard for iterative EKF schemes applied to nonlinear inverse problems governed by the Richards equation[40]. The inner criterion controls convergence of the local linearized update at each time step, while the outer criterion assesses global stability across the entire assimilation pass. In practice, relative changes below ${10}^{-3}$ already correspond to negligible adjustments of the estimated $ln{K}_s$ field and the associated hydraulic-head predictions. Together, the two criteria ensure that the estimated field has stabilized under repeated assimilation of the same dataset, mitigating the effect of the local linearization inherent in the EKF and promoting convergence toward a consistent solution.

4.2. Data-Driven Parameter Estimation for Flows in Unsaturated Soil Using the Ensemble Kalman Filter

The preceding sections have focused on forward modeling under the assumption that the hydraulic conductivity follows prescribed parametric forms. In practice, however, the van Genuchten–Mualem model, while widely used, imposes a rigid relationship between the hydraulic conductivity K and the pressure head $\varphi$ that may not hold in heterogeneous or structured soils. A purely parametric inversion that treats the VGM parameters as unknown still inherits this structural constraint, and the estimated K field is forced into a specific functional shape that may deviate from the true behavior.

An alternative is to abandon any parametric assumption and learn $K(\mathbf{x},\varphi )$ directly from data. This is appealing in its generality, but it introduces a severe ill-posedness: the unknown $K(\mathbf{x},\varphi )$ lives in a high-dimensional space indexed by both spatial position and pressure head, and the available observations are far too sparse to constrain it uniquely. Without additional structure, the inversion yields a highly non-unique solution with little physical interpretability.

A pragmatic compromise is to retain the physically interpretable saturation curve $\theta \left(\varphi \right)$—which can be measured independently in the laboratory—while relaxing the VGM constraint on the conductivity. This decomposition is not merely a matter of convenience. A direct Karhunen–Loève expansion of $K(\mathbf{x},\varphi )$ as a single field would require a two-dimensional covariance kernel indexed by both spatial location and pressure head, substantially increasing the number of modes needed to capture the full variability. Moreover, such an approach would discard the strong physical prior provided by the van Genuchten–Mualem model, which prescribes a well-established functional relationship between $K_r$ and $\varphi$. The spatial variability of $K_s\left(x\right)$ and the functional dependence of $K_r\left(\varphi \right)$ operate on different physical scales and carry very different prior knowledge: $K_s\left(x\right)$ is a stationary or quasi-stationary random field amenable to geostatistical modeling, while $K_r\left(\varphi \right)$ is a bounded, monotone function whose shape is well approximated by the van Genuchten–Mualem model. Merging them into a single object discards this structural distinction and forces the inversion to learn both aspects from the same limited data, which is neither efficient nor robust.

By separating the two factors, we retain the VGM prediction as a physically interpretable baseline for the relative permeability and restrict the KL expansion to the deviation from this baseline. The saturated conductivity $K_s\left(\mathbf{x}\right)$, which depends only on spatial position and is already well studied in hydrogeology, is estimated directly on the full grid using the ensemble approximation. In this way, the total number of degrees of freedom is kept moderate, the physical structure of unsaturated flow is preserved, and the data-driven correction is confined to the component that is least constrained by existing theory. Specifically, we decompose the hydraulic conductivity as

$$\begin{array}{c}\hfill \begin{array}{c}K(\mathbf{x},\varphi )=K_s\left(\mathbf{x}\right)K_r(\mathbf{x},\varphi ),\hfill \end{array}\end{array}$$

(4.24)

where $K_s\left(\mathbf{x}\right)$ is the saturated hydraulic conductivity and $K_r(\mathbf{x},\varphi )$ is the relative permeability function. The saturated conductivity $K_s\left(\mathbf{x}\right)$ is unknown and must be inferred from data. For the relative permeability, we treat the VGM prediction as a prior mean but allow spatially varying deviations. This is achieved by expanding $K_r$ in a Karhunen–Loève basis,

$$\begin{array}{c}\hfill \begin{array}{c}{K}_r(\mathbf{x},\varphi )=K_r^{\mathrm{VG}}\left(\varphi \right)+{\displaystyle \sum _{i=1}^m}\sqrt{{\lambda }_i}{\alpha }_i\left(\mathbf{x}\right){\psi }_i\left(\varphi \right),\hfill \end{array}\end{array}$$

(4.25)

where ${\{{\lambda }_i,{\psi }_i\left(\varphi \right)\}}_{i=1}^m$ are the leading eigenvalues and eigenfunctions of the covariance kernel of $K_r$ sampled from a prior ensemble, and ${\alpha }_i\left(x\right)$ are spatially correlated coefficient fields to be estimated. The number of modes m is chosen such that the expansion captures a prescribed fraction of the prior variance. This construction respects the VGM as a baseline while allowing data-driven corrections that are spatially heterogeneous and physically plausible.

The covariance kernel for the deviation of $K_r$ from the VGM baseline is prescribed as a stationary exponential form,

$$\begin{array}{c}\hfill \begin{array}{c}{C}_{K_r}(\varphi ,{\varphi }^{\prime })={\sigma }_{K_r}^{2}exp\left(-{\displaystyle \frac{|\varphi -{\varphi }^{\prime }|}{l_{K_r}}}\right),\hfill \end{array}\end{array}$$

(4.26)

where ${\sigma }_{K_r}^2$ denotes the prior variance and $l_{K_r}$ the correlation length in pressure-head space. This choice follows standard geostatistical practice for deriving KL expansions in hydrogeological parameter fields [6], and ensures positive definiteness and smoothness consistent with laboratory retention measurements.

The eigenvalues ${\lambda }_i$ and eigenfunctions ${\psi }_i\left(\varphi \right)$ are obtained by solving the eigendecomposition of a prescribed covariance kernel for the deviation of $K_r$ from the VGM baseline, either analytically from a chosen kernel $C_{K_r}(\varphi ,{\varphi }^{\prime })$ or empirically from a prior ensemble of plausible $K_r\left(\varphi \right)$ curves. This decomposition is performed before the assimilation, and the resulting ${\psi }_i\left(\varphi \right)$ are held fixed thereafter; they serve as a fixed functional basis on the $\varphi$-axis along which the spatially varying correction ${\alpha }_i\left(\mathbf{x}\right)$ acts. The temporal variation of $K_r(\mathbf{x},\varphi (\mathbf{x},t))$ thus arises solely from the transient $\varphi (\mathbf{x},t)$ passed through these fixed ${\psi }_i$ and the fixed VGM function $K_r^{\mathrm{VG}}\left(\varphi \right)$.

Unlike the EKF formulation in Section 4.1, where the state dimension must be kept small because the Jacobian matrix and the full covariance matrix must be computed and stored explicitly, the EnKF circumvents these requirements by approximating the state covariance through an ensemble of realizations. Consequently, $ln{K}_s\left(\mathbf{x}\right)$ can be estimated directly on the full computational grid without a separate KL expansion. This removes the need to impose a low-dimensional representation on the saturated conductivity field. Instead, the flexibility gained from the ensemble approximation is allocated to the relative permeability $K_r(\mathbf{x},\varphi )$, which is allowed to deviate from the VGM baseline through a KL expansion. In this way, the saturated conductivity retains its full spatial degrees of freedom, while the unsaturated conductivity benefits from a data-driven correction that is both parsimonious and physically constrained.

The relative permeability $K_r(\mathbf{x},\varphi )$ is required to be non-decreasing in $\varphi$ and to satisfy $0<K_r(\mathbf{x},\varphi )\le 1$, with $K_r(x,0)=1$. To enforce this physically motivated constraint within the ensemble Kalman filter, the Karhunen–Loève basis $\left\{{\psi }_i\left(\varphi \right)\right\}$ is constructed from a prior ensemble restricted to monotone $K_r$ curves consistent with the van Genuchten retention model. In addition, after each analysis step, ensemble members for which the reconstructed $K_r(\mathbf{x},\varphi )$ violates monotonicity or falls outside $(0,1]$ are projected back onto the nearest admissible curve by monotone isotonic correction and clipping. The saturated conductivity is handled in log-space and clipped to a physically reasonable range. These procedures ensure that all realizations entering the forward Richards solver remain thermodynamically consistent.

The unknown parameter field now consists of $ln{K}_s\left(\mathbf{x}\right)$ and the m coefficient fields ${\left\{{\alpha }_i\left(\mathbf{x}\right)\right\}}_{i=1}^m$. The total number of degrees of freedom is $(1+m)N_x$, where $N_x$ is the number of spatial grid points. Even for moderate m and a modest spatial discretization, the state dimension easily exceeds several thousand. This rules out methods that require explicit computation of the Jacobian, such as the extended Kalman filter, which would demand $O\left(N_s\right)$ forward solves per assimilation step and the storage of an $N_s\times N_s$ covariance matrix.

The ensemble Kalman filter offers a natural alternative. It approximates the state covariance using an ensemble of $N_e$ realizations, requiring only $N_e$ forward solves per step—typically two orders of magnitude fewer than the state dimension. No tangent-linear or adjoint model is needed, making the method agnostic to the complexity of the forward solver. The ensemble representation also captures nonlinear propagation of uncertainty without linearization, which is critical for the strongly nonlinear Richards equation. For these reasons, the EnKF is adopted as the inference engine in this work.

4.2.1. State Vector and Observation Operator

Let the computational domain be discretized into $N_x$ grid cells with grid points ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_{N_x}$. The state vector at assimilation step k is defined as

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{m}}_k={\left[ln{K}_s\left({\mathbf{x}}_1\right),\dots ,ln{K}_s\left({\mathbf{x}}_{N_x}\right),{\alpha }_1\left({\mathbf{x}}_1\right),\dots ,{\alpha }_1\left({\mathbf{x}}_{N_x}\right),\dots ,{\alpha }_m\left({\mathbf{x}}_1\right),\dots ,{\alpha }_m\left({\mathbf{x}}_{N_x}\right)\right]}^T\in {\mathbb{R}}^{N_s},\hfill \end{array}\end{array}$$

(4.27)

where $N_s=(1+m)N_x$. The logarithmic transformation ensures positivity of $K_s$ and improves the Gaussianity of the prior.

Given a state realization $\mathbf{m}$, the hydraulic conductivity field is reconstructed as

$$\begin{array}{c}\hfill \begin{array}{c}K(x,\varphi )=exp\left(ln{K}_s\left(x\right)\right)\left[K_r^{\mathrm{VG}}\left(\varphi \right)+{\displaystyle \sum _{i=1}^m}\sqrt{{\lambda }_i}{\alpha }_i\left(x\right){\psi }_i\left(\varphi \right)\right],\hfill \end{array}\end{array}$$

(4.28)

and the water retention curve $\theta \left(\varphi \right)$ is taken from the independently calibrated VGM. Substituting these into the Richards equation and solving numerically yields the predicted pressure head field $\varphi (x,t)$. The observation operator $\mathcal{H}$ extracts the predicted equivalents of the measurements. Let $N_{\mathrm{obs}}^{\mathrm{loc}}$ denote the number of observation locations and $N_t$ the number of assimilation time steps. The total number of observations is $N_{\mathrm{obs}}=N_{\mathrm{obs}}^{\mathrm{loc}}\times N_t$. Each observation location $x_j^{\mathrm{obs}}$ coincides with a computational grid node, so that the observation operator $\mathcal{H}$ simply extracts the simulated pressure head at the corresponding grid point. The observation vector is then

$$\begin{array}{c}\hfill \begin{array}{c}\mathcal{H}\left(\mathbf{m}\right)={\left[\varphi (x_1^{\mathrm{obs}},t_1),\dots ,\varphi (x_{N_{\mathrm{obs}}^{\mathrm{loc}}}^{\mathrm{obs}},t_1),\dots ,\varphi (x_1^{\mathrm{obs}},t_{N_t}),\dots ,\varphi (x_{N_{\mathrm{obs}}^{\mathrm{loc}}}^{\mathrm{obs}},t_{N_t})\right]}^T\in {\mathbb{R}}^{N_{\mathrm{obs}}}.\hfill \end{array}\end{array}$$

(4.29)

Measurements may include pressure heads from tensiometers or water content readings from soil moisture sensors; in the latter case, the conversion $\theta =\theta \left(\varphi \right)$ is applied before comparison.

4.2.2. Iterative EnKF with Parameter Localization

The ensemble Kalman filter is implemented in a standard form with covariance localization to mitigate spurious long-range correlations arising from finite ensemble size [35].

Initialization. Draw $N_e$ samples ${\left\{{\mathbf{m}}_0^{\left(j\right)}\right\}}_{j=1}^{N_e}$ from the prior distribution. For $ln{K}_s$, a Gaussian random field with exponential covariance is used. For the KL coefficients, the prior is ${\alpha }_i\sim \mathcal{N}(0,{\lambda }_i)$, where ${\lambda }_i$ is the eigenvalue associated with ${\psi }_i$.

1. Forecast step. For each ensemble member j, solve the Richards equation over the assimilation window $[t_{k-1},t_k]$ using the reconstructed $K^{\left(j\right)}(x,\varphi )$ and the fixed $\theta \left(\varphi \right)$. Collect the predicted observations,

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{y}}_k^{\left(j\right)}={\mathcal{H}}_k\left({\mathbf{m}}_{k-1}^{\left(j\right)}\right)+{ϵ}_k^{\left(j\right)},{ϵ}_k^{\left(j\right)}\sim \mathcal{N}(0,{\mathbf{R}}_k),\hfill \end{array}\end{array}$$

(4.30)

where ${\mathbf{R}}_k$ is the observation error covariance and the perturbation avoids ensemble collapse.

2. Analysis step. Compute the ensemble mean and covariance,

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{P}}_{my}={\displaystyle \frac{1}{N_e-1}}{\displaystyle \sum _{j=1}^{N_e}}\left({\mathbf{m}}_{k-1}^{\left(j\right)}-{\overline{\mathbf{m}}}_k^{-}\right){\left({\mathbf{y}}_k^{\left(j\right)}-{\overline{\mathbf{y}}}_k^{-}\right)}^T,{\mathbf{P}}_{yy}={\displaystyle \frac{1}{N_e-1}}{\displaystyle \sum _{j=1}^{N_e}}\left({\mathbf{y}}_k^{\left(j\right)}-{\overline{\mathbf{y}}}_k^{-}\right){\left({\mathbf{y}}_k^{\left(j\right)}-{\overline{\mathbf{y}}}_k^{-}\right)}^T.\hfill \end{array}\end{array}$$

(4.31)

where

$$\begin{array}{c}\hfill \begin{array}{c}{\overline{\mathbf{m}}}_k^{-}={\displaystyle \frac{1}{N_e}}{\displaystyle \sum _{j=1}^{N_e}}{\mathbf{m}}_{k-1}^{\left(j\right)},{\overline{\mathbf{y}}}_k^{-}={\displaystyle \frac{1}{N_e}}{\displaystyle \sum _{j=1}^{N_e}}{\mathbf{y}}_k^{\left(j\right)}.\hfill \end{array}\end{array}$$

(4.32)

Apply covariance localization by Schur product with a compactly supported correlation matrix $\mathbf{C}$,

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{P}}_{my}^{\mathrm{loc}}=\mathbf{C}\circ {\mathbf{P}}_{my},{\mathbf{P}}_{yy}^{\mathrm{loc}}={\mathbf{C}}_{\mathrm{obs}}\circ {\mathbf{P}}_{yy}.\hfill \end{array}\end{array}$$

(4.33)

The matrix $\mathbf{C}$ in the Hadamard product $\mathbf{C}\circ {\mathbf{P}}_{my}$ is a compactly supported correlation matrix that implements covariance localization. Its entry ${\mathbf{C}}_{ij}$ encodes the presumed correlation between the i-th state variable and the j-th observation, and is defined as a function of their spatial separation

$${\mathbf{C}}_{ij}=\rho \left({\displaystyle \frac{\parallel {\mathbf{x}}_i^{\mathrm{state}}-{\mathbf{x}}_j^{\mathrm{obs}}\parallel }{L}}\right),$$

where ${\mathbf{x}}_i^{\mathrm{state}}$ is the spatial coordinate associated with the i-th component of the state vector, ${\mathbf{x}}_j^{\mathrm{obs}}$ is the location of the j-th observation, L is the localization radius, and $\rho (·)$ is a positive definite function with compact support[41,42,43]. Similarly, ${\mathbf{C}}_{\mathrm{obs}}$ is defined for observation–observation pairs and enters the Hadamard product with ${\mathbf{P}}_{yy}$.

The purpose of localization is to suppress spurious long-range correlations that arise from sampling error when the ensemble size is much smaller than the state dimension. Without localization, the sample cross-covariance ${\mathbf{P}}_{my}$ contains non-negligible entries between distant locations, causing the analysis to update parameters using observations that are physically unrelated. The Hadamard product attenuates these entries proportionally to distance, retaining only those correlations that are physically plausible. This operation is standard in ensemble-based data assimilation for high-dimensional geophysical problems [42,43].

The Kalman gain is then given by

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{K}}_k={\mathbf{P}}_{my}^{\mathrm{loc}}{\left({\mathbf{P}}_{yy}^{\mathrm{loc}}+{\mathbf{R}}_k\right)}^{-1}.\hfill \end{array}\end{array}$$

(4.34)

3. Update each ensemble member.

$$\begin{array}{c}\hfill \begin{array}{c}{\mathbf{m}}_k^{\left(j\right)}={\mathbf{m}}_{k-1}^{\left(j\right)}+{\mathbf{K}}_k\left({\mathbf{y}}_k^{\mathrm{obs}}-{\mathbf{y}}_k^{\left(j\right)}\right).\hfill \end{array}\end{array}$$

(4.35)

4. Inflation. To prevent underestimation of spread, multiplicative inflation is applied after each analysis:

$${\mathbf{m}}_k^{\left(j\right)}\leftarrow {\overline{\mathbf{m}}}_k+\gamma \left({\mathbf{m}}_k^{\left(j\right)}-{\overline{\mathbf{m}}}_k\right),$$

with $\gamma \in [1,1.05]$ tuned heuristically[37,44].

5. Iteration. Advance to the next assimilation window and repeat until all observations have been assimilated.

The assimilation proceeds sequentially over disjoint time windows $[t_{k-1},t_k]$, $k=1,\dots ,N_t$. Within each window, the forward model is integrated from the initial condition inherited from the previous analysis step, using the most recently updated parameter fields. At the end of the window, the EnKF analysis updates both the parameter fields and the state variables, which then serve as the initial condition for the next window. This procedure assumes that the hydraulic parameters are time-invariant over the assimilation period, which is standard for unsaturated zone inverse problems.

The EnKF simultaneously updates the $ln{K}_s$ field and the KL coefficient fields, thereby learning a spatially heterogeneous relative permeability that can deviate from the VGM prediction wherever the data demand it. The KL truncation level m controls the flexibility of the model and is determined prior to inversion based on the prior spectrum.

4.2.3. Synthetic Experiments

To validate the proposed EnKF-based inversion strategy, synthetic experiments are conducted. A reference parameter field, denoted by $K_{\mathrm{true}}(\mathbf{x},\varphi )$, is first prescribed. The forward Richards solver is then invoked to generate the corresponding pressure head field ${\varphi }_{\mathrm{true}}(\mathbf{x},t)$. Observations are extracted at selected locations and time instants, and corrupted with additive Gaussian noise to mimic measurement errors,

$$y_j^{\mathrm{obs}}={\varphi }_{\mathrm{true}}(x_j^{\mathrm{obs}},t_k)+{ϵ}_j,{ϵ}_j\sim \mathcal{N}(0,{\sigma }_{\mathrm{obs}}^2).$$

The EnKF procedure described in Section 4.2.2 is then applied to estimate $ln{K}_s\left(\mathbf{x}\right)$ and ${\left\{{\alpha }_i\left(x\right)\right\}}_{i=1}^m$ from these synthetic observations. The quality of the inversion is assessed by comparing the estimated fields against the reference via root-mean-square error (RMSE) and correlation coefficient. This workflow allows a controlled evaluation of the inversion performance under various scenarios, including different observation densities, noise levels, and prior specifications.

5. Conclusions

We have established a unified framework linking the successive linear estimator (SLE) and the extended Kalman filter (EKF) for hydraulic conductivity estimation in groundwater flow. Starting from the least-squares formulation of the Kalman filter, we showed that the SLE update—traditionally derived from geostatistical principles—can be reinterpreted as an EKF correction step applied to a static parameter estimation problem, with the co-kriging estimate supplying the initial prior and the sensitivity matrix serving as the linearized observation operator. This equivalence places the SLE within the broader context of recursive Bayesian estimation and provides a rigorous foundation for its extension beyond the saturated steady-state setting.

Building on this equivalence, we extended the EKF framework to unsaturated flow governed by the transient Richards equation. Two complementary strategies were developed. The first strategy employs a low-dimensional basis expansion for the saturated conductivity field and applies an iterated EKF with inner–outer loop convergence criteria to assimilate spatiotemporal head observations. The second strategy adopts the ensemble Kalman filter (EnKF) to circumvent the explicit Jacobian computation required by the EKF, enabling direct estimation of the saturated conductivity on the full computational grid while representing the relative permeability through a KL expansion that allows data-driven deviations from the van Genuchten–Mualem baseline. Covariance localization and multiplicative inflation are incorporated to mitigate sampling errors inherent in finite-size ensembles.

Several methodological issues remain open. First, the choice between the EKF-based and EnKF-based strategies involves a trade-off between computational cost, scalability, and the ability to handle strong nonlinearity. The EKF route is more efficient for low-dimensional parameterizations but requires an adjoint solver; the EnKF route avoids the adjoint but demands careful tuning of localization radius, inflation factor, and ensemble size. Second, the KL expansion for the relative permeability relies on a prescribed covariance kernel whose parameters (variance and correlation length) are themselves uncertain and may benefit from hierarchical Bayesian treatment. Third, the current framework treats the van Genuchten–Mualem parameters as fixed; joint estimation of these parameters alongside the conductivity field warrants further investigation. Fourth, pressure head observations alone may not fully constrain the unsaturated conductivity, suggesting that additional data types—such as flux measurements or soil moisture profiles—should be incorporated in future work.

Future work will focus on numerical evaluation of the proposed framework using both synthetic and field-scale datasets, including a systematic comparison of the two strategies in terms of accuracy, robustness, and computational efficiency.