An Efficient and Robust Estimation of Spatio‐Temporally Distributed Parameters in Dynamic Models by an Ensemble Kalman Filter

The accuracy of Earth system models is compromised by unknown and/or unresolved dynamics, making the quantification of systematic model errors essential. While a model parameter estimation, which allows parameters to change spatio‐temporally, shows promise in quantifying and mitigating systematic model errors, the estimation of the spatio‐temporally distributed model parameters has been practically challenging. Here we present an efficient and practical method to estimate time‐varying parameters in high‐dimensional spaces. In our proposed method, Hybrid Offline and Online Parameter Estimation with ensemble Kalman filtering (HOOPE‐EnKF), model parameters estimated by EnKF are constrained by results of offline batch optimization, in which the posterior distribution of model parameters is obtained by comparing simulated and observed climatological variables. HOOPE‐EnKF outperforms the original EnKF in synthetic experiments using a two‐scale Lorenz96 model and a simple global general circulation model. One advantage of HOOPE‐EnKF over traditional EnKFs is that its performance is not greatly affected by inflation factors for model parameters, thus eliminating the need for extensive tuning of inflation factors. We thoroughly discuss the potential of HOOPE‐EnKF as a practical method for improving parameterizations of process‐based models and prediction in real‐world applications such as numerical weather prediction.


Introduction
Quantification of systematic model errors is of paramount importance in monitoring, understanding, and predicting Earth systems.The accuracy of process-based models is compromised by unknown and/or unresolved dynamics.Due to limited computational resources, many sub-grid scale parameterizations are required to solve the evolution of complex Earth systems.However, these sub-grid scale parameterizations inevitably introduce structural and parametric errors into the simulations.Although data assimilation has contributed to the prediction of Earth systems by integrating imperfect models and observations, it is essential to quantify and model the errors of process-based models to maximize the potential of data assimilation to accurately estimate the states of Earth systems.
Data-driven estimation of model errors within a data assimilation framework has been intensively investigated.For instance, Pathiraja and van Leeuwen (2022) proposed constructing a conditional model error probability density using kernel density estimation and incorporating it into an ensemble transform Kalman filter.Amemiya et al. (2023) applied a recurrent neural network to learn model errors from analysis increments estimated by local ensemble transform Kalman filter (LETKF).Tomizawa and Sawada (2021) demonstrated that emulating LETKF analysis through reservoir computing can mitigate the impact of model parameter misspecification on predictions (see also Brajard et al. (2020), Dueben and Bauer (2018), Weyn et al. (2019), and Weyn et al (2021)).All these approaches do not require knowledge of unresolved dynamics and can operate with partially observed systems.Although their real-world applications have been limited thus far, these methods are promising approaches to accurately estimate and predict the states of Earth systems by explicitly considering systematic model errors.
In the present paper, we propose estimating systematic model errors by spatio-temporally distributed model parameters.In our approach, these parameters, originally formulated as global fixed parameters in a process-based model, become spatio-temporally distributed, addressing the model bias introduced by their fixation.While this time-varying parameter estimation is particularly beneficial when the model parameters are known to be nonconstant (e.g., land cover in a hydrological model; see Pathiraja et al. 2018aPathiraja et al. , 2018b)), we will demonstrate and discuss the value of this approach even when the source of systematic model errors remains unclear.
In sequential data assimilation, model parameters can be estimated using an augmented state vector method, in which a parameter vector is concatenated to a state vector and then the augmented vector is updated in the analysis step (Moradkhani et al. 2005a).This simultaneous update of state variables and parameters has been successfully applied in various Earth system science domains, such as atmosphere (e.g., Kotsuki et al. 2018), land (e.g., Kurtz et al. 2016), and hydrology (e.g., Vrugt et al. 2013).However, there are few applications involving time-varying parameters in Earth system sciences, although some previous studies successfully estimated them for relatively low dimensional models (e.g., Ruiz et al. 2013a;Pathiraja et al. 2018aPathiraja et al. , 2018b;;Sawada and Hanazaki 2020).The primary challenge in estimating time-varying model parameters through sequential data assimilation is maintaining the appropriate background variance of the estimated parameters.Whenever model parameters are adjusted by data assimilation, the estimated variance of the updated parameters is always smaller than that of the background.To track rapid changes in model parameters, it is necessary to maintain a relatively large spread of adjusted parameters using covariance inflation methods, which is practically difficult since ones lack knowledge of the error growth in model parameters that are originally assumed to be time-invariant.There are many works to stabilize filters to choose the appropriate inflation for parameters.Moradkhani et al. (2005b) simply added Gaussian noise, whose variance is proportional to the background variance, to the adjusted parameters of each ensemble member.Ruiz et al. (2013a) proposed an adaptive inflation method for updating ensemble parameters in LETKF.Kotsuki et al. (2018) suggested maintaining the initial predefined ensemble spread using the relaxation to prior spread method (Whitaker and Hamill 2012).
As an alternative approach to stabilize the inflation of estimated model parameters, Sawada (2022) proposed constraining model parameters estimated by sequential data assimilation to the result of offline batch optimization in which the posterior distribution of model parameters is obtained by comparing simulated and observed climatological (long-term mean) variables.The proposed method, Hybrid Offline and Online Parameter Estimation with Particle Filtering (HOOPE-PF), was successfully applied to lowdimensional systems.The performance of HOOPE-PF is insensitive to hyperparameters of an ensemble inflation method, making it a practical and efficient approach for estimating time-varying parameters.Although HOOPE-PF requires performing computationally expensive offline batch optimization, a growing number of studies have recently used machine learning-based surrogate models to accelerate this process (e.g., Dunbar et al. 2021;Zhang et al. 2020;Sawada 2020;Teixeira Parente et al. 2019;Duan et al. 2017).This makes it promising to combine online data assimilation and offline batch optimization for estimating of time-varying model parameters.However, HOOPE-PF is difficult to apply to high-dimensional Earth system models since it relies on particle filtering.In this study, we aim at applying the same concept to the widely used ensemble Kalman filter (EnKF).Although variants of EnKF have been extensively used to estimate high-dimensional unknown parameters (e.g., Kang et al. 2011;Bellsky et al. 2014;Pulido et al. 2016;Katzfuss et al. 2020;Ruckstuhl and Janjić 2020), the idea of combing offline batch optimization and EnKF to estimate spatio-temporally distributed parameters has yet to be examined.We demonstrate that our proposed EnKF-based method (i.e., HOOPE-EnKF) is useful for the estimation of systematic model errors in a two-scale Lorenz96 model and a simple atmospheric general circulation model.

