{'fraction_non_alphanumeric': 0.034942417555844395, 'fraction_numerical': 0.010367862648490161, 'mean_word_length': 4.548213986248532, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Combining data assimilation and machine learning to infer unresolved scale parametrisation

In recent years, machine learning (ML) has been proposed to devise data-driven parametrisations of unresolved processes in dynamical numerical models. In most cases, the ML training leverages high-resolution simulations to provide a dense, noiseless target state. Our goal is to go beyond the use of high-resolution simulations and train ML-based parametrisation using direct data, in the realistic scenario of noisy and sparse observations. The algorithm proposed in this work is a two-step process. First, data assimilation (DA) techniques are applied to estimate the full state of the system from a truncated model. The unresolved part of the truncated model is viewed as a model error in the DA system. In a second step, ML is used to emulate the unresolved part, a predictor of model error given the state of the system. Finally, the ML-based parametrisation model is added to the physical core truncated model to produce a hybrid model.

In recent years, machine learning (ML) has been proposed to devise data-driven parametrisations of unresolved processes in dynamical numerical models. In most cases, the ML training leverages highresolution simulations to provide a dense, noiseless target state. Our goal is to go beyond the use of high-resolution simulations and train ML-based parametrisation using direct data, in the realistic scenario of noisy and sparse observations.
The algorithm proposed in this work is a two-step process. First, data assimilation (DA) techniques are applied to estimate the full state of the system from a truncated model. The unresolved part of the truncated model is viewed as a model error in the DA system. In a second step, ML is used to emulate the unresolved part, a predictor of model error given the state of the system. Finally, the ML-based parametrisation model is added to the physical core truncated model to produce a hybrid model.
The DA component of the proposed method relies on an ensemble Kalman filter while the ML parametrisation is represented by a neural network. The approach is applied to the two-scale Lorenz model and to MAOOAM, a reduced-order coupled ocean-atmosphere model. We show that in both cases the hybrid model yields forecasts with better skill than the truncated model. Moreover, the attractor of the system is significantly better represented by the hybrid model than by the truncated model.

Introduction
The Earth climate system is one example of a natural system that is reasonably well represented through known physical laws and that has been intensively observed for decades (see, e.g., [1]). Physical laws, in the form of ordinary (ODEs) or partial differential equations (PDEs), are implemented through numerical models providing the time evolution of the system's state. Although weather and climate predictions have constantly improved, and will likely continue to do so, uncertainties will ineluctably remain. Those usually fall into two major classes: (i) the internal variability driven by the sensitivity to the initial conditions, and (ii) the model errors. The former has to do with the amplification of the initial condition error and arises even in perfect models -it is mitigated by using data assimilation (DA) [2]. The latter is present even if one would perfectly know the initial conditions and has to do with the incorrect and/or incomplete representation of the laws governing the system. The two sources of errors are inevitably entangled and it is difficult to separate them in practice.
Machine learning (ML) was recently shown to be effective in reducing model error, in particular that originating from unresolved scales. This has been achieved using two approaches. The first consists in learning a subgrid parameterisation of a model from existing physicsbased expensive parametrisation schemes [3,4], or from the differences between high-and low-resolution simulations [5][6][7][8]. In those approaches, the unresolved part of the model is represented by a ML process while the core of the model is derived from ODEs. The second approach is to emulate the entire model using observations. With spatially dense and noise-free data, this approach has been based on sparse regression [9], echo state networks [10,11], recurrent neural networks [12], residual neural network [13] or convolutional neural networks [14,15]. The challenging problem of partial and/or noisy observations has been addressed using dedicated NN architecture [16] or in combination with data assimilation methods [17][18][19][20][21].
This work presents a new method to obtain a data-driven parameterisation of a model's unresolved scale. In particular, we aim at producing a hybrid model combining the physics-based core (encoding the best of our knowledge of the resolved scales physics) with the data-driven parameterisation. By leveraging the use of DA, our method efficiently handles noisy and sparse observations.

Objectives and definitions (a) Statement of the problem
We consider an autonomous chaotic dynamical system, seen as our reference "truth", represented by the ODE with z(t) ∈ R Nz being the system's state at time t. From an arbitrary state z(t) on the system's attractor the model can be integrated forward for one time step δt, to get: where z δt = z(t + δt).
Let us formally define a projection operator Π : R Nz → R Nx such as x = Π(z), with x being the projection of the full state into a reduced dimension state: Nx < Nz. For example, Π can be a subsampling operator retaining only a subset of z or a downscaling operator from a high-resolution state to a lower resolution. From Eq. 2.2, we also define x δt = Π(z δt ).
Consider now a scale-truncated model, M r , that provides an imperfect description of Eq. (2.2) in the reduced space R Nx , where the superscript r stands for "resolved". When initialised from the projection of the truth into the reduced space (i.e. no initial condition error: x r = x), the difference between the 1-time step predictions from Eqs. (2.2) and (2.3) defines the model error due to the neglected scales in the resolved model The objective of this work is to complement the resolved model (2.3) by an empirical representation of the unresolved scales based on a neural network (NN) trained on incomplete and noisy data. We will hereafter denote the NN-based unresolved scale representation by g(x, θ). The vector θ is the set of trainable parameters of the NN and is determined by minimising the loss function where p Z is the invariant probability density function (assumed to exist) on the attractor. We construct the hybrid model M h , parametrised by θ, such that Our objective is to optimally estimate the parameters θ so that the hybrid model is the most accurate representation of the true underlying dynamics. Apart from trivial cases, the loss function in Eq. (2.5) cannot be computed: x δt is unknown and p Z is generally intractable. Assuming ergodicity of the true dynamics, we can however estimate it by a Monte-Carlo approach such that: where z 0:K−1 = {z 0 , z 1 , · · · , z K−1 }, z k = z(t k ), is a set of samples of p Z , typically extracted from a time series of modelled state variables and x k = Π(z k ). Such samples are usually not independent (due to the underlying dynamics being deterministic), and only provide a biased approximation of p Z . Furthermore, the need to sample the whole attractor implies treating time series significantly longer than the decorrelation time (i.e. K very large in general) .

(b) Framework of the study
The loss function, Eq. (2.7), cannot be minimised directly because some of its key entries are unavailable, in particular, obviously, the true process Eq. (2.2) and the time series z 0:K−1 . The available terms in the loss function are: The truncated model. The truncated model, our best available knowledge about the true process, is usually very complex (high dimensional, nonlinear and with diagnostic variables). Hence, we shall assume that its gradient, ∇xM r , cannot be computed analytically. We will thus focus on developing a (model) adjoint-free approach that is more flexible and suitable to high dimensional nonlinear scenarios where deriving and maintaining an adjoint model is a difficult and costly task [22]. When the gradient can be computed, other efficient methods exist [21].

Observations.
Observations are incomplete and noisy and are obtained through: where x k is the true state in the reduced space, y k ∈ R Ny and H k ∈ R Ny×Nx are the observation vector and operator respectively at t k , while o k is the observation error, assumed to be uncorrelated in time and normally distributed with mean 0 and a variance-covariance matrix (σ o ) 2 I Ny , where I Ny is the identity matrix of size Ny × Ny and σ o is the standard deviation. For simplicity, the observation error standard deviation is taken constant and the observation rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 0000000 operator linear. Both assumptions can be relaxed without major drawbacks even if it can induce practical difficulties. The ideal, most favourable, situation in which the full system's state is observed in the reduced space with no error is referred as the "perfect observation" case: y k = x k . For convenience, we further assume that observations are available regularly at multiples of the model time step such that ∆t = t k+1 − t k = Ncδt, Nc ∈ N * . This also accounts for the fact that, in general, the observation sampling period is longer than the integration time step of the numerical model.

Method (a) Loss function approximation
Let assume that an estimation of x k (k ∈ {0, · · · , K}) is available at observation times, every ∆t = Ncδt, so that where M (Nc) is the composition of the model Nc times, and similarly for the truncated model, x r k+1 = M r(Nc) (x k , δt). Since observations are not available at each time step (Nc > 1), the model error m is not known at each time step neither, and the loss function Eq. (2.7) can not be exactly computed. In the following, we will present two key simplifying assumptions that will lead to a tractable approximation of the loss function.
The first consists in assuming ∆t to be short enough so that the state evolved by the truncated model is independent of the model error due to the unresolved scale so that the model error is an additive term to the truncated model forecast after a ∆t time integration. The second in that the variability of m is small within ∆t. The combination of these two hypotheses is known as the linear superposition assumption, and can be formalised as: Note that the same approximation was made in a similar setting by [8].
The optimisation can now be performed using the approximate loss function: The modified loss function, Eq. (3.3) can be computed without knowing the full state z k but only its projection x k . L a can be minimised using a gradient descent algorithm as long as the gradient of g(x k , θ) can be computed, which is standard for any neural network library.

(b) Description of the algorithm
In order to minimise the loss function defined in Eq. (3.3), a sequence of x 0:K has to be available. Two cases are considered: the first is the aforementioned "perfect observations" case, in which we have a complete and noise-free sequence of a state variable x k in the reduced space. This ideal situation will set the upper-bound performance in the algorithm evaluation that follows. The second case is the more realistic case of noisy (but yet unbiased) and possibly incomplete observations. Here, a complete sequence x 0:K is obtained by processing the incomplete and noisy observations using DA [2]: the observations y k are combined with the forecast from the truncated model x r k in order to provide the analysed state vector x a k . The DA method used in this work is the Finite-Size Ensemble Kalman Filter (EnKF-N) [23] implemented in the DAPPER framework [24]. Even if the proposed algorithm is general and suitable for any DA algorithm, the EnKF-N has been chosen because of its efficiency. In particular, the inflation factor, a needed fix to mitigate the impact of sampling errors in the ensemble-based DA methods, is automatically estimated (thus avoiding long tuning). This inflation factor accounts also implicitly and partially for the effect of the model error.
The correction x a k − x r k made by DA is called analysis increment and was used to estimate model error due to unresolved scales in sequential DA in [25,26]. An analysis increment is composed of a correction both for the model error (which is what the hybrid model aims at estimating) and for the initial conditions error (which cannot be represented by the hybrid model). There is also an additional uncertainty due to the observation errors. The initial conditions and the observations errors are two sources of uncertainties called the data uncertainties. If they are too large, relatively to the model error, they could bias the estimation of the loss function given in Eq. (3.3) causing the hybrid model to overfit on these data and to lack generalisation to other initial conditions. To mitigate this problem, the time series x a 0:K estimated by DA is filtered using a simple low-pass filter (a rolling mean) producing a smoothed time series x s 0:K . This filter is expected to correct for data uncertainty provided that the observation error is uncorrelated in time, and thus contains high-temporal frequencies. On the other hand, this filtering could remove the fastest scales of the unresolved part of the model. Although it has been assumed in the linear superposition assumption that these high-frequencies can be neglected, we do not know a priori to which extent, which can lead to possibly hamper the forecast skill. The scale separation between the model error acting on long time scales and the initial errors acting on faster scales has been used in previous studies either to estimate the model error in DA [18] or to improve the forecast of a NN model [11]. Finally, note that the filtering can be adapted separately for the fast and the slow variables contained in x k . This is, for instance, the case in coupled atmosphere-ocean models, and is addressed in section 5. In both perfect and imperfect observations cases, the loss function L a (θ) is minimised using a standard NN training procedure: if g(x k , θ) is represented by a NN, x k as the inputs and θ as weight, the problem is equivalent to a supervised regression problem in which the targeted output of the neural network is m The algorithm is summarised in Algorithm 1.

(c) Numerical experiment protocol
The performance of the algorithm is evaluated on twin experiments: a full model M is used to produce true states, synthetic observations, and to assess the forecast skill of the hybrid model. A series of true states x 0:K is produced using the true model after projection in the reduced space. This series is used as perfect observations to obtain the so-called "perfect observationderived hybrid model". Then, synthetic observations are generated using the definition in Eq. (2.8). Here, it is assumed that the observation operator and the observation error statistics are perfectly known. This assumption is not needed for the functioning of the proposed algorithm and could be relaxed at a future stage. Nonetheless, the assumption is done here to simplify the interpretation of the results so the focus is on the truncated model error. Furthermore, the observation error statistics can also be estimated within the DA process itself (see [40] for a review of such methods), and thus be integrated into our combined ML-DA method.
The observations generated using Eq. (2.8) are used in Algorithm 1 to produce the so-called "DA-derived hybrid model". Note that the perfect observation-derived hybrid model benefits from the optimal information at a given time step (the state is completely observed without noise). Given the conditions stated in section 2(b), it is expected to be the best possible model and represents the benchmark against which we will test the DA-derived hybrid model.

(d) Evaluation metrics
The skill of each model is evaluated by running an ensemble of N f = 20 forecasts for N f different initial conditions and for a time period τ : x f(l) (τ ) (l = 1, · · · , N f ). The N f different initial conditions are chosen on the attractor of the true model (by running the model long enough before setting the initial conditions), independent of the values of the time series x 0:K used to perform the DA and the NN training. This ensemble of forecast constitutes the so-called test dataset, as it is not used to optimise nor to tune the algorithm. If tuning the algorithm is needed, the original time series x 0:K can be split in a training part (used to optimise the NN) and a validation part (used to tune the NN algorithm).
As evaluation metric we will use the relative root mean square error (R-RMSE) n (τ )) is the forecast from the hybrid or the truncated model (resp. the true model) at time t for the n-th component of x and for the sample l corresponding to a simulation for one particular initial condition. Vn(x t ) is the n-th component of the variance over the time dimension of the true model forecast time series x t . Note that if the truncated and/or the hybrid models have the same variability as the true one (i.e. the same variance in time), the R-RMSE converges to 1.

