{'fraction_non_alphanumeric': 0.035401606851376846, 'fraction_numerical': 0.017149640684289744, 'mean_word_length': 4.693904646952324, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Observation Error Covariance Specification in Dynamical Systems for Data assimilation using Recurrent Neural Networks

Data assimilation techniques are widely used to predict complex dynamical systems with uncertainties, based on time-series observation data. Error covariance matrices modelling is an important element in data assimilation algorithms which can considerably impact the forecasting accuracy. The estimation of these covariances, which usually relies on empirical assumptions and physical constraints, is often imprecise and computationally expensive especially for systems of large dimension. In this work, we propose a data-driven approach based on long short term memory (LSTM) recurrent neural networks (RNN) to improve both the accuracy and the efficiency of observation covariance specification in data assimilation for dynamical systems. Learning the covariance matrix from observed/simulated time-series data, the proposed approach does not require any knowledge or assumption about prior error distribution, unlike classical posterior tuning methods. We have compared the novel approach with two state-of-the-art covariance tuning algorithms, namely DI01 and D05, first in a Lorenz dynamical system and then in a 2D shallow water twin experiments framework with different covariance parameterization using ensemble assimilation. This novel method shows significant advantages in observation covariance specification, assimilation accuracy and computational efficiency.

overcome these difficulties, we often rely on calibration (e.g., least square) methods based on some generic correlation kernels, often with homogeneous and isotropic characteristics [14]. Balanced operators can be employed for multivariate systems [15]. In terms of correlation kernels, the family of Matérn functions, including the Exponential kernel (Matérn 1/2), the Balgovind kernel (Matérn 3/2, also known as second order auto-regressive (SOAR) function) and the Gaussian kernel (Matérn 5/2), is often prioritized for covariance computing owing to its smoothness and capability to capture spatial correlations in physical processes [10,16,17]. Other stationary covariance models involve, for instance, convolution formulation [18] or diffusion-based operators [19], both contribute to an efficient storage of the covariance matrices. However, limited by homogeneous and isotropic assumption, it remains cumbersome to represent complex spatial correlation (often multidimensional and multivariate) using these one dimensional kernels.
In this study, we develop and test a novel data-driven approach based on recurrent neural networks (RNN) to improve both the accuracy and the efficiency of observation covariance specification in dynamical data assimilation problems. The novel approach is tested and compared with two state-of-the-art covariance tuning algorithms in two different digital experiments with parametric and non-parametric covariance estimation respectively.
The paper is organized as follows. In Sect. 2, we introduce the related work for error covariance specification. The problem statement and the contribution of this paper is described in Sect. 3. Data assimilation techniques and the ensemble methods are introduced briefly in Sect. 4. We then describe traditional posterior covariance tuning algorithms DI01 and D05 in Sect. 5. The novel LSTM-based method is introduced in Sect. 6, followed by the comparison in the Lorenz (Sect. 7) and the shallow water twin experiments (Sect. 8). We close the paper with a discussion in Sect. 9.

Related work
To gain a clearer insight into covariance evolution, some ensemble-based methods such as [1] (NMC) and [20] (EnKF), have been developed to provide a non-parametric covariance estimation. These methods depend on the propagation of an ensemble of simulated trajectories, initialized either at different forecasting time steps (NMC) or by adding some artificially set perturbations to the current state (EnKF). These methods are more appropriate for modelling the background matrix compared rather than the observation matrix. The latter, independent from the numerical simulations, can not be represented by the propagation of artificially added noises. The Particle-Aided Unscented Kalman Filter [21] can estimate systems with high nonlinearity with a real-time updating of the background matrix. However, the observation matrix can not be estimated directly via the Particle-Aided Unscented Kalman Filter. In practice, the observation matrix is often set to be diagonal or spatially isotropic for the sake of simplicity (e.g. [22]). However, it is shown in the work of [10] that well-specified correlated observation covariances can significantly improve the performance of DA algorithms.
Several methods of data-derived posterior diagnosis have also been developed based on the analysis of innovation quantities which consist of the difference between the observations and the projected background/analyzed state in the observation space. As a strong contributor to this topic, the meteorology community developed several well-known posterior diagnoses and their improved versions [23,24,25] to adjust the background/observation ratio, the correlation scale length or the full covariance structure in the observation space (both the observation matrix and the projected background matrix). Some iterative processes [26,27] based on the fixedpoint theory have also been proposed for error covariance tuning. Recent works of [28] and [29] have proved the convergence of so-called "Desroziers iterative methods" [24] (also known as D05) in the ideal case. In brief, they have mathematically proved that, by using a semi positive definite matrix as an initial guess, D05 iterative method converges on the exact * time-invariant (at least over a sufficiently long time period) observation error covariance when the background matrix and the transformation operator (which maps the state variables to real-time observations) are perfectly known a priori. On the other hand, it is also mentioned by [29] that a regularization step is necessary in practice for applying D05 and the convergence of the regularized iterations remains an open question [29,3]. To deal with time-varying systems, lag-innovation statistics are used for error covariance estimation [30]. The essential idea is to build a secondary Kalman-filtering process for adjusting error covariances using time-shifted innovation vectors. For more details of the innovation-based methods, we refer to the overview of [13] which also covers some other estimation methods, such as the family of likelihood-based approaches and expectation-maximum(EM) methods. * Here, by the term "exact", we refer to the covariance truly corresponding to the remaining errors present in the observation space

Problem statement and contribution
Our work lies in a similar condition of [24] and [29] where both the state forward model and the transformation operator are presumed to be well known. As the main difficulty concerns the non-synchronous time-variant observations in dynamical systems (which prevents empirical estimation), in this work we propose the use of recurrent neural networks (RNN) [31] for the specification of the observation matrix across the underlying dynamics of the observed quantities. RNN has been widely adapted for the prediction/reconstruction of dynamical systems, especially in natural language processing (NLP) [32] and image/video processing [33] due to its convincing capacity of dealing with time series. More recently, RNN has also made their way to other engineering fields such as biomedical applications and computational fluid mechanics [34]. In general, the combination of deep learning and data assimilation methods [35,36] has been widely adapted and analyzed in a variety of industrial applications, including air pollution [37] and ocean-atmosphere modelling [38]. A convolutional neural network (CNN) for covariance estimation has also been suggested in the work of [39]. In this study, we propose a novel methodology for LSTM-based covariance estimation which can be easily integrated into any DA schema for dynamical systems. Here, we first construct a set of training covariance matrices, being either parametric or non-parametric, within a certain range defined a priori. For each matrix in the training set, we then simulate a dynamic trajectory of the observation vector relying on the knowledge of the forward model where the noises at each time step are generated following a centered Gaussian distribution characterized by the error covariance. These trajectories are later used as input variables to train the long-short-term-memory (LSTM) RNN regression model where the time invariant observation matrices stand for the learning target. For the online evaluation, only the historical observation data is needed to predict the error covariances. Compared to traditional posterior tuning methods [24,40] which require several implementations of DA algorithms, the proposed machine learning (ML) method can be much more computationally efficient for real-time covariance estimation. Moreover, no prior knowledge concerning either the background or the observation matrix is necessary for the proposed ML approach unlike most of the traditional methods. For example, DI01 [23] requires precise knowledge of correlation structures for both background and observation matrices while D05 [24] make use of the perfect knowledge of the background covariance.
In order to make a comprehensive comparison with traditional methods, two different twin experiment frameworks are implemented in this paper, using respectively the Lorenz96 and the 2D shallow water models. The Lorenz system, characterized by only three state variables, is associated with a non-parametric covariance modelling while we use an isotropic correlation kernel to parameterize the observation matrix in the shallow water dynamics. In both cases, we compare the performance of the proposed LSTM-based method against the state-of-the-art tuning approaches D05 and DI01 in terms of both the covariance specification and the posterior DA accuracy. An ensemble DA schema is used for estimating the time variant background matrix for each of these methods.

Principle of data assimilation
The objective of data assimilation algorithms is to approach the estimation of system's states x to its true values x true , also known as the true state, by taking advantage of two sources of information: the prior estimation or forecast x b , which is also called the background state, and the measurement or observation y. DA algorithms aim to find an optimally weighted compromise between x b and y by minimizing the lost function J defined as: where H denotes the transformation operator from the state space to observation space. B and R are respectively the background and the observation error covariance matrices, i.e.
Errors b , y are supposed to be centered Gaussian following: In Eq. 1, the left side strives for incorporating the prior information x b , and the right side penalizes the difference between the observation y and the state variables after having been mapped to the observation space H(x). Both terms are weighted by the corresponding inverse of error covariance matrix (B −1 , R −1 ) to reflect confidences for each of them.
The optimization problem of Eq. 1, so called three-dimensional variational (3D-Var) formulation, is a general representation of variational assimilation which does not take into account model errors. The output of Eq. 1 is denoted as x a , i.e.
x a = argmin If H is the linear observation operator represented by H, Eq. 6 can be solved via BLUE (Best Linearized Unbiased Estimator) formulation: where A = Cov(x a − x true ) is the analyzed error covariance, and K is the Kalman gain matrix described by In the rest of this paper, we define H as a linear transformation operator. Nevertheless, it is usually more challenging to find the minimum of Eq. 1 when H is non-linear, even more challenging when states are highdimensional. The resolution for the minimization often involves gradient descent algorithms (such as "L-BFGS-B" [41] or adjoint-based [42] numerical techniques. DA algorithms could be applied to dynamical systems thanks to sequential applications expressed by a transition operator M t k →t k+1 (from discrete time t k to t k+1 ), where x b,t k+1 thus depends on the knowledge of M t k →t k+1 and the DA correcting state x a,t k , i.e., Obviously, the more accurate x a,t k is, the more reliable x b,t k+1 would be.
To leverage the information embedded in the background state and observations, covariance matrices modeling is a pivotal point in DA as they influence not only how prior errors spread but may also change the DA results [26].

Ensemble methods
Ensemble data assimilation (EnDA) [43,44] methods have shown a strong performance in dealing with non-linear chaotic DA systems by creating an ensemble with size M of the system state depicted as {x The latter is used to represent both the prior and the posterior probability distribution of the state variables. The system states of the ensemble evolve under M t k →t k+1 and DA is applied to each of these ensemble states at every assimilation windows. Furthermore, instead of evolving the system to obtain the B matrix, which is a time and computationally expensive process when a large number of states is available, B is estimated as a sample covariance: b,t k , and the estimation becomes more reliable with the increases of M . For applications in this study, EnDA, with a sufficiently large number of examples, is used to estimate x b,t k and B b,t k so that we can focus on the comparison of R matrix modelings.

Observation error covariances specification
For the estimation of R, under the assumption that the system model is stationary, a wide variety of methods have been explored, for example, the DI01 [23] method which adjusts accordingly the ratio between T r(B) and T r(R) and the D05 [24] approach which estimates the full observation space iteratively. However, these methods, based on posterior innovation quantities (i.e., y − H(x a )) which requires several applications of DA algorithms, can be computationally expensive. Moreover, these tuning methods, especially the D05 which estimates the full matrix, are not suitable for different matrix parameterizations. In this paper, working with time-series observation data, we use LSTM to predict the corresponding R matrix under similar assumptions of DI01 and D05. The two classical methods, introduced in Sect. 5, are implemented to compare the results with the proposed machine learning approach.

Desroziers & Ivanov (DI01) tuning algorithm
Because B and R determine the weight of background and observation information in the loss function ( Eq. 1), the knowledge of Tr(B) and Tr(R) is crucial to DA accuracy. DI01 [23] tuning algorithm, relying on the diagnosis of innovation quantities, has been widely adopted in meteorology [45,28] and geoscience [46]. Consecutive works have been carried out to improve its performance and feasibility in problems of large dimension [47]. Without modifying error correlation structures, DI01 adjusts the prior error amplitudes by applying an iterative fixedpoint procedure.
As demonstrated in [48,23], when B and R are perfectly specified, where H is a linearized observation operator. Based on Eq. 13 and 14 it is possible to iteratively correct the magnitudes of B and R, following using the two indicators where q is the current iteration. Acting as scaling coefficients, the sequences {s b,q } and {s o,q } modify the error variance magnitude in the iterative process. It is worth reminding that both the analyzed state x a and the gain matrix K q are obtained using B q , R q which depend on s b,q and s o,q . When the correlation patterns of both B and R are well known, DI01 is equivalent to a maximum-likelihood parameter tuning, as pointed out in [47,28].
Unlike other posterior covariance diagnosis/computations, such as [24,26], the estimation of the full matrices is not needed in DI01. Instead, only the estimation of two scalar values (J b , J o ) is required, which significantly reduce the computational cost. As a consequence, this method could be more appropriate for online covariance tuning.

Desroziers iterative method (D05) in the observation space
The Desroziers diagnosis (D05) [24], subject to prior and posterior state-observation residuals has been widely applied in engineering problems, including numerical weather prediction [28] and hydrology [3]. The work of [24] proved that when B and R are well known a priori, the expectation of the analysis state should satisfy: The difference between the left side and the right side of Eq.18, can be used as a validation indicator of the R matrix estimation where ||.|| F denotes the Frobenius norm. Applying this method, time variant observation/background data can contribute to the estimation of the R matrix because the expectation in Eq. 18 could be evaluated using residuals at different time steps. When the B matrix is well known, an iterative process has been introduced to estimate the R matrix: based on the fixed-point theory [29]. The current analysis state x a,q is obtained using the specification of R q while x b , B, y remains invariant. As proved in [28] and [29], under the assumption of sufficient observation data and well known B matrix, the iterative process of Eq.20 converges to the exact observation error covariance. However, as shown in [29], the intermediate matrices R q could be non-symmetric and possibly contain negative or complex eigenvalues, which is cumbersome for DA algorithms to deal with.
In practice a posterior regularization at each iteration step is often required to ensure the positive definiteness of R q [29] where the first step of the regularization could be symmetrizing the estimated R q matrix, i.e., The spectrum of R q now contains only real numbers but they are not necessarily positive. The hybrid method [2] is a standard approach in ensemble-based DA methods to ensure the positive definiteness, which consists of combining a prior defined covariance matrix C with the one obtained from empirical estimation. We thus obtain the formulation of the regularized observation matrix R r,n : following Eq.21 with µ ∈ (0, 1). The matrix C is often set as a diagonal matrix since it helps to enhance the matrix conditioning. In this work, we choose to set so that T r(R q ) will not be modified due to the post regularization. As mentioned in the discussion of [29,3], the convergence of regularized observation matrices remains an open question. Therefore, a small iteration number is often assigned for D05 tuning in industrial problems (e.g., q = 2 in [49] and [3]). Since the right side of Eq. 20 can be estimated using residual quantities at different time steps, D05 is often used to deal with time series observation data (e.g., [49,3]) when assuming the R matrix is time-invariant.

Introduction of RNN and LSTM
LSTM, first introduced in [50], is a kind of RNN [31], capable of solving long term dependency problems [51] that traditional RNN could not handle. As with other recurrent neural networks, LSTM has a chain-like structure. This structure is created by repeating the same module shown on the left side in Fig. 1. In addition, the repeating module comprises four neural networks instead of only one. The specific structure of the repeating module is on the right side in Fig. 1. An essentail part of LSTM is the cell state C t k−1 which is the long-term memory storing information about past behaviors. LSTM uses three gates with each composed out of a sigmoid layer neural network (single layer neural network with sigmoid activation function at the output layer) and a pointwise multiplication operation, to protect and control information of the cell state as shown in Fig. 1.
The first gate is the forget gate following: where the recurrent variable h t k−1 summarizing all the information about past behaviors, x t k resuming information about current behaviors, and W f and b f are weights and bias respectively parameterizing the sigma layer neural network. The forget gate decides what kind of information is going to be ignored in C t k−1 .
The second gate is the input gate, and it determines which new information is added into C t k−1 . This new informationC t k , as conforming toC C t k is then multiplied by weight coefficients i t k which is acquired by the input gate, i.e., i t k is applied to decide which new information would be employed to update C t k−1 .
Hence, the state cell C t k at the current time step t k can be attained using Finally, the acquisition of h t k requires the participation of the output gate and a tanh activation function: first, the tanh activation function tanh is used to create a cell state candidate information tanh(C t k ). tanh(C t k ) is then multiplied by some weight coefficients following to decide which information of tanh(C t k ) would contribute to the obtainment of h t k . Among them, o t k is generated by the output gate with neural network parameters W o and b o , i.e., 6.2 LSTM for R matrix estimation using time series observation data The tuning methods presented in Sect. 5 have been applied in various engineering applications with significant improvement of covariance specification and DA accuracy. However, these methods which require several applications of DA algorithms can be computationally expensive for high-dimensional problems. Another important drawback stands for the requirement of precise knowledge on either the correlation patterns of B and R (for DI01) or the full B matrix (for D05). In this study, we aim to build a data-driven surrogate model for efficient online R matrix specification using LSTM. Unlike DI01 or D05, no specific knowledge about the error covariances or the state/observation dynamical systems other than the transformation operator H and the forward model M t k →t k+1 (which is also indispensable for standard DA algorithms including variational methods and Kalman-filters) is required.
Based on an initial state x t k } as the system input and the R matrices as output, LSTM networks are used to learn the error distribution across the underlying observation dynamics. More precisely, a real function g R (.) : Φ R −→ R m×m is predefined where Φ R is an empirically estimated real space which defines the range of a set of parameters, such as marginal error variance, correlation scale length [18] for computing the R matrices. The generated observation matrices {R [iter] } are set to be symmetric positive definite (SPD) thanks to the function g R (.). Both g R (.) and Φ R vary for different applications. Generated states {x [iter] g,t k }, t k ∈ {0, · · · t T } with t T depicting the index of the final time step, could be attained by evolving the system respecting M t k →t k+1 : The observation {y [iter] t k }, t k ∈ {0, · · · t T }, is gained by mapping {x [iter] g,t k } through H and at the same time combining with random Gaussian noises: Algorithm 1: Simulated observations y generating process  t k }. The general process can be demonstrated in Fig. 2 where for each application, y [iter] is simulated by evolving the system knowing the system dynamics and the observation error covariance generator, and then y [iter] and R [iter] are applied to train LSTM so that LSTM can predict R [iter] when receiving unseen y [iter] .
Following the principle of many to one LSTM in Fig. 3, the input features of LSTM consist of observations t T }, while the output is R [iter] matrix. Different from classical covariance tuning algorithms, the LSTM network only makes use of historical observation data, requiring neither the Fig. 3: Many to one LSTM training and R prediction process background states nor the error covariance matrix. The advantage of using LSTM is more salient when the observation dimension is large, for example, millions or even billions, while such dimension is not uncommon in real-world applications [7,2].
To estimate R, LSTM is first trained to learn: -either related variables which can be used to constitute the symmetry observation error covariance matrix (i.e., input variables of the g R (.) function) in a parametric modeling; -or elements of the R matrix (e.g., variables in the upper triangle and those in the diagonal of the covariance matrix) in a non-parametric modeling. The whole process for R estimation using LSTM is described in Algorithm 2.  In Algorithm 2, it is suggested that LSTM training process is consisted of training and validation processes: the training process is comprised of the forward prediction and the backward neural network weight parameters Fig. 4: Online LSTM prediction of the observation matrix compared to traditional covariance tuning approaches, followed by online data assimilation updating processes; and validation process is used to predict desirable outputs or objectives and then calculate validation loss between predicted output and prior true output values. N epoch indicates the number of times that the entire example data set is passed forward and backward through the LSTM during the training process. Early stopping, which terminates the training process when the validation loss reaches the minimum and is always the minimum value after N patience epoch epochs, is applied to reduce the LSTM training time.
It is important to note that the offline data generation and LSTM training processes need to be carried out individually for different DA applications.