Overview of HOOPE-EnKF
The schematic of HOOPE-EnKF is shown in Figure 1.As explained in the previous section, the fundamental idea is that the analysis ensemble of model parameters in EnKF incorporates information from a predefined posterior distribution of time-invariant model parameters: (| ̂(0: )) ≅ (  , ) , where  ̂ is a time-invariant model parameter that can accurately simulate the climatology of the systems and (0: ) is observation from time 0 to T. In this study, we approximate (| ̂(0: )) as a Gaussian distribution with mean   and covariance matrix  (all parameters are assumed to be independent so that C is a diagonal matrix).In the next sub-section, we will explain how to obtain (| ̂(0: )) by offline batch optimization.In Section 2.3, our sequential data assimilation method, LETKF, will be formulated.In the section 2.4, we will introduce HOOPE-EnKF as a simple extension of LETKF.

Offline batch optimization
A discrete state-space dynamic system is written as: where () are the state variables at time t,  ̂ is the time-invariant model parameters which are the originally defined parameters in process-based models (see Section 1), () is the external forcing at time t, and () is the noise process representing model error at time t.() denotes the dynamic model.  () is the observation at time t,  is the observation operator, and () is the noise process representing observation error.
In this study, a simulated climatological index,   , is calculated from the long-term timeseries of the simulated observation   (0: ) where   is the simulated climatological index and () is the forward map.Note that   should not be a function of the system's initial conditions since the targeted climatological index for optimizing time-invariant model parameters should not be flowdependent.The purpose of the offline batch optimization in this study is to obtain the posterior distribution of  ̂, or (| ̂(0: )), based on the simulated climatological index   and the observed climatological index   .
The fully Bayesian inference using Equation 3 is impractical because (. ) is computationally expensive.To efficiently estimate (| ̂(0: )), a surrogate model of (. ) was developed.First, ensemble members of model parameters were generated from a parameter space.Second, the long-term integration of the dynamic model with these ensemble members was performed in parallel.Third, the statistical surrogate model was constructed from the data of this ensemble simulation by Gaussian process regression (Rasmussen and Williams 2006).There are two ways to replace the full process-based model integration with the statistical surrogate model.
where  () (. ) is the computationally cheap surrogate model of (. ) .We sampled  ̂ and performed ensemble simulation with the sampled parameters.Then, we fit the sampled relationship between  ̂ and our target metrices using Gaussian process regression with the Matern kernel.While the long-term integration of the dynamic model with input data such as initial and boundary conditions and external forcings is needed to obtain the model's climatology and evaluate ( ̂), the trained surrogate model directly estimates the relationship between model behaviors and parameters without performing the long-term integration.Equation (4a) approximates the simulated climatology by Gaussian process regression while we directly approximate the square difference between simulated and observed climatological index in Equation (4b).Equation (4b) is more suitable if the dimension of   is high.When we fit the square difference between simulation and observation, our dependent variable in the regression is a scaler, which makes the construction of a surrogate model simpler.We adopted (4a) in the two-scale Lorenz 96 model experiment, while (4b) is adopted in the experiment with a simple atmospheric model.
We applied the Metropolis-Hastings algorithm (Hastings, 1970) to draw samples from (| ̂(0: )).We assumed a prior distribution as a trimmed uniform distribution with predefined maximum and minimum values of parameters.Please refer to Sawada (2022) for details of this sampler (see also section 3 for the implementation of the sampler).
Although Sawada (2022)  It should be noted that there are many other practically promising approaches to obtain (  , ) .For instance, iterative ensemble smoothers (Evensen, 2018) and their extensions (e.g., Cleary et al. 2021) are effective in obtaining (  , ) with a relatively smaller number of iterations.If a dynamic model has already been extensively calibrated using a targeted climatological index, the expert knowledge of modelers can be beneficial to directly specify (  , ) (e.g., using their default settings as   ).

