Trajectory learning for ensemble forecasts via the continuous ranked probability score: a Lorenz '96 case study

TL;DR

This study explores trajectory learning for ensemble forecasts using CRPS as a loss function, demonstrating improved accuracy and calibration in the Lorenz '96 model, especially for short-term predictions.

Contribution

It introduces a CRPS-based trajectory learning method for stochastic parametrizations that outperform traditional derivative-fitting approaches in ensemble forecasting.

Findings

CRPS-based parametrizations are more accurate and sharp.

The approach simplifies calibration of ensemble models.

Outperforms derivative-fitting methods in short-term forecasts.

Abstract

This paper demonstrates the feasibility of trajectory learning for ensemble forecasts by employing the continuous ranked probability score (CRPS) as a loss function. Using the two-scale Lorenz '96 system as a case study, we develop and train both additive and multiplicative stochastic parametrizations to generate ensemble predictions. Results indicate that CRPS-based trajectory learning produces parametrizations that are both accurate and sharp. The resulting parametrizations are straightforward to calibrate and outperform derivative-fitting-based parametrizations in short-term forecasts. This approach is particularly promising for data assimilation applications due to its accuracy over short lead times.

Tables

Table 1. Table 1: Parameter values for the L96 experiments

Description	Parameter	Value
Number of $X$ variables	$K$	8
Number of $Y$ variables per $X$ variable	$J$	32
Coupling constant	$h$	1
Forcing term	$F$	20
Spatial scale ratio	$b$	10
Temporal scale ratio	$c$	4 or 10

Table 2. Table 2: Parameters for the stochastic processes used in the local model.

$c=4$	$c=10$
$S^2=4.61$	$S^2=4.04$
$\varphi =0.974$	$\varphi =0.977$

Table 3. Table 3: Summary of the tested model configurations.

Name	Model type	Model form	Noise
poly_gauss	Derivative fitting	Local	White, independent, additive
poly_ou	Derivative fitting	Local	Correlated, independent, additive
svd_gauss	Derivative fitting	Global	White, independent, additive
svd_ou	Derivative fitting	Global	Correlated, independent, additive
crps_ou2	Trajectory fitting, $N_t=2$	Global	Correlated, coupled, additive
crps_ou8	Trajectory fitting, $N_t=8$	Global	Correlated, coupled, additive
crps_ou32	Trajectory fitting, $N_t=32$	Global	Correlated, coupled, additive
crps_mult2	Trajectory fitting, $N_t=2$	Global	Correlated, coupled, multiplicative
crps_mult8	Trajectory fitting, $N_t=8$	Global	Correlated, coupled, multiplicative

Table 4. Table 4: Kolmogorov–Smirnov (KS) and Hellinger ( H 2 \mathrm{H}^{2} ) distances, computed from the empirical PDFs of the X X -variables over long simulation times. The PDFs are shown in Fig. 7 .

	$c=4$		$c=10$
Model	$\mathrm{KS}$	$H^2$	$\mathrm{KS}$	$H^2$
poly_gauss	0.136	0.258	0.238	0.466
poly_ou	0.204	0.390	0.269	0.513
svd_gauss	0.0693	0.160	0.160	0.377
svd_ou	0.240	0.449	0.330	0.624
crps_ou2	0.173	0.383	0.241	0.515
crps_ou8	0.0504	0.109	0.129	0.303
crps_ou32	0.0506	0.108	0.122	0.299
crps_mult2	0.0838	0.199	0.138	0.313

Table 5. Table 5: Parameters for the AR(1)-process used in the global model, per POD mode and per test case.

	$c=4$	$c=10$
${\xi }_1$	0.990	0.997
${\xi }_2$	0.991	0.997
${\xi }_3$	0.990	0.996
${\xi }_4$	0.991	0.995
${\xi }_5$	0.988	0.993
${\xi }_6$	0.988	0.992
${\xi }_7$	0.989	0.992
${\xi }_8$	0.985	0.993

Table 6. Table 6: Overview of the training data for the CRPS-based models. Shown are the trajectory length N t N_{t} and corresponding time interval in MTU, the batch size N b N_{b} , and the final losses per batch observed during training for the two considered test cases.

	Training parameters			Final loss
	$N_t$	MTU	$N_b$	$c=4$	$c=10$
Additive noise	2	0.010	200	0.290	0.402
	4	0.020	100	0.492	0.694
	8	0.040	50	0.721	1.08
	16	0.080	25	1.01	1.52
	32	0.16	12	1.55	2.32
	64	0.32	6	2.97	4.55
	128	0.64	3	5.07	8.51
	200	1.0	2	14.3	14.3
Multiplicative noise	2	0.010	200	0.0367	0.0707
	4	0.020	100	0.101	0.195
	8	0.040	50	0.275	0.410
	16	0.080	25	0.562	0.730

Equations

$$L_c(C(q))+M(C(q),\theta )=C(L(q)),$$