Application to the two-scale Lorenz model (a) Description of the model
The two-scale Lorenz model [27], hereafter L2S, is given by the following set of ODEs: where n = 0, · · · , Nx − 1 (Nx = 36) and m = 0, · · · , Nu − 1 (Nu = 360). The indices n are periodic, e.g., x 0 = x Nx . The values chosen for the parameters are the same as in [28]: the time scale ratio c is set to 10, the space-scale ratio b = 10, the coupling h = 1 and the forcing F = 10. Time t is expressed in model time unit, denoted MTU hereafter.
This set of two-scale ODEs, considered as the true model, is integrated using a fourth-order Runge-Kutta scheme with a time step of 0.005 MTU. The ODEs describing the evolution of x only represent the truncated model and are obtained by setting the coupling c to 0. It is integrated using a fourth-order Runge-Kutta scheme with a time step of 0.01 MTU.
(b) Setup of the reference experiment A so-called "reference experiment" is defined in this section. The true model is integrated over approximately 1500 MTU after a spinup of 3 MTU to produce the true state, on which observational noise is added. The EnKF-N is used to assimilate these observations using a large number N = 50 > Nx of ensemble members to reduce sampling errors. A noise is added to the state vector after each forecast to approximately account for the model error arisen from the model being truncated. It helps to avoid filter divergence and can be seen as additive inflation. This step is necessary given that, due to model error, the forecast ensemble would be otherwise underdispersive. The noise is assumed Gaussian with zero mean and standard deviation σ m = 0.06 optimised after tuning experiments (not shown here). In this reference experiment, the analysis obtained from the DA is not filtered, yet the sensitivity to the filtering of the DA analysis is studied in section 4(d).
The last step of the algorithm is to train a neural network to emulate the unresolved part of the model on the 1500 MTU time series produced by DA. The NN architecture is composed of convolutional layers (denoted "conv." in Table 1) where the non-linear activation function is a hyperbolic tangent (denoted "tanh"). Some additional parameters have been added, mainly to regularise the training: a batchnorm layer at the input layer, which normalises the training batch, and a L2-regularisation term on the parameters of the last layer. The parameters of the NN are optimised using the "RMSprop" [29] optimiser over 100 epochs. For each epoch, batches of 33 training examples are used to optimise the weights, until all the examples are consumed: this is a standard stochastic minimisation procedure [30]. Full details on the reference experiment are given in Table 1, in the column labelled L2S.