Local Ensemble Transform Kalman Filter (LETKF)
In LETKF and any other ensemble Kalman filters, the analysis probability distribution is entirely determined by its first and second moments which can be obtained from the minimum point and its associated Hessian matrix of the following cost function J: where  ̅  is the background ensemble mean of state variables,   is the background error covariance of . is the observation error covariance.The first and second terms in Equation 5 are called background and observation terms, respectively.
To estimate parameters , state augmentation is effective as discussed in Section 1.For the augmented vector, the cost function can be rewritten as: where  ̅  is the background ensemble mean of model parameters,   is the background error covariance of  , and   and   are the background error covariances showing the correlation between state variables and model parameters, which is the key for transferring information from observable variables to model parameters.
Note that these model parameters, , in Equation 6 are no longer time-invariant ones, so that they are updated in the analysis step.
When we redefine the augmented state vector   = (   ), the LETKF's analysis update is as follows (see Hunt et al. 2007 for the complete description): where   and  ̅  are the observation and the ensemble mean of the simulated observation, respectively.k is the ensemble size.The ith columns of   and   are  () −  ̅  and   () −  ̅   , respectively. () is the simulated observation of the ith ensemble member, and   () is the augmented state vector   of the ith background ensemble member. ̅   and  ̅   are the background and analysis ensemble means of the augmented state vector   , respectively.The perturbations of the analysis ensemble members can be found by: The ith column of   is   () −  ̅   , where   () is the augmented state vector of the ith analysis ensemble member.From   , we can immediately obtain the analysis ensemble members of both model state variables and parameters.
As we will show in the following sections, HOOPE-EnKF modifies Equation 6 and does not necessarily change the method for minimizing the cost function.Our proposed method and the results presented in this paper are applicable to other flavors of EnKF.

Filtering (HOOPE-EnKF)
As discussed in Section 2.1, the analysis ensemble of model parameters in HOOPE-EnKF incorporates information from a posterior distribution of time-invariant model parameters: (| ̂(0: )) ≅ (  , ).Since the background term in Equation 6 can be recognized as a regularization term, it is straightforward to incorporate (  , ) as an additional regularization term: Since (  , ) represent the "climatology" of the parameters, we call the third term of the right-hand side of Equation 12the climatological term.To solve Equation 12with the standard EnKF, we will convert this new form to the standard form of the cost function (Equation 6).We identified two approaches to accomplish this task.

Pseudo observation (PSO) approach
As a first approach, we combined the observation and climatological terms in Equation 12to get a more compact form: This form essentially has the same structure as Equations 5 and 6 by defining the "augmented observations"   = (    c ) , the associated observation operator ), and the new observation error covariance   = (     ).The implication of Equation 13 is that we can solve Equation 12with the standard EnKF by assimilating   into the system with variance  .We call this approach "pseudo observation" since we used   as a pseudo observation of model parameters.
It is important to discuss how to implement inflation and localization in the PSO approach.
In this study, we used prior multiplicative inflation which involves multiplying the background covariances with inflation factors.We simply applied the factors to the background perturbations in the PSO approach.In LETKF, observation localization is used.In this study, we let pseudo-observations of   to impact only the associate parameter variables.See also Section 3 for details of the implementation of inflation and localization.

Regression to climatology (RTC) approach
The second idea is to combine the background and climatological terms in Equation 12to form a general quadratic form Due to the presence of   and   , it is relatively complicated to find the analytic forms of the tilde variables on the right-hand side of ( 14).The detailed derivation is given in Appendix A and the result is reproduced here In EnKF, we approximate the pdf associated with the updated background term ( 14) by ensemble samples.In principle, we can draw new samples from the updated covariance (15b) by taking an appropriate square root of this matrix.However, this strategy is very ineffective in practice due to the high computation cost in estimating such a square root.
Note that the computational cost of RTC should be comparable to that of PSO, otherwise we will stick to the PSO approach due to its fast computation.Now note that the updated moments  ̃ and  ̃ does not contain any terms related to  , this suggests an approximation that treats  and  independently.Thus, our strategy is to draw samples of  from the updated marginal distribution ( ̃,  ̃), which is relatively easy since  ̃ is usually a diagonal matrix, while approximating the remaining moments.
Recall that in the PSO method, since the background term is not modified, it is straightforward to use the same samples { () } representing ( ̅  ,   ) as in the standard case of Equation 6.However, in the current method, since we update ( ̅  ,   ) to ( ̃,  ̃) , it is necessary to update { () } to new samples { ̃() } that conform with the updated pdf.We propose using optimal transport to move { () } drawn from ( ̅  ,   ) to { ̃() } that faithfully represent ( ̃,  ̃) .Because ( ̅  ,   ) and ( ̃,  ̃) are the Gaussian distributions, the optimal transport map has an analytical form (Peyré and Cuturi, 2019): where all matrix square roots denote symmetric positive-definite square roots.
The remaining task now is how we draw new samples { ̃() } that approximates moments related to  in (15a) and (15b).Since, it is important to maintain dynamical balance in each member, we simply keep all { () } intact by setting  ̃() =  () .The transformation ( 16) when ,   are diagonal matrices has the following form These updated perturbations enable us to check how our approximation alters the analytic form of the updated background term where all quantities with the suffice app denote the approximated moments.Comparing the approximated form (18) and the exact form (15), it is easy to see that our approximation retains the original uncertainty of the state   −   ( +   ) −1   →   while slightly increases the cross-correlations ( +   ) −1   →  1/2 ( +   ) −1/2   .This entails that a certain inflation has been introduced to analysis increments and analysis perturbations of parameters due to this increase.
This approximation leads to an approximated form of Equation 12 under the standard form of Equation 6: It is possible to replace the first moment  ̅  in ( 19) by the exact one  ̃ in (15a) by simply adding   [ +   ] −1 (  −  ̅  ) to all { () } .However, this may spoil dynamical balance in each member, therefore, we do not apply this kind of bias correction in this study.We call this approach the regression to climatology (RTC) since the prior pdf ( ̅  ,   ) is replaced with its regressed version ( ̃,  ̃).
Compared to the PSO method, it is subtle to apply prior multiplicative inflation to Equation 19.When the multiplicative inflation factor for parameters is defined as   , the correct inflated covariance of parameter,  ̃ , should not be taken as    ̃.It should be derived from Equation 16: which entails: Where  ̃ is the inflated version of the background ensemble mean of parameters after the optimal transportation.As the inflated background perturbations, we can obtain: The RTC approach appears more complicated and requires higher computational costs than the PSO approach.However, the computation can be significantly reduced when all parameters are assumed to be uncorrelated.In this case, the inflated background perturbation for each parameter takes a simple form: where   2 ,   2 are climatological and background error variances of the parameter θ , respectively.We can consider the factor     2 /(  2 +     2 ) as the effective inflation factor in the RTC method.The ensemble spread is inflated when   > thus explain the actual effect of   .We nudge background perturbations toward the climatological pdf (  , ).In the limit   → ∞, we replace all background samples by their climatological counterparts, ensuring the stable performance of our method.