Twin experiment principle
In order to overcome the drawback that, in a realistic experiment, x true is usually unknown and y is often mixed with noises, twin experiment, in which a prototypical test case is selected to simulate real situations, is applied so as to provide x true for comparison. In this experiment, a mapping is applied to some sampling true trajectory x true,t k at some points in space and time and arbitrary random noises are added to obtain simulated raw measurements y t k . DA is then implemented starting from initial background state x b,t 0 representing the prior information that could be obtained about corresponding state x true,t 0 , along with initial raw measurement y t 0 . The output state is then compared against x true , verifying the distance of these two states and minimizing it to evaluate and improve the performance of DA. In this section, we use a twin experiment to evaluate the performance of applying DA to a simple Lorenz system in which raw measurement error covariance is estimated/adjusted using respectively DI01, D05 and LSTM.

Experiment set up
The Lorenz system, first studied by Edward Lorenz, is a system of ordinary differential equations. For certain parameter values and initial conditions, the Lorenz system is notable for having chaotic solutions, in particular the Lorenz attractor, toward which a system tends to evolve. The Lorenz 96 system [52] has been widely used as a prototypical test case to compare the performance of DA algorithms [34,53,35]. Here we build a twin experiment framework with a simple three dimensional Lorenz system in which the state vector is denoted as (2) ]. The studied Lorenz system can be characterized as: (2) .
The initial values of the true state x true,t 0 are set to be [0, 1, 1.05] while the initial background state x b,t 0 is generated by combining x true,t 0 with a centered Gaussian noise b,t 0 : Then both of true states x true = {x true,t 0 , · · · , x true,t T } and background states x b,t k = {x b,t 0 , · · · , x b,t T } of the Lorenz system evolve by conforming respectively to the Lorenz equation in Eq. 32 until t=1s with total T = 1000 time steps.  Subsequently, observations y = {y t 0 , · · · , y t T } can be acquired by mapping x true through a linear observation operator and adding noises respecting multinormal distribution N (0, R) where R is randomly generated following the process described in Sect. 7.3. EnDA is then applied in this twin experiment to update the background ensemble using available observations. More precisely, every 10 time steps, EnDA is applied on x b,t k along with y t k to obtain the analysis states x a,t k . To simulate future background states before the next assimilation step, we add some artificial noises b,t k+1 to x b,t k+1 respecting multivariate distributions, i.e., The model error covariance matrix Q is supposed to be time-invariant for all generated trajectories.

