## 2 Methodology
### The Nested Ensemble Filtering Approach for Parameter Estimation and Uncertainty Quantification  
In the state and parameter estimation process through EnKF, one of the main problems is overshooting, which means that the parameters or states show abnormal values (e.g., negative, complex values) so that the data assimilation process cannot continue anymore, or wrong parameter and state values are generated. Therefore, to overcome this problem, a nested ensemble filtering (NEF) approach will be proposed, in which the SIR procedures will be integrated into the EnKF updating process to eliminate abnormal values [29; 30]. The full description of the nested ensemble filtering processes is illustrated in Figure 1.  
1. Model state initialization: Initialize \(N_{x}\)-dimensional model state variables for \(ne\) samples: \(x_{t,i}\), \(i=1,2,\ldots,ne\), \(x\in R^{N_{x}}\).
2. Parameter sampling: Sample \(N_{\theta}\)-dimensional model parameters for \(ne\) samples: \(\theta_{t,i}\), \(i=1,2,\ldots,ne\), \(\theta\in R^{N_{\theta}}\).
3. Sample weight assignment: Assign the particle weights uniformly: \(w_{t,i}=1/ne\).
4. Model state forecast step: Propagate the \(ne\) state variables and model parameters forward in time using model operator \(f\) \[x_{t+1,i}^{f}=f(x_{t,i}^{a},u_{t,i},\theta_{t+1,i}^{f})+\omega_{t+1,i},\] \[\omega_{t+1}\sim N(0,\sum_{t}^{m}),\ i=1,2,\ldots,\ ne.\]Observation simulation: Use the observation operator \(h\) to propagate the model state forecast:  
\[y^{f}_{t+1,i}=h(x^{f}_{t+1,i},\theta^{f}_{t+1,i})+v_{t+1,i},\]  
\[v_{t+1,i}\sim N(0,\sum\limits_{t+1}^{y}),\;i=1,\;2,\;\dots,\;ne.\]  
6. Parameters and states updating: Update the parameters and states via the EnKF updating equations \[x^{a}_{t+1,i}=x^{f}_{t+1,i}+K_{xy}[y_{t+1}+\varepsilon_{t+1,i}-y^{f}_{t+1,i}]\] \[\theta^{a}_{t+1,i}=\theta^{f}_{t+1,i}+K_{yy}[y_{t+1}+\varepsilon_{t+1,i}-y^ {f}_{t+1,i}].\]  
7. Estimate the likelihood:  
\[L(y_{t+1}|x^{a}_{t+1,i},\theta^{a}_{t+1,i})=\frac{1}{\sqrt{2\pi R_{t+1}}}\exp( -\frac{1}{2R_{t}}[y_{t+1}-h(x^{a}_{t+1,i},\theta^{a}_{t+1,i})]^{2})\]  
\[p(y_{t+1}|x^{a}_{t+1,i},\theta^{a}_{t+1,i})=\frac{L(y_{t+1}|x^{a}_{t+1,i}, \theta^{a}_{t+1,i})}{\sum\limits_{i=1}^{ne}L(y_{t+1}|x^{a}_{t+1,i},\theta^{a }_{t+1,i})}=p(y_{t+1}-h(x^{a}_{t+1,i},\theta^{a}_{t+1,i})|R_{t+1}).\]  
8. Obtain the updated weight for the analyzed ensemble values:  
\[w^{a}_{t+1,i}=\frac{w^{f}_{t+1,i}\cdot p(y_{t+1}-h(x^{a}_{t+1,i},\theta^{a}_{ t+1,i})|R_{t+1})}{\sum\limits_{i=1}^{ne}w^{f}_{t+1,i}\cdot p(y_{t+1}-h(x^{a}_{t+1,i},\theta^{a}_{t+1,i})|R_{t+1})}.\]  
9. Resampling: Apply the resampling procedure proposed by [PERSON] et al. [28] to eliminate the abnormal samples and replace the analyzed \(x^{a}_{t+1,i}\) and \(\theta^{a}_{t+1,i}\).  
10. Parameter perturbation: Take parameter evolution to the next stage by adding small stochastic error around the sample:  
\[\theta^{f}_{t+2,i}=\theta^{a}_{t+1,i}+\varepsilon_{t+1,i},\]  
\[\varepsilon_{t+1,i}\sim N(0,\eta S(\theta^{a}_{t+1,i}))\]  
where \(\eta\) is a hyper-parameter which determines the radius around each sample being explored and \(S(\theta^{a}_{t+1,i})\) is the standard deviation of the analyzed ensemble values.  
11. Check the stopping criterion: If measurement data are still available in the next stage, \(t=t+1\), return to Step 3. Otherwise, stop.  
Figure 1: The flow chart of the nested ensemble filtering method.