Two-scale Lorenz 96 model
We tested HOOPE-EnKF using the two-scale Lorenz 96 model: The dimension of the large-scale system,   , was set to 9, while that of the small-scale system,   , was set to 20.The external forcing,  , was set to 14.The time scale separation, , was set to 0.7.We set ℎ  and ℎ  to -2 and 1, respectively.These settings are identical to those used by Pathiraja and van Leeuwen (2022).
In our synthetic experiment, the true dynamics are represented by Equations 25-27.We assumed that our available forecast model is the single-scale Lorenz96 model: where .We used the 7200 MTUs timeseries of the two-scale Lorenz96 model as the truth.Observations of the resolved large-scale variables were generated by adding the Gaussian white noise to the true variables.We assumed that [] can be observed at  = 1,2,5,6 every 0.05 MTU.Note that Pathiraja and van Leeuwen (2022) also assumed that the system can be observed at the same grid points.The standard deviation of the Gaussian white noise was set to 0.1, so that we assumed that observation error is 0.1, which is much larger than that assumed by Pathiraja and van Leeuwen (2022).
In the offline batch optimization, we assumed that [] is a global time-invariant parameter,  ̂.First, we generated 100 ensemble members of  ̂ from  ̂= 0 to  ̂= 30 with constant intervals and then ran the single-scale Lorenz96 model for 28800 MTUs.
Second, we calculated a climatological index from the ensemble simulations,   in Equation 3. We used autocorrelation of [] for   .We calculated three autocorrelations with time lags of 2, 3, and 4 observation intervals, which correspond to 0.1, 0.15, and 0.2 MTU, respectively.Note that the dimension of   is 3.The autocorrelations were calculated for the last 2000 MTUs.Third, we smoothed the sampled relationship between  ̂ and   by the Gaussian process regression to construct the surrogate model shown in Equation 4. Fourth, by comparing   estimated by the surrogate model with   (i.e., autocorrelation estimated from observation), we drew the samples of  ̂ by the Metropolis-Hasting algorithm.Because of the small observation error, the effect of added noises to the truth on autocorrelation was negligible.Discarding the first 100,000 samples as the spin-up period, the remaining 400,000 samples are used to estimate (  , ).Note that the offline batch optimization gives a single sample mean and variance although (  , ) includes information of parameters in all grid points.
These single mean and variance are applied to all elements of the vector   and the diagonal part of C. For details on the implementation of the offline batch optimization, see Appendix B.
We implemented three types of LETKF.First, we applied the vanilla LETKF without using the climatology of parameters, (  , ).In this case, we simply solved Equation 6 with the augmented state vector to estimate both state variables and spatio-temporally distributed parameters.This experiment is referred to as the NOHOOPE experiment.Our second and third experiments are HOOPE-EnKF with the PSO approach (Section 2.4.1) and the RTC approach (Section 2.4.2),respectively.They are called the HOOPE-EnKF-PSO and HOOPE-EnKF-RTC experiments.We implemented three LETKFs with ensemble sizes of 10, 20, and 40.
We performed prior multiplicative inflation for all three LETKFs.For the NOHOOPE and HOOPE-EnKF-PSO experiments, the ensemble perturbations were inflated as follows: where   ,   ≥ 1 are multiplicative inflation factors for state variables and parameters, respectively.For the HOOPE-EnKF-RTC experiment, Equation 29 was used to inflate the ensemble perturbations of state variables, while Equations 23-24 were used for inflating parameter perturbations.Please note that we used different multiplicative inflation factors for state variables and parameters.We conducted three LETKFs (i.e.NOHOOPE, HOOPE-EnKF-PSO, and HOOPE-EnKF-RTC) with 31 × 31 = 961 combinations of   and   , with ranges [1.05,2.55]and [1.05,7.05],respectively.In addition, we performed three LETKFs with adaptive inflation factors.We used the Gaussian approach proposed by Miyoshi (2011) to adaptively estimate the multiplicative inflation factors.In this case, different multiplicative inflation factors estimated adaptively were used for state variables and parameters in different grid points.
For the observation localization of LETKF in this study, the inverse of observation covariance matrix  −1 was multiplied by the localization function (): where  is the distance from an observation location to an analyzed grid point.The use of tapering functions has been widely adopted in LETKF applications since Miyoshi and Amane (2006) proposed them. was set to 3.0.The choice of  does not qualitatively change our conclusions (not shown).Note that we consider parameters [] to be located at grid point k, so that the localization is applied in the same manner as for  [𝑘].
As discussed in Section 2.4.1, we allowed pseudo observations to impact only the associate parameter variables.It means that () is set to 0 if the locations of pseudo observations and analyzed parameters are not identical.
We used root-mean-square error (RMSE) and correlation coefficient (R) as the evaluation metrics.From the 7200 MTUs of the nature run, the simulation of the first 2500 MTUs was discarded as a spin-up of data assimilation.In the remaining 4700 MTUs, we compared the analysis mean of state variables and model parameters with the synthetic truth.