(c) Results
The terminology of the experiments described in the following is recalled in table 2. In Figure 1, both the true model and the DA-derived hybrid model (based on the reference experiment described in section 4(b)) are initialised from an initial condition on the attractor, chosen to be independent of the training set x 0:K . The true and the hybrid model are run over 5 MTU, and their difference is displayed. It can be noticed that both runs are very close until 2 MTU and that the hybrid model has predictive skill until 3-4 MTU for this particular set of initial conditions. Note that the Lyapunov time of the truncated model is 0.72 [33], meaning that the hybrid model provides accurate forecasts until 3 Lyapunov times.
In Figure 2, the R-RMSE is averaged over 20 members corresponding 20 initial conditions and also across all the 36 components of x. R-RMSE is displayed as a function of time. Several densities of observations have been considered: if Ny = 36 the full state is observed in the reduced space at each observation time. If Ny < 36, H k is a sub-sampling operator that draws randomly Ny values from the state following a uniform distribution changing the observation locations at each time step.
Results shown in Figure 2 (left panel) confirms that the DA-derived hybrid model has a predictive skill, significantly better than the truncated model until 4 MTU. The effect of reduced observation density is minor: the skill of the various hybrid models' forecasts is very similar. This shows the algorithm efficiency in handling sparse data to accurately train a NN model. This is a key strength of our method; most of the other approaches that parametrise a part of the model using ML assume dense observations, e.g. [3,4,7] (similarly to our perfect observation case).      Table 1.
probably due to the violation of the linear superposition assumption for high values of ∆t so that the coupling between the resolved part and the unresolved part of the model, as well as the effect of non-linearity of the unresolved part, are no longer negligible.
The right panel of Figure  is damped by the filter, but also by the increase of the observational noise that tends to add randomness on all scales, including the small ones. In this case, except for very strong noise, filtering does not seem to improve drastically the forecast skill. Notably, even a strong noise in the data has only a very small influence on the forecast skill of the hybrid model: when the noise on the observation is multiplied by a factor 20, the error in the forecast at t 0 + 2MTU is only multiplied by a factor 1.3.
From the results on the 2-scales Lorenz model, we conclude that the algorithm is robust against varying data spatial density, but is sensitive to their temporal distribution. Also, filtering the analyses obtained from DA may appear useful for slow processes but can deteriorate the results by filtering significant fast processes.