DA with LSTM-based covariance estimation
The Sect. 7.2 exhibits the process of generating artificial training data for the LSTM model. More details about the R matrix generation, the outputs of the LSTM model within the training data are revealed in this section. R are parameterized by four real coefficients r 0 , r 1 , r 2 and v R which determine respectively the three correlation coefficients and the error amplitude, i.e., In this study, v R is generated uniformly between 0 and 100, i.e., v R ∼ U(0, 100) with U denoting the uniform probability distribution. The correlation coefficients {r 0 , r 1 , r 2 } ∈ (−1, 1) 3 are obtained via randomly generated SPD matrices † . The LSTM network is thus trained to learn the R matrix by trying to build a function mapping time series observation data y t k to r 0 , r 1 , r 2 and v R . The specific structure of the LSTM model which consists of a LSTM input layer, a hidden layers with 200 neurons, and an output layer comprising four neurons applied to obtain R, is shown in Table 1. In this Lorenz twin experiment, two LSTM networks, respectively named as LSTM1000 and LSTM200 are designed with different input size of times series data. LSTM1000 is trained on a total of 1000 time steps for predicting the R matrix while LSTM200 makes use of only the first 200 time steps to simulate a realistic application where the time-invariant R matrix is estimated using historical data for improving future DA performance. The evaluation of both LSTM1000 and LSTM200, in terms of DA accuracy, is made using the full test dataset with 1000 times steps. By leveraging the LSTM model along with available observations, we can still perform DA algorithms even though R is not explicitly given. The results are then compared with the ones obtained using the exact R matrix. † https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make spd matrix.html

