## Abstract

There is a growing interest in developing stochastic schemes for the description of processes that are poorly represented in atmospheric and climate models, in order to increase their variability and reduce the impact of model errors. The use of such noise could however have adverse effects by modifying in undesired ways a certain number of moments of their probability distributions. In this work, the impact of developing a stochastic scheme (based on stochastic averaging) for the ocean is explored in the context of a low-order coupled (deterministic) ocean–atmosphere system. After briefly analysing its variability, its ability in predicting the oceanic flow generated by the coupled system is investigated. Different phases in the error dynamics are found: for short lead times, an initial overdispersion of the ensemble forecast is present while the ensemble mean follows a dynamics reminiscent of the combined amplification of initial condition and model errors for deterministic systems; for longer lead times, a reliable diffusive ensemble spread is observed. These different phases are also found for ensemble-oriented skill measures like the Brier score and the rank histogram. The implications of these features on building stochastic models are then briefly discussed.

## 1. Introduction

Stochastic forcings are currently used in atmospheric and climate models to describe part of the natural variability not accounted for by the usual deterministic parametrizations. This approach allows for introducing a variability at scales at which uncertainties seem to be inherently present (e.g. development of clouds), and at the same time compensating for the under dispersion of ensemble forecasts (e.g. [1]). Beyond the philosophical question raised by the use of such forcings, the potential improvements push us in exploring this approach, in particular for ensemble forecasts.

However, this approach could lead to side effects such as those pointed out in [2], in which it was shown that the use of stochastic white noises could correct part of the model errors present in weather and climate models such as the bias but could have an adverse effect for other moments. This could also have an impact on the predictability of the system, as suggested by the fact that adding noise increases the errors present initially in the forecasts [3].

In this paper, we explore the diverse effects of the presence of stochastic white noises on the variability and predictability of an ocean forced by the atmosphere. The central approach will be to reduce a coupled ocean–atmosphere system into a single ocean model forced stochastically; and once this new simplified model is built, to test its ability to reproduce the correct variability and to forecast the ocean generated by the full deterministic system. The stochastic modelling will follow closely the approach developed in [4–6], known as *stochastic averaging*. It provides a systematic approach for reducing a coupled slow–fast system into a stochastic system describing the dynamics of only the slow component. This approach is obviously very much suitable in the problem of modelling forcings of an ocean as sketched above. The system on which this approach will be tested is a low-order coupled model developed in [7,8], displaying properties with strong connections with the coupled dynamics between the upper layer of the ocean and the atmosphere at mid-latitudes. This model will be briefly presented in §2 and more information can be gathered in [8]. The reduction of the coupled system into a stochastic ocean through stochastic averaging will then be discussed in §3. Finally, the results of this reduction on its variability and predictability will be presented in §4 and a brief discussion of the results provided in §5.

## 2. The coupled low-order model

The classical strategy of building low-order models starting from a set of partial differential equations is to develop each field in a series of orthogonal functions, and truncate this series to the contributions of the dominant orthogonal modes. With the aid of an appropriate scalar product, the partial differential equations can then be projected on these different modes and a finite set of ordinary differential equations can be deduced describing the dynamics of their coefficients. For high truncations, the number of equations is low and the dynamics of such a system can be studied quite extensively. This programme has been followed on many occasions and has led to prominent achievements such as the understanding of the property of sensitivity to initial conditions (e.g. [9–12]).

In [7,8], two versions of a coupled ocean–atmosphere model have been proposed, the first one with 16 variables and the second one with 24. In both versions, the atmospheric model is based on a two-layer quasi-geostrophic flow defined on a β-plane, assuming a dry air, the hydrostatic equilibrium and a flow close to geostrophic balance. These assumptions allow for a reduction of the description of the atmospheric dynamics to one (well-known) pronostic equation for the potential vorticity (e.g. [13]). Although highly simplified, this description was very successful at the early stages of atmospheric modelling and forecasting. These equations are defined over a squared domain, representing an idealized (channel) zone at mid-latitudes. The Fourier series of the barotropic and baroclinic components of the quasi-geostrophic flow are truncated to the first six dominant modes as in [14] and the first 10 modes as in [15,16], compatible with the usual (ad hoc) boundary conditions of a squared domain (periodic in the *x*-direction and no meridional flow at the domain boundaries *y*=0,*L*_{y}). These modes are
2.1