Application to a low-order coupled ocean-atmosphere model (a) Description of the model
We consider here the Modular Arbitrary-Order-Ocean-Atmosphere Model (MAOOAM) introduced by [34]. MAOOAM has 3 layers (2 for the atmosphere and 1 for the ocean) and is a reduced-order quasi-geostrophic model resolved in the spectral space. Its state is composed of na modes of the atmospheric barotropic streamfunction ψ a,i and the atmospheric temperature anomaly θ a,i respectively, plus no modes of the oceanic streamfunction ψ o,j and the oceanic temperature anomaly θ o,j respectively. The total number of variables is Nx = 2na + 2no. We consider two versions of MAOOAM: the true model with dimension Nz = 56 (na = 20,no = 8) corresponding to the configuration "atm. 2x-2y oc. 2x-4y" in [34] and the truncated model with Nx = 36 (na = 20,no = 8) corresponding to the configuration atm. "2x-4y oc. 2x-4y" in [34]. The truncated model is missing 20 high-order atmospheric variables (10 for the streamfunction and 10 for the temperature anomaly). Thus the truncated model does not resolve the atmosphere-ocean coupling related to these high order atmospheric modes.
The true model is used to generate synthetic observations. The forecast skill and the long-term properties of the truncated and the hybrid models will be evaluated by inspecting 3 crucial model variables, called key variables, that are ψ o,2 , θ o,2 and ψ a,1 (i.e. the second components of ocean streamfunction and temperature and the first component of the atmospheric streamfunction). They account for 42%, 51%, and 18% respectively of the variability of a reanalysis of 2-dimensional fields [35], and have been already used in previous studies (e.g. [36]). MAOOAM has also been recently used to study coupled data assimilation methods [37,38]. Unsurprisingly, in MAOOAM the ocean variables are considered the slow ones while the atmospheric variables are the fast ones.