Simple atmospheric general circulation model
To demonstrate the effectiveness of HOOPE-EnKF in addressing realistic problems within Earth system sciences, we employed the SPEEDY model (Molteni 2003) We focused on estimating parameters within the convection scheme of the SPEEDY model.This scheme is a simplified version of the mass-flux scheme developed by Tiedke (1993).For a detailed description, please refer to Molteni and Kucharski.The convection scheme is activated in conditionally unstable regions where saturation moist static energy decreases with height and where humidity exceeds a prescribed threshold in both the planetary boundary layer and the layer one level above it.The humidity criteria can be described in the following equation: where   and  −1 are specific humidity in the lowest and second lowest vertical layers, respectively.In the original SPEEDY model, the thresholds,  ℎ1 and  ℎ2 , are defined as functions of saturation specific humidity: where    and  −1  are saturation specific humidity in the lowest and second lowest vertical layers, respectively.  is a parameter.Ruiz et al. (2013b) found that the model's behavior is notably sensitive to   .They assumed that   was unknown and varied over time.They successfully estimated this time-varying   from conventional observations using LETKF.
In the present study, we directly estimated  ℎ1 and  ℎ2 without the knowledge that these thresholds are dependent on saturation specific humidity.Given that saturation specific humidity is a function of temperature and pressure, our target parameter is different in the different model grids and its dimension amounts to 9216 (96 grids × 48 grids × 2 levels), in contrast to Ruiz et al. (2013b) who estimated a global parameter.
Furthermore,  ℎ1 and  ℎ2 temporally change since saturation specific humidity changes.We tested if this spatio-temporally distributed parameter can be estimated under the similar observation settings to Ruiz et al. (2013b).
We performed an observation system simulation experiment.Observations were taken at the grids shown in Figure C1  In the offline batch optimization, we assumed that both  ℎ1 and  ℎ2 represent a global time-invariant parameter, denoted as  ℎ ̂ .First, we generated 100 ensemble members of  ℎ ̂ , ranging from  ℎ ̂= 1.0 to  ℎ ̂= 20.0 with constant intervals and then executed the SPEEDY model.Note that we no longer applied the equation ( 33), and the convection was activated when humidity exceeds a fixed threshold.Second, we calculated a climatological index from the ensemble simulations, represented as   in Equation 3.This index was derived from the year-averaged specific humidity at the observed grids.An average of 6-hourly simulated specific humidity from January 1982 to December 1982 was taken at the observed grids depicted in Figure C1 and then compared with the averaged observed specific humidity.Note that the dimension of   is identical to the number of observable grids and it is much larger than that in the twoscale Lorenz 96 model experiment.Therefore, we decided to construct the surrogate model to directly estimate square differences between simulation and observation.Third, we smoothed the sampled relationship between  ℎ ̂ and   by the Gaussian process regression to construct the surrogate model shown in Equation 4. Fourth, by comparing   estimated by the surrogate model with   , we drew the posterior samples of  ℎ ̂ by the Metropolis-Hasting algorithm.The observation error of   was designated as 1.0 [gkg -1 ] mirroring the original observation error for specific humidity.After discarding the first 100,000 samples as the spin-up period, the remaining 400,000 samples were used to estimate (  , ).For details on the implementation of the offline batch optimization, see Appendix B.
As we mentioned in section 3.1, we performed three flavors of LETKF: NOHOOPE, HOOPE-EnKF-PSO, and HOOPE-EnKF-RTC.Wind components, temperature, and surface pressure were analyzed every 6 hours. ℎ1 and  ℎ2 were concatenated to the state vector and analyzed every 2 days.We found that frequently updating  ℎ1 and  ℎ2 led to unstable estimations.This instability arises because observable state variables are sensitive to  ℎ1 and  ℎ2 only when the model and/or the nature run exhibit convective activity at each model grid.Consequently, we updated them less frequently.Note that we refrained from adjusting specific humidity to increase the sensitivity of  ℎ1 and  ℎ2 to discrepancies between simulation and observation of specific humidity, ensuring accurate parameter estimation.We adopted the same prior multiplicative inflation that we used in the two-scale Lorenz 96 experiments (see section 3.1) although our focus was solely on experiments with adaptive inflation.The localization function was defined as: where  and  are horizontal and vertical distances between observation and analyzed grid points, respectively.The horizontal and vertical localization scales,  ℎ and   , were set to 1000 km and 0.1 (sigma level), respectively.Please refer to section 3.1 for our strategy of localization in the parameter space.The ensemble size was set to 30.