where non-dimensional space variables are introduced, *x*′=*x*/*L* and *y*′=*y*/*L*, where *L* is a characteristic length fixed to *πL*=5000 km, *n*=2*L*_{y}/*L*_{x} is the aspect ratio between the meridional and zonal domain lengths, and *L*_{y}=*πL*. The barotropic, *ψ*=(*ψ*^{u}+*ψ*^{l})/2 (where *ψ*^{u} and *ψ*^{l} are the upper and lower layers streamfunction fields, respectively), and the baroclinic, *θ*=(*ψ*^{u}−*ψ*^{l})/2, parts of the flow are then expressed as
2.2with *K*=6 [14] or 10 [15,16]. In the following, we consider only the case *K*=10. These horizontal fields and time are also appropriately non-dimensionalized (*ψ*′ = *ψ*/(*L*^{2}*f*_{0}), *θ*′ = *θ*/(*L*^{2}*f*_{0}), *t*′=*tf*_{0}, where *f*_{0} is the Coriolis parameter at a latitude of 45^{°}). The pronostic quasi-geostrophic equations (for the barotropic and baroclinic modes) can then be projected on each mode *F*_{i} using the usual scalar product
2.3and one gets a set of 12 or 20 ordinary differential equations for *K*=6 or 10, respectively, listed in [14,15]. A few corrections of the model coefficients were also provided in [16]. Apart from *L* and the aspect ratio *n*, a few other important parameters can be identified, namely the static stability *σ*_{0}, the surface friction coefficient *k*=*k*_{d}/(2*f*_{0}), the internal friction coefficient between the two layers, *k*′=*k*′_{d}/*f*_{0}, the relaxation coefficient towards thermal equilibrium, *h*′′=*h*′_{d}/*f*_{0}, and finally the constant thermal coefficient forcing, *θ**, providing energy to the atmosphere. *θ** only affects the dynamics of the first mode *F*_{1} of the baroclinic component of the system since this forcing is assumed to be of the form, *θ***F*_{1}. The presence of the ocean is felt through the addition of a coupling reflecting that friction is occurring on a moving ocean surface. A term of the form *k*_{d}∇^{2}*Ψ*, where *Ψ* is the streamfunction of the (single layer) oceanic flow, is thus added to the right-hand side of the vorticity equation of the lower atmospheric layer, leading to a surface friction term equal to −*k*_{d}∇^{2}(*ψ*^{l}−*Ψ*). This quantity is also properly non-dimensionalized and projected on the different atmospheric modes using the inner product (2.3). The detailed expressions are given in [8].

The ocean model choice also relies on our strategy of low-order modelling, allowing for discarding many processes that are currently not essential for our purpose (namely the analysis of stochastic physics in multiple time-scale ocean–atmosphere systems), but retaining minimal physics allowing for the development of a coupled ocean–atmosphere dynamics. One therefore focuses on an ocean model based on the reduced gravity quasi-geostrophic shallow water model [17], a single oceanic layer of uniform density in which only momentum transport is allowed. Thermal, salinity and biological effects are neglected. Additional fields (temperature, salinity, etc.) could be incorporated in the future to evaluate their impacts on the coupled dynamics and its predictability. The vorticity equation is also non-dimensionalized as for the atmospheric component, and two additional parameters are present, *γ*=−*L*^{2}/*L*^{2}_{R}, where *L*_{R} is the Rossby radius of deformation, and *r*′ a (non-dimensional) friction coefficient at the bottom of the fluid layer [8].

The atmosphere affects the ocean through a forcing of the form, curl_{z}*τ*, the curl of the wind stress. Assuming that the wind stress is given by (*τ*_{x},*τ*_{y})=*C*(*u*−*U*,*v*−*V*), where *u* and *v* are the horizontal components of the geostrophic wind, −∂*ψ*^{l}/∂*y* and ∂*ψ*^{l}/∂*x*, respectively, and *U* and *V* , the corresponding quantities in the ocean, one gets
2.4Here, the wind stress is proportional to the relative velocity between the flow in the ocean layer and the surface wind. The main coupling parameter is *δ*=*C*/(*ρhf*_{0}), where *ρ* is the water density, *h* the depth of the fluid layer and *C* a friction coefficient whose expression is given below.

