## 2 Mathematical Framework
### gPCE-ES and gPCE-MDA  
In this section, we provide the mathematical framework for the Bayesian inverse modelling approach based on gPCE. We aim at estimating the state variable and model parameter distributions where few information is available on such inputs. This is pursued by establishing a Bayesian inverse modelling approach with the model outcome constrained by spatio-temporal observations of measurable quantities. Let us assume that \(\mathbf{d}\in\mathbb{R}^{m}\) is the vector of noisy empirical measurements of the true observable vector \(\mathbf{d}^{t}\in\mathbb{R}^{m}\):  
\[\mathbf{d}=\mathbf{d}^{t}+\mathbf{e}, \tag{24}\]  
where \(\mathbf{d}^{t}=\mathcal{H}(\mathbf{u}^{t})\) with \(\mathcal{H}\) a measurement operator mapping the true state variables \(\mathbf{u}^{t}\in\mathbb{R}^{\ell}\) into the true observable vector. The vector \(\mathbf{e}\in\mathbb{R}^{m}\) is the additive observational error with probability distribution function \(\rho(\mathbf{e})=\prod_{i=1}^{m}\rho(e_{i})\), that is, assuming \(e_{i}\) with \(i=1,...,m\) mutually independent and identically distributed variables. Of course, \(\mathbf{u}^{t}\) is not available because of the uncertainties of the mathematical modelling process. Through inverse modelling we aim at providing the best estimate of the \"true\" state/parameters variables conditioned by (i) available mathematical models, (ii) prior parameters knowledge, and (iii) measurement data. Here, we employ the ES technique and its iterative version in form of ES-MDA [14; 49; 50]. The linear unbiased smoother estimate can be written in a matrix form as:  
\[X^{u}=X^{f}+K(D-HX^{f}), \tag{25}\]  
where \(X^{f}=(U^{f},P^{p})^{T}\in\mathbb{R}^{(\ell+n)\times n_{MC}}\) is the augmented matrix of the forecast ensemble states \(\mathbf{u}^{f}\) and prior parameters \(\mathbf{p}^{p}\), with \(U^{f}\) obtained as the collection of the outcome of \(n_{MC}\) MC simulations and \(P^{p}\) the corresponding prior parameters for every simulations assuming no model errors. The augmented matrix \(X^{u}=(U^{u},P^{u})^{T}\in\mathbb{R}^{(\ell+n)\times n_{MC}}\) is the updated solution by correction of \(X^{f}\) with assimilated measurements. Matrix \(D\in\mathbb{R}^{m\times n_{MC}}\) contains the perturbed measurements matrix with \(\mathbf{e}\sim N(0,R)\) and \(R\in\mathbb{R}^{m\times m}\) the error covariance matrix, while \(H\) is a linearization of the measurement operator \(\mathcal{H}\). The innovation \((D-HX^{f})\) is weighted over the so-called Kalman gain \(K\) calculated as:  
\[K=C^{f}H^{T}(HC^{f}H^{T}+R)^{-1}, \tag{26}\]  
where \(C^{f}\) is the covariance matrix of the forecast augmented state \(X^{f}\). In the form of gPCE-ES, the MC simulation of the forward operator are approximated by straightforward polynomial evaluations to cast the matrix \(X^{f}\) by using the surrogate solution provided by Equation (12). Thus, a large number of MC samples can be adopted for the forecast ensemble and reduce sampling errors [27]. It is sometimes difficult to achieve a proper constraining of the model through the available data with the single-update of ES, especially when the forward model outcome have a strongly non-linear dependence on the uncertain parameters [14; 51]. For this reason, Reference [50] introduced the MDA technique, that was later combined with ES [14], in order to improve the results of a single assimilation for history-matching problems. The ES-MDA is widely recognized as a tool to refine the conditioning of the model especially when non-linearities are significant, even if there is no rigorous proof about its convergence [52]. The idea behind ES-MDA is to repeat the assimilation of the same measurement data for \(n_{MDA}\) times. In order to avoid an over confidence in the available data, the covariance matrix \(R\) is inflated by a coefficient \(\alpha_{MDA}\geq 1\).