Two-scale Lorenz 96 model
Figure 2 shows the non-parametric posterior distribution of  ̂ , which can accurately simulate autocorrelations of the system.The dashed red line shows the fitted Gaussian distribution.Although the original sampled distribution is skewed, we directly used sample mean and variance for the Gaussian distribution of the climatological parameter (  , ).HOOPE-EnKF-PSO and HOOPE-EnKF-RTC substantially stabilize the filters.They accurately estimate both state variables and model parameters within a wide range of inflation factors.This is the empirical evidence that the performance of our proposed method is not greatly sensitive to the choice of covariance inflation.We confirm that they can work even with extremely large inflation factors for model parameters (e.g., 500) (not shown), since when   → ∞, we replace all background samples by their climatological counterparts, ensuring the stable performance of our methods, as discussed in Section 2.4.2.HOOPE-EnKF-RTC is stabler than HOOPE-EnKF-PSO, as we can achieve low RMSE and high R with wider ranges of inflation factors especially when the ensemble size is large.The optimal inflation factors for the state variables are relatively large, representing the large model errors induced by the parameter fixation.Even if it were possible to extensively tune the inflation factors, the performance of NOHOOPE is significantly worse than HOOPE-EnKF-RTC and HOOPE-EnKF-PSO (Tables 1 and 2).
Note that it is practically difficult to obtain the optimal performance of NOHOOPE due to the high sensitivity of the performance to inflation factors, while Figures 5 and 6 imply that HOOPE-EnKF-PSO and HOOPE-EnKF-RTC can easily achieve performance similar to the optimal one shown in Tables 1 and 2. It is impractical to extensively tune the inflation factors for maximizing the skill in estimating model parameters since model parameters cannot be observed.Therefore, the insensitivity of the performance to the inflation factors in HOOPE-EnKF-PSO and HOOPE-EnKF-RTC is an important advantage.
To avoid tuning the inflation factors, adaptive inflation was applied.Tables 1 and 2 show that NOHOOPE with adaptive inflation does not work with the small ensemble size.
When the ensemble size is large, the adaptive inflation is somehow effective for  To assess three LETKFs in estimating state variables, we used specific humidity at the lowest level.No notable differences were found in the other state variables (not shown).

Simple atmospheric general circulation model
Figure 10a shows RMSE between analysis and synthetic true specific humidity.While two variants of HOOPE-EnKF stably maintain RMSE below 1.0 [gkg -1 ], NOHOOPE occasionally gets higher RMSE after March 1982.
This poor performance of NOHOOPE in estimating low level specific humidity can be attributed to its inaccurate estimation of  ℎ1 and  ℎ2 .Figure 10b shows RMSE between analysis and synthetic true  ℎ1 .The results for  ℎ2 closely mirror those of  ℎ1 (not shown).After March 1982, the parameter estimation of NOHOOPE gets unstable, leading to extremely large RMSE.Conversely, both HOOPE-EnKF-PSO and HOOPE-EnKF-RTC offer accurate and consistent estimations of  ℎ1 .Figure 11 shows the global distribution of   ×    (synthetic true convection thresholds) and estimated  ℎ1 .Both HOOPE-EnKF-PSO and HOOPE-EnKF-RTC reproduce high  ℎ1 along the intertropical convergence zone.In contrast, NOHOOPE produces extremely large positive and even negative  ℎ1 .Figure 12 shows that both HOOPE-EnKF-PSO and HOOPE-EnKF-RTC accurately estimate the meridional shift of  ℎ1 .
Overall, the difference of accuracies between HOOPE-EnKF-PSO and HOOPE-EnKF-RTC is negligible.Their estimated  ℎ1 in mid-latitude and polar regions is consistently greater than the synthetic truth.Given that the error in specific humidity in the tropics is significantly more pronounced than in other regions, offline batch optimization tends to prioritize parameters that minimize the error in the tropics, overlooking the comparatively minor errors elsewhere.Consequently, our estimated posterior distribution of  ℎ ̂ presents considerably higher values than the synthetic true  ℎ1 in mid-latitude and polar regions.HOOPE-EnKF thus inhibits  ℎ1 from reaching such low values suitable for mid-latitude and polar regions.This is the side effect of introducing global climatological parameters to stabilize the filters.