Results
To evaluate the LSTM performance of predicting R, we first compare the LSTM predicted R and the predefined true R matrix in the test set and analyze the impact of predicted R on DA accuracy. Since the LSTM prediction of R matrix is non-parametric in these twin experiments, we compare each element of the predicted matrix (i.e.,r 0 , r 1 , r 2 and v R ) against the ground truth. As for the DA accuracy, we calculate the difference between (2),t } (obtained via Eq. 35) refined in DA using R estimated via DI01 (q = 2), D05 (q = 3) or LSTM and the error-free true states x true,t = {x true,(0),t , x true,(1),t , x true,(2),t }. Both D05 and DI01 initialize with a random R matrix, generated through g R (.) with the same range of parameters Φ R as for the LSTM training. Fig. 5 shows the elements of the estimated R matrix (i.e.,r 0 , r 1 , r 2 and v R ) obtained by LSTM1000 and LSTM200, both trained on 103486 Lorenz system observation samples, and predicts on a test dataset of 10000 samples. The result of the D05 approach, with a calibration of r 0 , r 1 , r 2 and v R against the true value, has also been displayed in Fig. 5 (i-l) in the same test dataset for comparison purposes. The blue circles are the LSTM/D05 prediction results, while the red line is the true value of the corresponding parameter. We observe that both LSTM1000 and LSTM200 prediction results fit very well to the predefined true value of each parameter, compared to the D05 tuning approach, especially for the error amplitude v R which is of most importance in covariance specification. Based on these experimental results, we can conclude that the proposed LSTM approach is capable of predicting the observation matrix, in terms of both correlation coefficients and error amplitude, when time series observation data y t k are given. Fig. 6 and Table 2 illustrate the averaged DA performance along with R attained in different ways with 10000 observations in the test dataset. We remind that for each observation sample {y t k }, t k ∈ {0, · · · , t T }, 100 background trajectories are generated to perform EnDA. Among these algorithms, DI01 uses two iterations to correct the magnitudes of B and R by conforming to Eq.15. Each D05 iteration calculates the innovation quantities (y − H(x a )) and (y − H(x b )) every 10 time steps through a DA procedure, and then applies Eq.20 to attain the updated R for the corresponding Lorenz system. Fig. 6 displays the evolution of the mean square error (MSE) std mse,(i),t between DA refined x b,t (following Eq. 35) and the true states x true,t , i.e.,   Table 2: DA performance of the Lorenz system evaluated in mse for the three state variables based on R varied by the predefinition, the algorithm LSTM1000, LSTM200, DI01 and D05 where i ∈ {0, 1, 2}, N = 1988 is the number of examples, and M = 100 is the size of the ensemble DA. Table 2 shows the averaged (against time) absolute error to interpret the difference between assimilated and true states, where t T = 1000 is the total time steps that the Lorenz system has evolved. It should be noted that the DI01 approach, which exclusively adjusts the T r(B)/T r(R) ratio without modifying the correlation structure, is over-performed by more refined covariance tuning/specification methods such as D05 and LSTM as shown in Fig. 6. Furthermore, Table 2 shows that lstm1000: std mse is smaller than d05: std mse , which is consistent with the results shown in Fig. 5. Thus, we can conclude that LSTM1000 is sound at predicting R, contributing to a better DA performance, compared to DI01 and D05. Such a conclusion is further proved in Table 2 in which mse based on LSTM1000 has values, for all three parameters displayed, close to those based on true R, which is predefined and used throughout the Lorenz system sample generations.
What we have not exhibited in Fig.6 are the curves of lstm200: std mse and true: std mse , representing the MSE based on LSTM200 predicted R and manually predefined true R respectively, as they have been almost totally overlapped with lstm1000: std mse . Such a fact permits both the same conclusion to be reached with LSTM1000, and that LSTM200 which makes use of the observation data of only the first 200 time steps, is sound at predicting R to contribute to a good DA performance with future observations. This conclusion is supported in Table 2 where mse based on LSTM200 predicted R has almost the same values as that based upon manually predefined true R.