$$\mathrm{Loss}(\theta ,q)={\left|L_c(C(q)+M(C(q),\theta )-C(L(q))\right|}^p.$$

$$\mathrm{Loss}(\theta ,q)=\sum _{j=1}^{N_t}{\left|z_{\theta }(t_j)-C(q(t_j))\right|}^p$$

$${\pi }_{\mathcal{A}}^{t_1}(q^{\prime })={\int }_{\mathcal{A}}\delta (q^{\prime }-{\psi }_{[t_0,t_1]}(q)){\pi }_{\mathcal{A}}^{t_0}\text{d}q.$$

$${\pi }_{\mathrm{CG}}^{t_1}(z)={\int }_{\mathcal{B}}{\pi }_{\mathcal{A}}^{t_1}(q)\text{d}q,$$

$${\phi }_{[t_0,t_1]}(z,\theta ):=z+L_{c,[t_0,t_1]}(z)+M_{[t_0,t_1]}(z,\theta ),$$

$$\begin{array}{rl}\mathrm{Loss}(\theta ,q)& =\sum _{j=1}^{N_t}\mathrm{CRPS}(z_{\theta },C(q),t_j),\\ \mathrm{CRPS}(z_{\theta },q,t_j)& =\frac{1}{M}\sum _{i=1}^M\left|z_{\theta }^i(t_j)-C(q(t_j))\right|-\frac{1}{2M^2}\sum _{j=1}^M\sum _{k=1}^M\left|z_{\theta }^i(t_j)-z_{\theta }^k(t_j)\right|,\end{array}$$

$$\frac{\text{d}{X}_k}{\text{d}t}$$

$$\frac{\text{d}{Y}_j}{\text{d}t}$$

$$\frac{\text{d}{X}_k}{\text{d}t}=-X_{k-1}(X_{k-2}-X_{k+1})-X_k+F+M(X_k,t);k=1,\dots ,K.$$

$$U(X,t)=-X_{k-1}(X_{k-2}-X_{k+1})-X_k+F-\left(\frac{X_k(t+\Delta t)-X_k(t)}{\Delta t}\right),$$

$$M(X_k,t)=M_{\det }(X_k)+r(t),$$

$$M_{\det }(X_k)$$

$$r(t+\Delta t)=\varphi r(t)+\sqrt{\Delta t{S}^2(1-{\varphi }^2)}\text{penalty}10000\epsilon (t).$$

$$U(t)={\xi }_0+\sum _{i=1}^8{\xi }_i{\lambda }_i{a}_i(t).$$

$$M(t)={\xi }_0+\sum _{i=1}{\xi }_i{\lambda }_i{r}_i(t),$$

$$M(t)=\sum _i{\xi }_i{c}_i(X){\mathbf{r}}_i(t),$$

$$\mathbf{r}(t+\Delta t)=\mathbf{r}(t)+\mathbf{A}(\mathbf{\mu }-\mathbf{r}(t))+\mathbf{B}\mathbf{\epsilon }\sqrt{\Delta t}$$

$$\text{Loss per batch}=\frac{1}{N_b}\frac{1}{N_t}\sum _{i=1}^{N_b}\sum _{j=1}^{N_t}\mathrm{CRPS}(X(\theta ),C(q),t_j),$$

$$\begin{array}{rl}\mathrm{Loss}(\theta ,q,\alpha )& =\sum _{j=1}^{N_t}\mathrm{CRPS}(z_{\theta },C(q),t_j,\alpha ),\\ \mathrm{CRPS}(z_{\theta },C(q),t_j,\alpha )& =\frac{1}{M}\sum _{i=1}^M\left|z_{\theta }^i(t_j)-C(q(t_j))\right|-\frac{\alpha }{2M^2}\sum _{j=1}^M\sum _{k=1}^M\left|z_{\theta }^i(t_j)-z_{\theta }^k(t_j)\right|,\end{array}$$

$$\mathrm{MSE}(x,y,t)=\frac{{\sum }_{i=1}^M{\left(x_i(t)-y(t)\right)}^2}{M}.$$

$$\mathbb{E}{\left(\frac{1}{M}\sum _{i=1}^M{x}_i-y\right)}^2$$

$$\mathbb{E}\frac{1}{M}\sum _{i=1}^M{\left(x_i-\frac{1}{M}\sum _{j=1}^M{x}_j\right)}^2$$

$$\mathrm{CRPS}(F,y_{\mathrm{obs}})={\int }_{-\infty }^{\infty }(F(s)-{\mathbb{I}}_{\{s\ge y_{\mathrm{obs}}\}}{)}^{2}ds.$$

$$\mathrm{CRPS}(x,y_{\mathrm{obs}})=\frac{1}{M}\sum _{i=1}^M\left|x_i-y_{\mathrm{obs}}\right|-\frac{1}{2M^2}\sum _{i=1}^M\sum _{j=1}^M\left|x_i-x_j\right|.$$

$$\mathrm{KS}(P,P^{\prime })=\underset{i}{\max }|P_i-P_i^{\prime }|,$$

$${\mathrm{H}}^2(p,p^{\prime })=\frac{1}{\sqrt{2}}\sqrt{\sum _{i=1}^M{\left(\sqrt{p_i}-\sqrt{p_i^{\prime }}\right)}^2},$$

Full text

1 Introduction

Addressing uncertainty stemming from unresolved physical effects or uncertain measurements constitutes a core challenge in geophysical fluid dynamics. Incomplete representations of physical processes, finite truncation of computational models used for dynamic simulations, and uncertain initial conditions all contribute to numerical prediction uncertainty. Stochastic parametrizations provide a means to simultaneously compensate for model error and predict the inherent uncertainty via ensemble forecasts, based on available observational data. The modeling objective thus shifts to calibrating a noise model or parametrization. In this work, we explore the calibration of stochastic models for ensemble forecasting via trajectory learning based on the continuous ranked probability score (CRPS).

Stochastic methods have been effectively utilised in numerical weather prediction for several decades [22, 5]. In recent years, data-driven stochastic models have been actively developed. The two-scale Lorenz ’96 system [32], a simplified representation of atmospheric dynamics with interactions between slow and fast scales, commonly serves as an initial test for novel calibration strategies [4]. Neglecting the fast variables is analogous to coarsening in fluid simulations, and hence introduces error and uncertainty in the evolution of the slow variables. Consequently, data from fast scales can inform stochastic parametrizations. For instance, [45] and [2] employed explicitly computed sub-grid scale data to tune parametrizations, fitting them as polynomials with a stochastic component for residuals, modeled as correlated or uncorrelated noise. Ensemble predictions with stochastic forcing were found to outperform perturbed parameter ensembles, with correlated additive noise yielding reliable forecasts. [10] modeled the sub-grid effects using Markov chains conditional on the resolved variables, achieving good agreement in the autocorrelation of the resolved variables and in their long-time distribution. Generative adversarial networks were used in [16] to construct maps between resolved variables and sub-grid scale forcings, producing accurate spatiotemporal correlations.

The development of data-driven stochastic models extends beyond low-dimensional dynamical systems, with significant theoretical and methodological advancements made in fluid models. Examples include approaches based on statistical turbulence theory [15] or the assumption of spatiotemporal coherent error patterns [5]. Modern approaches include developing structure-preserving stochastic models that exploit the geometry of geophysical fluid-dynamical models to construct specific multiplicative noise. For example, the framework of Stochastic Advection by Lie Transport (SALT) [26] is well suited for advection-dominated flows. This approach has yielded efficient uncertainty quantification for the two-dimensional Euler equations [9, 13, 47] and (thermal) quasi-geostrophic equations [8, 25] in a reduced-order setting, with applications in data assimilation [7].

A pragmatic low-rank approach to account for uncertainty and model error in fluid-dynamical systems involves separating spatiotemporal forcing into fixed spatial profiles and corresponding time series. These time series are subsequently modeled as stochastic processes. Proper orthogonal decomposition modes (POD) modes, also referred to as empirical orthogonal functions (EOFs) in atmospheric science [44, 46, 20], efficiently capture coherent structures observed in data. For example, [1] compared data-driven methods for ocean models and found that stochastic processes with memory effects performed well in both short- and long-time forecasts. Global basis methods have also been employed to capture flow statistics in coarsened turbulence simulations, resulting in stochastic parametrizations that reproduce energy spectra, total energy, or enstrophy [24, 11, 14, 12].

In this work, we aim to produce a forecast with an CRPS and thereby tackle the model calibration challenge. While this challenge does not coincide with parameter estimation in stochastic differential equations (SDEs), they are related. Methods based upon maximum likelyhood estimation [39, 38], generalised method of moments [21], and the efficient method of moments [17] are among the approaches outlined in the review [37] that bear similarity to the CRPS learning approach adopted in this paper. However, parameter estimation in SDEs differs from producing a skillful forecast. The CRPS approach used in this study, can be viewed as both a data assimilation technique for ensemble forecasting and as a model calibration strategy [27].

Compensating for model error and unresolved dynamics remains a longstanding challenge in computational fluid dynamics, with the resulting models generally referred to as ‘closure models’. Machine learning has proven to be an efficient approach for model calibration (see [41] and references therein). A robust method is trajectory fitting, where the discrepancy between predicted solution trajectories and reference trajectories is minimised [36]. However, many learned closure models are deterministic. Learning stochastic models is a logical next step necessary for probing prediction uncertainty via ensemble methods [6], and currently forms an active research area. A particular inspiration for the present study is the recent work by [30], who trained a stochastic model based on the almost fair CRPS to avoid problems with score degeneracy that can arise when all but one ensemble member has the same value as the reference. This approach is employed in the Artificial Intelligence Forecasting System (AIFS) developed at the ECMWF in test cases at two resolutions, where the trained stochastic model demonstrated improved variability compared to train deterministic models. The stochastic model was optimised for short-range forecasts but also showed good performance for longer range forecasts. Similarly, in [29] neural global circulation models for the atmosphere are trained using a stochastic loss based on the CRPS, giving strong practical evidence suggesting that minimising the CPRS to predict short-term weather is an effective way to tune parameterizations for climate.

Contributions and paper outline

This work investigates the feasibility of calibrating stochastic parametrizations for ensemble forecasts through trajectory learning, employing the CRPS as a loss function. This choice of loss function directly incorporates forecast reliability and skill into the training procedure. A simple linear stochastic model form is used to ensure interpretability and serves as a proof of concept for CRPS-based trajectory learning. Specifically, a multidimensional coupled Ornstein–Uhlenbeck process is adopted in additive and multiplicative settings to account for unresolved dynamics and estimate uncertainty. The presented approach and results inform future studies on stochastic parametrization in fluid models based on trajectory learning of ensemble forecasts.

The paper is structured as follows. The noise calibration methodology is introduced in Section 2. The governing equations are provided in Section 3, along with the adopted parametrizations and calibration details. The model performance is assessed in Section 4, after which the paper is concluded in Section 5.

2 Noise calibration methodology

This section introduces the noise calibration methodology. We first summarize the derivation of a general deterministic sub-grid scale measurement procedure in Section 2.1. This is then extended to a probabilistic procedure in Section 2.2, where the use of the CRPS as a loss function is motivated.

2.1 Deterministic sub-grid scale measurement procedure

We consider two dynamical systems. The first is the ‘true’ dynamics evolving according to $\dot{q}=L(q)$. The second is the ‘coarse’ dynamics described by $\dot{z}=L_{\mathrm{c}}(z)$. The true dynamics is integrated in time using the ‘high-fidelity solver’, whereas the coarse dynamics is integrated using the ‘low-fidelity solver’. Corresponding solution snapshots are referred to as high-fidelity and low-fidelity snapshots, respectively. We introduce the solution spaces $\mathcal{V}$ for the true solution and $\mathcal{W}$ for the coarse solution. This formulation allows $\mathcal{V}$ and $\mathcal{W}$ to be either function spaces or spaces of solutions on discrete meshes. In the following, we treat these as the latter, where elements of $\mathcal{V}$ are high-resolution numerical solutions and elements of $\mathcal{W}$ are low-resolution numerical solutions. We further introduce a coarsening operator $C:\mathcal{V}\to\mathcal{W}$, which is generally assumed to discard information and is thus not invertible.

The true dynamics will generally be either unknown or too computationally expensive for practical forecasting. On the other hand, the coarse dynamics provides a computationally feasible alternative but suffers from model incompleteness and other sources of error. Therefore, a forcing term will be included in the coarse dynamics to compensate for these errors. Generally, we refer to the discrepancy between $L_{c}$ and $L$ as ‘unresolved scales’ leading to ‘unresolved dynamics’.

A perfect model is one where the modeled coarse dynamics equal the coarsened true dynamics. This requires defining a model correction term, or parametrization, $M(C(q),\theta)$, such that the complete coarse dynamics $L_{c}+M$ satisfies

[TABLE]

where $\theta$ is a tunable set of parameters.

It is generally not possible to find closed-form expressions for this perfect model, as it depends on properties of the governing dynamics as well as on adopted numerical methods. Nonetheless, approximations to this model can be found using observational data. A direct approach based on (1) is derivative fitting, where the parameters $\theta$ are chosen to minimize the loss function

[TABLE]

That is, the objective is to match the instantaneous coarse-grained dynamics over a given training data set. Here, $p$ determines how strongly outliers are penalised. Naturally, regularization terms can be included in the loss function without loss of generality, but they will not be further considered here.

An alternative approach is trajectory learning, where the model aims to minimize the distance between the low-fidelity dynamics and observed coarse-grained high-fidelity dynamics over specified integration intervals. We denote a reference trajectory by $q(t_{j}),\penalty 10000\ j=1,\ldots,N_{t}$, generated by the dynamics $\dot{q}=L(q)$, and a model trajectory by $z_{\theta}(t_{j}),\penalty 10000\ j=1,\ldots,N_{t}$, generated by the low-fidelity dynamics including the model as $\dot{z}=L_{c}(z)+M(z,\theta)$. Then, the model parameters $\theta$ are chosen such that the loss function

[TABLE]

is minimised over trajectories within a training data set. Here, $N_{t}$ denotes the length of each trajectory. Trajectory fitting can thus account for long-term dynamics (temporal non-locality) and error accumulation when the trajectory length is appropriately chosen.

A concise comparison between derivative fitting and trajectory learning is provided in [36, Table 1]. In brief, the difference between derivative fitting and trajectory learning decreases when the trajectory length is shortens and the time instances $t_{j}$ become sufficiently close. However, when data are sparse or noisy, approximations of the instantaneous dynamics based on available measurements become less accurate and derivative fitting loses reliability. Simultaneously, derivative fitting requires time derivatives of training data, which may not always be readily available. Conversely, trajectory fitting requires either a discretised adjoint equation or a differentiable numerical solver. Furthermore, if the trajectory length is too long or the time instances are spaced too far apart, chaotic behavior and perturbation sensitivity of the dynamical system may decrease the effectiveness of trajectory learning.

2.2 Probabilistic sub-grid scale measurement procedure

The assumption that the coarsening operator $C$ discards information has considerable consequences for model derivation. Namely, the preimage of a coarse configuration $z$ through $C$ is a distribution of high-fidelity solutions $\{q_{i}\}_{i}\in\mathcal{V}$. Since each $q_{i}$ satisfies $C(q_{i})=z$, it is impossible to determine which high-fidelity solution produced the particular coarse configuration $z$. Hence, there is uncertainty in the ‘perfect’ evolution of a coarse-grained model. This motivates a probabilistic approach to sub-grid model derivation, aiming to match distributions of solutions rather than a single solution. Stochastic modeling is well-suited for this purpose.

This intuition can be made more precise as follows. We let $\pi_{\mathrm{full}}$ denote the distribution of fully resolved configurations, where $\pi_{\mathrm{full}}^{t_{0}}(q)$ describes the probability of a configuration $q$ being the realization of the dynamics at time $t_{0}$. The true dynamics are defined via a high-fidelity flow map $\psi_{[t_{0},t_{1}]}(q):=q+L_{[t_{0},t_{1}]}(q)$, evolving the dynamics from $t=t_{0}$ to $t=t_{1}$ using the high-fidelity solver.

Given a coarse configuration $z_{0}$ at $t_{0}$, we denote the set of corresponding full configurations by $\mathcal{A}=\{q\in\mathcal{V}\mid C(q)=z_{0}\}$ and its distribution by $\pi_{\mathcal{A}}^{t_{0}}$. The high-fidelity flow map is deterministic, hence the evolution of this distribution under the high-fidelity dynamics can be expressed via a transition kernel $\tau(q^{\prime}|q)=\delta(q^{\prime}-\psi_{[t_{0},t_{1}]}(q))$. At time $t_{1}$, the distribution is expressed as

[TABLE]

The configurations that make up $\pi_{\mathcal{A}}^{t_{1}}$ are elements of $\mathcal{V}$, i.e., these are fully resolved configurations. The corresponding distribution of coarse-grained fields is obtained by marginalizing,

[TABLE]

where $\mathcal{B}=\{q\in\mathcal{V}\mid C(q)=z\}$. Intuitively, $\pi_{\mathrm{CG}}^{t_{1}}(z)$ is the distribution of attainable coarse-grained fields at time $t_{1}$, given the coarse configuration $z_{0}$ at $t_{0}$ and the uncertainty induced by the coarsening operator $C$.

Hence, from this perspective, the sub-grid scale model and the low-fidelity solver should define a transition kernel $\tau_{\mathrm{model}}$ that satisfies $\tau_{\mathrm{model}}(z^{\prime}|z)=\pi_{\mathrm{CG}}^{t_{1}}(z^{\prime})$. In this work, the transition kernel is defined implicitly via the stochastic forward flow map

[TABLE]

where $L_{c,[t_{0},t_{1}]}$ denotes low-fidelity integration from time $t_{0}$ to $t_{1}$. The sub-grid scale model is denoted by $M$ and depends on parameters $\theta$. Subsequently, an ensemble prediction obtained with the flow map (6) will serve as a Monte Carlo approximation to the distribution $\pi_{\mathrm{CG}}^{t_{1}}$.

The distribution $\pi_{\mathrm{CG}}^{t_{1}}$ of coarse-grained field will generally not be known. Monte Carlo approximations can be computed using the high-fidelity flow map $\psi_{[t_{0},t_{1}]}$, yet this map can be unknown or may be prohibitively expensive to compute. Instead, a single observed trajectory can be treated as a time-varying delta distribution to facilitate probabilistic trajectory learning. To this end, we consider an ensemble forecast $\{z_{\theta}^{i}(t_{j})\}_{i=1,\ldots,M},\penalty 10000\ j=1,\ldots,N_{t}$, generated by the flow map (6). Here, the superscript $i$ denotes the $i^{\mathrm{th}}$ ensemble member, $M$ is the ensemble size, and $N_{t}$ is the trajectory length. As before, we let $q(t_{j}),\penalty 10000\ j=1,\ldots,N_{t},$ be the reference trajectory generated by the high-fidelity dynamics. The continuous ranked probability score (CRPS) then serves as a probabilistic extension to the mean absolute error, enabling trajectory learning for ensemble predictions. The objective becomes to choose model parameters $\theta$ such that the loss function

[TABLE]

is minimised over a set of training trajectories.

The first term in the CRPS represents the ensemble mean absolute error, while the second term quantifies the ensemble spread and is subtracted from the mean error. On its own, a larger spread would reduce the CRPS. However, in practice, an increase in spread usually leads to a larger mean error, which causes the CRPS to increase. Hence, a small CRPS implies both an accurate ensemble mean absolute error and a sharp ensemble forecast (small spread). Furthermore, the CRPS is suitable for ensemble forecasts of nonlinear systems, as no assumptions are made on the distribution of the samples. As a proper scoring rule, the CRPS is minimised only when the forecast distribution matches the observed distribution. This ensures that training with the aforementioned loss function favors both ensemble accuracy and sharpness. Moreover, the absence of squared errors in the loss function makes the training less sensitive to large outliers, thereby increasing training robustness in chaotic systems. For additional details regarding the CRPS, we refer to [3, 18, 19, 23, 35]. An illustration of CRPS-based trajectory learning of ensemble forecasts is provided in Fig. 1.

3 Governing equations and models

In this section, we detail the governing equations and stochastic parametrizations. We differentiate between ‘true’ reference dynamics, given by the deterministic two-scale Lorenz ’96 (L96) model, integrated at high temporal resolution. A computationally cheaper ‘coarsened’ forward model comprises a single-scale L96 model, integrated at low temporal resolution. Thus, the latter exhibits deficiencies due to reduced temporal resolution and the omission of fast-scale dynamics. These deficiencies are parametrised by means of stochastic models, using synthetic data generated from the high-fidelity reference dynamics.

3.1 Governing equations

True dynamics

The two-scale L96 system describes the interaction of variables across multiple spatial and temporal scales. The large (slow) scales are denoted by $X$, and these are coupled to the small (fast) scales, denoted by $Y$. The governing equations read

[TABLE]

Periodic boundary conditions are imposed on the variables, which translates to $X_{k+K}=X_{k}$ and $Y_{j+JK}=Y_{j}$. The two-scale L96 system (eqs. 8 and 9) is a low-dimensional deterministic representation of atmospheric dynamics at a latitude circle. The slow scale evolution of $X$ is evolved with the nonlinear term $X_{k-1}(X_{k-2}-X_{k+1})$, which mimics advection in a fluid flow. The term $-X_{k}$ represents a linear damping, and $F$ denotes a constant forcing, which for sufficiently large values results in chaotic dynamics. In the two-scale extension, each slow variable $X_{k}$ is associated with a set of fast variables $Y_{j,k}$. The fast variables in turn interact nonlinearly with themselves. This system was studied by Lorenz and Emanuel [34] to determine error growth in ensemble forecasting. For a discussion of how to interpret the parameters we refer to [33].

The two parameter sets considered in this study are identical to those presented by [2] and are summarised in Table 1. Time is expressed in terms of model time units (MTU), where one MTU corresponds to approximately 5 atmospheric days.

Coarse dynamics

The distinction between slow and fast variables renders the L96 system a well-suited testbed for stochastic parametrizations. The coarse dynamics consist of only the slow variables with an additional term $M$ to parametrize the effect of the unresolved fast variables,

[TABLE]

Here, $M$ denotes the parametrization of both the missing unresolved scales $Y$ and adverse numerical effects arising from integrating the coarsened system using a large time step. The parametrization will introduce stochasticity to account for uncertainty.

Data gathering

A high-fidelity data set (referred to as the ‘truth’) is generated by integrating the true dynamics using the three-stage third-order strong stability preserving (SPPRK3) time-stepping scheme [42] with fixed time-step size $\Delta t=0.001$ MTU. The coarse dynamics are simulated with a step size $\Delta t=0.005$ MTU. The stochastic integrals arising from the parametrizations are interpreted in the Stratonovich sense and are integrated using a stochastic SSPRK3 scheme [48], which has strong order of convergence 1/2 in the stochastic setting and reverts to deterministic order 3 in the absence of noise [40].

The high-fidelity solution is initialised by assigning to all variables an independent sample from a standard normal distribution, after which the solution is integrated for 1500 MTU to remove any transitory behavior. The obtained solution at this time serves as the starting point for the data measurement procedure. Subsequently, the truth is simulated for an additional 500 MTU to generate training and verification data.

3.2 Derivative-fitting-based models

Sub-grid scale models can be derived from direct measurements of the sub-grid tendency [2]. This procedure measures the discrepancy between the true and coarse dynamics. Consequently, models derived from this data are based on derivative fitting.

The sub-grid tendency $U$ is approximated as follows,

[TABLE]

and is computed exclusively from the snapshots of the true dynamics. The measurements represent the difference between the evolution of the coarse dynamics and the true dynamics, with the latter approximated by an Euler forward step of size $\Delta t$. Here, $\Delta t=0.005$ in accordance with the coarse time step.

Treating the $X$ variables as ‘spatial’, the measurements are stored a space-time array $\mathbb{U}$. This serves as input for two distinct approaches to stochastic parametrization, referred to as the local model and the global model. Both forecast models are of the form

[TABLE]

where $M_{\mathrm{det}}$ is a deterministic parametrization and $r(X,t)$ denotes a stochastic term, which may be space-dependent.

3.2.1 Local model

The first forecast model employs a third-degree least-squares polynomial fit of $U$ in terms of $X$ and a stochastic model for the residuals. That is, $M_{\mathrm{det}}(X_{k})=b_{0}+b_{1}X_{k}+b_{2}X_{k}^{2}+b_{3}X_{k}^{3}$, as presented by [2]. This parametrization is local, as it is applied to each $X_{k}$ separately, and the model term does not depend on the neighboring values.

For the first and second test cases, we obtain respectively

[TABLE]

The stochastic components of the parametrization follow from the residuals of the polynomial fit. Specifically, two additive stochastic processes are considered, based on the sample variance $S^{2}$ of the residuals and their autocorrelation $\phi$. The first process consists of uncorrelated, normally distributed noise, where the increments over a single time step are drawn from $\mathcal{N}(0,S^{2}\Delta t)$. The second process is an Ornstein–Uhlenbeck (OU) process that takes into account measured autocorrelation of the residuals. It is discretised as a zero-mean first-order autoregressive (AR(1)) process, which evolves as

[TABLE]

Here, $\varepsilon(t)\sim\mathcal{N}(0,1)$ and the scaling term ensures that the resulting samples have the same variance as the uncorrelated model. The value of $\phi$ is set to the autocorrelation of the residuals after $\Delta t=0.005$ MTU. All parameter values are summarised in Table 2.

3.2.2 Global model

The second forecast model is based on global forcing, where the stochastic parametrization is defined as a sum of spatial profiles multiplied by scalar stochastic processes. This approach is motivated by turbulence, which is ultimately a nonlocal phenomenon. The nonlinear advection term in the L96 system ensures that the evolution of each slow variable is intrinsically linked to the state of its neighbors, thereby motivating parametrizations that utilize global basis functions. Here, the spatial profiles are selected as POD modes to capture coherent patterns within the measurements. Alternative approaches can be used analogously, for instance by employing a Fourier basis, but these are not considered in this study.

The spatial profiles are computed from the sub-grid tendency matrix $\mathbb{U}$ in a direct manner. The time-mean is subtracted from $\mathbb{U}$, yielding the average spatial profile $\xi_{0}$ and an anomaly matrix $\mathbb{A}$ with zero-mean rows. The singular value decomposition $\mathbb{A}=B\Sigma V^{T}$ then yields spatial profiles as columns of $B$, the respective squared eigenvalues on the diagonal of $\Sigma$, and the corresponding time series as columns of $V$. The POD modes will be denoted by $\xi_{i}(x),i=1,\ldots,8$, and are of unit norm. The corresponding time series are denoted by $a_{i}$ and have unit variance. Scaling is accounted for via the eigenvalues $\lambda_{i}$. The decomposition of the original signal $U(X,t)$ then takes the form

[TABLE]

The spatial profiles are depicted in Fig. 2, with $k$ treated as the spatial variable. It is observed that spatially coherent patterns can be extracted from the sub-grid scale data.

The time series $a_{i}(t)$ are replaced by independent scalar stochastic processes $r_{i}(t)$ to define the forecast model, analogous to the local model. This results in a sub-grid scale model of the form

[TABLE]

where uncorrelated Gaussian noise, $r_{i}(t)\sim\mathcal{N}(0,\Delta t)$, is the first approach taken. The second approach uses correlated noise, implemented by replacing the process parameter $\phi$ in (13) by the measured autocorrelation of $a_{i}(t)$. An overview of the model parameters is given in Appendix A.

3.3 CRPS-based forecast model

The CRPS-based forecast models are defined in terms of the POD modes with coefficients determined by a learned stochastic process. Here, we consider both additive and multiplicative model forms. Generally, these models can be expressed as

[TABLE]

where $\mathbf{r}_{i}$ is the $i^{\mathrm{th}}$ component of coupled OU processes. The function $c_{i}(X)$ introduces a possibly nonlinear scaling based on the current configuration $X$.

The adopted additive noise model follows from

[TABLE]

and setting $c_{i}(X)=1$. Here $\mathbf{r}\in\mathbb{R}^{K}$ denotes the value of the stochastic process, $\boldsymbol{\mu}\in\mathbb{R}^{K}$ is the mean, $\mathbf{A}\in\mathbb{R}^{K\times K}$ defines the coupling between the different components, $\mathbf{B}\in\mathbb{R}^{K\times q}$ determines the noise scaling, and $\boldsymbol{\varepsilon}\in\mathbb{R}^{q}$ are $q$ independent Brownian motions, that is, $\varepsilon_{i}(t)\sim N(0,1)$. In what follows, we choose $q=K$. The coupled OU process offers greater flexibility compared to the previously used independent scalar OU processes. Namely, interaction between the modes is encoded in $\mathbf{A}$, and incorporating the mean $\boldsymbol{\mu}$ removes the need to include the mean profile $\xi_{0}$ in the model.

The multiplicative noise model is presently defined by setting $c_{i}(X)=\left(a_{i}+b_{i}(\xi_{i}\cdot X)^{2}\right)$ to modulate the noise with the energy of the solution expressed in terms of the POD modes. The parameters $a_{i}$ and $b_{i}$ will be tuned during the training procedure, and allow the learned model to incorporate state-dependent forcing in a straightforward manner.

Modeling the error of a nonlinear dynamical system with a linear stochastic process, such as the coupled OU process (17), is a deliberate simplification. Although the stochastic process itself is linear, its influence on the forecast is nonlinear owing to its integration within the nonlinear system dynamics. This provides a dynamically consistent approach to incorporate the impact of unresolved scales [49]. In particular, the additive coupled OU process provides a tractable and interpretable method for introducing temporally correlated noise. The multiplicative model represents a step towards a more realistic state-dependent parametrization.

In the present study, all model variables are obtained through an optimization procedure detailed in the next subsection. Specifically, no restrictions are imposed on the learnable parameters; the mean $\boldsymbol{\mu}$, the coupling matrix $\mathbf{A}$, the noise scaling $\mathbf{B}$, and the parameters $a_{i}$ and $b_{i}$ are determined without manual intervention. This approach is deemed feasible given the low dimensionality of the noise model, which results in only 160 tunable model parameters in the multiplicative case.

3.3.1 Model training

The model parameters are obtained via an iterative optimization procedure. The high-fidelity data, spanning 500 MTU and stored at intervals of 0.005 MTU, are divided into calibration data and validation sets in an 80/20 split. Denoting the model parameters collectively by $\theta$, and the obtained ensemble prediction by $X(\theta)$, the loss is defined as

[TABLE]

where $C(q)$ is the coarsened reference configuration. The batch size is denoted by $N_{b}$ and depends on the chosen trajectory length $N_{t}$, such that $N_{b}N_{t}$ maintains a comparable magnitude when $N_{t}$ is varied. Furthermore, a batch comprises $N_{b}$ ensembles, each consisting of 20 ensemble members during training. The CRPS in (18) is computed per $X$-variable separately and subsequently summed over all variables to obtain the loss per batch.

The model is trained during 400 epochs, after which the loss per batch was found to have stabilised. Each epoch consists of 125 batches. The Adam optimizer [28] is used with a fixed learning rate set at $10^{-2}$ for training the additive noise model and $10^{-4}$ for the multiplicative noise model. A series of trajectory lengths is used, ranging from $N_{t}=2$ (0.01 MTU) to $N_{t}=200$ (1 MTU). A full description of the batch size, trajectory length, and final measured loss per batch is provided in Appendix A.

The loss per batch during training is presented in Fig. 3 for different values of $N_{t}$. We observe that a short trajectory length for the additive noise model yields only a marginal decrease in the loss. Since each training trajectory starts from a single perfect initial condition, the error and spread after two time steps are small. Hence, for such short time intervals, the initial loss shows is barely improved upon during training. This effect diminishes for longer trajectories, and a decrease of the loss during training is realised until the trajectories become too long. Unstable training behavior is observed for $N_{t}=128$ and $N_{t}=200$, which is attributed to the chaotic dynamics of the system. Specifically, long model trajectories will at some point diverge from the reference trajectories. Minor changes in model parameters will induce significant changes in trajectory due to high system sensitivity, and too much effort is put in correcting long lead-time dynamics. This can lead to unstable gradients and inefficient optimization. This issue could be mitigated by weighting the loss function such that long lead-time predictions contribute less to the loss [36]; however, this is not further investigated here.

The evolution of the loss for the multiplicative noise model differs from that of the additive noise model. We observe that a short trajectory length ($N_{t}=2$) leads to a significant, rapid decrease of the loss, owing to the increased model parameters when compared to the additive noise model. Furthermore, training becomes less effective even at modest trajectory lengths, suggesting higher sensitivity to the model parameters as the trajectory length increases.

4 Model performance

The prediction quality of each model is assessed at short lead times (i.e., ‘weather forecasts’ of the L96 system) and at long lead times (i.e., the ‘climate’ of the system).

The weather forecasts are carried out via ensemble predictions for 2 MTU, where each ensemble comprises 50 members. The slow variables sampled from 100 high-fidelity solution snapshots are used as initial conditions. These snapshots, each separated by 10 MTU, lie outside of the training/validation dataset. We consider perfectly known initial conditions in Section 4.1 and perturbed initial conditions in Section 4.3. The prediction ensembles are assessed via the mean squared error (MSE), the spread-error relationship, and the CRPS. The predicted climate is assessed via a single long-time ensemble prediction, computed over 3000 MTU. The resulting distributions of the slow variables are compared via the Kolmogorov–Smirnov distance and the Hellinger distance. All error assessment metrics are further detailed in Appendix B. A summary of the models tested in this section is given in Table 3.

4.1 Short-time predictions

The mean squared error (MSE) for short lead times is shown in Fig. 4, for both test cases. Qualitatively, the results are consistent across both test cases. The models with independent correlated noise (‘poly_ou’, ‘svd_ou’) are found to perform poorly. This is attributed to poor estimation of the correlation parameters in the autoregressive model. Namely, a visual comparison of the autocorrelation of the fitted autoregressive processes and the measured residuals and POD time series (Fig. 11 in Appendix A) suggests that the fitted temporal correlation parameters are too high. The uncorrelated derivative-fitting-based models (‘poly_gauss’, ‘svd_gauss’) initially exhibit a similar growth of the MSE, but the global model saturates at a lower value than the local model, suggesting that incorporating spatial correlation patterns improves prediction accuracy.

A further increase in average accuracy is observed for the CRPS-based models. Increasing the trajectory length generally leads to parametrizations with improved accuracy, up to a certain limit. For example, the additive noise model trained at $N_{t}=8$ is more accurate than that trained at $N_{t}=2$, and is approximately as accurate as that trained at $N_{t}=32$. The multiplicative models exhibit an even smaller MSE regardless of the trajectory length used in training, suggesting that the additional degrees of freedom allow for increased prediction accuracy.

The average CRPS is shown in Fig. 5. The derivative-fitting-based models with correlated noise show rapid growth of the CRPS, owing to the large mean error of these ensembles. However, this becomes less pronounced when $c=10$, where we observe that the local models (‘poly_gauss’, ‘poly_ou’) perform approximately equally well. The global model with uncorrelated noise exhibits similar evolution of the CRPS as the additive noise model trained at $N_{t}=2$, and generally performs well. Overall, the lowest CRPS is obtained for the additive noise models trained with $N_{t}=8$ and $N_{t}=32$. For $c=4$, the multiplicative noise models show a similar evolution to these additive noise models. In contrast, for $c=10$, the multiplicative noise models yield a higher CRPS on average.

Although the multiplicative models display a smaller MSE (Fig. 4) and smaller spread (Fig. 6) than the other methods, this is not clearly reflected in the CRPS. The former metrics are measured as 2-norms, whereas the CRPS can be regarded as a measure of error and spread in the 1-norm. This suggests that the error and spread of the multiplicative model contains numerous small values that contribute little to the 2-norm but accumulate in their contribution to the 1-norm. In turn, this implies that the corresponding predictions are both accurate and sharp, in line with the training objective specified by the CRPS loss function.

A visual comparison between the ensemble error and spread is provided in Fig. 6. A necessary condition for statistical consistency in ensembles is that the mean ensemble variance estimate equals the squared ensemble mean error [31]. Here, the ensemble error is computed following (22). We observe that all additive noise models are overdispersive, since the ensemble spread is larger than the ensemble error, which is maintained for all lead times. The forecast ensembles obtained with the multiplicative model are found to be underdispersed, implying an overconfident forecast based on the measured error. This is maintained until a lead time of approximately 1.5 MTU for $c=4$ and 1 MTU for $c=10$.

For the trained models, it could be desirable to diminish the weight of the self-spread term in the CRPS loss to possibly favor statistical consistency. This adaptation requires an additional hyperparameter or a term penalizing the discrepancy between the ensemble spread and MSE. Such a modification might diminish the ensemble sharpness during training without sacrificing accuracy. We investigate the use of a hyperparameter to scale the self-spread term in Section 4.4.

4.2 Long-time predictions

The climate prediction using the stochastic parametrizations is shown in Fig. 7, defined via the empirical PDFs of the slow variables over a long simulation time (3000 MTU). Here, the ‘crps_mult8’ is not included in the results, as it exhibited numerical instabilities during extended simulation. Generally, all stochastic parametrizations produce qualitatively similar symmetric PDFs with an excessively large standard deviation. An improvement of the CRPS-based models is found compared to the derivative-fitting-based models, which becomes clearer when increasing the trajectory length. This is summarised in Table 4, indicating that the additive noise models learned using longer trajectory lengths outperform the other tested models. However, the results are not as accurate as those presented in literature (see, e.g., [16], [10]), suggesting that combining CRPS-based learning with other state-dependent model forms might be beneficial for long-time prediction accuracy.

4.3 Short-time predictions with perturbed initial conditions

The ensemble predictions at short lead times presented in Section 4.1 are repeated with perturbed initial conditions to assess the robustness of the stochastic parametrizations. The slow variables constituting the initial conditions for each ensemble are perturbed by a random amount, achieved by sampling $r\sim\mathcal{N}(0,0.1\hskip 0.85358ptS^{2})$. Here $S^{2}$ is the sample variance of the slow variables in the reference data, measured as $S^{2}=6.45$ for $c=4$ and $S^{2}=5.08$ for $c=10$. The perturbations are applied independently per variable and per initial condition, however, they are the same for each model configuration to ensure a fair comparison.

The resulting ensemble error and spread are shown in Fig. 8. The resulting spread values are qualitatively similar to those from the perfect initial condition experiments (Fig. 6). The multiplicative noise model exhibits an increased error growth until approximately 0.75 MTU when compared to the other models. This increased error is attributed to the state-dependent forcing losing accuracy when perturbations are introduced to the initial conditions. Overall, the qualitative similarity observed in the error of the additive models suggests that the learned parametrizations are robust against initial condition uncertainty.

4.4 Sensitivity to self-spread penalty

The previous subsections highlight that the multiplicative noise models generate sharp ensemble forecasts, particularly visible at shorter lead times. This is a consequence of using the CRPS as the loss function, which balances accuracy and ensemble spread, favouring slightly dispersed forecasts that still reproduce the observed value during training. From a statistical perspective, however, forecast reliability requires that the ensemble spread equals the ensemble error [43, 31]. We therefore investigate the role of the spread term in the CRPS loss function and quantify its influence on short-term forecasts.

To this end, we modify the loss function (7) so that the self-spread penalty is scaled by a factor $\alpha$,

[TABLE]

where $\alpha\in[0,1]$. We henceforth refer to (19) as the scaled CRPS loss function. The addition of $\alpha$ allows for training models that minimize the ensemble mean error while reducing the importance of the ensemble spread. The standard CRPS loss function (7) is recovered when $\alpha=1$. The forecast spread term is progressively diminished when $\alpha$ decreases, and vanishes entirely when $\alpha=0$, leaving only the ensemble mean error in the loss function. The value of $\alpha$ cannot exceed 1 to ensure nonnegativity of the scaled CRPS.

The effects of scaling the CRPS loss function are assessed by repeating the model training and forecasting procedures for a sequence of values of $\alpha$. A single model form is selected for the sake of clarity. Here, we adopt the multiplicative noise model and use a trajectory length $N_{t}=8$ to train the model, corresponding to the ‘crps_mult8’ setting considered in the previous subsections. This model was found to produce accurate yet sharp ensemble forecasts, and is therefore suited to demonstrate the ability to change the ensemble spread via choice of the loss function.

The evolution of the training loss is shown in Fig. 9, demonstrating consistent behaviour across the two test cases. Starting from the same initial loss value, the loss decreases at a similar rate for each values of $\alpha$. Smaller final losses are achieved for smaller $\alpha$. This is expected, since a reduced weight on the self-spread term in the scaled CRPS imposes fewer or less stringent constraints during training, thereby facilitating easier optimisation of the loss function.

A comparison between the ensemble error and ensemble spread is provided in Fig. 10. A clear effect on the ensemble spread is observed in both test cases, with the forecast sharpness increasing as $\alpha$ decreases. This occurs because, when $\alpha=1$, the CRPS formulation favours a larger spread to minimise the loss. At the same time, the ensemble error increases with decreasing $\alpha$. This suggests that retaining the self-spread term in the CRPS is beneficial for reducing mean error, despite not directly contributing to the error itself. Larger values of $\alpha$ regularise training by preventing a collapse of ensemble members, thereby reducing overfitting and generally yielding more stable predictions.

For the same reason, the spread observed for $c=10$ and small values of $\alpha$ initially remains small but later exceeds the spread obtained with $\alpha=1$. This behaviour is also attributed to the regularising property of the self-spread term in the loss function. Namely, models trained with small $\alpha$ are more sensitive to perturbations, and their forecast spread is therefore amplified more at later simulation times. These results demonstrate that the ensemble spread and model robustness can be systematically controlled through the choice of loss function, specifically via the weighting of the self-spread term.

5 Conclusion

In this paper, we investigated the calibration of stochastic parametrizations for ensemble forecasting using trajectory learning based on the continuous ranked probability score (CRPS). Using the two-scale Lorenz ’96 system as a testbed, we demonstrated that this approach yields models that provide both accurate and sharp short-term forecasts. Both additive and multiplicative noise models, leveraging global spatial patterns from POD modes, were found to outperform models based on derivative fitting. The improved accuracy at short lead times suggests that CRPS-based trajectory fitting is promising approach when combined with data assimilation, where accurate short-term forecasts are essential. Moreover, the self-spread term in the CRPS was found to have a stabilising effect on the model. Ultimately, this work establishes CRPS-based trajectory learning as a promising methodology for developing data-driven stochastic parametrizations for complex dynamical systems, representing a first step towards its application in geophysical flows.

The stochastic models presented here serve as a proof-of-concept that trajectory learning can be applied to yield accurate ensemble forecasts. A natural next step is to develop additional state-dependent parametrizations using tools such as neural networks, which could capture more complex nonlinear error dynamics. However, many challenges remain. As previously observed in deterministic trajectory fitting, the length of the training trajectory remains a critical hyperparameter. Performance generally improved with longer trajectories, but became unstable at certain lengths due to the chaotic system dynamics. For complex parametrizations, such as those based on neural networks, this sensitivity might become even more pronounced. Future work could explore strategies to mitigate these adverse effects, for instance, by assigning greater weight on short lead times within the loss function to emphasize initial accuracy.

Acknowledgements

SE gratefully acknowledges the Royal Society of Arts and Sciences in Gothenburg (KVVS) for funding the research visit to Imperial College London that started this work. This work was supported in part by the Swedish Research Council (VR) through Grant No. 2022-03453 and by the European Research Council (ERC) Synergy grant “Stochastic Transport in Upper Ocean Dynamics” (STUOD) – DLV-85640.

Appendix A Model and training details

This appendix provides supplementary details for the models presented in the text.

The parameters for the derivative-fitting-based models with correlated noise, introduced in Section 3.2, are summarised in Table 5. These models use increments generated by an Ornstein–Uhlenbeck process, discretised as a first-order autoregressive (AR(1)) model, for the stochastic model component. The specific autocorrelation coefficients per POD mode are provided. Fig. 11 provides a visual comparison between the measured autocorrelation from the data and the autocorrelation function of the AR(1)-processes with fitted parameters. As discussed in Section 4.1, this figure highlights a potential overestimation of the temporal correlation.

The training parameters used for the CRPS-based models, as described in Section 3.3.1, are given in Table 6. A series of trajectory lengths $N_{t}$ were used to investigate the effect on model performance. The corresponding batch sizes $N_{b}$ are chosen to maintain a comparable computational load per epoch. For each value of $N_{t}$, the final loss per batch achieved during training for the additive and multiplicative models is included.

Appendix B Forecast assessment metrics

In the following definitions, we denote the truth by $y$ and denote by $x$ an ensemble forecast of size $M$, consisting of ensemble members $x_{i},\penalty 10000\ i=1,\ldots,M$. All $x_{j}$ and $y$ are regarded as independent identically distributed random variables with mean $\mu$ and variance $\sigma^{2}$.

B.1 Short lead time diagnostics

A general diagnostic of model accuracy is the mean square error (MSE),

[TABLE]

If many ensemble predictions are made, for example by using many different initial conditions, one may instead consider taking the mean value of the MSE. 2. 2.

The spread-error relationship assesses the reliability of the forecast ensemble [43, 31]. A statistically consistent ensemble requires the ensemble spread ensemble and squared mean error to be equal. This is based on the identities of the ensemble mean error and the ensemble variance,

[TABLE]

This metric should ideally be applied after correcting the bias of the forecast, although this is not always feasible in practice. This relation can still be used when the mean forecast error is small compared to the random component of the forecast error. 3. 3.

The Continuous Ranked Probability Score (CRPS), is a proper probabilistic forecast score used to assess the quality of a ensemble forecast. Given a empirical Cumulative Distribution Function (CDF) $F$, one can define the CRPS with respect to the observation $y_{\mathrm{obs}}$ by integrating the CDF as follows

[TABLE]

Here, lower scores indicate better forecast skill. For a discrete forecast, $\{x_{1},...,x_{M}\}$ the CRPS score is approximated as follows,

[TABLE]

Note that $1/2M^{2}$ can be replaced by $1/2M(M-1)$ when used as a estimator. The above CRPS definition is for a single point, and typically is subsequently averaged over points in which there is observations/data, either in time or space. The first term in the CRPS describes the accuracy to the data, and the second term penalizes self-spread as to enforce a sharper forecast.

B.2 Long lead time diagnostics

Diagnostics for long-time simulations generally concern the distribution of predicted variables and statistics of measured time series. Following We denote by $P$ the empirical forecast cumulative distribution function (CDF) and by $P^{\prime}$ the verification (truth) CDF. The empirical probability distribution functions (PDFs) are respectively denoted by $p$ and $p^{\prime}$. The PDFs are approximated by histograms with $M$ bins, where $p_{i}$ ($p^{\prime}_{i}$) denotes the $i^{\mathrm{th}}$ bin of $p$ ($p^{\prime}$). We compare the long-time distributions via the following diagnostics.

The Kolmogorov–Smirnov (KS) distance

[TABLE]

which measures the largest mismatch between two distributions based on their CDFs. It is sensitive to shifts of distributions or spikes, since these yield large local discrepancies. 2. 2.

The Hellinger distance

[TABLE]

which is suited to measure the overlap between distributions defined by their probability density. Compared to the KS distance, it is more sensitive to shape differences in the distributions.