Discussion and Conclusions
We proposed combining offline batch optimization and EnKF toward an (2016) successfully reduced the dimension of targeted parameters, the advantage of HOOPE-EnKF over this method lies in the absence of a need for a priori knowledge of the functional form.Instead, we directly estimate high-dimensional parameters.Kang et al. (2011) and Bellsky et al. (2014) emphasized the significance of localizing the optimization problem of high-dimensional parameters to accurately estimate spatiotemporally distributed parameters.However, Kang et al. (2011) could not achieve sufficient skill when they assumed that their targeted parameters varied over time.
Although Bellsky et al. (2014) successfully solved a problem similar to ours in the twoscale Lorenz 96 model, they posited two constraints: all model grids had to be observed, and localization scales needed to be narrowly set.It implies a necessity to avoid using remote observations.In contrast, our approach is not bound by these constraints.When remote observations are insensitive to parameters, the parameters are forced to align with climatology, ensuring the stability of the filter.Therefore, we can take a risk to use remote observations to constrain parameters, and it is unnecessary for us to observe all model grids.Katzfuss et al. (2020) rigorously derived new filters and smoothers in which highdimensional state variables were estimated by EnKF, while parameters were separately estimated by non-parametric Bayesian inference techniques such as Gibbs sampling and particle filtering.Although their method holds promise for joint state-parameter estimation, it is best suited for scenarios where the number of unknown parameters is not too large, which differs from the context of our study.
We demonstrated that the estimation of spatio-temporally distributed model parameters  12) and method to solve it can be shared with variational methods.
The primary disadvantage of HOOPE-EnKF is its inability to predict systematic model We proposed two variants of HOOPE-EnKF: HOOPE-EnKF-PSO and HOOPE-EnKF-RTC.Both variants perform significantly better than the vanilla LETKF, and their performances are less sensitive to inflation factors.Apparently, the advantage of HOOPE-EnKF-PSO is its ease of implementation.HOOPE-EnKF-PSO only requires generating pseudo-observations of model parameters to be assimilated and no modification of a data assimilation code is necessary.Although HOOPE-EnKF-RTC is more complex than HOOPE-EnKF-PSO, the additional computational cost is comparable to HOOPE-EnKF-PSO in our proposed implementation.HOOPE-EnKF-RTC is better than HOOPE-EnKF-PSO in our experiment of the two-scale Lorenz96 model.This is likely because, in HOOPE-EnKF-RTC, we conditioned ensemble members of model parameters before assimilating observations, which contribute to eliminating spurious correlations between observation and model parameters.However, no notable differences are found in the more realistic experiment by the simple global atmospheric SPEEDY model.Our future work should focus on testing both methods in real-world applications, such as numerical weather prediction.
variance of these observed climatological indices was recognized as   .In the SPEEDY model experiment, we directly used the prescribed observation error in specific humidity as   .

C. Observation network in the SPEEDY model experiment
The simulated observation network in the SPEEDY model experiment can be found in Figure C1.
F[k] is a model parameter.To make Equation 28 represent the true dynamics, [] =  + [] , which implies that F[k] should be spatio-temporally distributed parameters.The objective of our experiment is to estimate [] =  + [] from sparsely distributed observations without any knowledge of S and Equations 26-27.Both the two-scale and single-scale Lorenz96 models were numerically solved by a fourth-order Runge-Kutta method with a time step of ∆ = 0.0005.Following Pathiraja and van Leeuwen (2022) and Amemiya et al. (2023), we used model time unit (MTU), of Appendix C, spanning every vertical level of the model grid.We observed wind components, temperature, specific humidity, and surface pressure with respective observation errors of 1.0 [ms -1 ], 1.0 [K], 1.0 [gkg -1 ], and 1.0 [hPa].These observation errors served as the standard deviations of the Gaussian white noise added to the nature run.This setting of observation errors is identical to Ruiz et al. (2013b).For the nature run, we ran the original SPEEDY model with   = 0.5 from January 1982 to March 1983.Observations are collected at 6-hour intervals, matching the frequency of data assimilation cycles.

Figures 3 -
Figures 3-6 show the performance of three LETKF flavors as a function of the two inflation factors.NOHOOPE has no hope to simultaneously estimate state variables and parameters accurately, even when the ensemble size is extremely large.NOHOOPE can accurately estimate state variables only when it chooses large inflation factors for state variables and small inflation factors for model parameters.This setting of inflation factors prevents the temporal variability of the parameters from being estimated, and makes them converge to a constant value.By specifying large inflation factors for state variables, NOHOOPE can resolve the model error introduced by a constant model parameter.
NOHOOPE to simultaneously estimate state variables and parameters accurately, which was difficult when the inflation factors were fixed.HOOPE-EnKF-PSO and HOOPE-EnKF-RTC outperform NOHOOPE with adaptive inflation.Note that adaptively inflated NOHOOPE with 40 ensemble members does not significantly outperform HOOPE-EnKF-PSO and HOOPE-EnKF-RTC with 10 ensemble members.The results of adaptive inflation indicate that HOOPE-EnKF-RTC is superior to HOOPE-EnKF-PSO due probably to the wider range of inflation factors that yield accurate estimation in HOOPE-EnKF-RTC (Figures 5 and 6).Figures 7 and 8 visualize the estimated [] by three LETKF flavors with adaptive inflation.HOOPE-EnKF-PSO and HOOPE-EnKF-RTC track the rapid changes in [] more stably and accurately than NOHOOPE.However, the skill of HOOPE-EnKF to estimate extremely small and large values of [] is marginal.This is a side effect of introducing climatological parameters to stabilize the filters.