(b) Experimental setup
We will express time in real time units (minutes, hours, days, ...) but, in practice, the model time is non-dimensional. Consequently, the dimensioned time values presented hereafter are not round numbers.
Given the diverse time scales and amplitudes of the MAOOAM variables, the noise parameters are all scaled on a climatological standard deviation of high frequencies σ hf ∈ R Nx , which is defined as the temporal standard deviation of the state vector after filtering out slow variations of a period longer than 1 month. This high-pass filter is carried out by subtracting the 1-month running average.
The parameters of the experiments are presented in Table 1, in the column labelled MAOOAM. The true model is integrated over approximately 62 years after a spinup of 30, 000 years, in the same configuration as in [34]. In all experiments with MAOOAM, the state is fully observed every 27 hours (∆t = 27 hours) (corresponding to Nc = 1, 000). A small modification was made to the observations from Eq. (2.8) to account for the fact that observations of the ocean are not at the same scale as those of the atmosphere: before being assimilated, instantaneous ocean observations are averaged over a 55 days rolling period centred at the analysis times. The EnKF-N is used as DA algorithm. The noise on the model forecast is added only to the atmospheric variables with a standard deviation of σ m = 10 −3 σ hf 1:20 . As mentioned in section 3(b), the analysis obtained from the DA is filtered. The slow processes are expected to occur mainly in the ocean, so only the ocean components of the state vector x a 0:K are filtered to produce x s 0:K . Differently from the L2S model experiments, filtering the analysis has proven necessary to train the hybrid model using MAOOAM.
The NN-architecture is a simple 3 layers multi-layer perceptrons; see Table 1 for full details. As opposed to the L2S model, the state vector has no locality properties (because it is defined in the spectral space), so the convolutional layers are not applicable (see the discussion about locality in [20]). The training of the NN is performed in the same way as for the L2S experiments.