Following [18], we use the following set of modes for the description of the streamfunction of the horizontal velocity field in the ocean:
2.5in order to get the free-slip boundary conditions in the domain over which the flow is defined at *x*=0,2*π*/*n* and *y*=0,*π*. The non-dimensional ocean streamfunction field is then defined as
2.6In addition, a specific inner product is adopted for the oceanic model as in [18],
2.7Using the same procedure as for the equations of the atmospheric component, the dynamical vorticity equation is projected using (2.7) on the different modes *ϕ*_{i}. The resulting set of ordinary differential equations of the ocean model are explicitly given in [8] together with the wind stress forcing terms. For convenience, we rewrite these equations in a more synthetic form as follows:
2.8where
2.9where are the coefficients of the non-dimensional streamfunction of the lower layer of the atmosphere. We refer to the study of Vannitsem & De Cruz [8] for the details on the coefficients. Note that the forcing *f*(*i*) does only explicitly depend on the atmospheric variables, but the latter are in turn influenced by the ocean flow field (*A*_{i}’s). These forcings are therefore *a priori* also implicitly dependent on the *A*_{i}’s. This aspect will be investigated further when modelling this forcing using stochastic processes in §3.

The estimation of the main physical parameters is made as follows. For the atmosphere, the parameter *k* is related to the surface drag felt by the lower layer of the two-layer QG model. This is estimated based on the Ekman layer theory [17], p. 115 as
2.10after dividing by *f*_{0}, and where *H* and *d* are the thickness of the lower atmospheric layer and the thickness of the Ekman surface layer, respectively. Typically, *H* is of the order of 5000 m and *d* of the order of 100–1000 m. This implies that *k* is falling in a range of [0.01,0.1]. Here, the value is fixed to *k*=0.02 (and the other dissipation parameters are fixed to *h*′′=*k*′=2*k*).

For parameter *δ*, one can use the estimate of [19]. The dimensional forcing coefficient is given by
2.11where *ρ*_{a} and *ρ*_{o} are the densities of the air and of the sea water, respectively, *h* is the thickness of the ocean layer and *C*_{D} the surface friction coefficient. With *C*_{D}≈0.001, *h*≈20–500 m, |*V* |≈5–10 m s^{−1}, *ρ*_{a}≈1 kg m^{−3} and *ρ*_{o}≈1000 kg m^{−3}, one gets values (once normalized by *f*_{0}) in the range [0.0001,0.01]. Note that *C* in (2.4) is equivalent to *C*=|*V* |*ρ*_{a}*C*_{D}. *δ* is fixed to 0.0019305 in the following.

For the evaluation of the thermal forcing *θ**, the same approach as in [14,15] is adopted, through the use of the thermal wind relation
2.12where *ΔT** is the climatological temperature difference between the meridional limits at *y*=0 and *y*=*L*_{y}, and *R* is the ideal gas constant. For realistic climatological differences, *θ** should be of the order of 0.1–0.2. All non-dimensional parameter values that are used in this work (except when explicitly stated in the text or figure captions) are provided in table 1.

Figure 1 displays the Lyapunov spectrum characterizing the instability properties of the flow (the 24 Lyapunov exponents *σ*_{i} are ranked in decreasing order and given in day^{−1}), as obtained with the parameters *θ**=0.14, *n*=1.5 and *δ*=0.0019305 after an integration of about 10^{6} days. This spectrum is computed using a standard algorithm [20]. Three positive exponents are obtained indicating the chaotic nature of the dynamics of the coupled system. The fourth one is very close to 0 with an amplitude 100 times smaller than the third and the fifth exponents.

## 3. The reduced deterministic and stochastic models

Stochastic modelling has long been recognized as a way to represent sufficiently random processes acting on fast time scales and affecting the slow time-scale variables of interest [21]. The reduction of the fast components of a coupled fast–slow system to stochastic forcing is a relatively active subject nowadays and has been explored and developed by different groups [4,5,22–25]. In this paper, we follow the approach proposed by Arnold *et al.* [5], and subsequently applied in [6,22] in the context of different climate models. This section will closely follow the description of the stochastic averaging approach settled in [5,6].

The approach consists in assuming that there is a time-scale separation between the slow and fast time scales. Let us consider that the system under investigation is described by a system of equations
3.1where **x** and **y** are, say, the ocean and atmospheric variables, respectively. The fast time-scale processes are described by the *y*-variables. Applying the stochastic averaging consists in averaging out the *y*-variables of the first set of equations for a given **x** through
3.2where *p*(**y**|**x**) is the conditional probability of having **y** given **x**.

One then gets a first approximation for the dynamics of the slow system as
3.3where is the slow variable, and this approximation was referred to as (A) approximation in [5,6]. In the following, we will refer to it as the Deterministic Model (DetM). Note that this approximation is valid for infinite separation time scales, and of course this limit is never reached in real applications. There are therefore corrections to this approximation that can be expressed as a white noise as discussed in detail in [5]. The full stochastic representation of the slow system is given by the following Ito stochastic differential equation:
3.4where *W* is a vector of Wiener processes, **G**(**x**) the noise-induced drift term and ** σ**(

**x**), the covariance function of the Wiener processes. This full approximation was named (N+) by Arnold

*et al.*[5] and is referred in the following as the Stochastic model (StM). The noise-induced drift is present only when the stochastic forcing effectively depends on the slow variables themselves,

**G**(

**x**)=(1/4)∂

_{x}

**(**

*σ***x**). The covariance matrix

*σ*^{2}(

**x**) is estimated through 3.5where the entries of the matrix

**C**

_{S′S′}(

*τ*;

**x**) are 3.6where , and the overbar refers to the averaging (3.2). The indices

*i*,

*j*are spanning the phase space dimension of the slow variables,

**x**.

In the context of the analysis of the two-scale system presented in §2, we will further assume that
3.7simplifying the development of the stochastic model. This choice is motivated by the typical fast decorrelation found in chaotic systems, and its validity can be checked *a posteriori* by looking at the actual variation of [**C**_{S′S′}(*τ*;**x**)]_{ij}.

In the coupled model described in §2, the slow component is the ocean (**x**=(*A*_{1},…,*A*_{4})), forced by the fast atmospheric component (**y**=(*ψ*_{1},…,*ψ*_{10},*θ*_{1},…,*θ*_{10})) through wind forcing (2.9). This wind forcing will be modelled based on the approach just outlined above through stochastic forcings. Here, the Deterministic approximation (3.3) describing the dynamics of the slow variables is given by (2.8) in which the forcings (2.9) are averaged in time. *S*_{i}′ appearing in the stochastic terms of (3.6) are the fluctuations around the mean values, .

The dependence of the covariance matrix given in (3.7) on the oceanic modes, together with the decorrelation time scale *ξ*_{i,j}, should be first evaluated. The parameter *ξ*_{i,j} has been estimated through the analysis of several correlation functions of the *f*_{i}′ given in equation (2.9) based on a model integration of 10^{6} days. This parameter has been found to vary from 0.5 to 2 days depending on the specific coefficients [**C**_{S′S′}(*τ*;**x**)]_{ij} investigated, for the model parameters of the coupled system used in figure 1. Such values are typical of the decorrelation time scale that can be found at synoptic scales in the atmosphere.

In order to figure out what is the impact of the choice of the parameter *ξ*_{i,j} on the properties of the stochastic model, we have developed two stochastic model. A first version (StM1) uses a rough estimate of *ξ*_{i,j}=*ξ*=0.6 days for all the coefficients, and a second version (StM2) based on a separate tuning of these parameters for the four different entries [**C**_{S′S′}(*τ*;**x**)]_{ii}. These values are *ξ*_{1,1}=*ξ*_{1}=1.5, *ξ*_{2,2}=*ξ*_{2}=0.7, *ξ*_{3,3}=*ξ*_{3}=1.7, *ξ*_{4,4}=*ξ*_{4}=0.7 days. For the dominant covariance entries (as shown in figure 2) with *i*≠*j*, the parameters are fixed such that *ξ*_{1,3}=*ξ*_{1}=1.5, *ξ*_{2,4}=0.7 days. For the other entries of smaller amplitudes, the parameters have been arbitrarily fixed to *ξ*_{1,2}=*ξ*_{1,4}=*ξ*_{1}, *ξ*_{2,3}=*ξ*_{2}, and *ξ*_{3,4}=*ξ*_{3} days.

The dependence of the covariance matrix on the variables *A*_{i} has also been explored using an integration of the coupled model of 10^{6} days. Figure 2*a*,*b* displays this dependence in *A*_{2} and *A*_{4} for the dominant entries of the covariance matrix, respectively (averaged over all the other configurations of *A*_{i} with *i*≠2 or 4, respectively). There is no substantial dependence of the entries of the covariance matrix as a function of *A*_{i} (although, as discussed in §2, we could expect some dependences through the influence of the ocean on the atmosphere). One therefore considers that [**C**_{S′S′}(0;**x**)]_{ij} is independent of **x**. The drift term can then be assumed to be equal to 0 and a constant covariance matrix *σ*^{2} will be used. The square root matrix ** σ** is then obtained through an SVD decomposition. The entries of

**for**

*σ**θ**=0.14 and

*δ*=0.0019305 are provided in table 2, together with the mean forcings .

It is worth emphasizing that the stochastic forcing adopted here is an additive white noise. This feature differs to the usual findings that the noise representing the sub-grid scale processes is multiplicative (e.g. [26–28]), but it clearly provides a good representation of the impact of the atmosphere on the ocean in the present idealized context. A similar finding was also found in the context of sub-grid scale modelling of a quasi-geostrophic model in [29], indicating that additive noises could be of a broader interest for sub-grid scale or boundary forcing modelling. Finally, the mean forcing is quite weak compared with the standard deviation of the fluctuations. This already suggests that the averaged system is weakly forced by the atmosphere, and the deterministic approximation should therefore be stable. This feature will be discussed in §4.

## 4. Dynamics and predictability

The stochastic model is integrated using the Heun method as the coupled system (e.g. [30]), taking into account the stochastic nature of the forcing by multiplying the stochastic square root matrix *σ* by , where *N*(0,1) is a a random draw from a Gaussian distribution of mean 0 and variance 1 and *Δt*, the time step of the integration.

Figure 3 displays the temporal evolution of (a) *A*_{1} and (b) *A*_{2} for the coupled system, DetM, (red) and the stochastic ocean model, StM2, (green). The temporal evolution of both systems displays strong similarities that precludes the possibility to distinguish which trajectory is coming from the coupled system or the stochastic model. Note that for the deterministic system DetM, the solution converges towards a stationary solution (*A*_{1}=22883, *A*_{2}=−3202, *A*_{3}=−19052, *A*_{4}=20574 m^{2} *s*^{−1}). The stochastic trajectory partly displayed in figure 3 is oscillating around this stationary solution. Drawing on this encouraging result, one can now look at more involved quantities. The power spectra of the mode *A*_{1} generated by both systems are shown in figure 4. Very similar spectral structures are found with a power-law (red) decay from very large up to small time scales. This feature is reminiscent of the behaviour already discovered by Hasselman [21] and afterward revisited by Arnold [4]. Interestingly, the red spectrum is present up to scale of the order of 6000 days (for the series investigated).

Obviously, the coupled model does not incorporate additional effects associated with the presence of a deep ocean component. This affects the structure of the spectrum on very long time scales (centennial), and one cannot therefore expect to match real ocean spectra on such time scales.

Let us now inspect in details the probability density of the different variables generated by the two versions of the stochastic model (StM1 and StM2, see §3) when compared with the full coupled system. Figure 5 displays the probability density of the different ocean variables, *A*_{1},…,*A*_{4}, obtained with the full coupled system and with the two stochastic model versions. StM1 is displaying a too large variability for variables *A*_{1} and *A*_{3}. For the second model version StM2, the probability densities are much closer to the reference system ones. This makes of this second model version a very good candidate for the description of the ocean component of the coupled system.

Let us now turn to the error dynamics. In order to analyse the use of a stochastic model in a forecasting context, we set up the following experiment. First assume that all the variables of the forecasting system are affected by an initial condition error drawn from a Gaussian distribution with mean 0 and variance 10^{−12} (in adimensional units). For the ocean variables, this amounts to considering a standard deviation of the error of the order of 1000th of the climatological standard deviations of the model variables. These errors are then superimposed on the chaotic trajectory generated by the coupled system. A set of states are then selected in this series and are used as initial conditions of the stochastic ocean model and of its deterministic approximation.

For reference, we first plot in figure 6*a*, the mean squared error (m.s.e.) for the oceanic modes in the case where the forecasting model is considered as perfect (equivalent to the full coupled system) and only affected by initial random errors with variance 10^{−12}. The curves have been obtained using 5000 realizations starting from different initial conditions on the reference system’s attractor of this twin experiment. The three curves refer to different values of the coupling parameter, *δ*=0.0004826,0.0019305,0.009653. For all these curves, the error first slowly decreases during a period whose duration is directly related to *δ*. Then after about 10 days, the error rapidly amplifies (due to the large errors present in the atmosphere after a few days). A third phase then emerges during which the error follows a diffusive dynamics (the m.s.e. is proportional to time) for a long period lasting up to 1000 days. Finally, the saturation occurs associated with the finite size of the attractor. The third phase of diffusion is the regime for which extended forecasts are obviously possible.

Let us now consider the stochastic and the deterministic approximation described in §3 as forecasting models. Both models are integrated in time starting from the same initial conditions sampled on the attractor of the reference system, and affected by the Gaussian initial error mentioned earlier. In order to emphasize the importance of using stochastic forcings, a set of 100 ensemble members are used for each initial state. For each member a different realization of the stochastic forcing is used. Figure 6*b* displays the m.s.e. of the ensemble mean (averaged over the 100 members) of stochastic model version StM1, green, together with the m.s.e. of the ensemble mean of the deterministic approximation model version, dark blue. Clearly, both models provide the same result for the ensemble mean dynamics. It is also interesting to note that at short time scales (up to a few days) the m.s.e. of the ensemble mean is increasing following a *t*^{2} evolution (after a short exponential-like behaviour associated with the presence of initial condition errors). This evolution is reminiscent of the model error dynamics of deterministic systems as discussed in [3,12], which takes over the initial condition errors after a certain period. For long times, the m.s.e. of the ensemble mean saturates at a value about twice smaller than the reference m.s.e. of randomly chosen states on the attractor of the reference system. The latter feature is related to the fact that two different trajectories (of the ‘perfect’ coupled system, *z* and *w*) are becoming completely decorrelated for long lead times, implying that the m.s.e. 〈(*z*−*w*)^{2}〉=〈(*z*−〈*z*〉)^{2}〉+〈(*w*−〈*w*〉)^{2}〉=2〈(*z*−〈*z*〉)^{2}〉 with 〈*w*〉=〈*z*〉. When comparing this quantity to the evolution of the m.s.e. of the ensemble mean of a perfect ensemble, one should therefore get a factor of 2 asymptotically. In this experiment, StM and DetM are not perfect model versions of the coupled system, and one should expect some deviation to the latter relation; however, it seems to be well respected owing to the fact that the climatological averages of the model versions are close to the one of the reference system.

On the same panel, the averaged variances around the mean of the ensemble forecasts (referred as the spread in the figures) are displayed for the StM1 and the DetM models. For the latter, the variance is decreasing reflecting the stable nature of the dynamics generated by the deterministic approximation. For the former, the variance is increasing diffusively (linear in time) and after about a few days displays an amplitude around 1.7 larger than the m.s.e. of the ensemble mean. As mentioned in [31], the equality of the m.s.e. of the ensemble mean and of the ensemble variance is the signature of the use of a correct variability for the ensemble, which is then qualified as reliable. For this stochastic model version (StM1), the ensemble variance is too large when compared with the m.s.e. of the ensemble mean for all lead times, leading to an overdispersed ensemble forecast. This suggests that the stochastic forcing imposed in this first version is too large.

For the second version of the stochastic model (StM2) displayed in figure 6*c*, the results are better since both curves (m.s.e. of the ensemble mean and variance around the ensemble mean) are very close to each other after about a few days, reflecting the use of a reliable ensemble for intermediate and long lead times. For short lead times, the overdispersion persists. A similar picture is found when the initial condition error amplitude is increased (but to a lesser extend), as illustrated in figure 6*d*. In this case, the role of the initial condition error on the ensemble mean error dynamics persists during a longer period before the rapid amplification associated with the modelling error.

The feature at short lead times is related to the specific behaviour of the ensemble mean, displaying the typical error evolution of a purely deterministic system. In order to understand this behaviour, it is worth recalling the dynamics of the ensemble mean of a system. Following [2], one can write the dynamics of the ensemble mean, **m**, for a system with quadratic nonlinearities as
4.1where *V* _{j,k}=〈*δx*_{j}*δx*_{k}〉, and *δx*_{i} represents the error relative to the ensemble mean along direction *i*. If the initial error variance is small, the dynamics of the ensemble mean will mainly follow the dynamics defined by the first term on the right-hand side of equation (4.1), namely a purely deterministic dynamics. The stochastic noise will enter into play for longer lead times through the second and higher order moments of the error distribution.

To summarize, the stochastic forcing is here overdispersive for short times but necessary for a correct evaluation of the skill of the ocean at long lead times. It is also worth pointing out that the variance of the ensemble follows a linear evolution, again suggesting the diffusive nature of the dynamics of the model (and of the ocean fields generated by the coupled system after the transient impact of the atmosphere). This behaviour is typical of the error dynamics of stochastic systems with additive forcings as discussed in [32].

In order to evaluate the quality of the ensemble forecasts, two other important quantities have been computed, the Brier score and the rank histogram (e.g. [33]). The Brier score represents the m.s.e. between the predicted probability of occurrence of a certain event and its actual realization (0 or 1). Figure 7 displays the results for the four variables of the ocean, separately, for the two stochastic model versions discussed in §3 and for the deterministic model version. The event considered here is the crossing of a threshold of the variable’s amplitudes equal to the climatological mean + 1 s.d. A clear improvement is obtained with the stochastic model when compared with the DetM model. It is even better when the decorrelation parameters *ξ*_{i,j} are well tuned. The brier score is saturating (as the m.s.e.) after about 1000 days indicating the potential for long range forecasts. As for the error dynamics, two regimes emerge before and after about a few days, typical time scale before which the ocean is feeling the presence of a (determinstic) coloured atmosphere [8]. The picture is similar for other thresholds.

The rank histogram provides an information on the reliability of an ensemble forecast. It is built by computing the rank of the actual observation (our reference trajectory) among each ensemble forecast (and at each lead time) and to accumulate their frequency of occurrence as a function of their rank within the ensemble. For a reliable ensemble with an infinite number of realizations, this frequency should be equal 1/(*N*+1), where *N* is the number of ensemble members. This further means that the observation can be considered as an equally likely member of the ensemble.

Figure 8 displays the rank histogram for the model StM2 as obtained with a number of 5000 different initial conditions on the attractor of the reference system (i.e. realizations) and a number of ensemble members *N*=100. Figure 8*a*–*d* corresponds to a lead time of 1.01 days, figure 8*e*–*h* to a lead time of 11.1 days and figure 8*i*–*l* to a lead time of 112 days. These lead times were chosen arbitrarily within the different regimes already mentioned for the error dynamics in figures 6 and 7. At shortlead times, the ensemble is overdispersive as reflected by a histogram weakly populated on both right and left sides. This is in agreement with the overdispersion pointed out with the analysis of the m.s.e. in figure 6. Interestingly, for variables *A*_{1} and *A*_{3}, a double peak structure is emerging, indicating that observations overpopulate intermediate ranges of values within the ensemble. This reflects the drastic difference at such a time scale between the actual (deterministic or coloured) forcing of the ocean in the coupled system and the Gaussian white noise approximation used for stochastic modelling. At longer lead times the histogram is flattening but still keeping a complicated structure up to about 10 days. The flat rank histogram for long lead times also reveals the quality of the stochastic modelling associated with the model version StM2 for long range forecasts of the ocean.

## 5. Discussion

Stochastic forcings play an increasing role in modelling unresolved processes in atmospheric and climate models. Adding noise is however far from obvious from a theoretical point of view, since important assumptions should be made from the beginning such as the nature of the noise (e.g. Gaussian white, see [4]) and its correspondence to the real dynamics, and the assumption of a clear separation of scales between the resolved and unresolved processes. These simplifications could have a strong impact on the statistical and dynamical properties of the system under investigation, such as the moments of the distribution generated by the stochastic model or the skill of the model in predicting the evolution of the underlying dynamics.

These aspects have been explored here through the stochastic modelling (based on stochastic averaging [4–6]) of the ocean in the context of a (low-order) coupled ocean–atmosphere system. The results were surprisingly good, except for the short forecasting time scales. The variability is well reproduced and the long-term probabilistic forecasts are reliable, with a diffusive dynamics of the error for long forecasting time scales. At short lead times, the picture is different with an error dynamics of the ensemble mean reminiscent of the one of deterministic systems in the presence of model errors, while the ensemble variance is behaving diffusively (≈*t*). This discrepancy is inherited from the fact that the ensemble mean is behaving primarily in a deterministic fashion (as the reference ocean generated by the coupled model), while the ensemble spread is diffusive.

These two typical time scales at short and long lead times were also isolated in the analysis of the evolution of the Brier score and of the rank histogram. In both cases, the deterministic (or coloured) nature of the actual atmospheric forcing has a clear impact on the skill for short lead times.

The specific structure of the stochastic forcing (table 2) also suggests that the dominant contribution to the dynamics of the ocean layer is the natural variability of the wind forcing of this idealized system. It is therefore expected that once bifurcation analyses are performed on such ocean models for which the bifurcation parameter is the (mean) strength of the wind forcing, one could get potential solutions far from the real underlying dynamics.

The problem of model reduction investigated here can be extended in different ways, first by exploring other approaches to build the stochastic forcings (e.g. [2]) and compare their respective advantages, second by applying stochastic averaging when part of the atmosphere is resolved. In the latter case, the separation of scales will not be met anymore as in [22]. Finally, this work can be extended by applying this approach to build a realistic forcing for a state-of-the-art ocean model such as the one of [34], and investigating the intrinsic limits of skill of such models.

## Funding statement

This work is partially supported by the Belgian Federal Science Policy Office under contracts SD/CA/04A and BR/121/A2/STOCHCLIM.

## Acknowledgements

The constructive comments of the reviewers are gratefully acknowledged.

## Footnotes

One contribution of 14 to a Theme Issue ‘Stochastic modelling and energy-efficient computing for weather and climate prediction’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.