Figure 9
Figure 9 shows the non-parametric posterior distribution of  ℎ ̂, which can accurately simulate the climatology of specific humidity.The dashed red line shows the fitted Gaussian distribution.Notably, the original sampled distribution exhibits significant skewness.The sensitivity of the global distribution of specific humidity to  ℎ ̂ becomes negligible when  ℎ ̂ is larger than 15 [gkg -1 ] since convection rarely occurs in this range.The pronounced skewness of the distribution arises from this nonlinear behavior of the parameter's sensitivity.As we did in section 4.1, we directly used sample mean and variance for the Gaussian distribution of the climatological parameter (  , ).
efficient and practical method to estimate time-varying parameters in high-dimensional spaces.In our newly proposed method, HOOPE-EnKF, model parameters estimated by EnKF are constrained by the results of offline batch optimization, in which the posterior distribution of model parameters is obtained by comparing simulated and observed climatological variables.HOOPE-EnKF outperforms the original LETKF in the synthetic experiments using the two-scale Lorenz96 model and the SPEEDY atmospheric model.One advantage of HOOPE-EnKF over the traditional EnKFs is that its performance is not greatly affected by inflation factors for model parameters, thus eliminating the need for extensive tuning of inflation factors.The inflation of the online estimated error of model parameters has been extensively discussed in previous works (e.g., Kotsuki et al. 2018; Ruiz et al. 2013a; Moradkhani et al. 2005b), and HOOPE-EnKF contributes to this body of published literature.Despite numerous studies related to EnKFs that address the estimation of highdimensional parameters, including real-world weather forecasting (e.g., Ruckstuhl and Janjić 2020), our work offers a unique contribution to this established literature.Pulido et al. (2016) suggested that sub-grid scale contributions of the two-scale Lorenz 96 model were initially parameterized by variables at a larger scale.Then, a functional form of this parameterization was estimated online by ETKF.While the approach byPulido et al.
errors, whilePathiraja and van Leeuwen (2022) andAmemiya et al. (2023) successfully constructed the statistical models to directly predict the error term induced by unresolved processes.Although HOOPE-EnKF accurately estimates the contribution of unresolved physics, it still assumes a persistent model for the evolution of model parameters during the forecast phase.Consequently, the prediction error may increase for lead times longer than the timescale of change in time-varying parameters.A reasonable countermeasure to this limitation is to nudge the time-varying parameters towards (  , ) during a forecast step.It is promising to analyze the spatio-temporal distribution of parameters in specific parameterizations and deepen the understanding of their uncertainty, as demonstrated by our SPEEDY model experiment.By understanding the source of uncertainty in a specific parameterization from spatio-temporally distributed model parameters, it becomes straightforward to predict and reduce the systematic model error.

Figure C1 .
Figure C1.Red squares show the location of the simulated observations in the SPEEDY model experiment.

Figure 1 .
Figure 1.Schematics of the (a) original sequential data assimilation and (b) Hybrid Online and Offline Parameter Estimation (HOOPE) concept.Black lines are timeseries of the true model parameter.Blue lines are the estimation of the online data assimilation algorithm and red arrows show the adjustments by the data assimilation steps.The gray area in (b) shows the posterior distribution of model parameters calculated by the offline batch optimization.When estimated parameters substantially deviate from the truth due to erroneousobservations and/or sampling error with an insufficient ensemble size, there is a risk of filter degeneracy.HOOPE-EnKF is expected to prevent filter degeneracy in this situation by efficiently move ensemble members to the probabilistic distribution of time-invariant parameters.See also Section 2. This figure is adopted fromSawada (2022).

Figure 2 .
Figure 2. Posterior distribution of a time-invariant model parameter from (blue) samples by the Metropolis-Hastings algorithm and (red) a fitted Gaussian distribution in the two-scale Lorenz 96 experiment.

Figure 3 .
Figure 3. RMSE of state variables in (a-c) NOHOOPE, (d-f) HOOPE-EnKF-PSO, and (g-i) HOOPE-EnKF-RTC with the ensemble size of (a,d,g) 10, (b,e,h) 20, and (c,f,i) 40.Horizontal axis shows inflation factors for state variables, while vertical axis shows inflation factors for model parameters.White tiles indicate that the model diverges.

Figure 4 .
Figure 4. Same as Figure 3 but for R of state variables.
Amemiya et al. (2023)em does not need to be observed, whereas the observation network ofAmemiya et al. (2023)included observations in all grid points.HOOPE-EnKF does not rely on the assumption of homogeneity in sub-grid scale physics, which is a requirement for Pathiraja and van Leeuwen (2022) andAmemiya et al. (2023).When a model parameter is not sensitive to any observations assimilated within a data assimilation window, HOOPE-EnKF ensures that this model parameter converges to (  , ).This property of HOOPE-EnKF may significantly contribute to the stability of the filter avoiding the negative impact of spurious correlation in numerical weather prediction, where the sensitivity of observation to model parameters varies spatio-

Table 1 .
Summary of RMSE of model parameters in three LETKF variants."Optimal" means the performance by inflation factors which minimize RMSE.

Table 2 .
Summary of R of model parameters in three LETKF variants."Optimal" means the performance by inflation factors which maximize R.