Experiment set up
For further evaluating the performance of error covariance estimation using LSTM when incorporated with predefined correlation kernels, we also set up a twin experiment framework with a simplified 2D shallow water dynamical model which is frequently used for testing data assimilation algorithms ( e.g., [42,26]). A cylinder of water is positioned in the middle of the study field with size 20mm × 20mm and released at the initial time step t k = t 0 s (i.e., with no initial speed), leading to a non-linear wave-propagation. The dynamics of the water level h (in mm), as well as the horizontal and vertical velocity field (respectively denoted as u and v in 0.1m/s), is given by the non-conservative shallow water equations, In Eq.39, b = 0.1 is the viscous drag coefficient, while g is the constant earth gravity. These equations are discretized in a 20 × 20 regular grid, solved by first-order finite difference method with a time discretization δ t = 10 −4 s, from t 0 = 0s to t 1000 = 0.1s. This resolution is considered as the reference (i.e. the true state x true ) when performing DA algorithms. The state variables in this DA modeling are the combination of the velocity fields {u} 20×20 and {v} 20×20 , leading to the state dimension of 800. The evolution of the reference state (x true,t k ), together with the error-free model equivalent observations (i.e., H(x true )), is illustrated in Fig. 8. Spatially correlated prior observation errors are generated artificially and combined with the transformation operator to simulate real-time observations. More precisely, the observations are generated from the model equivalent H(x true ) separately for the fields u and v. H is defined as a sparse matrix to imply the fact that measurements in real-world applications are sparser than true states due to the interference existing in the former situations as well as the limited performances of sensors. As shown in Fig. 7, the spatial observations at time t k is defined as the average of u t k and v t k in a 2 × 2 cells area with an observation error y t k , y u,i,j,t k = 1 4 (u true,2i,2j,t k + u true,2i+1,2j,t k + u true,2i,2j+1,t k + u true,2i+1,2j+1,t k ) + y u,i,j,t k , and identical for y v,i,j,t k . Therefore, the dimension of the observation vector y = [y u,t k , y v,t k ] is 200. In this experiment, we suppose that the observation error y u,i,j,t k and y v,i,j,t k , respectively of the velocity fields u and v, follow the same Gaussian distribution N (0, R). Thus, the observation error covariance in this shallow water system can be fully characterized by a 100 × 100 R matrix after the observations (originally in a 2D grid) being converted to a 1D vector. Here we adopt a different parameterization of the R matrix thanks to an isotropic correlation function ψ R (.), where D = D 0 , ..., D 99 , representing the error variances in the 2D (10 × 10) velocity field. Each element of D is generated individually following an uniform distribution, which produces only positive elements to guarantee the positive definiteness of R. ψ R (.) is the second-order auto-aggressive (also known as Balgovind) function,  where L R is the correlation scale length, fixed as L R = 10 in this application. r denotes the correlation scale length in the 2D space and is also generated uniformly with r ∼ U (1,5). Being part of Matern kernels, the SOAR function is often used in DA for prior error covariance modeling [26,6] thanks to its smoothness and good conditioning. The simulation of via the same discretization of Eq.39 (except the initial conditions) is used as background states at time t k in the DA modeling. Similar to the Lorenz experiment(i.e., Eq. 35), {x b,t k } is acquired by combining {x a,t k } with randomly generated Gaussian errors, while {x a,t k } is obtained every 100 time steps (i.e., 0.01s) from ensemble DA with time series observation data {y t k } and the estimated observation error covariance R.