(c) Results
The forecast skill metrics are presented in Table 3 for the truncated model as well as for the perfect observation-derived and the DA-derived hybrid models. Given the different time scales involved, the forecast lead time of the key atmospheric variable ψ a,1 is 1 day whereas the forecast lead time of the two key oceanic variables ψ o,2 and θ o,2 is 2 years. It can be seen that both perfect observation-derived and DA-derived hybrid models have superior skill to the truncated model. The improvement is larger for the ocean, with a factor of 2 to 3, and is similar for both hybrid models. Recall that the true model has here the same oceanic variables as the truncated model, so there is no difference in the representation of the pure oceanic processes. The improvement is thus fully due to an enhanced representation of the atmosphere-ocean coupling processes, the hybrid model better representing the interplay between the unresolved fast atmospheric variables and the slow oceanic variables.
The atmospheric key variable is improved to a lesser extent by the two hybrid models, and the perfect observations-derived model is significantly better than the DA derived model. This proves the limited capability of the hybrid model to represent a fast process, a situation further exacerbated in the case of the DA-derived hybrid model, when only noisy and partial data are at disposal. This result was indeed expected given the assumptions made on the unresolved term of the model in section 3(a) when a slow variation of the unresolved term was assumed. The fast processes are also less accurately represented because the sampling rate of the observations (27 hours) is well beyond the atmospheric time scale, and because of the presence of the observation errors (when applied).
In Figure 3, the attractors of the different models are displayed in the phase space defined by two key variable: ψ o,2 and ψ a,1 . A significant difference can immediately be seen between the attractors of the truncated and the true models: the truncated model visits areas of the phase space that are not admitted in the real dynamics. Remarkably, these discrepancies are reduced by the hybrid models both derived from perfect observation and from DA. Some states seem to remain out of the true model attractor, however, but much fewer.
Quantitative characterisation of the attractors is presented in Table 4, which provides the 3 quartiles (including the median) for each key variable. For the truncated model and the hybrid models, the difference of the quartiles is given relatively to the true model. For all the oceanic variables the distribution of the values of both hybrid models is significantly closer to the true distribution than for the truncated model. For the key atmospheric variable ψ a,1 , only the hybrid model derived from perfect observations shows an improvement. It confirms the conclusion made on the forecast skill that the hybrid model represents well the slow process, in particular oceanic variables in this case, and that the fast processes are not fully retrieved, in particular in case of the DA-derived hybrid model.    on realistic assumptions that only noisy and incomplete observations are available at a lower frequency than the model integration time step.

Conclusion
With a two-scale low-order chaotic system [27], we showed that the hybrid model forecast skill is sensitive to the observation frequency but very robust against high observational noise and sparse spatial distribution. This is probably due to the fact that increased observation frequencies challenge the validity of the linear superposition assumption more than large observational noise (see the discussion in Appendix of [39]). We then applied the method to the low-order coupled ocean-atmosphere model MAOOAM [34] which contains multiple temporal scales. Forecast skill and global statistics were significantly improved by the hybrid model compared with the truncated model encouraging further studies in high-dimensional and more realistic scenarios. Notably, the hybrid model derived from noisy observations has comparable forecast skill on the oceanic variables to that of the hybrid model derived from perfect observations.
In view of operational systems, it should be noted that the proposed algorithm relies on two existing data assimilation and neural networks training techniques that both scale well in highdimension (see, e.g., [41] and [42]). In principle, the present algorithm can be applied to larger and more realistic problems. In particular, the fact that the method does not rely on the adjoint of the truncated model is an advantage in terms of code maintenance. However, we foresee some practical challenges: for instance, the computational architecture and the data types used for physics-based numerical models and for machine learning algorithms can be very different (e.g. multi-core supercomputers and graphics processing units). Training and running hybrid models efficiently imposes heavy requirements on both the hardware and software and may come with an overhead even if some tools are very promising [43].
The approach presented here can also accommodate the additional representation of the remaining model error (i.e. the model error of the hybrid model): it could either be done within the numerical model by parameterising the model error [21], or by training stochastic neural networks [44].
Data Accessibility. Data used in all the experiments are available at ftp://ftp.nersc.no/reddaml/.
The instructions to download the data, to run all the Lorenz 2 scale application and to reproduce figures and tables of the article can be found at https://github.com/brajard/reddaml/releases/tag/v1.0.