Inference of stochastic parameterizations for model error treatment using nested ensemble Kalman filters
Guillermo Scheffler Juan Ruiz Manuel Pulido







   



#1(<ref>)






plainnat








[enumerate]label*=*.
.1cm
.1cm
16.3cm
-1.5truecm












    =[2]d #1/d #2
            
            
            #1#1#1[-1mm]0mm6mm2mm0mm !2mm0mm[-1mm]0mm6mm2mm0mm ?2mm0mm
            
            
            #1ℝ
            ℋℰℳ
            
            
            
            ./../=1
            ==10000
            =10000
            
            1,2]Guillermo Scheffler1,3,4]Juan Ruiz4,5,6]Manuel Pulido[1]Centro de Investigaciones del Mar y la Atmosfera, CONICET-UBA, Buenos Aires, Argentina[2]Department of Mathematics, Universidad Nacional del Nordeste,Corrientes,Argentina[3]Department of Atmospheric and Oceanic Sciences, Universidad de Buenos Aires, Buenos Aires, Argentina[4]UMI-IFAECI, Buenos Aires, Argentina[5]Data Assimilation Research Centre, Department of Meteorology, University of Reading, UK[6]Department of Physics, Universidad Nacional del Nordeste and CONICET, Corrientes, ArgentinaInference of stochastic parameterizations for model error treatment using nested ensemble Kalman filters
                [
                December 30, 2023
            ========================================================================================================
            
            
            Stochastic parameterizations are increasingly being used to represent the uncertainty associated with model errors in ensemble forecasting and data assimilation.  One of the challenges associated with the use of these parameterizations is the optimization of the properties of the stochastic forcings within their formulation. In this work a hierarchical data assimilation approach based on two nested ensemble Kalman filters is proposed for inferring parameters associated with a stochastic parameterization. The proposed technique is based on the Rao-Blackwellization of the parameter estimation problem. The technique consists in using an ensemble of ensemble Kalman filters, each of them using a different set of stochastic parameter values. We show the ability of the technique to infer parameters related to the covariance structure of stochastic representations of model error  in the Lorenz-96 dynamical system.  The evaluation is conducted with stochastic twin experiments and imperfect model experiments with unresolved physics in the forecast model. The proposed technique performs successfully under different model error covariance structures. The technique is proposed to be applied offline as part of an a priori optimization of the data assimilation system and could in principle be extended to the estimation of other hyperparameters of a data assimilation system.
            
            
            
            
            
            
            § INTRODUCTION
            
            
            Model error treatment has  become a key ingredient for data assimilation systems. The most significant sources of the so-called model errors in numerical geophysical models are a consequence of simplifications in the representation of the dynamics, such as errors associated to the discretization of large-scale dynamics and unresolved or under-resolved physical processes represented by parameterizations. When combining observations and forecasts in data assimilation systems, a proper account of both observational and forecast model errors is crucial for a successful state estimation. In data assimilation, model errors are usually separated into two components: systematic model error and random model errors. Systematic model errors are considered as the mean model error over a sufficiently long time window, while random model errors are the departures from said mean.
             
            The conventional approach for systematic model error treatment involves the estimation of a forecast bias term, which is augmented to the state vector during the data assimilation process <cit.>. In addition to explicit bias estimation, a random error term can also be incorporated to the state variables before the data assimilation step. The probability distribution of these errors can be inferred, for example, using innovation statistics from previous assimilation cycles  <cit.>. 
            
            In the context of the ensemble Kalman filter, neglecting random model errors usually results in an underestimation of the forecast error covariances. <cit.> gives a review of operational approaches to deal with these issues. The strategies typically involve a representation of model uncertainties either as an ad-hoc inflation of the forecast errors covariance matrix, multi-model and multi-parameterization ensembles or the use of stochastic parameterizations. The role of covariance inflation to compensate for model errors has been widely studied <cit.>. The inflation can be incorporated either as a random noise added to the different states of the ensemble, known as additive inflation <cit.>, or by amplifying the ensemble members deviations from the mean state, referred to as multiplicative inflation <cit.>. The latter approach however assumes that the model errors have the same structure as the dynamically evolved internal errors <cit.>. While this hypothesis does not hold for most of the atmospheric numerical models, the multiplicative inflation approach has been rather successfully used for this purpose as well as for dealing with sampling errors due to the small ensemble sizes used operationally.  The combination of parametric model error treatments with additive and multiplicative inflation was examined by <cit.>. A significant improvement in analysis error was found when compared to using each of them separately.
            
            Explicit representation of random model errors typically require an estimate of their spatio-temporal covariance structure, which can be either constructed arbitrarily, or parameterized as a function of free parameters. Several algorithms have been proposed to infer these type of parameters.  <cit.>  used innovation statistics to estimate horizontal decorrelation length scale of model error and vertical covariances on a three-level quasi-geostrophic model. An iterative implementation of the expectation-maximization algorithm combined with an ensemble Kalman filter was successfully applied by <cit.> to infer different forms of covariance matrices of an additive Gaussian model error in a nonlinear state-space model.  On a similar model scenario, <cit.> proposed a Bayesian framework to explicitly account for the marginal posterior distribution of parameters, by  either an exhaustive grid-based exploration of the parameter space or using a particle filter approximation. This scheme can be combined with the ensemble Kalman filter to produce sequential estimations of state and parameters. Maximum likelihood estimators based on expectation-maximization and Newton-Raphson minimization were presented by <cit.>. These methods allow for the simultaneous estimation of deterministic parameters and parameters associated to stochastic processes representing model error. 
            
            The use of stochastic parameterizations represents a promising approach for model error treatment in ensemble forecasting and data assimilation. These parameterizations introduce stochastic processes directly to the model evolution, either as a random perturbation added at every time step of the model integration, or by stochastically perturbing the physical parameterizations tendencies. Stochastic parameterizations may account for a physically consistent representation of subgrid processes <cit.> and so, a state-dependent representation of model errors. The use of stochastic parameterizations impacts positively on ensemble prediction skill as shown by <cit.> and <cit.>. Stochastic parameterizations are particularly successful in triggering noise-induced transitions <cit.>. When the spatial and temporal correlation structures are properly tuned, stochastic parameterizations may potentially provide a more consistent representation of the interactions between the resolved dynamics and the subgrid parameterizations <cit.>. Within a mesoscale data assimilation system, <cit.> showed that the use of a stochastic backscatter scheme consistently outperformed the multiplicative covariance inflation scheme and the multiphysics ensemble approach. However, stochastic parameterizations require a careful tuning of the properties of the stochastic forcings in order to account for the model uncertainties. The parameters that characterize the covariance matrix of the stochastic process are referred to as stochastic parameters from now on.  Whereas estimation of deterministic parameters of the dynamical model is straightforward within the ensemble Kalman filter using state augmentation <cit.>, stochastic parameters cannot be estimated in this way. Previous studies on the use of the augmented state approach in the ensemble Kalman filter showed that the lack of correlation between the mean of the ensemble of state variables and the stochastic parameters may lead to unreliable estimations <cit.>.
            
            To overcome the lack of sensitivity of the forecast mean state to stochastic parameters, here we propose the application of a hierarchical Bayesian  framework based on two nested ensemble Kalman filters. The model state is estimated with an ensemble Kalman filter in an inner cycle, as in any conventional ensemble Kalman filter implementation. However, the filtering process is applied over different independent ensembles, each of them integrated using a different set of stochastic parameters. Hence, an ensemble of ensemble Kalman filters is assimilated in an outer filter cycle. The outer cycle is used for the estimation of the stochastic parameters. Stochastic parameters have a direct impact on the forecast error covariances, thus playing a critical role on the state analysis quality on the inner cycle. In order to increase the sensitivity of the analysis to stochastic parameters, we propose the use of a longer assimilation window in the outer cycle, composed of several internal cycles. In this way, the information in the outer filter for parameter estimation includes an ensemble of trajectories of analysed model states which is expected to be sensitive to the stochastic parameters. 
            
            Whereas few stochastic parameters can be tuned using an exhaustive exploration of the parameter space, the proposed technique can be used to simultaneously estimate several stochastic parameters. Its main advantage lays, indeed, in the ability to explicitly estimate multiple elements of the covariance matrix used in the stochastic parameterization, at a computational cost that is comparable to the state-of-the art expectation-maximization algorithms <cit.>  and schemes based on nesting sequential Monte Carlo algorithms <cit.>. It should be remarked that the proposed methodology based on nested ensemble Kalman filters is intended to be applied offline, as an optimization tool for both the stochastic dynamical model and the data assimilation system. 
            
            In Section 2 the stochastic parameter estimation based on the nested ensemble Kalman filters is introduced under a Bayesian framework. The dynamical model used in the experiments is described on Section 3. Section 4 describes the experimental setup and the different covariance matrix structures, which are evaluated in the experiments. Results from stochastic twin experiments and from imperfect model experiments are shown in Section 5. We conclude with a brief summary and discussions in Section 6.
              
            
            
            § METHODOLOGY
            
            
            To avoid the limitations of the conventional augmented state approach for ensemble based stochastic parameter estimation <cit.>, the proposed hierarchical inference technique is conducted using an ensemble of N_J  data assimilation systems, each of them using a different set of stochastic parameter values. The technique involves two nested data assimilation cycles. Firstly, the state estimation is performed independently by each data assimilation system with a set of fixed stochastic parameter values, using the ensemble Kalman filter. This step is referred to as inner cycle. Secondly, parameter estimation is performed using the resulting mean forecast states of the inner cycles as a priori states. Thus, the mean a priori states are in turn the members of an ensemble that describes the density of the state conditioned to the parameter values. The Kalman filter equations are applied to this stochastic parameter ensemble to update their values. This step will be referred as outer cycle. This procedure is repeated sequentially. 
            
            The assimilation window of the parameter estimation cycle (i.e. the outer cycle) is composed by K state estimation cycles (inner cycles). Since stochastic parameters are assumed to change slowly with time, parameters are assumed static within the outer cycle and are denoted as θ_l, which means that the stochastic parameters will be updated only every K state estimation cycles to increase the identifiability of the stochastic parameters <cit.>.
            
            At a given time,  𝐱_l, k denotes the model state at the k-th state estimation cycle during the l-th parameter estimation cycle and θ_l denotes the parameters. For simplicity, we denote any quantity at time (l,0) with a single subindex l, so for example 𝐱_l,0=𝐱_l-1,K=𝐱_l. A schematic representation of time indexes is shown in Fig. <ref>.
            
            
            
            
            
            
            
            Given an initial prior joint density of the model state variables and parameters p(𝐱_0,θ_0), the inference goal is to sequentially obtain the posterior density of the augmented state at time l given K observation sets distributed in time from (l,1) to (l,K), this is p(𝐱_l+1,θ_l+1| 𝐲_l,1:K). Model state is then estimated every time at which observations are available, while the augmented state is estimated every K state assimilation cycles. Using the definition of conditional density, the joint parameter–state inference is given by 
            
                p(𝐱_l+1,θ_l+1| 𝐲_l,1:K)=p(𝐱_l+1| θ_l+1,𝐲_l,1:K) p(θ_l+1|𝐲_l,1:K).
            
            The posterior density of the augmented state at the final time of a set of observations is given by the posterior density of the state conditioned to the value of the parameters and the observations multiplied by the posterior density of the parameters given the observations. The first term in the RHS of (<ref>) is determined through the sequential process
            
                p(𝐱_l,k| θ_l+1,𝐲_l,1:k)= 
                p(𝐲_l,k| 𝐱_l,k,θ_l+1,𝐲_l,1:k-1) p(𝐱_l,k| θ_l+1,𝐲_l,1:k-1)/p(𝐲_l,k|𝐲_l,1:k-1)
             
            from k=1 to k=K.  In practice, (<ref>) corresponds to the application of a Kalman filter for the model state given a certain set of parameters and observations at time k.
            
            The parameter posterior density, the second term in the RHS of (<ref>), is rewritten through a sequential conditioning process as
            
                p(θ_l+1|𝐲_l,1:K)∝  p(θ_l+1|𝐲_l-1,1:K) ∏^K_k=1 p(𝐲_l,k|𝐲_l,1:k-1,θ_l+1),
            
            where p(𝐲_l,1|𝐲_l,1:0,θ_l+1) ≐ p(𝐲_l,1|θ_l+1) and  p(θ_l+1 | 𝐲_l-1,1:K) is the prior density of the parameters, given the previous observations:  𝐲_l-1,1:K. Note that we leave implicit the conditioning to observations from previous outer cycles. As mentioned before, we have assumed parameters are static within an outer cycle. Alternatively, a random walk or any Markov process could also be considered, in that case p(θ_l+1 |   y_l-1,1:K ) = ∫ p(θ_l+1|θ_l) p(θ_l|𝐲_l-1,1:K) dθ_l. The parameters are assumed to be updated immediately after the k-th state assimilation cycle.
            
            
            Since the parameters are not directly observed, their likelihood to a new observation 𝐲_l,k is taken into account through the marginalization of the model state,
            
                p(𝐲_l,k|𝐲_l,k-1,θ_l+1)=∫ p(𝐲_l,k|𝐱_l,k,θ_l+1) p(𝐱_l,k|𝐲_l,1:k-1,θ_l+1) d𝐱_l,k.
            
            
            Equation (<ref>) involves a prediction density given the previous estimated state and the observation likelihood given the parameters, θ_l+1, which are then integrated along the model state. This marginalization is over the full model state and it is likely to produce errors if it is conducted directly through Monte Carlo sampling (e.g. using samples of a particle filter or an ensemble Kalman filter). Instead of using a direct sampling from the joint density, we take a Gaussian assumption for both the prediction density and the observation likelihood. Under this assumption, the sufficient statistics ℰ(𝐱_l,k| θ_l+1,𝐲_l,1:k-1), is used to assimilate the observations in (<ref>),  where ℰ(·) denotes the expectation operator. As it is known from the Rao-Blackwell theorem, the use of sufficient statistics in the estimator is expected to improve the inference of parameters. The Gaussian assumption is of course also taken for the ensemble Kalman filter that we use, so this assumption is coherent with the filter and does not imply an additional hypothesis.
            
            Following  the derivation of the sequential marginalized observation likelihood from <cit.>, the observation likelihood conditioned on  θ_l+1 is therefore
            
                p(𝐲_l,k|𝐱_l,k,θ_l+1) ∝exp[(𝐲_l,k - ℋ(𝐱_l,k))^T𝐑^-1 (𝐲_l,k - ℋ_l,k(𝐱_l,k))],
            
            where ℋ is the observation operator that transforms from model to ovservational space. Observational errors are assumed additive and Gaussian with covariance 𝐑. In coherence with the assumption taken in the ensemble Kalman filter, we now assume that the forecast density can be represented approximately by a Gaussian density, namely, 
            
            
                p(𝐱_l,k|𝐲_l,1:k-1,θ_l+1) ∝
                exp[(𝐱_l,k - 𝐱^f_l,k(θ_l+1))^T𝐏_l,k(θ_l+1)^-1 (𝐱_l,k - 𝐱^f_l,k(θ_l+1))],
            
            where  𝐱^f_l,k(θ_l+1) ≐ℰ (𝐱_l,k|𝐲_l,1:k-1,θ_l+1)≐∫𝐱_l,k p(𝐱_l,k|𝐲_l,1:k-1,θ_l+1) d 𝐱_l,k  is the mean forecast state conditioned on θ_l+1 and 𝐏_l,k(θ_l+1) is the forecast error covariance given θ_l+1. In what follows, The dependencies on θ_l+1 are dropped to reduce notation and defining 𝐇_l,k as the linearized observation operator at 𝐱^f_l,k, the resulting approximated observation likelihood in the time interval (l,1) to (l,K) is
            
                ∏_k=1^K p(𝐲_l,k|𝐲_l,k-1,θ_l+1)∝
                ∏_k=1^K exp[(𝐲_l,k - ℋ_l,k(𝐱^f_l,k))^T (𝐇_l,k𝐏_l,k𝐇_l,k^T+𝐑 )^-1 (𝐲_l,k - ℋ_l,k(𝐱^f_l,k))],
            
            
            which is equivalent to the approximated observation likelihood obtained in <cit.>.  That work also constrains the values of the statistical parameters within a time interval where K observation sets are available. However,  a point estimation of the parameters is conducted there through maximization of the logarithm of the approximated observation likelihood. On the other hand, here we assume a Bayesian framework, see  (<ref>), in which an inference of the density of the parameters conditioned to the set of observations is obtained given some prior knowledge of the parameters. Our procedure resembles the Rao-Blacwellized particle filter <cit.>, where N_J ensemble Kalman filters are conducted in order to marginalize the parameters posterior distribution.
            
            
            In this work, the parameters to be estimated are assumed to be associated with an additive Gaussian model error. We also assume that parameters follow a Gaussian distribution. While the latter hypothesis is not warranted, it allows us to treat the state-parameter estimation problem by using two nested ensemble Kalman filters.
            
            
            Let us consider an ensemble of initial parameters θ_0^(j), with j=1,⋯,N_J, sampled from  p(θ_0). Each parameter θ_0^(j) is associated with an ensemble of N_I model states {𝐱_0^(j,i), i=1,⋯,N_I}. Therefore, a set of N_J ensembles is initialized and each of them represents different values of the parameters θ. The updates to the  ensemble state members are determined in the inner cycles with ensemble Kalman filters. The filters should be identical to the data assimilation system for which parameters are being estimated, e.g. same physical parameterizations, number of ensemble members. Note that the the N_J ensembles are assumed to evolve independently, hence the state update neglects any correlation between the ensembles. The mean state of the j-th ensemble is given by
            
                𝐱_l,k^a(j)=𝐱_l,k^f(j) + 𝐏^(j)_l,k𝐇^T (𝐇_l,k𝐏^(j)_l,k𝐇_l,k^T + 𝐑)^-1  [𝐲_l,k-  ℋ_l,k(𝐱_l,k^f(j))],
            
            where 𝐱^a(j), 𝐱^f(j) and 𝐏^(j) denote the analysed state, forecasted state and the forecast covariance matrix of the j-th ensemble respectively.
            
            
            Next, the parameter posterior density conditioned to the observations, p(θ_l+1|𝐲_l,1:K), is inferred using (<ref>). In this sense, Equation (<ref>) can be interpreted  as the serial assimilation of observations along the l-th state assimilation window. Under the already taken assumptions, the likelihood function and the forecast density are assumed Gaussian in  (<ref>). Therefore, in the “outer” cycle we apply the ensemble Kalman filter to infer p(θ_l+1|𝐲_l,1:K). Because the initial condition of the hidden state is not known with complete certainty, the observational error for the inference of the parameter posterior density increases with the forecast error covariance matrix of the hidden state (see Eq. <ref>) . Thus, an increased“observational” error covariance matrix for the assimilation of parameters  is obtained in (<ref>), 𝐇𝐏_l,k(θ_l) 𝐇^T+𝐑. 
            
            
            The application of the Kalman filter to the ensemble representing the parameter distribution (<ref>) with (<ref>) results in the analyzed mean parameters given by
            
                θ_l+1^a = θ^f_l+1 + 
                ∑_k=1^K 𝐏^θ𝐱_l,k𝐇^T[𝐇  (𝐏_l,k + 𝐏^𝐱𝐱_l,k) 𝐇^T + 𝐑]^-1 [𝐲_l,k-  ℋ(𝐱^f_l,k)],
            
            where 𝐱^f_l,k is the average of ensemble forecast mean states at time (l,k) over the N_J ensembles. The covariance matrices 𝐏_l, 𝐏^𝐱𝐱_l,k and 𝐏^θ𝐱_l  are defined empirically as follows
            
            
            
                𝐏_l,k= 1N_J∑_j=1^N_J𝐏^(j)_l,k
            
             is the sample forecast state covariance averaged among the N_J ensembles,
            
            
                𝐏^𝐱𝐱_l,k = 1/N_J-1∑_j=1^N_J (𝐱^f(j)_l,k-𝐱^f_l,k) (𝐱^f(j)_l,k-𝐱^f_l,k)^T
            
            is the sample covariance of the N_J ensembles, and similarly,
            
            
                𝐏^θ𝐱_l,k = 1/N_J-1∑_j^N_J  (θ^f(j)_l+1-θ^f_l+1) (𝐱^f(j)_l,k-𝐱^f_l,k)^T
            
            is the parameter-state covariance matrix.
            
            Note that the forecast state error covariance for the parameter estimation in the outer cycle is  the sum of the mean covariance of the N_J ensembles and the covariance of the outer ensemble. Equation (<ref>) shows that parameters are estimated using ensemble mean states of the inner cycle as individual state members in the outer cycle. This formulation defines a parameter-state covariance matrix that is able to transfer model state innovations to statistical parameters. This is a key difference between the nested ensemble Kalman filters and the  standard state augmentation approach for parameter estimation, since the impact of stochastic parameters is accounted from an ensemble mean, and not over individual members.
            
            
            
            
             §.§ Implementation details
            
            
            The implementation of the inner cycles and the outer cycles are both based on the Ensemble Transform Kalman Filter <cit.>. For the outer cycles, the parameter update equation (<ref>) is implemented as an asynchronous ETKF <cit.>.  For this purpose, an aggregated vector is constructed by column-wise concatenating observations 𝐲_l,1:k in a single observation vector 𝐲^∗_l. A similar concatenation is performed with the ensemble members in the state and in the observational spaces. The aggregated observational error covariance matrices  𝐑^∗_l and mean covariance 𝐏^∗_l are constructed with the k-th diagonal block 𝐑_l,k and 𝐏_l,k respectively.  
            The main steps of the nested ensemble transform Kalman filters in this work are implemented as follows:
            
                  
            
              *  Given N_J parameters θ^(1:N_J)_0 and N_J independent ensembles of N_I state members, 𝐱^(1:N_J,1:N_I)_0,  and (L × K) observations y_1:L,1:K.	
            	
              *  State estimation: For each assimilation cycle l, (with l=1,⋯,L) do:
            	
              *  For each inner assimilation cycle k, (with k=1,⋯,K) do:
            		
              *  Calculate the ensembles of analysed states 𝐱_l,k^a(j,i)  performing N_J-EnKFs independently.
            	 		
              *  Store  ℋ_l,k(𝐱̅_l,k^f(j)) for each ensemble and the average of the forecast error covariance matrix in the observational space over the N_J ensembles 𝐇𝐏_l,k𝐇^T.
            		
              *  Parameter estimation: Obtain the ensemble of estimated parameters:
            	
              *  Concatenate the K mean predicted observations ℋ_l,k(  𝐱̅_l,1:K^f(j)) to construct an (n_x× K)-dimensional ensemble of N_J members.
            		
              *  Construct the agreggated observation vector 𝐲^∗_l=[𝐲_l,1,⋯,𝐲_l,K]^T and the tangent linear observation operator 	𝐇_l^∗=[𝐇_l,1, 𝐇_l,2, ⋯,𝐇_l,K]^T
              *  Construct the block diagonal extended observational error covariance matrix 𝐑_l^∗, whose k-th diagonal block is 𝐑^∗_l,k=𝐇_l,k^T𝐏̅_l,k𝐇_l,k+𝐑.
            		
              *  Obtain the updated parameter ensemble mean and perturbation ETKF using the aggregated matrices calculated in steps <ref>-<ref>
            
            Since the experiments were conducted in a low-dimensional dynamical system (see Section <ref>), the use of covariance localization is not explored in this work. Localization may become mandatory for systems in which the number of state-space dimensions exceeds the number of ensemble members, as occurs in numerical forecast models. Additionally, we have assumed that sampling errors of the filter can be partially accounted by the stochastic parameterization so that multiplicative covariance inflation of state variables is not included in most of the experiments. As has been described in <cit.> and <cit.> assuming a persistence model for the parameters can result in the collapse of the parameter ensemble spread and divergence of the parameter estimation. However, in our experiments, such mechanisms were not needed to increase parameter ensemble spread. 
            
            
            
            
            § DESCRIPTION OF THE EXPERIMENTS
            
            
             §.§ The Lorenz-96 dynamical model
            
            The two-scale Lorenz-96 dynamical model has been extensively used as a testbed model for the development of data assimilation schemes due to its reduced computational cost and its ability to mimic specific properties of the atmospheric predictability <cit.>. It represents the dynamics of a cyclical set of large-scale variables over a circle of latitude, each coupled to a set of high-frequency small-scale variables.  Each model equation contains terms that represent non-linear advection, dissipation and external forcings. The small-scale variables are coupled to the large-scale variables through an additive forcing term. 
            
            The set of equations of the two-scale Lorenz-96 dynamical model is given by large-scale variable equations,
            
            
                dx_ndt=-x_n-1(x_n-2-x_n+1)-x_n+F-hcb∑_m=M(n-1)+1^Mn y_m,
            
            and small scale variable equations,
            
                dy_mdt=-cby_m+1(x_m+2-x_m-1)-cy_m+hcb x_1+int[m-1/M],
            
            where n=1,⋯,N and m=1,⋯,MN. 
            
            Both sets of variables have cyclic boundaries conditions: x_n+N=x_n and y_m+MN=y_m.  In this work, the coupling and scale parameters are set to the standard values of h=1, b=10 and c=10 as in <cit.> and <cit.>. The number of large-scale variables was set to N=8, each coupled to M=32 small-scale variables y, so that MN=256. To achieve a chaotic solution, the external forcing is set to F=20 for all the experiments. 
            
            In the imperfect model experiments, the small-scale variables can be interpreted as unknown physical processes which cannot be explicitly resolved in numerical models, so that only the dynamics of large-scale variables are represented by the model. The effect of the small-scale variables is introduced as a parametrization that is a function of the resolved large-scale variables only. This mimics in a very simple way model errors associated with the parametrization of unresolved processes in realistic atmospheric or oceanic numerical models.  The Lorenz-96 system results particularly suitable for proof-of-concept experiments involving subgrid model error representation and parameterizations <cit.>. The truncated version of the model can be expressed as
            
            
                dx_ndt=-x_n-1(x_n-2-x_n+1)-x_n-U(x_n),
            
            where U(x_n) represents the parameterization of small-scale processes. The forcing term F is also assumed to be part of the parameterization U.  In this work, the parameterization is of the form
            
            
                U(x_n)=a_0+a_1 x_n+e_n(t).
            
            The first two terms represent a deterministic forcing that is a function of only the resolved variable x_n. The coefficients a_0 and a_1 can be estimated via a least-squares fitting using an integration of the complete system (Eqs. <ref>-<ref>) as in <cit.>, or inferred via data assimilation using only noisy observations of the resolved variables of the full system <cit.>.
            
            The processes that cannot be accounted by a deterministic function of the state variables, are included as a state-independent red-noise stochastic forcing, discretized as the realization of a zero mean first-order autoregressive process (AR(1)),
            
                𝐞(t)=ϕ 𝐞(t-Δ t) + (1-ϕ^2)^1/2 η.
            
            
            The coefficient ϕ represents the lag-1 autocorrelation of 𝐞(t), and Δ t is the model integration timestep. The vector η∼𝒩(0,Σ) represents a random draw from a zero-mean Gaussian distribution with covariance Σ. The adequacy of a stochastic parametrization for model error representation in the Lorenz-96 model was proved by <cit.>, <cit.> and by <cit.>.
            
            
            
            
            
             §.§ Experimental setup 
            
            
            We first evaluate the nested ensemble Kalman filters using twin experiments. In these experiments, the “true" integration  consists of an integration of the truncated model (<ref>)-(<ref>) with a stochastic forcing generated using a prescribed covariance structure. The same model and covariance structure are then used as forecast model, but with uncertain parameters. Since stochastic processes are present in both “true" and forecast models, it is not possible to replicate the true integration using the forecast model, even when using identical initial conditions and parameter configuration. Hence, these experiments are referred to as “stochastic twin experiments".  
            
            
            The truncated Lorenz-96 dynamical system was integrated with a fourth-order Runge-Kutta  scheme. The system was initialized after a spinup of 1460 dimensionless model time units, which is roughly equivalent to 20 years of atmospheric evolution.  The nature run was generated integrating the “true model" for 250 model time units (i.e. 50000Δ t) with a timestep of Δ t=0.005. The coefficients a_0 and a_1 were set to a_0=19.169 and a_1=-0.813. These values were estimated via least square fitting using an integration of the two-scale Lorenz-96 system with F=20. The structure and specific values of stochastic parameters used in the experiments are specified below. Synthetic observations are then generated by perturbing the nature run with zero-mean Gaussian uncorrelated noise of variance 𝐑=σ_R^2 𝐈 and σ_R^2=1, where 𝐈 is the identity matrix.  All the variables are observed simultaneously, with a frequency of δ t=10Δ t=0.05. 
            
            In the experiments, ensembles of N_I=30 members are used for the inner cycles, whereas for the outer cycle,   N_J=15 independent ensembles are considered. Initial conditions for the states of the ensembles were randomly chosen between uncorrelated states from the true model integration.  The number of inner cycles within each outer cycle is set to K=5.  This value was chosen to balance the parameter convergence speed, precision and computational cost, through preliminary sensitivity experiments.
            
            
            Imperfect model experiments are conducted using the two-scale Lorenz-96 model as the true state evolution. In this case, the model was integrated during 250 model time units, with a time step of Δ t =0.001.  Observations are generated using the same observational error and operator as in the stochastic twin experiments.  Note that only the large-scale variables are observed in these experiments. The truncated Lorenz-96 model is used as forecast model, with the same integration scheme as in the stochastic twin experiments. 
            
            To evaluate the sensitivity of the estimations to observational sampling errors, each of the proposed assimilation experiments is repeated 10 times, with different realizations of observational error and stochastic forcing, and changing the ensemble of initial states and parameters.  The computation of verification scores excludes the first 200 state assimilation cycles to avoid the effect of the filter spinup.
            
            
            
            
             §.§ Model error covariance structure
            
            
             An explicit estimation of the stochastic processes covariance matrix might result intractable for geophysical models. In practice, several assumptions and simplifications can be considered in  an attempt to replicate the structure of said covariances. In this work, the nested ensemble Kalman filters are used to estimate parameters related to different structures of the covariance matrix Σ.  The proposed parameterizations of the covariance matrix represent different hypothesis of the behavior of the model error which are usually assumed in practice (i.e. Gaussian errors, spatially symmetric covariances, isotropy). In this work, parameters for the following covariance matrix structures Σ are estimated:
            
            
            
            
              IIsotropic non-correlated: In this case the covariance matrix is expressed as Σ= σ^2 𝐈 where 𝐈 is the identity matrix. This model assumes that the model error variance is the same for all the resolved variables and that model errors for different model variables are uncorrelated. In this case, the standard deviation σ is the only parameter to be estimated.
            
            
              IIIsotropic exponential covariance:  The covariance is parameterized as Σ_i,j=σ^2 e^-ρ d_i,j, where d_i,j is the minimum distance between variables x_i and x_j, indicating an exponential spatial decrease of covariances. This approach assumes again that the variance of the model error is the same for all the variables, but it incorporates an a priori spatial covariance structure.  Smaller values of ρ are associated with longer model error spatial correlations.  In this case, the parameters to estimate are the standard deviation σ and the spatial scale parameter ρ.
            
            
              IIIHorizontally symmetric homogeneous covariance matrix  Σ: All the variables are assumed to have the same spatial covariance structure and that covariances are horizontally symmetric (namely Σ_n,n-i=Σ_n,n+i). The stochastic parameters to estimate in this case are the variance σ^2 and the model error neighbouring covariances (in our model only 5 parameters). 
            
            
              IVNon-isotropic non-correlated: This covariance structure ignores spatial correlations but assumes that the stochastic forcing associated to each variable has a different standard deviation. In this case  the covariance matrix is represented as Σ=diag(σ_1^2,σ_2^2,⋯,σ_N^2) representing a spatially heterogeneous model error distribution and the parameters {σ_1^2,σ_2^2,⋯,σ_N^2} are estimated independently.
            
            
            The temporal autocorrelation parameter in the experiments was fixed to  ϕ=0.984, as in <cit.>, representing a persistent stochastic forcing. Although the autocorrelation parameter could be included for estimation, it has been found that the optimal solution is not unique. There is a wide range of optimal combinations, in terms of root mean squared error and ignorance skill scores, between the stochastic forcing amplitude and the autocorrelation time parameter <cit.>.
            
            
            
            
            
            
            § RESULTS FROM STOCHASTIC TWIN  EXPERIMENTS
            
            
             §.§ Isotropic non-correlated stochastic noise
            
            
            The nature run for the first experiment uses the isotropic non-correlated covariance structure (case I) with σ_nat=2. This configuration of the stochastic forcing mimics the 2-scales Lorenz-96 model <cit.>. The only parameter to estimate is the standard deviation σ. The initial values for parameter σ were drawn from a 𝒩(1.5,0.5^2) distribution.
            
            Estimation results for parameter σ are shown in Fig. <ref>a. In all the experiments, the parameter values converge rapidly during the first 10 (50) outer (inner) assimilation cycles. Independently of the values used in the parameter prior distribution at time l=0, estimated parameters converge toward a narrow range of values after approximately 250 parameter assimilation cycles  (1250 state assimilation cycles). Though different experiments do not converge to an identical parameter value, final estimations in the different experiments have a small relative standard deviation of less than 1.2%.  Figure <ref>b shows the parameter ensemble evolution during the first 100 parameter assimilation cycles for one of the experiments shown in Fig. <ref>a. After the short assimilation spinup period, parameter updates are small and the parameter ensemble spread remains relatively stable throughout the duration of these experiments.
            
            
            
            
            
            The mean estimation of σ averaged among different experiments is σ^a=2.2, which is slightly larger than the parameter used in the nature integration σ_nat=2. To analyze the validity of the inferred parameter, we conducted an exhaustive sampling of the parameter space. For this purpose, we performed data assimilation experiments in which only the state variables were assimilated and different fixed values of the parameter σ were used during the entire assimilation experiment. Values of σ were evenly distributed, covering the range [1.25,3.25] (with Δ_σ =0.05). Each data assimilation experiment consisted of 2300 assimilation cycles, excluding an initial spinup of 200 assimilation cycles. To avoid sampling issues, the experiments were repeated 25 times for each  parameter value using different observational errors and ensemble of initial conditions for the first assimilation cycle. The mean over space, time and different experiment realizations of the analyzed state RMSE is shown in Fig. <ref>. A clear global minimum is found in the experiments. The cost function has an overall convex geometry, with increased sensitivity towards smaller values of σ. The optimal standard deviation parameter found through exhaustive parameter evaluation was σ_ex=2.15, which is also larger than σ_nat.  The  discrepancy found between σ_ex and σ_nat is expected to be a consequence of the usage of a finite ensemble size without using multiplicative or additive covariance inflation in the data assimilation process.  Hence, the larger stochastic noise amplitude attempts to correct sampling errors due to the finite ensemble size. As shown in Fig. <ref>, the values of σ estimated with the nested ensemble Kalman filters is rather coherent with this cost function.
            
            
            
            
            
             §.§ Parameterized spatial correlations
            
            
            In these experiments, stochastic parameterizations of the true model and the forecast model  use the isotropic double exponential covariance structure (Case II) with a decaying function.
            
            The parameters to estimate in these experiments are σ and ρ. In the nature integration, the standard deviation of the process is set to σ_nat=2,  and the decorrelation scale parameter is set to  ρ_nat=0.3. The latter leads to a moderate decaying rate, i.e. the covariance between the most distant variables is Q_i,i+4≈ 0.3σ^2. Initial values for parameter σ are again drawn from a 𝒩(1.5,0.5^2) distribution, while a 𝒩(0.5,0.15^2) distribution was used for parameter ρ.
            
            Results from simultaneous estimations of σ and ρ are shown in Fig. <ref>. Assuming complete ignorance of the parameter values used in the nature integration, on average, the estimations converge to parameter values σ^a=2.12 and ρ^a=0.29. The estimated values for σ are on average at least 5% larger than the value used in the nature integration. It is worth reminding that no inflation is being added to the state ensemble, so the variance overestimation may also be associated to the requirement of additional covariance inflation to alleviate the effect of sampling errors.  
            
            The optimal parameter combination obtained through exhaustive search that minimizes RMSE is σ_ex=2.15 and ρ_ex=0.36. Therefore, the optimal values estimated with the nested filters are slightly biased towards lower values of both parameters.  Smaller values of ρ are associated with larger covariances between distant variables, hence the experiment with exhaustive evaluation of the parameter space suggests weaker spatial correlations and an inflated variance. The RMSE cost function is again convex and asymmetric, especially for ρ (Fig. <ref>). However, all the estimates lay close to the set of parameter values that produce the minimum RMSE (blue dots in Fig. <ref>). The RMSE associated with these estimations is at most 1% larger than the minimum RMSE.
             
            
            The ETKF implementation of the nested ensemble Kalman filters requires the inversion of matrix 𝐑^∗=𝐇𝐏̅𝐇^T + 𝐑 in every outer cycle. A significant reduction of the computational cost is obtained if  𝐏̅_l,k can be assumed diagonal for R^∗ computation. Results for different repetitions of the experiment considering a diagonal 𝐏̅_l,k are marked with stars in Fig. <ref>. Slight differences in the estimated parameters are found when this assumption is considered. On average, the parameter σ is approximately 1% larger than when using off-diagonal elements of 𝐏̅_l,k, while the difference in ρ is almost negligible (i.e. about 0.5% smaller). The effect on the analyzed state RMSE is rather small (<0.01%). The practical tweak of assuming a diagonal matrix 𝐏̅_l,k does not degrade significantly the quality of estimation, while reducing its computational cost. Note that if a non-square root ensemble Kalman filter was used in the outer cycle, there would be no computational benefits in assuming  𝐏̅_l,k to be diagonal, since such schemes would require computation of (𝐇𝐏̃^𝐱𝐱𝐇+𝐇𝐏̅𝐇 + R)^-1.
            
            
            
            
            
            
            
            
            
            
            
            
             §.§ Non-isotropic variance estimation
            
            
            In state-of-the-art geophysical models, model error is usually non isotropic since each physical variable at each location might be affected differently by model errors. It is interesting to study if the proposed technique can retrieve the structure of Σ when removing the isotropic assumption (i.e. using covariance model III). With this purpose, an idealized experiment was conducted in which the stochastic parameterization is driven by an uncorrelated zero-mean Gaussian process, with variances σ^2_1=σ^2 _4=2.5^2 and the rest of the variances set to 1.5^2. The number of parameters to be estimated in this case is 8. The initial values of parameters were sampled from  𝒩(2;0.5^2).   
            
            
            
             
            
            Figures <ref>a,b show the estimation of the parameters σ^2_1 and σ^2_2 as a function of time for independent experiment repetitions.  In most experiments, estimated parameters are larger than the the values used in the nature integration, with a net difference of up to +0.5.  Similar results were found for the remaining parameters (Fig. <ref>c). The overestimation is expected to compensate for the limited ensemble size in the state ensemble. In spite of the noise, estimations on most experiments converged to the same parameter range after a spin-up period of around 300 outer cycles.
            
            Under this experimental setting, the number of parameters to estimate is equivalent to the number of state variables. This is the experiment with the largest number of estimated parameters, therefore it is important to evaluate the impact of the parameter ensemble size  upon the quality of the estimations. Table <ref> shows the temporal standard deviation of the estimated parameters over the last 300 parameter assimilation cycles for experiments using different number of ensemble members in the outer cycle. Values were averaged over 25 experiments with different observational error realizations and initial parameter ensembles. The variability of the estimated parameters has some dependence with the number of ensembles used in the outer cycle. In particular, an improved convergence is found when the number of ensembles is increased.  In the case with N_J=5 ensembles, as expected the parameter estimations show large variations among experiments. This suggests that for larger dimensional systems, the use of covariance localization in the parameter space may become mandatory.  Note, however, that the RMSE of these experiments in average  does not improve significantly when using more than N_J=15 ensembles for a parameter space with 8 degrees of freedom.
            
            
            
            § RESULTS FROM IMPERFECT MODEL EXPERIMENTS
            
            
            In these experiments, the two-scale Lorenz-96 model is used to generate the nature run. In the first set of experiments, the stochastic forcing used in the truncated model has covariance structure I, then, the only parameter to estimate is σ. Figure <ref>a shows estimation results for independent experiments with different realizations of observational errors, stochastic forcing and initial conditions. In most cases, convergence is achieved during the first 300 parameter assimilation cycles. The mean estimated value is σ^a=1.92. These estimations are compared with exhaustive evaluation of the parameter space (Fig. <ref>b). The  minimum value found  through exhaustive search corresponds to σ_ex=1.95, which is very close to the mean estimated value and also within the range of the estimated parameters in the different realizations of the estimation experiments.
            
            Unlike the experiments in the previous section, the optimal structure of the covariance matrix Σ for the imperfect model scenario is not known and the covariance structure I may result in a suboptimal representation of the model error. Thus, we inferred empirically, and offline, the characteristics of the covariance matrix that best fits the truncated model to the two-scale Lorenz-96 system. For these diagnostics, the true state evolution is assumed known including the evolution of the small-scale variables, contrary to the data assimilation experiments in which we assume that we only know a set of noisy observations of the large-scale variables. A large integration of 10000 model time units of the two-scale Lorenz-96 system was conducted. Using the least-squares  deterministic parameters a_0 and a_1,  the covariance of the residuals are given by
            
                r(x_n,t)=[U_det(𝐱_n,t) - ℱ(𝐱_n,t)],
            
            where U_det  is the forcing estimated by the deterministic parametrization (first two terms in (<ref>)) and ℱ is the forcing obtained in the two-scale Lorenz-96 system (last two terms in Eq. <ref>)
            
                ℱ(x_n,t)=F-hcb∑_m=M(n-1)+1^Mn y_m
            
            
            The covariance of the residuals r(𝐱,â_0,â_1,t) can be seen as an approximation of the model error covariance matrix Σ of the truncated Lorenz-96 model when using only the deterministic part of the parameterization.  Figure <ref>a shows the covariance of the residuals. Model errors have a variance of σ^∗ 2=4.93± 0.07, while the covariances between neighboring variables are σ^∗_i,i±1≈ -0.55 and σ^∗_i,i± 2. ≈ 0.7. Similar model error covariance structures were found for other configurations of the two-scale Lorenz-96 system <cit.>. Inferring this type of intricate covariance structure is not straightforward.
            
            
            
            These results are compared with parameters estimations calculated via exhaustive sampling of the parameter space. For this case, 10 independent experiments were performed, using different observational errors. Since the computational cost grows exponentially in this methodology, the parameter space was explored with a spatial grid of Δ_σ=0.125 and a 5-dimensional guess given by the nested ensemble Kalman filters. While both cases have a similar variance (σ^∗ 2≈ 4.93 and σ^2_ex≈ 4.25), the optimal stochastic forcing covariance between distant variables is significantly larger than the ones estimated offline (Fig. <ref>b). However the offline estimated parameters are not expected to be optimal in an RMSE sense for a data assimilation system <cit.>. The covariances estimated with the nested ensemble Kalman filters are expected to be similar to the ones found through the costly exhaustive exploration.
            
            We evaluate the potential of the nested ensemble Kalman filters to uncover the covariance structure using covariance model III. This structure is flexible enough to represent the complex covariance associated with model error in the truncated Lorenz-96 equations. Results of 10 repetitions of the experiment  are shown in Fig. <ref>. The estimations are less precise than in the previous experiments and require around 400 parameter assimilation cycles to converge. The mean parameter values obtained with the nested ensemble Kalman filters are in general consistent with the offline estimations shown in Fig. <ref>a, but with pronounced differences in the magnitudes of the off-diagonal elements. However, estimated values are close to the parameter values that effectively minimize the analysis RMSE (dashed lines). The nested ensemble Kalman filters are able to accurately estimate the variance  σ^2_i and the first covariance σ_i,i± 1. It is also able to recover the sign of the upcoming covariances. Further experiments are needed to assess the possibility of estimating more distant covariances.  For estimating distant covariances, an increase of the ensemble size is of paramount importance to diminish the impact of spurius distant correlations due to undersampling.
            
            
            
            
            
            § DISCUSSION
            
            
            In this work we introduce a novel data assimilation technique to infer stochastic parameters that uses on a nested implementation of two ensemble Kalman filters, based on a hierarchical Bayesian framework. The estimation of stochastic parameters requires an ensemble of data assimilation systems that is identical to the system whose parameters are being estimated i.e. model configuration, resolution, number of ensemble members. In this way, the technique not only can be used to infer parameters for stochastic parameterizations,  but it can also be used to estimate other hyperparameters associated to the data assimilation system. Whereas, we implemented an ensemble transform Kalman filter <cit.> for both state and stochastic parameter estimation, the nested ensemble Kalman filters can be extended to different flavours of the ensemble Kalman filters and the use of other data assimilation schemes in the inner cycles such as a hybrid-variational one. For the outer cycle, the use of particle filters is also possible, especially if the distribution of the parameters is strongly non-Gaussian, resembling the Rao-Blackwellized particle filter <cit.>.
            
            The proposed technique is intended to be used offline, as a tool for optimizing the data assimilation system. While the computational cost of the proposed technique is relatively large, it proves to be remarkably more economic than manually tuning model parameters, particularly when estimating more than 4-5 parameters. The computational cost might be comparable with other state-of-the-art schemes like the expectation-maximization algorithm <cit.> and less expensive than SMC^2<cit.>.
            
            We proved that the technique is able to successfully estimate parameters on stochastic twin experiments with simple model error covariance structures like the double exponential function or the diagonal isotropic case.  The estimated parameters are close to the optimal parameter values found through an exhaustive exploration of the parameter space at a significantly lower computational cost. The technique is also robust for the simultaneous estimation of multiple stochastic parameters. Additionally, and more importantly, the nested ensemble Kalman filters were able to recover the structure of model error covariances in an experiment with missing subgrid dynamics without making any a priori assumptions on the covariance structure nor the missing physics. The experiments were performed on a low-dimensional chaotic model. Further experiments in high-dimensional systems are required for which  covariance localization both in the state and in the parameter space becomes necessary. Further research is also required to evaluate the potential of reconstructing model error covariances between variables of different types and scales.
            
            The stochastic parameterizations used in this work couple the model error representation to the model dynamics by incorporating the stochastic forcing directly on the model equations. This is an important difference with respect to other model error treatment schemes that incorporate background state perturbations in the instant prior to the assimilation. Additionally, in this work, the stochastic forcing is assumed to be state independent. In the Lorenz-96 model, the amplitude of stochastic perturbations may be partially controlled by its associated state variable <cit.>. The proposed scheme is expected to handle state-dependent stochastic parameterizations, as well as other like, parameterizations of stochastically perturbed tendencies <cit.>.
            
            
            The possibility of estimating other types of hyperparameters in the context of the ensemble Kalman filter is not explored in this work but it appears as a promising venue. The nested ensemble Kalman filters, in principle, could be applied to the estimation of parameters related to the observational error covariance matrix, as well as covariance localization length-scales for state estimation. Hybrid schemes, like the ensemble 4DVar <cit.>, could also benefit from the proposed technique for inferring the optimal covariances weighting coefficients.
            
            
            
            99
            
            [Aksoy et al.(2006)]aksoy06
            Aksoy A, Zhang F, Nielsen-Gammon J. 2006. Ensemble-Based Simultaneous State and Parameter Estimation in a Two-Dimensional Sea-Breeze Model. Mon. Weather Rev, 134, 2951–2970.
            
            [Annan et al.(2005)]annan05
            Annan JD, Hargreaves JC, Edwards NR, Marsh R. 2005. Parameter estimation in an intermediate complexity earth system model using an ensemble Kalman filter. Ocean Model., 8, 135–154.[Anderson and Anderson(1999)]anderson99
            Anderson JL, Anderson SL. 1999. A Monte Carlo Implementation of the Nonlinear Filtering Problem to Produce Ensemble Assimilations and Forecasts. Mon. Weather Rev., 127, 2741–2758.
            
            [Anderson(2001)]anderson01
            Anderson J. 2001. An Ensemble Adjustment Kalman Filter for Data Assimilation. Mon. Wea. Rev., 129, 2884–2903.
            
            [Anderson(2009)]anderson09
            Anderson JL. 2009. Spatially and temporally varying adaptive covariance inflation for ensemble filters. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 61, 72–83. [Arnold et al.(2013)]arnold13
            Arnold HM, Moroz IM, Palmer TN. 2013. Stochastic parametrizations and model uncertainty in the Lorenz 96 system.  Phil Trans R Soc A, 371, 20110479. [Baek et. al(2006)]baek06
            Baek SJ, Hunt BR, Kalnay E, Ott E, Szunyogh I. 2006. Local ensemble Kalman filtering in the presence of model bias. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 58, 293–306.[Birner and Williams(2008)]birner08
            Birner T, Williams PD. 2008. Sudden Stratospheric Warmings as Noise-Induced Transitions. J. Atmos. Sci., 65, 3337-3343[Buizza et al.(1999)]buizza99
            Buizza R, Milleer M, Palmer TN. 1999. Stochastic representation of model uncertainties in the ECMWF ensemble prediction system. Q. J. Royal Meteorol. Soc, 125(560), 2887–2908.[Carrasi et al.(2017)]carrassi17
            Carrassi A, Bocquet M, Hannart A, Ghil M. 2017. Estimating model evidence using data assimilation. ‎Q. J. Royal Meteorol. Soc, 143(703), 866–880.
            
            [Chopin et al.(2013)]chopin13
            Chopin N, Jacob PE, Papaspiliopoulos O. 2013. SMC2: an efficient algorithm for sequential analysis of state space models. J. R. Stat. Soc. Series B Stat. Methodol. , 75(3), 397–426.
            
            [Crommelin and Vanden-Eijnden(2008)]crommelin08
            Crommelin D, Vanden-Eijnden E. 2008. Subgrid-Scale Parameterization with Conditional Markov Chains. J. Atmos. Sci., 65, 2661–2675.[Christensen et al.(2015)]christensen15
            Christensen HM, Moroz IM, Palmer TN. 2015. Stochastic and perturbed parameter representations of model uncertainty in convection parameterization. J. Atmos. Sci., 72(6), 2525–2544.
            
            [Danforth et al.(2007)]danforth07
            Danforth CM, Kalnay E, Miyoshi T. 2007. Estimating and correcting global weather model error. Bull. Am. Meteorol. Soc., 88, 303–304. [Dee(1995)]dee95
            Dee DP. 1995. On-line estimation of error covariance parameters for atmospheric data assimilation. Mon. Weather Rev., 123, 1128–1145.
            
            [Dee and Da Silva(1998)]dee98
            Dee DP, Da Silva AM. 1998. Data assimilation in the presence of forecast bias. ‎Q. J. Royal Meteorol. Soc, 124(545), 269-295. 
            
            [DelSole and Yang(2010)]delsole10
            DelSole T, Yang X. 2010. State and parameter estimation in stochastic dynamical models. Phys. D, 239, 1781–1788. 
            [Doucet et al.(2000)]doucet00
            Doucet A, Godsill S, Andrieu C. 2000. On sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and computing,  10(3), 197–208.
            
            [Dreano et al.(2017)]dreano17
            Dreano, D., P. Tandeo, M. Pulido, T. Chonavel, and I. Hoteit, 2017: Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximisation algorithm. Q. J. R. Meteorol. Soc., 1877–1885.[Frei and Künsch(2012)]frei12
            Frei M,  Künsch HR. 2012. Sequential State and Observation Noise Covariance Estimation Using Combined Ensemble Kalman and Particle Filters. Mon. Weather Rev., 140(5), 1476–1495.[Griffith and Nichols(2000)]griffith00
            Griffith AK, Nichols NK. 2000. Adjoint methods in data assimilation for estimating model error. Flow, turbulence and combustion, 65(3), 469–488.
            
            [Grooms et al.(2015)]grooms15
            Grooms I, Lee Y, Majda AJ. 2015. Ensemble Filtering and Low-Resolution Model Error: Covariance Inflation, Stochastic Parameterization, and Model Numerics. Mon. Weather Rev., 143(10), 3912–3924.[Ha et al.(2015)]ha15
            Ha S, Berner J, Snyder C. 2015. A Comparison of Model Error Representations in Mesoscale Ensemble Data Assimilation. Mon. Weather Rev., 143(10), 3893–3911.[Hamill and Whitaker(2005)]hamill05
            Hamill TM, Whitaker JS. 2005. Accounting for the Error due to Unresolved Scales in Ensemble Data Assimilation: A Comparison of Different Approaches. Mon. Weather Rev., 133(11), 3132–3147. [Hamill and Whitaker(2011)]hamill11
            Hamill TM, Whitaker JS. 2011. What Constrains Spread Growth in Forecasts Initialized from Ensemble Kalman Filters? Mon. Weather Rev., 139, 117–131.[Harlim and Hunt(2007)]harlim07
            Harlim J, Hunt BR. 2007. Four-dimensional local ensemble transform Kalman filter: Numerical experiments with a global circulation model.Tellus, Ser. A Dyn. Meteorol. Oceanogr., 59, 731–748.[Houtekamer et al.(2009)]houtekamer09
            Houtekamer PL, Mitchell HL, Deng X. 2009. Model Error Representation in an Operational Ensemble Kalman Filter, Mon. Wea. Rev., 137, 2126–2143.[Hunt et al.(2007)]hunt07
            Hunt BR, Kostelich EJ, Szunyogh I. 2007. Efficient data assimilation for spatio-temporal chaos: A local ensemble transform Kalman filter.  Physica D,  77, 437–471.[Jazwinski(2007)]jazwinski07
            Jazwinski AH. 2007. Stochastic processes and filtering theory. Courier Corporation.
            
            [Koyama and Watanabe(2010)]koyama10
            Koyama H, Watanabe W. 2010. Reducing forecast errors due to model imperfectness using ensemble
            Kalman filtering.  Mon. Wea. Rev.,  138, 3316–3332.
            
            [Leutbecher et al.(2017)]leutbecher17
            Leutbecher M, Lock SJ, Ollinaho P, Lang S, Balsamo G, Bechtold P, and Coauthors. 2017. Stochastic representations of model uncertainties at ECMWF: state of the art and future vision. Q. J. R. Meteorol. Soc., 143, 2315–2339.[Li et al.(2009)]li09
            Li, H., E. Kalnay, T. Miyoshi, and C. M. Danforth, 2009: Accounting for Model Errors in Ensemble Data Assimilation. Mon. Weather Rev., 137, 3407–3419.[Lorenz(1996)]lorenz96
            Lorenz EN. 1996. Predictability — A problem partly solved. Pp. 1-18 in Proceedings of seminar on predictability: Volume 1. ECMWF, Reading, UK.
            
            [Mitchell and Houtekamer (2000)]mitchell00
            Mitchell HL, Houtekamer PL. 2000. An adaptive ensemble Kalman filter. Mon. Weather Rev., 128, 416-433.[Mitchell and Carrassi(2015)]mitchell15
            Mitchell L, Carrassi A. 2015. Accounting for model error due to unresolved scales within ensemble Kalman filtering. ‎Q. J. Royal Meteorol. Soc, 141(689), 1417-1428. [Miyoshi(2011)]miyoshi11
            Miyoshi T. 2011. The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter. Mon. Weather Rev., 139, 1519-1535[Navon(1998)]navon98
            Navon IM. 1998. Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography. Dynamics of Atmospheres and Oceans, 27(1-4), 55-79. [Orrell(2003)]orrell03
            Orrell D. 2003. Model Error and Predictability over Different Timescales in the Lorenz’96 Systems. Journal of the Atmospheric Sciences, 60(17), 2219–2228.
            
            [Palmer(2001)]palmer01
            Palmer TN. 2001. A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parameterization in weather and climate prediction models.  Q.J.R. Meteorol. Soc.,  127, 279–304.[Palmer et al.(2009)]palmer09
            Palmer, T. N., Buizza, R., Doblas-Reyes, F., Jung, T., Leutbecher, M., Shutts, G. J., and co-authors (2009). Stochastic parametrization and model uncertainty.  ECMWF Tech. Memo, 598, 1-42.
            
            [Palmer and Weisheimer(2011)]palmer11
            Palmer TN, Weisheimer A. 2011. Diagnosing the causes of bias in climate models - why is it so hard? Geophys. Astrophys. Fluid Dyn., 105, 351-365, [Pulido et al.(2016)]pulido16
            Pulido M, Scheffler G, Ruiz JJ, Lucini MM, Tandeo P. 2016. Estimation of the functional form of subgrid-scale parametrizations using ensemble-based data assimilation: a simple model experiment. Q. J. R. Meteorol. Soc., 142(701), 2974-2984.[Pulido and Rosso(2017)]pulido17
            Pulido M, Rosso OA. 2017. Model Selection: Using Information Measures from Ordinal Symbolic Analysis to Select Model Subgrid-Scale Parameterizations. J. Atmos. Sci., 74(10), 3253–3269.
            
            [Pulido et al.(2018)]pulido18
            Pulido M, Tandeo P, Bocquet M, Carrassi A, Lucini M. 2018. Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 70(1), 1442099.
            
            [Ruiz et al.(2013a)]ruiz13a
            Ruiz J, Pulido M. Miyoshi T. 2013. Estimating parameters with ensemble-based data assimilation. A review.  J. Meteorol. Soc. Japan.  91, 79–99.  
            
            [Ruiz et al.(2013b)]ruiz13b
            Ruiz J., M. Pulido and T. Miyoshi, 2013: Estimating parameters with ensemble-based data assimilation. Parameter covariance treatment.   J. Meteorol. Soc. Japan.  91, 453–469. 
            
            [Ruiz and Pulido(2015)]ruiz15
            Ruiz J, Pulido M, 2015: Parameter Estimation Using Ensemble-Based Data Assimilation in the Presence of Model Error. Mon. Weather Rev., 143, 1568–1582.[Santitissadeekorn and Jones(2015)]santiti15
            Santitissadeekorn N, Jones C. 2015. Two-Stage Filtering for Joint State-Parameter Estimation. Mon. Weather Rev., 143, 2028–2042.[Shutts(2005)]shutts05
            Shutts G. 2005. A kinetic energy backscatter algorithm for use in ensemble prediction systems. Q. J. R. Meteorol. Soc., 131, 3079–3102.[Smith(2001)]smith01
            Smith LA. 2001. Disentangling uncertainty and error: on the predictability of nonlinear systems. Pp. 31—64  In Nonlinear dynamics and statistics . Ed. A. I.Mees. Birkhauser, Boston, USA
            
            [Stroud and Bengtsson(2007)]stroud07
            Stroud JR, Bengtsson T. 2007. Sequential State and Variance Estimation within the Ensemble Kalman Filter. Mon. Weather Rev., 135, 3194–3208.[Stroud et al.(2018)]stroud17
            Stroud JR, Katzfuss M, Wikle CK. 2018. A Bayesian adaptive ensemble Kalman filter for sequential state and parameter estimation. Mon. Weather Rev.,  146, 373–386.
            
            [Wang et al.(2007)]wang07
            Wang XG, Snyder C, Hamill TM. 2007. On the theoretical equivalence of differently proposed ensemble–3DVAR hybrid analysis schemes. Mon. Wea. Rev., 135, 222–227
            
            [Whitaker and Hamill(2012)]whitaker12
            Whitaker JS, Hamill TM. 2012. Evaluating Methods to Account for System Errors in Ensemble Data Assimilation. Mon. Weather Rev., 140(9), 3078–3089.[Wilks(2005)]wilks05
            Wilks DS. 2005. Effects of stochastic parameterizations in the Lorenz ’96 system.  Q. J. R. Meteorol. Soc.,  131, 389–407.
             
            [Yang and Delsole(2009)]yang09
            Yang X, Delsole T. 2009. Using the ensemble Kalman filter to estimate multiplicative model parameters. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 61, 601–609.
             
            [Zupanski(1997)]zupanski97 
            Zupanski D, 1997. A General Weak Constraint Applicable to Operational 4DVAR Data Assimilation Systems. Mon. Weather Rev., 125, 2274-2292. 
            
            § INTRODUCTION
            
            
            
             Model error treatment has  become a key ingredient for data assimilation systems. The most significant sources of the so-called model errors in numerical geophysical models are a consequence of simplifications in the representation of the dynamics, such as errors associated to the discretization of large-scale dynamics and unresolved or under-resolved physical processes represented by parameterizations. When combining observations with model forecasts in data assimilation systems, a proper account of both observational and forecast model sources of? errors is crucial for a successful state estimation. In data assimilation, model errors are usually separated into two components: systematic model error and random model errors. Systematic model errors are considered as the mean model error over a sufficiently long time window, while random model errors are the departures from said mean.
             
            The conventional approach for systematic model error treatment involves the estimation of a forecast bias term, which is augmented to the state vector during the data assimilation process <cit.>. In addition to explicit bias estimation, a random error term can also be incorporated to the state variables before the data assimilation step. The covariance structureprobability distribution of these errors can be inferred, for example, using innovation statistics from previous assimilation cycles  <cit.>. 
            
            In the context of the ensemble Kalman filter, neglecting random model errors usually results in an underestimation of the forecast error covariances. <cit.> gives a review of operational approaches to deal with these issues. The strategies typically involve a representation of model uncertainties either as an ad-hoc inflation of the forecast errors covariance matrix, multi-model and multi-parameterization ensembles or the use of stochastic parameterizations. The role of covariance inflation to compensate for model errors has been widely studied <cit.>. The inflation can be incorporated either as a random noise added to the different states of the ensemble <cit.>, or by amplifying the ensemble members deviations from the mean state <cit.>.  The latter approach however implicitly assumes that the model errors have the same structure as the dynamically evolved internal errors <cit.>. While this hypothesis does not hold for most of the atmospheric numerical models, the multiplicative inflation approach has been rather successfully used for this purpose as well as for dealing with sampling errors due to the small ensemble sizes used operationally.  The combination of parametric model error treatments with additive and multiplicative inflation was examined by <cit.>. A significant improvement in analysis error was found when compared to using each of them separately.
            
            Explicit representation of random model errors typically require an estimate of their spatio-temporal covariance structure, which can be either constructed arbitrarily, or parameterized as a function of free parameters. Several algorithms have been proposed to infer these type of parameters.  <cit.>  used innovation statistics to estimate model error horizontal decorrelation length scale and vertical covariances on a three-level quasi-geostrophic model. An iterative implementation of the expectation-maximization algorithm combined with an ensemble Kalman filter was successfully applied by <cit.> to infer different forms of covariance matrices of an additive Gaussian model error in a nonlinear state-space model.  On a similar model scenario, <cit.> proposed a Bayesian framework to explicitly account for the marginal posterior distribution of parameters, by  either an exhaustive grid-based exploration of the parameter space or using a particle filter approximation. This scheme can be combined with the ensemble Kalman filter to produce sequential estimations of state and parameters. Maximum likelihood estimators based on expectation-maximization and Newton-Raphson minimization were presented by <cit.>. These methods allow for the simultaneous estimation of deterministic parameters and parameters associated to stochastic processes representing model error. 
            
            
            A complementary approach for model error treatment in ensemble forecasting and data assimilation consists in the use of stochastic parameterizations. These parameterizations introduce stochastic processes directly to the model evolution, either as a random perturbation added at every time step of the model integration, or by stochastically perturbing the physical parameterizations tendencies. Stochastic parameterizations may account for a physically consistent representation of subgrid processes <cit.> and so, a state-dependent representation of model errors. The use of stochastic parameterizations has shown to have a positive impact in ensemble prediction skill <cit.>. Stochastic parameterizations are particularly successful in triggering noise-induced transitions <cit.>. When the spatial and temporal correlation structures are properly tuned, stochastic parameterizations may potentially provide a more consistent representation of the interactions between the resolved dynamics and the subgrid parameterizations <cit.>.Within a mesoscale data assimilation system, <cit.> showed that the use of stochastic backscatter scheme consistently outperformed the multiplicative covariance inflation scheme and the multiphysics ensemble approach. However stochastic parameterizations require carefully tuning the properties of the stochastic forcings in order to account for the model uncertainties. The parameters that characterize the stochastic process covariance matrix will be referred as stochastic parameters from now on.  Whereas optimization of deterministic parameters of the dynamical model is straightforward  within the ensemble Kalman filter using state augmentation <cit.>, stochastic parameters cannot be estimated in this way. Previous studies on the use of augmented state ensemble Kalman filter showed that the lack of correlation between the mean of the ensemble of state variables and the stochastic parameters may lead to completely unreliable estimations <cit.>.
            
            To overcome the lack of sensitivity of the forecast mean state, here we propose the application of two nested ensemble Kalman filters. The model state is estimated in an inner ensemble Kalman filter cycle, as in any conventional ensemble Kalman filter implementation. However, the filtering process is applied over different independent ensembles, each of them integrated using a different set of stochastic parameters. Hence, an ensemble of ensemble Kalman filters is assimilated in an outer filter cycle. The outer cycle is used for the estimation of the stochastic parameters. Stochastic parameters have a direct impact on the forecast error covariances, thus playing a critical role on the state analysis quality on the inner cycle. In order to recover the sensitivity of the analysis to stochastic parameters, we propose the use of a longer assimilation window in the outer cycle, composed of several internal cycles. In this way, the information in the outer filter for parameter estimation includes an ensemble of trajectories of analysed model states which is expected to be sensitive to the stochastic parameters. 
            
            Whereas few stochastic parameters can be tuned using an exhaustive exploration of the parameter space, the proposed technique can be used to simultaneously estimate several stochastic parameters. Its main advantage lays, indeed, in the ability to explicitly estimate elements of the covariance matrix used in the stochastic parameterization, at a computational cost that is comparable to the state-of-the art expectation-maximization algorithms <cit.>  and schemes based on coupling sequential Monte Carlo algorithms<cit.>. It should be remarked that proposed methodology based on nested ensemble Kalman filters is intended to be applied offline, as an optimization tool for both the stochastic dynamical model and the data assimilation system. 
            
            In Section 2 the stochastic parameter estimation based on the nested ensemble Kalman filters is introduced under a Bayesian framework. The dynamical model used in the experiments is described on Section 3.  In Section 4 we explain the experimental setup and the different covariance matrix structures, which are evaluated in the experiments. Results from stochastic twin experiments and from imperfect model experiments are shown in Section 5. We conclude with a brief summary and discussions in Section 6. 
              
            
            
            
            
            
            
            § METHODOLOGY
            
            
            
            To avoid the limitations of the conventional augmented state approach for ensemble based stochastic parameter estimation <cit.>, the proposed hierarchical inference technique is conducted using an ensemble of N_J  data assimilation systems, each of them using a different set of stochastic parameters. The technique involves two nested data assimilation cycles. Firstly, the state estimation is performed independently by each data assimilation system with a set of fixed stochastic parameter values, using the ensemble Kalman filter. This step is referred to as inner cycle. Secondly, parameter estimation is performed using the inner cycles mean a priori states of each of the ensembles. The mean a priori states are in turn the members of an ensemble that describes the density of the state conditioned to the parameter values. The Kalman filter equations are applied to this stochastic parameter ensemble to update their values. This step will be referred as outer cycle. This procedure is repeated sequentially. 
            
            The assimilation window of the parameter estimation cycle (i.e. the outer cycle) is composed by K state estimation cycles (inner cycles). Since stochastic parameters are assumed to change slowly with time, parameters are assumed static within the innerouter cycles and are denoted as θ_l, which means that the stochastic parameters will be updated only every K state estimation cycles to increase the identifiability of the stochastic parameters <cit.>.
            
            At a given time,  𝐱_l, k denotes the model state at the k-th state estimation cycle during the l-th parameter estimation cycle and θ_l denotes the parameters. For simplicity, we denote any quantity at time (l,0) with a single subindex l, so for example 𝐱_l,0=𝐱_l-1,K=𝐱_l. A schematic representation of time indexes is shown in Fig. <ref>.
            
            
            
            
            
            
            
            Given an initial prior joint density of the model state variables and parameters p(𝐱_0,θ_0), the inference goal is to sequentially obtain the posterior density of the augmented state at time l given K observation sets distributed in time from (l,1) to (l,K), this is p(𝐱_l+1,θ_l+1| 𝐲_l,1:K). The augmented state is estimated every K time steps, while the model state is estimated every time at which observations are available.Model state is then estimated every time at which observations are available, while the augmented state is estimated every K state assimilation cycles. Using the definition of conditional density, the joint parameter–state inference is given by 
            
            p(𝐱_l+1,θ_l+1| 𝐲_l,1:K)=p(𝐱_l+1| θ_l+1,𝐲_l,1:K) p(θ_l+1|𝐲_l,1:K). 
            The posterior density of the augmented state at the final time of a set of observations is given by the posterior density of the state conditioned to the value of the parameters and the observations multiplied by the posterior density of the parameters given the observations. The first term in the RHS of augPost is determined through the sequential process
            
            p(𝐱_l,k| θ_l+1,𝐲_l,1:k)=p(𝐲_l,k| 𝐱_l,k,θ_l+1,𝐲_l,1:k-1) p(𝐱_l,k| θ_l+1,𝐲_l,1:k-1)/p(𝐲_l,k|𝐲_l,1:k-1) 
            from k=1 to k=K.  In practice, sttPost corresponds to the application of a Kalman filter for the model state given a certain set of parameters and observations at time k.
            
            The parameter posterior density, the second term in the RHS of augPost, is rewritten through a sequential conditioning process as
            
            p(θ_l+1|𝐲_l,1:K)∝  p(θ_l+1|𝐲_l-1,1:K) ∏^K_k=1 p(𝐲_l,k|𝐲_l,1:k-1,θ_l+1), 
            where p(𝐲_l,1|𝐲_l,1:0,θ_l+1) ≐ p(𝐲_l,1|θ_l+1) and  p(θ_l+1 | 𝐲_l-1,1:K) is the prior density of the parameters, given the previous observations:  𝐲_l-1,1:K. Note that we leave implicit the conditioning to observations from previous outer cycles. As mentioned before, we have assumed parameters are static within an outer cycle. Alternatively, a random walk or any Markov process could also be considered, in that case p(θ_l+1 |   y_l-1,1:K ) = ∫ p(θ_l+1|θ_l) p(θ_l|𝐲_l-1,1:K) dθ_l. The parameters are assumed to be updated immediately after the k-th state assimilation cycle.
            
            
            Since the parameters are not directly observed, their likelihood to thea new observation 𝐲_l,k is taken into account through the marginalization of the model state,
            
            p(𝐲_l,k|𝐲_l,k-1,θ_l+1)=∫ p(𝐲_l,k|𝐱_l,k,θ_l+1) p(𝐱_l,k|𝐲_l,1:k-1,θ_l+1) 𝐱_l,k. 
            
            Equation (<ref>) involves a prediction density given the previous estimated state and the observation likelihood given the parameters, θ_l+1, which are then integrated along the model state. This marginalization is over the full model state and it is likely to produce errors if it is conducted directly through Monte Carlo sampling (e.g. using samples of a particle filter or an ensemble Kalman filter). Instead of using a direct sampling from the joint density, we take a Gaussian assumption for both the prediction density and the observation likelihood. Under this assumption, the sufficient statistics (𝐱_l,k| θ_l+1,𝐲_l,1:k-1), is used to assimilate the observations in (<ref>),  where (·) denotes the expectation operator. As it is known from the Rao-Blackwell theorem, the use of sufficient statistics in the estimator is expected to improve the inference of parameters. The Gaussian assumption is of course also taken for the ensemble Kalman filter that we use, so this assumption is coherent with the filter and does not imply an additional hypothesis.
            
            Following  the derivation of the sequential marginalized observation likelihood from <cit.>, the observation likelihood conditioned on  θ_l+1 is therefore
            
                p(𝐲_l,k|𝐱_l,k,θ_l+1) ∝exp[(𝐲_l,k - ℋ(𝐱_l,k))^T𝐑^-1 (𝐲_l,k - ℋ_l,k(𝐱_l,k))],
            
            where ℋ is the observation operator that transforms from model to ovservational space. Observational errors are assumed additive and Gaussian with  observational error covariance matrix𝐑. In coherence with the assumption taken in the ensemble Kalman filter, we now assume that the forecast density can be represented approximately by a Gaussian density, namely, [[Agrego dependencias de θ en x^f, porque se estan mencionando solo para 𝐏; y tambien una dependencia temporal del H]]
                p(𝐱_l,k|𝐲_l,1:k-1,θ_l+1) ∝exp[(𝐱_l,k - 𝐱^f_l,k(θ_l+1))^T𝐏_l,k(θ_l+1)^-1 (𝐱_l,k - 𝐱^f_l,k(θ_l+1))],
            
            where  𝐱^f_l,k≐ (𝐱_l,k|y̌_l,1:k-1,θ_l+1)≐∫x̌_l,k p(x̌_l,k|y̌_l,1:k-1,θ_l+1) d x̌_l,k  is the mean forecast state conditioned on θ_l+1 and 𝐏_l,k is the forecast error covariance given θ_l+1. The dependencies on θ_l+1 are dropped from equations to reduce notation and definig 𝐇 as the linearized observation operator, tThe resulting approximated observation likelihood in the time interval (l,1) to (l,K) is
            
            
            
                ∏_k=1^K p(𝐲_l,k|𝐲_l,k-1,θ_l+1)∝∏_k=1^K exp[(𝐲_l,k - _l,k(𝐱^f_l,k))^T (𝐇_l,k𝐏_l,k𝐇_l,k^T+𝐑 )^-1 (𝐲_l,k - _l,k(𝐱^f_l,k))],
            
            which is equivalent to the approximated observation likelihood obtained in <cit.>.  That work also constrains the values of the statistical parameters within a time interval where K observation sets are available. However,  a point estimation of the parameters is conducted there through maximization of the logarithm of the approximated observation likelihood. On the other hand, here we assume a Bayesian framework, see  bayesPar, in which an inference of the density of the parameters conditioned to the set of observations is obtained given some prior knowledge of the parameters. Our procedure resembles the Rao-Blacwellized particle filter <cit.>, where N_J ensemble Kalman filters are conducted in order to marginalize the parameters posterior distribution.
            
            
            In this work, the parameters to be estimated are assumed to be associated with an additive Gaussian model error. We also assume that parameters follow a Gaussian distribution. While the latter hypothesis is not warranted, it allows us to treat the state-parameter estimation problem by using two nested ensemble Kalman filters.
            
            
            Let us consider an ensemble of initial parameters θ_0^(j), with j=1,⋯,N_J, sampled from  p(θ_0). Each parameter θ_0^(j) is associated with an ensemble of N_I model states {𝐱_0^(j,i), i=1,⋯,N_I}. Therefore, a set of N_J ensembles is initialized and each of them represents different values of the parameters θ. The updates to the  ensemble state members are determined in the inner cycles with ensemble Kalman filters. The filters should be identical to the data assimilation system for which parameters are being estimated, e.g. same physical parameterizations, number of ensemble members. Note that the the N_J ensembles are assumed to evolve independently, hence the state update neglects any correlation between the ensembles. The mean state of the j-th ensemble is given by
            
                𝐱_l,k^a(j)=𝐱_l,k^f(j) + 𝐏^(j)_l,k𝐇^T (𝐇_l,k𝐏^(j)_l,k𝐇_l,k^T + 𝐑)^-1  [𝐲_l,k-  ℋ_l,k(𝐱_l,k^f(j))],
            
            where 𝐱^a(j), 𝐱^f(j) and 𝐏^(j) denote the analysed state, forecasted state and the forecast covariance matrix of the j-th ensemble respectively.
            
            
            Next, the parameter posterior density conditioned to the observations, p(θ_l+1|y̌_l,1:K), is inferred using bayesPar. In this sense, Equation bayesPar can be interpreted  as the serial assimilation of observations along the l-th state assimilation window. Under the already taken assumptions, the likelihood function and the forecast density are assumed Gaussian in  bayesPar. Therefore, in the “outer” cycle we apply the ensemble Kalman filter to infer p(θ_l+1|y̌_l,1:K). Because the initial condition of the hidden state is not known with complete certainty, the observational error for the inference of the parameter posterior density increases with the forecast error covariance matrix of the hidden state (see Eq. <ref>) for the inference of the parameter posterior density, anThus, anincreased“observational” error covariance matrix for the assimilation of parameters  is obtained in obsLik, ȞP̌_l,k(θ_l) Ȟ^+Ř. 
            
            
            The application of the Kalman filter to the ensemble representing the parameter distribution bayesPar with obsLik results in the analyzed mean parameters given by
            
            
                θ_l+1^a = θ^f_l+1 + ∑_k=1^K 𝐏^θ𝐱_l,k𝐇^T[𝐇  (𝐏_l,k + 𝐏^𝐱𝐱_l,k) 𝐇^T + 𝐑]^-1 [𝐲_l,k-  ℋ(𝐱^f_l,k)],
             where 𝐱^f_l,k is the average of ensemble forecast mean states at time (l,k) over the N_J ensembles. The empirical covariance matrices 𝐏_l, 𝐏^𝐱𝐱_l,k and 𝐏^θ𝐱_l  are defined as follows
            
            
            
            
                𝐏_l,k= 1N_J∑_j=1^N_J𝐏^(j)_l,k =1(N_I-1) N_J∑_j=1^N_J∑_i=1^N_I (𝐱^f(j,i)_l,k - 𝐱^f(j)_l,k)  (𝐱^f(j,i)_l,k - 𝐱^f(j)_l,k)^T
            
             is the sample forecast state covariance averaged among the N_J ensembles,
            
            
            
                𝐏^𝐱𝐱_l,k = 1/N_J-1∑_j=1^N_J (𝐱^f(j)_l,k-𝐱^f_l,k) (𝐱^f(j)_l,k-𝐱^f_l,k)^T
            
            
            is the sample covariance of the N_J ensembles, and similarly,
            
            
                𝐏^θ𝐱_l,k = 1/N_J-1∑_j^N_J  (θ^f(j)_l+1-θ^f_l+1) (𝐱^f(j)_l,k-𝐱^f_l,k)^T
            
            is the parameter-state covariance matrix.
            
            Note that the forecast state error covariance for the parameter estimation in the outer cycle is  the sum of the mean covariance of the N_J ensembles and the covariance of the outer ensemble. Equation lparUpdate shows that parameters are estimated using ensemble state means as individual state members in the outer cycle. This formulation defines a parameter-state covariance matrix that is able to transfer model state innovations to parameters. This is a key difference between the nested ensemble Kalman filters and the  standard state augmentation for parameter estimation, since the impact of stochastic parameters is accounted from an ensemble mean, and not over individual members.
            
            
            
            
             §.§ Implementation details
            
            
            In this work both the inner cycles and the outer cycles are based on the Ensemble Transform Kalman Filter <cit.>. For the outer cycles, the parameter update equation lparUpdate is implemented as an asynchronous ETKF <cit.>.  For this purpose, an aggregated vector is constructed by column-wise concatenating observations 𝐲_l,1:k in a single observation vector 𝐲^∗_l. A similar concatenation is performed with the ensemble members in the state and in the observational spaces. The aggregated observational error covariance matrices  𝐑^∗_l and mean covariance 𝐏^∗_l are constructed with the k-th diagonal block 𝐑_l,k and 𝐏_l,k respectively.  
            The nested ensemble transform Kalman filters in this work are implemented as follows:
            
                  
            
              *  Given N_J parameters θ^(1:N_J)_0 and N_J independent ensembles of N_I state members, x̌^(1:N_J,1:N_I)_0,  and (L × K) observations y_1:L,1:K.	
            	
              *  State estimation: For each assimilation cycle l, (with l=1,⋯,L) do:
            	
              *  For each inner assimilation cycle k, (with k=1,⋯,K) do:
            		
              *  Calculate the ensembles of analysed states 𝐱_l,k^a(j,i)  performing N_J-EnKFs independently.
            	 		
              *  Store  ℋ_l,k(𝐱̅_l,k^f(j)) for each ensemble and the average of the forecast error covariance matrix in the observational space over the N_J ensembles 𝐇𝐏_l,k𝐇^.
            		
              *  Parameter estimation: Obtain the ensemble of estimated parameters:
            	
              *  Concatenate the K mean predicted observations ℋ_l,k(  𝐱̅_l,1:K^f(j)) to construct an (n_x× K)-dimensional ensemble of N_J members.
            		
              *  Construct the agreggated observation vector 𝐲^∗_l=[𝐲_l,1,⋯,𝐲_l,K]^T and the tangent linear observation operator 	𝐇_l^∗=[𝐇_l,1, 𝐇_l,2, ⋯,𝐇_l,K]^T
              *  Construct the block diagonal extended observational error covariance matrix 𝐑_l^∗, whose k-th diagonal block is 𝐑^∗_l,k=𝐇_l,k^T𝐏̅_l,k𝐇_l,k+𝐑.
            		
              *  Obtain the updated parameter ensemble mean and perturbation ETKF using the aggregated matrices calculated in steps <ref>-<ref>
            
            Since the experiments were conducted in a low-dimensional dynamical system (see Section <ref>), the use of covariance localization is not explored in this work. Localization may become mandatory for systems in which the number of state-space dimensions exceeds the number of ensemble members, as occurs in numerical forecast models. Additionally, we have assumed that sampling errors of the filter can be partially accounted by the stochastic parameterization so that multiplicative covariance inflation of state variables is not included in most of the experiments. As has been described in <cit.> and <cit.> assuming a persistence model for the parameters can result in the collapse of the parameter ensemble spread and divergence of the parameter estimation. However, in our experiments, such mechanisms were not needed to increase parameter ensemble spread. 
            
            § DESCRIPTION OF THE EXPERIMENTS
            
            
             §.§ The Lorenz-96 dynamical model
            
            The two-scale Lorenz-96 dynamical model has been extensively used as a testbed model for the development of data assimilation schemes due to its reduced computational cost and its ability to mimic specific properties of the atmospheric predictability <cit.>. It represents the dynamics of a cyclical set of large-scale variables over a circle of latitude, each coupled to a set of high-frequency small-scale variables.  Each model equation contains terms that represent non-linear advection, dissipation and external forcings. The small-scale variables are coupled to the large-scale variables through an additive forcing term. 
            
            The set of equations of the two-scale Lorenz-96 dynamical model is given by large-scale variable equations,
            
            
                dx_ndt=-x_n-1(x_n-2-x_n+1)-x_n+F-hcb∑_m=M(n-1)+1^Mn y_m,
            
            and small scale variable equations,
            
                dy_mdt=-cby_m+1(x_m+2-x_m-1)-cy_m+hcb x_1+int[m-1/M],
            
            where n=1,⋯,N and m=1,⋯,MN. 
            
            Both sets of variables have cyclic boundaries conditions: x_n+N=x_n and y_m+MN=y_m.  In this work, the coupling and scale parameters are set to the standard values of h=1, b=10 and c=10 as in <cit.> and <cit.>. The number of large-scale variables was set to N=8, each coupled to M=32 small-scale variables y, so that MN=256. To achieve a chaotic solution, the external forcing is set to F=20 for all the experiments. 
            
            In the imperfect model experiments, the small-scale variables can be interpreted as unknown physical processes which cannot be explicitly resolved in numerical models, so that only the dynamics of large-scale variables are represented by the model. The effect of the small-scale variables is introduced as a parametrization that is a function of the resolved large-scale variables only. This mimics in a very simple way model errors associated with the parametrization of unresolved processes in realistic atmospheric or oceanic numerical models.  The Lorenz-96 system results particularly suitable for proof-of-concept experiments involving subgrid model error representation and parameterizations <cit.>. The truncated version of the model can be expressed as
            
            
                dx_ndt=-x_n-1(x_n-2-x_n+1)-x_n-U(x_n),
            
            where U(x_n) represents the parameterization of small-scale processes. The forcing term F is also assumed to be part of the parameterization U.  In this work, the parameterization is of the form
            
            
                U(x_n)=a_0+a_1 x_n+e_n(t).
            
            The first two terms represent a deterministic forcing that is a function of only the resolved variable x_n. The coefficients a_0 and a_1 can be estimated via a least-squares fitting using an integration of the complete system (Eqs. <ref>-<ref>) as in <cit.>, or inferred via data assimilation using only noisy observations of the resolved variables of the full system <cit.>.
            
            The processes thath cannot be accounted by a deterministic function of the state variables, are included as a state-independent red-noise stochastic forcing, discretized as the realization of a zero mean first-order autoregressive process (AR(1)),
            
                𝐞(t)=ϕ 𝐞(t-Δ t) + (1-ϕ^2)^1/2 η.
            
            
            The coefficient ϕ represents the lag-1 autocorrelation of 𝐞(t), and Δ t is the model integration timestep. The vector η∼𝒩(0,Σ) represents a random draw from a zero-mean Gaussian distribution with covariance Σ. The adequacy of a stochastic parametrization for model error representation in the Lorenz-96 model was proved by <cit.>, <cit.> and by <cit.>.
            
            
            
            
            
             §.§ Experimental setup 
            
            
            We first evaluate the nested ensemble Kalman filters using twin experiments. In these experiments, the “true" integration is generated byconsists of an integration of the truncated model l96par-usto with a stochastic forcing generated using a prescribed covariance structure. The same model and covariance structures are then used as forecast model, but with uncertain parameters. Since stochastic processes are present in both “true" and forecast models, it is not possible to replicate the true integration using the forecast model, even when using identical initial conditions and parameter configuration. Hence, these experiments are referred to as “stochastic twin experiments".  
            
            
            The truncated Lorenz-96 dynamical system was integrated with a fourth-order Runge-Kutta  scheme. The system was initialized after a spinup of 1460 dimensionless model time units, which is roughly equivalent to 20 years of atmospheric evolution.  The nature run was generated integrating the “true model" for 250 model time units (i.e. 50000Δ t) with a timestep of Δ t=0.005. The coefficients a_0 and a_1 were set to a_0=19.169 and a_1=-0.813. These values were estimated via least square fitting using an integration of the two-scale Lorenz-96 system with F=20.  Synthetic observations are then generated by perturbing the nature run with zero-mean Gaussian uncorrelated noise of variance 𝐑=σ_R^2 𝐈 and σ_R^2=1, where 𝐈 is the identity matrix.  All the variables are observed simultaneously, with a frequency of δ t=10Δ t. 
            
            In the experiments, ensembles of N_I=30 members are used for the inner cycles, whereas for the outer cycle,   N_J=15 independent ensembles are considered. Initial conditions for the states of the ensembles were randomly chosen from the true model integration.  The number of inner cycles within each outer cycle is set to K=5.  This value was chosen to balance the parameter convergence speed, precision and computational cost, through preliminary sensitivity experiments.
            
            
            Imperfect model experiments are conducted using the two-scale Lorenz-96 model. In this case, the true integration corresponds to an integration of the two-scale Lorenz-96 model during 250 model time units, with a time step of Δ t =0.001.  Observations are generated as in the stochastic twin experiments.  Note that only the large-scale variables are observed in these experiments and that the truncated Lorenz-96 model is used as forecast model, with the same integration scheme as in the stochastic twin experiments. 
            
            To evaluate the sensitivity of the estimations to observational sampling errors, each of the proposed assimilation experiments is repeated 10 times, with different realizations of observational error and changing the ensemble of initial states and parameters and the seed of the random number generator in the stochastic parameterization.  The computation of verification scores excludes the first 200 state assimilation cycles to avoid the effect of the filter spinup.
            
            
            
            
             §.§ Model error covariance structure
            In this work, the nested ensemble Kalman filters are evaluated using four different stochastic forcing covariance matrix structures Σ. These proposed covariance matrix structures represent different hypothesis of the behavior of the model error which are usually used in practice and are used in both the stochastic twin experiments and in the imperfect model experiments of the upcoming sections. An explicit estimation of the stochastic processes covariance matrix might result intractable for geophysical models. In practice, several assumptions and simplifications can be considered in  an attempt to replicate the structure of said covariances. In this work, the nested ensemble Kalman filters are used to estimate parameters related to different structures of the covariance matrix Σ.  The proposed parameterizations of the covariance matrix represent different hypothesis of the behavior of the model error which are usually assumed in practice (i.e. Gaussian errors, spatially symmetric covariances, isotropy). In this work, parameters for the following covariance matrix structures Σ are estimated:
            
              IIsotropic non-correlated: In this case the covariance matrix is expressed as Σ= σ^2 𝐈 where 𝐈 is the identity matrix. This model assumes that the model error variance is the same for all the resolved variables and that model errors for different model variables are uncorrelated. In this case, the standard deviation σ is the only parameter to be estimated.
            
            
              * Isotropic exponential covariance:  The covariance is parameterized as Σ_i,j=σ^2 e^-ρ d_i,j, where d_i,j is the minimum distance between variables x_i and x_j, indicating an exponential spatial decrease of covariances. This approach assumes again that the variance of the model error is the same for all the variables, but it incorporates an a priori spatial covariance structure.  Smaller values of ρ are associated with longer model error spatial correlations.  In this case, the parameters to estimate are the standard deviation σ and the spatial scale parameter ρ.  No estoy convencido de que esta sea la covarianza de un SOAR. Por lo que vi en los papers, me faltaria un coeficiente (1+ρ d_i,j) para que sea exactamente eso. i.e. https://journals.ametsoc.org/doi/pdf/10.1175/1520-0493%281999%29127%3C2293%3ATDCFFN%3E2.0.CO%3B2 
            
            
              * Horizontally symmetric homogeneous covariance matrix  Σ: All the variables are assumed to have the same spatial covariance structure and that covariances are horizontally symmetric (namely Σ_n,n-i=Σ_n,n+i). The stochastic parameters to estimate in this case are the variance σ^2 and the model error neighbouring covariances (in our model only 5 parameters). 
            
            
              * Non-isotropic non-correlated: This covariance structure ignores spatial correlations but assumes that the stochastic forcing associated to each variable has a different standard deviation. In this case  the covariance matrix is represented as Σ=diag(σ_1^2,σ_2^2,⋯,σ_N^2) representing a spatially heterogeneous model error distribution and the parameters {σ_1^2,σ_2^2,⋯,σ_N^2} are estimated independently.
            
            
            The temporal autocorrelation parameter in the experiments was fixed to  ϕ=0.984, as in <cit.>, representing a persistent stochastic forcing. Although the autocorrelation parameter could be included for estimation, it has been found that the optimal solution is not unique. There is a wide range of optimal combinations, in terms of root mean squared error and ignorance skill scores, between the stochastic forcing amplitude and the autocorrelation time parameter <cit.>.
            
            
            
            
            
            
            § RESULTS WITH STOCHASTIC TWIN  EXPERIMENTS
            
            
             §.§ Isotropic non-correlated stochastic noise
            
            
            The nature run for the first experiment uses the isotropic non-correlated covariance structure (case I) with σ_nat=2. The only parameter to estimate is the standard deviation σ. The initial valuesensemble for parameter σwerewas drawn from a 𝒩(1.5,0.5^2) distribution.
            
            
            Estimation results for parameter σ are shown in Fig. <ref>a. In all the experiments, the parameter values converge rapidly during the first 10 parameter assimilation cycles10 (50) outer (inner) assimilation cycles?. Independently of the values used in the parameter prior distribution at time l=0, estimated parameters converge toward a narrow range of values after approximately 250 parameter assimilation cycles  (1250 state assimilation cycles). Though different experiments do not converge to an identical parameter value, final estimations in the different experiments have a small relative standard deviation of less than 1.2%.  Figure <ref>b shows the parameter ensemble evolution during the first 100 parameter assimilation cycles for one of the experiments shown in Fig. <ref>a. After the short assimilation spinup period, parameter updates are small and the parameter ensemble spread remains relatively stable throughout the duration of these experiments.
            
            
            
            
            
            The mean estimation of σ averaged among different experiments is σ^a=2.2, which is slightly larger than the parameter used in the nature integration σ_nat=2. To analyze the validity of the inferred parameter, we conducted an exhaustive sampling of the parameter space. For this purpose, we performed data assimilation experiments in which only the state variables were assimilated and different fixed values of the parameter σ were used during the entire assimilation experiment. Values of σ were regularlyevenly distributed, covering the range [1.25,3.25] (with Δ_σ =0.05). Each data assimilation experiment consisted of 2300 assimilation cycles, excluding an initial spinup of 200 assimilation cycles. To avoid sampling issues, the experiments were repeated 25 times for each  parameter value using different observational errors and ensemble of initial conditions for the first assimilation cycle. The mean over space, time and different experiment realizations of the analyzed state RMSE is shown in Fig. <ref>. A clear global minimum is found in the experiments. The cost function has an overall convex geometry, with increased sensitivity towards smaller values of σ. The optimal standard deviation parameter found through exhaustive parameter evaluation was σ_ex=2.15, which is also larger than σ_nat.  The  discrepancy found between σ_ex and σ_nat is expected to be a consequence of the usage of a finite ensemble without using multiplicative or additive covariance inflation in the data assimilation process.  Hence, the larger stochastic noise amplitude attempts to correct sampling errors due to a finite ensemble size. As shown in Fig. <ref>, the values of σ estimated with the nested ensemble Kalman filters is rather coherent with this cost function.
            
            
            
            
            
             §.§ Parameterized spatial correlations
            
            
            In these experiments, stochastic parameterizations of the true model and the forecast model  use the the isotropic double exponential covariance structure (Case II), with a decaying function.
            
            The parameters to estimate in these experiments are σ and ρ. In the nature integration, the standard deviation of the process is set to σ_nat=2,  and the decorrelation scale parameter is set to  ρ_nat=0.3. The latter leads to a moderate decaying rate, i.e. the covariance between the most distant variables is Q_i,i+4≈ 0.3σ^2. Initial values for parameter σ are again drawn from a 𝒩(1.5,0.5^2) distribution, while a 𝒩(0.5,0.15^2) distribution
            
            Results from simultaneous estimations of σ and ρ are shown in Fig. <ref>. Assuming complete ignorance of the parameter values used in the nature integration, on average, the estimations converge to parameter values σ^a=2.12 and ρ^a=0.29. The estimated values for σ are on average at least 5% larger than the value used in the nature integration. It is worth reminding that no inflation is being added to the state ensemble, so the variance overestimation may be associated to the requirement of additional covariance inflation to alleviate the effect of sampling errors.  
            
            The optimal parameter combination estimatedfound? through exhaustive search that minimizes RMSE is σ_ex=2.15 and ρ_ex=0.36. Therefore, the optimal values estimated with the nested filters are slightly biased towards lower values of both parameters.  Smaller values of ρ are associated with larger covariances between distant variables, hence the experiment with exhaustive evaluation of the parameter space suggests weaker spatial correlations and an inflated variance. The RMSE cost function is again convex and asymmetric, especially for ρ (Fig. <ref>). However, all the estimations lay close to the set of parameter values that produce the minimum RMSE (blue dots in Fig. <ref>). The RMSE associated with these estimations is at most 1% larger than the minimum RMSE.
             
            
            The ETKF implementation of the nested ensemble Kalman filters requires the inversion of matrix 𝐑^∗=𝐇𝐏̅𝐇^ + 𝐑 in every outer cycle. A significant reduction of the computational cost is obtained if  𝐏̅_l,k can be assumed diagonal for R^∗ computation. Results for different repetitions of the experiment considering a diagonal 𝐏̅_l,k are marked with stars on Fig. <ref>. Slight differences in the estimated parameters are found when this assumption is considered. On average, the parameter σ is approximately 1% larger than when using off-diagonal elements of 𝐏̅_l,k, while the difference in ρ is almost negligible (i.e. about 0.5% smaller). The effect on the analyzed state RMSE is rather small (<0.01%). The practical tweak of assuming a diagonal matrix 𝐏̅_l,k does not degrade significantly the quality of estimation, while reducing its computational cost. Note that if a non-square root ensemble Kalman filter was used in the outer cycle, there would be no computational benefits in assuming P̃̅̃_l,k to be diagonal, since such schemes would require computation of (𝐇𝐏̃^𝐱𝐱𝐇+𝐇𝐏̅𝐇 + R)^-1.
            
            
            
            
            
            
            
            
            
            
            
            
             §.§ Non-isotropic variance estimation
            
            
            In state-of-the-art geophysical models, model error is usually non isotropic since might be affected differently by model errors. It is interesting to study if the proposed technique can retrieve the structure of Σ when removing the isotropic assumption (i.e. using covariance model III). With this purpose, we propose an idealized experiment in which the stochastic parameterization is driven by an uncorrelated zero-mean Gaussian process, with variances σ^2_1=σ^2 _4=2.5^2 and the rest of the variances set to 1.5^2. The number of parameters to be estimated in this case is 8. The initial ensemble of parameters was chosensampled? from  𝒩(2;0.5^2).   
            
            
            
             
            
            Figures <ref>a,b show the evolution of parameters σ^2_1 and σ^2_2 for independent experiment repetitions.  In most experiments, estimated parameters are larger than the the values used in the nature integration, with an net difference of up to +0.5.  Similar results were found for the remaining parameters (Fig. <ref>c). The overestimation is expected to compensate for the limited ensemble size in the state ensemble. In spite of the noise, estimations on most experiments converged to the same parameter range after a spin-up period of around 300 outer cycles.
            
            InUnder? this experimental setting, the number of parameters to estimate is equivalent to the number of state variables,  and is the experiment with the largest number of estimated parameters. Because of this,This is the experiment with the largest number of estimated parameters, therefore it is important to evaluate the impact of the parameter ensemble size  upon the quality of the estimations. Table <ref> shows the temporal standard deviation of the estimated parameters over the last 300 parameter assimilation cycles for experiments using different number of ensemble members in the outer cycle. Values were averaged over 25 experiments with different observational error realizations and initial parameter ensembles. The variability of the estimated parameters has a large dependence with the number of ensembles used in the outer cycle. In particular, an improved convergence is found when the number of ensembles is increased.  In the case with N_J=5 ensembles, the parameter estimations show large variations among experiments. This suggests that for larger dimensional systems, the use of covariance localization in the parameter space may become mandatory.  Note, however, that the RMSE of these experiments in average  does not improve significantly when using more than N_J=15 ensembles, with the caveat that we are estimating only 8 parametersfor a parameter space with 8 degrees of freedom. [[Sugerencia de Juan, para no llamarlo "caveat"]]
            
            § RESULTS WITH THE TWO-SCALE LORENZ96 MODEL
            
            
            For the imperfect model experiments, we use the two-scale Lorenz-96 model to generate the nature run. In the first set of experiments, the stochastic forcing used in the truncated model has covariance structure I, then, the only parameter to estimate is σ. Figure <ref>a shows estimation results for independent experiment repetitions. In most cases, convergence is achieved during the first 300 parameter assimilation cycles. The mean estimated value is σ^a=1.92. These estimations are compared with exhaustive evaluation of the parameter space (Fig. <ref>b). The  minimum value found  through exhaustive search corresponds to σ_ex=1.95, which is very close to the mean estimated value and also within the range of the estimated parameters in the different realizations of the estimation experiments.
            
            Unlike the experiments in the previous section, the optimal structure of the covariance matrix Σ for the imperfect model scenario is not known and the covariance structure I may result in a suboptimal representation of the model error. Thus, we inferred empirically, and offline, the characteristics of the covariance matrix that best fits the truncated model to the two-scale Lorenz-96 system. For these diagnostics, the true state evolution is assumed known including the evolution of the small-scale variables, contrary to the data assimilation experiments in which we assume that we only know a set of noisy observations of the large-scale variables. A large integration of 10000 model time units of the two-scale Lorenz-96 system was conducted. Using the least-squares  deterministic parameters a_0 and a_1,  the covariance of the residuals are given by
            
                r(x_n,t)=[U_det(𝐱_n,t) - ℱ(𝐱_n,t)],
            
            where U_det  is the forcing estimated by the deterministic parametrization (first two terms in udet) and ℱ is the forcing obtained in the two-scale Lorenz-96 system (last two terms in Eq. <ref>)
            
                ℱ(x_n,t)=F-hcb∑_m=M(n-1)+1^Mn y_m
            
            
            The covariance of the residuals r(𝐱,â_0,â_1,t) can be seen as an approximation of the model error covariance matrix Σ of the truncated Lorenz-96 model when using only the deterministic part of the parameterization.  Figure <ref>a shows the covariance of the residuals. Model errors have a variance of σ^∗ 2=4.93± 0.07, while the covariances between neighboring variables are σ^∗_i,i±1≈ -0.55 and σ^∗_i,i± 2. ≈ 0.7. Similar model error covariance structures were found for other configurations of the two-scale Lorenz-96 system <cit.>. Inferring this type of intricate covariance structure is not straightforward.
            
            
            
            These results are compared with parameters estimations calculated via exhaustive sampling of the parameter space. For this case, 10 independent experiments were performed, using different observational errors. Due to the exponential growth of computational cost of this methodology, parameter space was explored with a spatial grid of Δ_σ=0.125 and a 5-dimensional guess given by the nested ensemble Kalman filters. While both cases have a similar variance (σ^∗ 2≈ 4.93 and σ^2_ex≈ 4.25), the optimal stochastic forcing covariance between distant variables is significantly larger than the ones estimated offline (Fig. <ref>b). However the offline estimated parameters are not expected to be optimal in an RMSE sense for a data assimilation system <cit.>. The covariances estimated with the nested ensemble Kalman filters should ideally be closer to the ones found through the costly exhaustive exploration.
            
            We evaluate the potential of the nested ensemble Kalman filters to uncover the covariance structure using covariance model III, since it is flexible enough to represent the complexe covariance associated with model error in the truncated Lorenz-96 equations. Estimation results of 10 repetitions of the experiment  are shown in Fig. <ref>. The estimations are less precise than in the previous experiments and require around 400 parameter assimilation cycles to converge. The mean parameter values obtained with the nested ensemble Kalman filters are in general consistent with the offline estimations shown in Fig. <ref>a, but with pronounced differences in the magnitudes of the off-diagonal elements. However, estimated values are closer to the parameter values that effectively minimize the analysis RMSE (dashed lines). The nested ensemble Kalman filters are able to accurately estimate the variance  σ^2_i and the first covariance σ_i,i± 1. It is also able to recover the sign of the upcoming covariances. Further experiments are needed to assess the possibility of estimating more distant covariances.
            
            
            
            
            
            § DISCUSSION
            
            
            In this work we introduce a novel data assimilation technique to infer stochastic parameters based on a nested implementation of the ensemble Kalman filter.  The estimation of stochastic parameters requires an ensemble Kalman filter where each member is associated to a data assimilation system that is identical to the system whose parameters are being estimated i.e. model configuration, resolution, number of ensemble members. In this way, the technique not only can be used to infer parameters for stochastic parameterizations,  but it also can be used to estimate other hyperparameters associated to the data assimilation system. Whereas, we implemented an ensemble transform Kalman filter <cit.> for both state and stochastic parameter estimation, the nested ensemble Kalman filters can be extended to different flavours of the ensemble Kalman filters and the use of other data assimilation schemes in the inner cycles such as hybrid-variational. For the outer cycle, the use of particle filters is also possible, especially if the distribution of the parameters is strongly non-Gaussian, resembling the Rao-Blackwellized particle filter <cit.>.
            
            It is worth noticing that the proposed technique is intended to be used offline, as a tool for a priori parameter estimation. While the computational cost of the proposed technique is relatively large, it proves to be remarkably more economic than manually tuning model parameters, particularly when the number of parameters to estimate increases. The computational cost might be comparable with other state of the art schemes like the expectation-maximization algorithm <cit.> and less expensive than SMC^2<cit.>.
            
            In this work we proved that the technique is able to successfully estimate parameters on stochastic twin experiments with simple model error covariances structures like the double exponential decay function or the diagonal isotropic case.  The estimated parameters are close to the optimal parameter values found through an exhaustive exploration of the parameter space, at a significantly lower computational cost. The technique is also robust for the simultaneous estimation of multiple stochastic parameters. Additionally, and more importantly, the nested ensemble Kalman filters are able to recover the structure of model error covariances in an experiment with missing subgrid dynamics, without making any a priori assumptions on the covariance structure nor the missing physics. The experiments were performed on a low-dimensional chaotic model, further experiments in high-dimensional systems are required for which  covariance localization both in the state and in the parameter space becomes necessary. Further research is also required to evaluate the potential of reconstructing model error covariances between different types of variables.
            
            The stochastic parameterizations used in this work couple the model error representation to the model dynamics, by incorporating the stochastic forcing directly on the model equations. This is an important difference with respect to other model error treatment schemes that incorporate background state perturbations in the instant prior to the assimilation. Additionally, in this work, the stochastic forcing is assumed to be state independent. In the Lorenz-96 model, the amplitude of stochastic perturbations may be largely influenced by its associated state variable <cit.>. The proposed scheme is expected to handle state-dependent stochastic parameterizations, as well as other like, stochastically perturbed tendencies parameterizations <cit.>.
            
            
            The possibility of estimating other types of other hyperparameters in the context of the ensemble Kalman filter or other assimilation methods remains unexplored. The nested ensemble Kalman filters, in principle, could be applied to the estimation of parameters related to the observational error covariance matrix, as well as covariance localization length-scales for state estimation. Hybrid schemes, like the ensemble 4DVar <cit.>, could benefit from the proposed technique for inferring the optimal covariances weighting coefficients.
            
            
            
            
            50[Aksoy et al.(2006)]aksoy06
            Aksoy, A., Zhang, F., Nielsen-Gammon, J. (2006). Ensemble-Based Simultaneous State and Parameter Estimation in a Two-Dimensional Sea-Breeze Model. Monthly Weather Review, 134(2003), 2951-2970.
            
            
            [Annan et al.(2005)]annan05
            Annan, J. D., J. C. Hargreaves, N. R. Edwards, and R. Marsh, 2005: Parameter estimation in an intermediate complexity earth system model using an ensemble Kalman filter. Ocean Model., 8, 135-154.[Anderson and Anderson(1999)]anderson99
            Anderson, J. L., and S. L. Anderson, 1999: A Monte Carlo Implementation of the Nonlinear Filtering Problem to Produce Ensemble Assimilations and Forecasts. Mon. Weather Rev., 127, 2741-2758.
            
            [Anderson(2001)]anderson01
            Anderson J. 2001. An Ensemble Adjustment Kalman Filter for Data Assimilation. Mon. Wea. Rev., 129, 2884–2903.
            
            [Anderson(2009)]anderson09
            Anderson, J. L. (2009). Spatially and temporally varying adaptive covariance inflation for ensemble filters. Tellus, Series A: Dynamic Meteorology and Oceanography, 61 A(1), 72-83. https://doi.org/10.1111/j.1600-0870.2008.00361.x
            
            [Arnold et al.(2013)]arnold13
            Arnold H.M., Moroz I.M., Palmer T.N., 2013: Stochastic parametrizations and model uncertainty in the Lorenz 96 system.  Phil Trans R Soc A, 371, 20110479 [Baek et. al(2006)]baek06
            Baek, S. J., B. R. Hunt, E. Kalnay, E. Oott, and I. Szunyogh, 2006: Local ensemble Kalman filtering in the presence of model bias. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 58, 293-306, doi:10.1111/j.1600-0870.2006.00178.x.
            
            [Birner and Williams(2008)]birner08
            Birner, T., and P. D. Williams, 2008: Sudden Stratospheric Warmings as Noise-Induced Transitions. J. Atmos. Sci., 65, 3337-3343[Buizza et al.(1999)]buizza99
            Buizza, R., Milleer, M., Palmer, T. N. (1999). Stochastic representation of model uncertainties in the ECMWF ensemble prediction system. Quarterly Journal of the Royal Meteorological Society, 125(560), 2887-2908. https://doi.org/10.1002/qj.49712556006
            
            [Carrasi et al.(2017)]carrassi17
            Carrassi, A., Bocquet, M., Hannart, A., Ghil, M. (2017). Estimating model evidence using data assimilation. Quarterly Journal of the Royal Meteorological Society, 143(703), 866-880.
            
            [Chopin et al.(2013)]chopin13
            Chopin, N., Jacob, P. E., Papaspiliopoulos, O. (2013). SMC2: an efficient algorithm for sequential analysis of state space models. J. R. Stat. Soc. Series B Stat. Methodol. , 75(3), 397-426.
            
            [Crommelin and Vanden-Eijnden(2008)]crommelin08
            Crommelin, D., Vanden-Eijnden, E. (2008). Subgrid-Scale Parameterization with Conditional Markov Chains. Journal of the Atmospheric Sciences, 65(2003), 2661-2675. https://doi.org/10.1175/2008JAS2566.1
            
            [Christensen et al.(2015)]christensen15
            Christensen, H. M., Moroz, I. M., Palmer, T. N. (2015). Stochastic and perturbed parameter representations of model uncertainty in convection parameterization. J. Atmos. Sci., 72(6), 2525-2544.
            
            [Danforth et al.(2007)]danforth07
            Danforth, C. M., E. Kalnay, and T. Miyoshi, 2007: Estimating and correcting global weather model error. Bull. Am. Meteorol. Soc., 88, 303-304, [Dee(1995)]dee95
            Dee, D. P., 1995: On-line estimation of error covariance parameters for atmospheric data assimilation. Mon. Weather Rev., 123, 1128-1145.
            
            [Dee and Da Silva(1998)]dee98
            Dee, D. P., Da Silva, A. M. (1998). Data assimilation in the presence of forecast bias. Quarterly Journal of the Royal Meteorological Society, 124(545), 269-295. 
            
            [DelSole and Yang(2010)]delsole10
            DelSole, T., and X. Yang, 2010: State and parameter estimation in stochastic dynamical models. Phys. D, 239, 1781-1788, 
            [Doucet et al.(2000)]doucet00
            Doucet, A., Godsill, S.,  Andrieu, C., 2000: On sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and computing,  10(3), 197-208.
            
            [Dreano et al.(2017)]dreano17
            Dreano, D., P. Tandeo, M. Pulido, T. Chonavel, and I. Hoteit, 2017: Estimating model error covariances in nonlinear state-space models using Kalman smoothing and the expectation-maximisation algorithm. Q. J. R. Meteorol. Soc., 1877-1885, doi:10.1002/qj.3048.
            
            
            [Frei and Künsch(2012)]frei12
            Frei, M.,  Künsch, H. R. (2012). Sequential State and Observation Noise Covariance Estimation Using Combined Ensemble Kalman and Particle Filters. Monthly Weather Review, 140(5), 1476-1495.[Griffith and Nichols(2000)]griffith00
            Griffith, A. K., Nichols, N. K. (2000). Adjoint methods in data assimilation for estimating model error. Flow, turbulence and combustion, 65(3), 469-488.
            
            [Grooms et al.(2015)]grooms15
            Grooms, I., Lee, Y., Majda, A. J. (2015). Ensemble Filtering and Low-Resolution Model Error: Covariance Inflation, Stochastic Parameterization, and Model Numerics. Monthly Weather Review, 143(10), 3912-3924. https://doi.org/10.1175/MWR-D-15-0032.1
            
            [Ha et al.(2015)]ha15
            Ha, S., Berner, J., Snyder, C. (2015). A Comparison of Model Error Representations in Mesoscale Ensemble Data Assimilation. Monthly Weather Review, 143(10), 3893-3911. https://doi.org/10.1175/MWR-D-14-00395.1
            
            [Hamill and Whitaker(2005)]hamill05
            Hamill, T. M., Whitaker, J. S. 2005: Accounting for the Error due to Unresolved Scales in Ensemble Data Assimilation: A Comparison of Different Approaches. Monthly Weather Review, 133(11), 3132–3147. http://doi.org/10.1175/MWR3020.1
            
            [Hamill and Whitaker(2011)]hamill11
            Hamill, T. M., and J. S. Whitaker, 2011: What Constrains Spread Growth in Forecasts Initialized from Ensemble Kalman Filters? Mon. Weather Rev., 139, 117-131, doi:10.1175/2010MWR3246.1. http://journals.ametsoc.org/doi/abs/10.1175/2010MWR3246.1 (Accessed October 5, 2017).
            
            [Harlim and Hunt(2007)]harlim07
            Harlim, J., and B. R. Hunt, 2007: Four-dimensional local ensemble transform Kalman filter: Numerical experiments with a global circulation model. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 59, 731-748[Houtekamer et al.(2009)]houtekamer09
            Houtekamer, P. L., Mitchell, H. L., Deng, X. (2009). Model Error Representation in an Operational Ensemble Kalman Filter. https://doi.org/10.1175/2008MWR2737.1
            
            [Hunt et al.(2007)]hunt07
            Hunt, B. R., E. J. Kostelich, and I. Szunyogh, 2007: Efficient data assimilation for spatio-temporal chaos: A local ensemble
            transform Kalman filter.  Physica D,  77, 437-471, doi:10.1016/j.physd.2006.11.008.
            
            [Jazwinski(2007)]jazwinski07
            Jazwinski, A. H. (2007). Stochastic processes and filtering theory. Courier Corporation.
            
            [Koyama and Watanabe(2010)]koyama10
            Koyama, H., and W. Watanabe, 2010: Reducing forecast
            errors due to model imperfectness using ensemble
            Kalman filtering.  Mon. Wea. Rev.,  138, 3316–3332.
            
            [Leutbecher et al.(2017)]leutbecher17
            Leutbecher, M., Lock, S.J., Ollinaho, P., Lang, S., Balsamo, G., Bechtold, P., and Coauthors, 2017: Stochastic representations of model uncertainties at ECMWF: state of the art and future vision. Q. J. R. Meteorol. Soc., 143, 2315-2339[Li et al.(2009)]li09
            Li, H., E. Kalnay, T. Miyoshi, and C. M. Danforth, 2009: Accounting for Model Errors in Ensemble Data Assimilation. Mon. Weather Rev., 137, 3407-3419, doi:10.1175/2009MWR2766.1. [Lorenz(1996)]lorenz96
            Lorenz, E. N. 1996: Predictability—A problem partly solved. Pp. 1-18 in Proceedings of seminar on predictability: Volume 1. ECMWF, Reading, UK.
            
            [Mitchell and Houtekamer (2000)]mitchell00
            Mitchell, H. L., and P. L. Houtekamer, 2000: An adaptive ensemble Kalman filter. Mon. Weather Rev., 128, 416-433[Mitchell and Carrassi(2015)]mitchell15
            Mitchell, L., Carrassi, A. (2015). Accounting for model error due to unresolved scales within ensemble Kalman filtering. Quarterly Journal of the Royal Meteorological Society, 141(689), 1417-1428. https://doi.org/10.1002/qj.2451
            
            [Miyoshi(2011)]miyoshi11
            Miyoshi, T., 2011: The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter. Mon. Weather Rev., 139, 1519-1535[Navon(1998)]navon98
            Navon, I. M. (1998). Practical and theoretical aspects of adjoint parameter estimation and identifiability in meteorology and oceanography. Dynamics of Atmospheres and Oceans, 27(1-4), 55-79. https://doi.org/10.1016/S0377-0265(97)00032-8
            
            
            [Orrell(2003)]orrell03
            Orrell, D. (2003). Model Error and Predictability over Different Timescales in the Lorenz’96 Systems. Journal of the Atmospheric Sciences, 60(17), 2219-2228.
            
            [Palmer(2001)]palmer01
            Palmer, T. N. 2001: A nonlinear dynamical perspective on model error: A proposal for non-local stochastic-dynamic parameterization in weather and climate prediction models.  Q.J.R. Meteorol. Soc.,  127, 279-304. doi: 10.1002/qj.49712757202
            
            [Palmer et al.(2009)]palmer09
            Palmer, T. N., Buizza, R., Doblas-Reyes, F., Jung, T., Leutbecher, M., Shutts, G. J., and co-authors (2009). Stochastic parametrization and model uncertainty.  ECMWF Tech. Memo, 598, 1-42.
            
            [Palmer and Weisheimer(2011)]palmer11
            Palmer, T. N., and A. Weisheimer, 2011: Diagnosing the causes of bias in climate models - why is it so hard? Geophys. Astrophys. Fluid Dyn., 105, 351-365, [Pulido et al.(2016)]pulido16
            Pulido, M., Scheffler, G., Ruiz, J. J., Lucini, M. M., and Tandeo, P. (2016). Estimation of the functional form of subgrid-scale parametrizations using ensemble-based data assimilation: a simple model experiment. Q. J. R. Meteorol. Soc., 142(701), 2974-2984. https://doi.org/10.1002/qj.2879
            
            [Pulido and Rosso(2017)]pulido17
            Pulido, M. and Rosso, O. A. (2017). Model Selection: Using Information Measures from Ordinal Symbolic Analysis to Select Model Subgrid-Scale Parameterizations. J. Atmos. Sci., 74(10), 3253-3269.
            
            [Pulido et al.(2018)]pulido18
            Pulido, M., Tandeo, P., Bocquet, M., Carrassi, A., Lucini, M. (2018). Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 70(1), 1442099.
            
            [Ruiz et al.(2013a)]ruiz13a
             Ruiz J., M. Pulido and T. Miyoshi, 2013: Estimating parameters with ensemble-based data assimilation. A review.  J. Meteorol. Soc. Japan.  91, 79-99.  
            
            [Ruiz et al.(2013b)]ruiz13b
            Ruiz J., M. Pulido and T. Miyoshi, 2013: Estimating parameters with ensemble-based data assimilation. Parameter covariance treatment.   J. Meteorol. Soc. Japan.  91, 453-469. 
            
            [Ruiz and Pulido(2015)]ruiz15
            Ruiz, J., and M. Pulido, 2015: Parameter Estimation Using Ensemble-Based Data Assimilation in the Presence of Model Error. Mon. Weather Rev., 143, 1568-1582, doi:10.1175/MWR-D-14-00017.1. [Santitissadeekorn and Jones(2015)]santiti15
            Santitissadeekorn, N., and C. Jones, 2015: Two-Stage Filtering for Joint State-Parameter Estimation. Mon. Weather Rev., 143, 2028-2042[Shutts(2005)]shutts05
            Shutts, G., 2005: A kinetic energy backscatter algorithm for use in ensemble prediction systems. Q. J. R. Meteorol. Soc., 131, 3079-3102, doi:10.1256/qj.04.106. [Smith(2001)]smith01
            Smith, L. A. 2001. Disentangling uncertainty and error: on the predictability of nonlinear systems. Pp. 31—64  In Nonlinear dynamics and statistics . Ed. A. I.Mees. Birkhauser, Boston, USA
            
            
            [Stroud and Bengtsson(2007)]stroud07
            Stroud, J. R., and T. Bengtsson, 2007: Sequential State and Variance Estimation within the Ensemble Kalman Filter. Mon. Weather Rev., 135, 3194-3208[Stroud et al.(2017)]stroud17
            Stroud, J. R., Katzfuss, M., and Wikle, C. K. (2017). A Bayesian adaptive ensemble Kalman filter for sequential state and parameter estimation. Monthly Weather Review, (2017).
            
            
            [Wang et al.(2007)]wang07
            Wang, X. G., Snyder, C,. Hamill, T. M. (2007) On the theoretical equivalence of differently proposed ensemble–3DVAR hybrid analysis schemes. Mon. Wea. Rev., 135, 222-227
            
            [Whitaker and Hamill(2012)]whitaker12
            Whitaker, J. S., Hamill, T. M. (2012). Evaluating Methods to Account for System Errors in Ensemble Data Assimilation. Monthly Weather Review, 140(9), 3078-3089. https://doi.org/10.1175/MWR-D-11-00276.1
            
            [Wilks(2005)]wilks05
             Wilks D. S., 2005: Effects of stochastic parameterizations in the Lorenz ’96 system.  Q. J. R. Meteorol. Soc.,  131, 389–407. doi: 10.1256/qj.04.03
             
            [Yang and Delsole(2009)]yang09
            Yang, X., and T. Delsole, 2009: Using the ensemble Kalman filter to estimate multiplicative model parameters. Tellus, Ser. A Dyn. Meteorol. Oceanogr., 61, 601-609
             
            [Zupanski(1997)]zupanski97 
            Zupanski, D., 1997: A General Weak Constraint Applicable to Operational 4DVAR Data Assimilation Systems. Mon. Weather Rev., 125, 2274-2292,
