{'fraction_non_alphanumeric': 0.030099904974702725, 'fraction_numerical': 0.010324370136374142, 'mean_word_length': 4.743067846607669, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Equation-free surrogate modeling of geophysical flows at the intersection of machine learning and data assimilation

There is a growing interest in developing data-driven reduced-order models for atmospheric and oceanic flows that are trained on data obtained either from high-resolution simulations or satellite observations. The data-driven models are non-intrusive in nature and offer significant computational savings compared to large-scale numerical models. These low-dimensional models can be utilized to reduce the computational burden of generating forecasts and estimating model uncertainty without losing the key information needed for data assimilation to produce accurate state estimates. This paper aims at exploring an equation-free surrogate modeling approach at the intersection of machine learning and data assimilation in Earth system modeling. With this objective, we introduce an end-to-end non-intrusive reduced-order modeling (NIROM) framework equipped with contributions in modal decomposition, time series prediction, optimal sensor placement, and sequential data assimilation. Specifically, we use proper orthogonal decomposition (POD) to identify the dominant structures of the flow, and a long short-term memory network to model the dynamics of the POD modes. The NIROM is integrated within the deterministic ensemble Kalman filter (DEnKF) to incorporate sparse and noisy observations at optimal sensor locations obtained through QR pivoting. The feasibility and the benefit of the proposed framework are demonstrated for the NOAA Optimum Interpolation Sea Surface Temperature (SST) V2 dataset. Our results indicate that the NIROM is stable for long-term forecasting and can model dynamics of SST with a reasonable level of accuracy. Furthermore, the prediction accuracy of the NIROM gets improved by one order of magnitude by the DEnKF algorithm. This work provides a way forward toward transitioning these methods to fuse information from Earth system models and observations to achieve accurate forecasts.


Introduction
The integration of models and observations has greatly transformed the predictions of the Earth system, including weather forecasts and climate projections [1,2]. Observations can be either used to estimate the initial condition for the prediction that is consistent with the present state of the Earth system through the process of data assimilation or to reduce model errors by advancing the representation of certain processes within a model. Data assimilation (DA) is a well-established method that involves combining information coming from the forecast model with available observations and are used extensively in numerical weather prediction (NWP) [3][4][5]. DA can also be viewed from the Bayesian perspective that involves fusing data (observations) with the prior knowledge (i.e., mathematical representation of physical processes; model output) to obtain an estimate of the distribution of the true state of the process [6,7].
DA methods are usually classified into two types, variational approaches and sequential approaches [8]. Variational DA is formulated as the constrained optimization problem defined by a cost function to minimize the discrepancy between the prior knowledge (i.e., the computational model) and observations. On the other hand, the sequential DA involves evolving the state of the system with background information until observations get available, and then updating the system's state. One of the key components of the DA cycle is the forecast model. The forecast models used within DA are based on solving the governing equations or the best approximation of physical processes numerically using spatial and temporal discretization on a computational grid [9]. Despite their success, the current forecast models have difficulty representing complex processes like turbulence, convection and clouds [10,11] leading to high uncertainty in the prediction [12]. Additionally, the computational cost of ensemble forecasting (using many realizations with perturbed initial conditions) with numerical models is huge. As a result, the forecast model is one of the major limiting factors in DA, and this has spurred interest in using data-driven methods for the Earth system modeling that can deliver both computational efficiency and better representation of physical processes derived from data [13,14].
Data-driven methods have been applied to a wide range of problems in Earth system modeling [15]. For example, machine learning approaches have been explored for parameterization of subgrid-scale processes [16][17][18], precipitation nowcasting [19], superresolution of wind and solar data [20], and weather forecasting [21]. Some of the studies with the data-driven weather prediction (DDWP) have already started showing promising results with similar if not superior performance compared to state-of-the-art NWP models [22][23][24]. These DDWP models leverage deep learning methods like a convolutional neural network, Fourier neural operator, and recurrent neural network that are trained on data obtained from reanalysis products [25]. Furthermore, data-driven models can be augmented by taking up the information from physical knowledge, or at least respecting the conservation properties [26][27][28][29]. Some of the work towards these directions include enforcing the physical laws or statistical constraint into the loss function [30,31], tailoring the architecture design to enforce certain symmetry or conservation laws [32][33][34], symbolic regression for equation discovery [35], and end-to-end learning strategy [36].
These studies indicate that ML has the potential to improve scientific knowledge (and hence models), especially when we cannot express our understanding in the form of mathematical equations in physics-based models [37,38]. DA could benefit a lot from ML, where a forecast model is replaced with a hybrid model that incorporates a neural network as a component of the physical model [39,40] or as a complete replacement with a data-driven model for forecasting and state estimation [41,42]. Another advantage of an ML-based emulator is the quick computation of backpropagation gradient with automatic differentiation that can be used as a replacement for expensive adjoint solvers within variational DA [43,44]. Similarly, some of the challenges with ML such as handling uncertain and sparsely sampled data can be mitigated with DA. Brajard et al. [45] proposed a two-step process to learn the parameterization of unresolved processes where a neural network was trained to learn the model error obtained from the analysis state of the truncated model. Some other works at the intersection of DA and ML include iterative application of neural network and DA to learn the chaotic dynamics [46], enforcing conservation of mass constraint in DA [47], and for building tangent-linear and adjoint models of parameterization schemes in variational DA [48].
In this work, we seek to build a data-driven surrogate model trained on the analysis data and integrate it with a sequential DA cycle. More specifically, we will focus on a class of data-driven methods that are based on projection-based reduced-order modeling (ROM). These methods first extract the recurrent spatial structures from the high-dimensional data using singular value decomposition or convolutional autoencoder [49,50]. The next step is to either solve the projected governing equations on the dominant modes [51,52] or learn the projected dynamics with data-driven methods [53][54][55]. The recurrent neural network is one of the most popular algorithms in modeling the dynamics of lower-dimensional latent space, and has been applied for a wide range of applications [56,57]. One of the advantages of using projection-based ROM as a surrogate model is that it allows for DA on a reduced-order space rather than on a full-state and this latent assimilation framework leads to substantial speed improvement in contrast to DA [58][59][60][61]. The key challenge with latent assimilation is to link the real-time observations in physical space with the latent observation space [62]. Our approach for latent observation space construction is by first identifying the near-optimal sensor locations using QR pivoting and then using the measurements at these discrete locations for reconstruction on a tailored basis [63].
The paper is organized as follows. In Section 2.1, the surrogate modeling with proper orthogonal decomposition (POD) for dimensionality reduction, and a long short-term memory (LSTM) network for modeling dynamics is discussed. The detailed discussion on the deterministic ensemble Kalman filter (DEnKF) is provided in Section 2.2. The computation of the latent observation space is described in Section 2.3. In Section 3, the performance of the proposed framework is demonstrated for the NOAA Optimum Interpolation (OI) Sea Surface Temperature (SST) V2 dataset. Finally, concluding remarks are provided in Section 4.

Methods
In this section, we introduce different components of the proposed framework including surrogate modeling, sequential data assimilation, and optimal sensor placement strategy. The surrogate model is based on the projection-based model order reduction which relies on proper orthogonal decomposition (POD) and a recurrent neural network to model the dynamics of the latent space. This surrogate model is integrated within the sequential data assimilation to correct the forecast using real-time observations. The latent observation space is obtained by first reconstructing the full state from point sensor measurements and then projecting it onto POD bases.

Surrogate Modeling
Our surrogate model is constructed using the projection-based reduced-order modeling where the high-dimensional system is first compressed to a low-dimensional system and the dynamics of the low-dimensional system is modeled. The compression is achieved by projecting the high-dimensional data onto a set of optimal linear basis functions obtained via POD [64,65]. POD provides the optimal linear basis functions as they minimize the error between the true data and its projection in the L 2 sense compared to all other linear basis functions of the same dimension. Given the N s snapshots of the data for a state variable x ∈ R n , we can form the matrix A as follows where x (k) corresponds to an individual snapshot in time. Then, we perform the reduced singular value decomposition (SVD) as follows where U ∈ R n×Ns is a matrix with orthonormal columns which corresponds to the left-singular vectors, V ∈ R Ns×Ns is matrix with orthonormal columns representing the right-singular vectors, and Σ ∈ R Ns×Ns is a matrix with nonnegative diagonal entries, called singular values, and are arranged such that σ 1 ≥ σ 2 ≥ · · · ≥ σ Ns ≥ 0. For the dimensionality reduction task, only the first N r columns of U and V (denoted as U and V ) are retained along with the upper-left N r × N r sub-matrix of Σ (denoted as Σ). The reduced-order approximation of the matrix A can be written as follows The total L 2 error between the snapshot data matrix A and its reduced-order approximation A satisfies the following equalities [66] A where · 2 is the matrix-2 norm. Equation 4 states that across all possible matrices B ∈ R n×Ns of rank N r or less, the matrix A provides the best approximation in the L 2 sense and the error between A and its N r -order approximation A equals σ Nr+1 . From here on, the N r columns of U are called as the POD modes or basis functions and we denote them as Φ = {φ k } Nr k=1 . Once the POD modes are obtained, the compressed latent space for a single state vector x can be written as where α is the reduced-order approximation of the full state vector x and is also referred to as the POD modal coefficients. The reconstruction of the state vector is then computed as where x is the optimal reconstruction of full state vector x. The number of retained modes, i.e., N r , is decided based on the variance of the data captured by the retained modes. The singular values σ i give a measure of the quality of information that is retailed in N r -order approximation of the matrix A. The amount of the energy retained by POD modes can be calculated using the quantity called the relative information content (RIC) as follows The second step of the surrogate modeling is to model the evolution of the latent variables. If the exact equations governing the physical system are known, intrusive approaches like Galerkin projection can be applied to build the ROM. The ROM is obtained by substituting the low-rank approximation into the full order model and then taking the inner product with test basis functions to yield a system of N r ordinary differential equations (ODEs). However, for geophysical systems, the exact governing equations are unavailable or insufficient for the desired purpose (for example due to coarse grid resolution). On the other hand, the reasonably accurate and complete observational data for the evolution of the state of the system has been collected for many decades. This situation makes the equation-free techniques for predicting the future state of the dynamical system very attractive for multiscale and chaotic dynamical system modeling [67][68][69]. Recently, machine learning methods based on recurrent neural network have been applied successfully in many studies to build the non-intrusive ROM [56].
We apply the long short-term memory (LSTM) neural network to model the evolution of the latent variables. LSTM has been successfully applied in many studies dealing with modeling high-dimensional spatio-temporal chaotic dynamics of physical systems [54,70,71]. LSTM is a type of recurrent neural network (RNN) that can capture the long-term dependencies in the evolution of time series data [72]. RNNs contain loops that allow them to persist information from one time step to another and can be expressed as where h (t) ∈ R d h is the hidden state at time t, o (t) ∈ R d h is the input vector at time t, f h→h is the hidden to hidden mapping, and f h→o is the hidden to output mapping. The output of the model is the forecastõ (t+1) at time step t + 1.
One of the limitations of RNNs is vanishing (or exploding) gradient to capture the long-term dependencies. This problem occurs because the gradient is multiplied with the weight matrix repetitively during backpropagation through time (BPTT) [73]. The LSTM mitigates this issue by employing a memory cell composed of gating mechanism that decides which information to memorized or forgotten. The equations that implicitly define the mapping from hidden state of the previous time step (i.e., h (t−1) ) and input vector at the current time step (i.e., o (t) ) to the forecast hidden state (i.e., h (t) ) can be written as follows where g o ∈ R d h are the forget gate, input gate, and output gate, respectively. The The l is referred to as the lookback window which governs how much amount of the old temporal information is required to forecast the future state of the system accurately.
When we utilize the LSTM network for constructing a surrogate model, the reduced-order state of the system at a future time step, i.e., α (k+1) , is learned as the function of a short history of l past temporally consecutive reduced-order states as follows where F (· ; θ) is the nonlinear function parameterized by a set of parameters θ, and α is the low-dimensional approximation of the full state vector, i.e., the POD modal coefficients given in Eq. 5. Since the surrogate model does not use any governing equations of the system, it is also referred to as the non-intrusive reduced-order model (NIROM).

Data Assimilation
We consider the dynamical system whose evolution can be represented as where x (k) ∈ R n is the state of the system at discrete time t k , and M t k →t k+1 : R n → R n is the nonlinear model operator that defines the evolution of the system over the interval [t k , t k+1 ]. The term w (k+1) denotes the model error that takes into account any type of uncertainty in the model that can be attributed to boundary conditions, imperfect models, etc. Let z (k) ∈ R m be observations of the state vector obtained from sparse, noisy measurements and can be written as where q(·) is a nonlinear function that maps R n → R m , and v (k+1) ∈ R m is the measurement noise. We assume that the measurement noise is a white Gaussian noise with zero mean and the covariance matrix R (k+1) , i.e., v (k+1) ∼ N (0, R (k+1) ). Additionally, the noise vectors w (k+1) and v (k+1) are assumed to be uncorrelated to each other at all time steps.
The sequential DA can be considered as a Bayesian inference framework that estimates the state x (k+1) of the system given the observations up to time t k+1 , i.e., z (1) , . . . , z (k+1) . When we utilize observations to estimate the state of the system, we say that the data are assimilated into the model, and use the notation x (k+1) to denote an analyzed state estimate of the system at time t k+1 . When all the observations before (but not including) time t k+1 are applied for estimating the state of the system, then we call it the forecast estimate and denote it as x (k+1) f . The ensemble Kalman filter (EnKF) [74] follows the Monte Carlo estimation method to approximate the covariance matrix in the Kalman filter equations [75]. Instead of modeling the exact evolution of a probability density function under nonlinear dynamics, ensemble methods maintain an empirical approximation to the target distribution in the form of a set of ensemble members X (k) (i) for i = 1 . . . N . We begin by initializing the state of the system for different ensemble members , where x (0) represents the best-known state estimate at time t 0 , and P (0) is the initial covariance error matrix.
The propagation of the state for each ensemble member over the time interval [t k , t k+1 ] can be written as The prior state and the covariance matrix are approximated using the sample mean and error covariance matrix P (k+1) f as follows where the superscript T denotes the transpose, and A (k+1) f (i) is the anomalies between the forecast estimate for the ith ensemble and the sample mean. Once the observations get available at time t k+1 , the forecast state estimate is assimilated using the Kalman filter analysis equation as follows Unlike the EnKF algorithm, the DEnKF does not employ any perturbed observations. The Kalman gain matrix is computed using its square root version (without storing or computing P (k+1) f explicitly) as follows where Q ∈ R m×n is the Jacobian of the observation operator (i.e., Q kl = ∂q k ∂x l ), and the matrix A The anomalies for all ensemble members are then updated separately with half the Kalman gain as shown below The state for all ensemble members is updated by adding ensemble anomalies to analysis state estimate and can be written as where λ is the inflation factor to account for modeling errors. Inflation and covariance localization approaches are usually used in the EnKF framework to mitigate small number of ensembles [76][77][78]. The above ensembles are used as initial ensembles for the next assimilation cycle and the procedure is continued.
The DA procedure described so far corresponds to the full state vector of the system. However, in this study, the forward model is replaced with the surrogate model, and therefore, we perform the data assimilation in the latent space. The similar ideas have also been used in other studies that deals with the data assimilation for reduced-order models [43,59,61]. Once the LSTM network is trained, it is used to forecast the future state of the POD modal coefficients in an auto-regressive manner [79]. The evolution of the reduced-order model over the interval [t k , t k+1 ] with the LSTM network is given as follows where α (k+1) f is the low-dimensional forecast estimate of the system, and α (k) is the low-dimensional analyzed state of the system at time t k . Once the observations at time t k+1 gets available, they need to be processed to obtain the latent observations α (k+1) . Usually, the observations are available at very few sparse locations, i.e., m << n, and therefore, they need to be mapped from observation-space R m to state-space R n through some reconstruction technique. Peyron et al. [61] used the simple interpolation to learn this mapping while applying data assimilation to the indoor air quality problem. We utilize the POD reconstruction augmented with QR pivoting for learning the map from observations to full state, and its further details are provided in Section 2.3. The analysis equation in the latent space can be written as follows The Kalman gain matrix in Eq.29 is calculated as follows where Q ∈ R Nr×Nr is the Jacobian of the latent observation operator, and the matrix A (k+1) f ∈ R Nr×N is formed in latent space similar to Eq.25.

Reconstruction from Discrete Sensor Locations
As discussed previously, we need to reconstruct the full state of the system from limited number of discrete sensor locations. The low-rank approximation methods based on POD modes are one of the most popular method for the reconstruction task, where the sensor data is used to estimate the POD modal coefficients [63,80]. Additionally, the QR decomposition with column pivoting is used for the near-optimal sensor placement in contrast to random sensor placement. We are interested in estimating the full state of the system given the sensor measurements s ∈ R m at discrete locations. For this case, we have s i = x j for some 1 ≤ i ≤ m and 1 ≤ j ≤ n. We can write this as follows where Θ ∈ R m×n is constructed by taking m rows of n × n identity matrix. Therefore, each row of the matrix Θ will consist of all zeros except for the corresponding observation location, where it will have the value of one. We use the first m POD modes obtained from singular value decomposition of the matrix A given in Eq. 1 for the reconstruction task and is denoted as Ψ ∈ R n×m (i.e., the first m columns of U ). The Ψ is also referred to as the tailored basis functions in the literature. The approximation of the full state is given by where a ∈ R m . Once the sensor measurements gets available, the POD modal coefficients for the new sample can be calculated as follows Once the POD modal coefficients for the new sample is determined, the full state of the system can be reconstructed using Eq. 32. The latent observations for data assimilation are computed by projecting the reconstructed data onto the POD basis functions Φ as follows So far, we did not discuss the choice of the matrix Θ. The choice of the sensor locations can have a significant impact on the accuracy of the full state reconstruction. The sensor locations can be chosen either randomly or based on some heuristics or intuition of the physical problem. There are several techniques like optimal experimental design [81] and Bayesian criteria [82] that can be used to determine the optimal sensor locations for moderately sized problems. The QR factorization is a powerful and robust data-driven method to determine the near-optimal sensor locations solely based on the data [63,80]. In this method, we perform the QR decomposition with column pivoting of the matrix Ψ T . The QR decomposition calculates a column permutation matrix C T ∈ R n×n , an orthogonal matrix Q ∈ R m×m , and an upper triangular matrix R ∈ R m×n such that Ψ T C T = QR. The greedy approximation of the optimal sensor locations is obtained from the first m rows of matrix C, i.e., by setting Θ = C, where C corresponds to first m rows of the matrix C [63,83]. We note here that one can also use the QR decomposition method for oversampled case, where the number of sensor exceeds the number of tailored basis functions used for reconstruction.

Results and Discussion
We first describe the NOAA Optimum Interpolation (OI) Sea Surface Temperature (SST) V2 dataset used in this study for building a surrogate model and then integrating it within the data assimilation cycle. Then we present our numerical experimental results where we analyze the performance of our surrogate model and surrogate model assisted data assimilation framework.

Dataset and Preprocessing
We use the NOAA OI SST V2 analysis dataset for building our surrogate model. The analysis uses in situ (ship and buoy) and satellite SSTs plus SSTs simulated by the sea-ice cover. Before computing the analysis, the satellite data is adjusted for biases using the method of Reynolds [84] and Reynolds and Marsico [85]. This dataset consists of the weekly average sea surface temperature on a 1 • latitude × 1 • longitude global grid (180 × 360). The SST dataset exhibits a strong periodic structure due to seasonal fluctuations. Despite this seasonal periodicity, complex ocean dynamics lead to rich flow physics in this dataset. This dataset has been used in number of recent studies on flow reconstruction [86], geophysical emulation [87], and dynamic mode decomposition [88]. Here, we use the data from October 1981 to December 2000 (1000 snapshots) for building a surrogate model and the data from January 2001 to June 2018 (914 snapshots) for comparing the performance of the surrogate model for forecasting.

Numerical Experiments
In this subsection, we detail the numerical experiment design to assess the performance of NIROM and the integration of NIROM with DA, i.e., NIROM-DA. First, we use the data from October 1981 to December 2000 (1000 snapshots) to identify the POD basis functions. The masking operation is used to remove the data that correspond to the land area prior to its utilization for surrogate modeling. The temporal mean is also subtracted from the data and the POD is carried out for the unsteady component of the SST field. The first four POD modes capture approximately 90% of the variance of the data and are enough to capture the long-term and seasonal trends in the SST. Therefore, we fix N r = 4 in this study. The first four POD bases computed with POD are shown in Fig. 1. The next step is to train the LSTM network to learn the dynamics of a reduced-order system and a lookback window of l = 4 is used. Once the LSTM is trained, the model is deployed in auto-regressive manner for dynamical system forecasting. In the auto-regressive deployment, the initial condition for l time steps is used to predict the forecast state of the system at (l + 1)th time step. Then the state of the system from (2) − (l + 1) is used to determine the forecast state at (l + 2)th time step. The LSTM network architecture used in this study employs skip-connection after every hidden layer and is displayed in Fig. 2. The input and output data for training the LSTM network is scaled between [−1, 1] to accelerate the training. The LSTM network utilizes three layers of stacked cells, 80 neurons per LSTM cell, and a ReLU activation function. The LSTM network is trained using an Adam optimizer with the learning rate of 1 × 10 −3 and the weights of the network are initialized with the random normal initializer. The hyperparameters of the network are determined based on previous studies [89], and considerations of computational efficiency.
Once the surrogate model for SST is built, it is integrated into the DA cycle. The performance of the NIROM-DA framework is analyzed using a twin experiment for the period of January 2001 to June 2018. This data was not utilized in any stage of the NIROM construction. This ensures that there is no overlap in the training and online deployment of the surrogate model. For our twin experiment, the observations are generated by adding noise drawn from the Gaussian distribution with zero mean and the covariance matrix R k , i.e., v k ∼ N (0, R k ). We use R k = σ 2 b I, where σ b is the standard deviation of measurement noise and is set at σ b = 1. The observations are assumed to be collected only at discrete locations and the sensor locations are determined using the QR pivoting. We utilize only 300 discrete sensors for the reconstruction and this corresponds to less than 1% of the full state of the system. This number is also very Figure 1: The first four leading POD modes extracted from the SST dataset. The temporal mean was subtracted from the data and the POD modes are computed for he unsteady component of the SST field.
close to the optimal rank truncation threshold determined for this dataset from previous studies [63]. The number of tailored basis functions is also set equal to the number of sensors to avoid the intractable computation associated with the oversampled case. As discussed in Section 2.3, the tailored basis functions are just the POD modes, and we can measure the variance of the data captured by these modes using their singular values. The first 300 POD modes capture approximately more than 99.5% of the variance of the data and this leads to the improved reconstruction of small-scale features. Fig. 3 depicts the near-optimal sensor locations determined through QR pivoting.  The initialization of ensembles is very important for the sequential DA [90]. We experimented with three methods for initialization of ensembles, the first method is using random snapshots collected over some time to initialize the ensembles, the second one is adding random perturbation to the full state and then projecting it onto POD basis functions, and the third method is to add a random perturbation to the reduced-order state of the system. The first method is not suitable for our problem because of the seasonality in the dataset and this leads to the phase difference between the initial condition of ensembles. The second method did not lead to sufficient variability in the initial condition of the ensembles. The third approach ensures that the initial condition for all ensembles is sufficiently different without having any phase difference. Since we are doing both the forecast and data assimilation in latent space, the approach of adding a random perturbation to the true latent space is more suitable for our problem. Once all the ensembles are initialized, the forecasting is started using the NIROM in an auto-regressive manner. The observations are assimilated every third week to correct the initial condition for the future forecast. As discussed previously, the latent observations are obtained by first reconstructing the full state using the information at discrete sensor locations and tailored basis functions. The full state is then projected onto POD modes to compute latent observations for assimilation. The number of ensembles is set at 40 and the inflation factor of 1.5 is used to account for model error. Fig. 4 shows the evolution of POD modal coefficients for the forecasting period, i.e., 2001-2018. The POD modal coefficients can be interpreted as the contribution of different spatial frequencies (i.e., basis functions) in the evolving flow. One can observe that the first three modal coefficients are responsible for capturing the seasonal dynamics, and hence do not exhibit strong chaotic behavior. The fourth modal coefficient is significantly more chaotic than the first three modal coefficients due to the stochastic nature of small-scale fluctuations. The prediction from the NIROM is quite accurate in the initial period and the error gradually increases with time. This is a well-known limitation of the auto-regressive deployment of the LSTM for modeling the evolution of dynamical systems [87,91] and can be attributed to the error accumulation over time. The difference is particularly considerable for the first modal coefficient where there is a large phase difference between the true dynamics and predicted dynamics near the final time. Similar errors are also observed for the second and third modal coefficients along with a very inaccurate prediction for the fourth modal coefficient. The NIROM-DA framework can provide an analysis state that is very close to the true modal coefficients. The phase difference between the true and analyzed modal coefficients is very small leading to accurate modeling of seasonal dynamics. The analysis state of the fourth modal coefficient is also very close to the true state, meaning that the small-scale fluctuations are captured accurately with the NIROM-DA framework.
Next, we analyze the performance of NIROM and NIROM-DA frameworks in the reconstruction of the temperature field. The temperature field is reconstructed using the predicted modal coefficients and the POD basis functions with Eq. 6. We note here that the performance of our surrogate model strategy is limited by the number of POD modes used and the energy captured by those POD modes. Thus, we can at the most recover the true projection (TP) of the analysis temperature field (also referred to as the FOM) onto the POD modes. One can also interpret TP as the filtered version of the analysis temperature field. Thus, for our analysis, the root mean squared error (RMSE) is computed between the TP temperature field and the reconstructed temperature field from NIROM and NIROM-DA frameworks. Fig. 5 displays the evolution of RMSE in degrees Celsius for the forecasting period. Overall the RMSE for NIROM-DA is one order of magnitude less than using only NIROM for the prediction. The probability density of the RMSE for NIROM and NIROM-DA approaches is displayed in Fig. 6. The RMSE for NIROM lies mostly between 0.5 • Celsius to 1.5 • Celsius, while the RMSE for NIROM-DA is mostly centered around 0.2 • Celsius. This shows that a pure data-driven surrogate model can be built by exploiting the observational data collected over many decades and it can provide a sufficient level of accuracy in forecasting. Furthermore, the prediction from the forecast model can be improved using online sensor measurement in a computationally efficient manner using latent assimilation.
The quantitative findings for the NIROM and NIROM-DA frameworks are further supported by the visualization of the temperature field at two different times. Fig. 7 displays the temperature field on September 14, 2009 and the temperature field on June 21, 2018 is shown in Fig. 8. As we are retaining only four POD modes for the surrogate model, some of the small-scale features of the exact temperature field (i.e., FOM) are not captured by TP. While these small-scale features can be captured by retaining a large number of POD modes, training a recurrent neural network to learn the dynamics of lower-energy POD modes will require additional treatment. One potential solution for this can be the non-linear proper orthogonal decomposition (NLPOD) where an autoencoder is used to compress a large number of POD modal coefficients [92]. From Fig. 7, it is seen that even though larger structures in the flow are captured by NIROM, the error is higher for NIROM compared to NIROM-DA, especially in the northern hemisphere region. Similar   observations are also noted in Fig. 8, where the error for NIROM is larger than the NIROM-DA in both the northern and southern hemisphere regions.

Concluding Remarks
In this paper, we propose a novel framework to construct a non-intrusive reduced order model (NIROM) for geophysical flow emulation and integrate this low-dimensional model within the deterministic ensemble Kalman filter (DEnKF) algorithm to assimilate sparse and noisy observations in the latent space. The NIROM is based on a proper-orthogonal decomposition to identify the dominant modes from the data and then uses a recurrent neural network to learn the dynamics of low-dimensional latent space. This surrogate modeling strategy exploits the archival of observational data without relying on any kind of governing equations. One of the critical components of the proposed framework is the reconstruction of full state form sparse and noisy discrete sensor measurement. This is achieved using the least square estimation to map the information at near-optimal sensor locations determined through QR pivoting to the full state of the system. The latent observations for assimilation are generated by projecting this full state onto the POD modes.
We demonstrate the performance of our framework for the NOAA Optimum Interpolation Sea Surface Temperature (SST) analysis dataset. The NIROM is able to achieve a stable long-range forecast along with predicting larger structures in the temperature field with a sufficient level of accuracy. Once the NIROM is integrated within the DEnKF, the prediction is improved quantitatively by almost one order of magnitude. The results show that the NIROM can be readily coupled with the DA and latent space trajectory can be corrected every time sparse and noisy observations get available. With this framework, the computational saving is achieved through the replacement of a forward numerical solver with a data-driven model, and assimilation in latent space in contrast to the full state of the system. The proposed framework is extremely flexible and other algorithms can be easily accommodated. For example, one can use algorithms like convolutional autoencoder for dimensionality reduction, or shallow decoder for the reconstruction of full state from sparse discrete measurements. Furthermore, the NIROM can be easily differentiated with automatic differentiation and can be utilized in variational DA settings.
Finally, we point out that there is a number of future research directions that one can take. One such future direction is to improve the performance of the surrogate model by incorporating more modes to capture the low-energy content and frameworks like non-linear proper orthogonal decomposition [92] can be exploited for that. Another future research direction is to enhance the robustness of a long short-term memory network by reducing error accumulation during its auto-regressive deployment [91]. Our future work will also include exploring multi-fidelity data assimilation where a few ensembles of high-fidelity numerical solvers will be complemented with a large ensemble of data-driven models for forecasting [93] and employ a shallow decoder for full state reconstruction [94].
Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.