DA with LSTM estimated R
As with the Lorenz experiment, simulated observations {y t k }, generated in the same process with that in the Lorenz system, are used as input training data for the LSTM model, while the D vector and the correlation scale r served as training output. The specific structure of this LSTM network is shown in Table 3. This model has the same structure as the one of the Lorenz system shown in Table 1, except the input and ouput dimension. Besides, two types of LSTM are also proposed as what have been realized in Lorenz system: LSTM1000 employs the whole 1000 time steps observation data as LSTM training and prediction inputs while LSTM200 makes use of only the first 200 time steps of observation data as the LSTM inputs. After the LSTM is well trained on the training set of 173000 generated observation trajectories in this experiment, R can be gained even only observation time series data {y t k },is acknowledged. Similar to the Lorenz system, EnDA is performed here with an ensemble of 100 state trajectories initialized from the same initial state x t 0 for each observation series. EnDA takes place every 200 time steps with the R matrix estimated through different methods. Fig. 9 and Fig. 10 illustrate respectively the predictions results of LSTM1000 and LSTM200 against the true value on 10000 test examples in the test dataset, which demonstrate that the trained LSTM exhibits a good performance at predicting related values applied to compose R, including both marginal error variances (i.e., the elements of the D vector) and the correlation scale length r. The prediction results of LSTM200 (Fig. 10) are almost as accurate as the ones obtained via LSTM1000 (Fig. 9). Training and evaluating the LSTM network using the first 200 time steps is sufficient to obtain an accurate estimation of the R matrix. The so composed R matrices, based on the prediction results from LSTM1000, are shown in Fig. 11 in comparison with the true R matrix and that obtained via D05( q = 3) after regularization. In order to estimate the high dimensional R matrix in this application, D05 makes use of the DA residuals every 10 time steps, which is different from the final DA algorithm, in each of the three iterations. Four different sample R matrices in the test dataset are displayed in Fig. 11. We observe that both LSTM and D05 manage to acquire a similar covariance structure as the true R matrix while the considerable advantage of the LSTM approach can still be noticed. These results confirm the finding of Fig. 9 and Fig. 10 that LSTM can perform well in the parametric prediction of the R matrix even in high dimensional systems. To assess the performance of DA accuracy, the metric mse,(i) (cf., Eq.38) estimated using a set of 53 observation time series, is displayed in Table 4. Since it would take too much space to display the total mse,(i) i=0,1,··· ,800 of all state variables, we select only some of them demonstrated in Table 4, together with the averaged value¯ mse of all 800 cell coordinates in the field of u and v. The displayed results further prove what we have concluded from Fig. 9 to Fig. 11, that is, the LSTM-based approaches own an advantage in terms of the DA accuracy compared to D05 (q = 3) and DI01 (q = 2) tuning algorithms . Furthermore, the performance of LSTM200 is very close to LSTM1000 with an even slightly smaller average MSE, probably due to the sampling randomness. This result further confirms our analysis of Fig. 10 Table 4: DA performance of the shallow water system evaluated in mse based on R varied by the predefinition, the algorithm LSTM1000, LSTM200, DI01 and D05. The execution time only consider the process of R matrix specification without the DA procedure afterwards computational time (of a laptop CPU) of online covariance tuning/specification is also shown in Table 4 where only the evaluation time of the trained LSTM model is taken into account. The training of LSTM can be totally performed offline. As for D05 and D01, we exclude the final DA step (i.e., the computational time is estimated for q − 1 iterations). As shown in Table 4, the LSTM approach is also considerably faster than traditional tuning methods, which allows a near-real-time application in dynamical systems.

Discussion
The precision of DA reconstruction/prediction depends heavily on the specification of both the background and the observation error correlation. The latter is often challenging to estimate in real-world applications because of the dynamic nature of the observation data. Furthermore, the observation matrix R can not be empirically estimated from an ensemble of simulated trajectories, unlike the background error covariance. In this paper, we review in detail some well-known observation covariance tuning algorithms [23,54], based on time-variant posterior innovation quantities. These methods, being widely adopted in geoscience, rely on some specific prior assumptions such as knowledge of the correlation structure [23] or the background matrix [24]. This is difficult to fulfill in some domains where very little knowledge about the prior error is available.
In this study we have proposed a novel machine learning approach based on LSTM neural networks to predict the R matrix using time series observation data as model input. Similar to the work of [23] and [24], R is assumed to be time-invariant, at least over a sufficiently long time period. Both the Kalman-and variationaltype assimilation methods can benefit from the method proposed in this paper for improving the assimilation accuracy. The proposed data-driven approach also contributes to tackling one of the major bottlenecks of DA: it is time consuming and computationally expensive to update covariance matrix, by mapping raw sensor observations to observation error covariance matrix. In both the Lorenz96 and the shallow water models presented in this paper, the LSTM-based approach displays significant strength, compared to classical posterior tuning methods DI01 and D05, in terms of: (i)estimation accuracy of the observation covariance R; (ii) reconstruction and prediction accuracy of the DA schema using the estimated R matrix; (iii) computational efficiency of the online covariance estimation; (iv) flexibility of different model parameterization. It is worth mentioning that an important limitation of the proposed LSTM-based method is the specification of Φ R which defines the range of parameters for training.
Since we assume that the observation matrix is time-invariant, the proposed approach could only deal with fixed sensor placement for dynamical systems, which is also the case of DI01 and D05 tuning algorithms. The possibility of time-variant sensor placement warrants further investigation. As pointed out by [55], DL model can be stolen or reverse engineered by model inversion or model extraction attack. Despite the fact that all data used in the current study is generated from toy models, it is important to ensure the data privacy when applying the model to real applications. Future research should also consider applying the new method to a broader range of real-world problems, including NWP, hydrology, and object tracking, where the offline data simulation could be more computationally expensive compared to the two test models presented in this paper. To this end, future studies could also investigate the combination of model reduction methods, such as domain localization [56], proper orthogonal decomposition, information-based data compression [57], auto-encoder neural networks [58], and the current covariance estimation method. More precisely, the data assimilation can be performed in the compressed low dimensional space (e.g., obtained from POD or auto-encoder). The LSTM-based covariance specification algorithm developed in this work can be used to estimate the observation error covariance matrices in the low dimensional space for improving the accuracy of reduced-order data assimilation approaches.

Code availability
Code for the proposed LSTM-based covariance specification, together with DI01, D05 methods, for both Lorenz and shallow water models is available at https://github.com/scheng1992/LSTM_Covariance