A generic approach to explicit simulation of uncertainty in the NEMO ocean model

. In this paper, a generic implementation approach is presented, with the aim of transforming a deterministic ocean model (like NEMO) into a probabilistic model. With this approach, several kinds of stochastic parameterizations are implemented to simulate the non-deterministic effect of unresolved processes, unresolved scales and unresolved diversity. The method is illustrated with three applications, showing that uncertainties can produce a major effect in the circulation model, in the ecosystem model, and in the sea ice model. These examples show that uncertainties can produce an important effect in the simulations, strongly modifying the dynamical behaviour of these three components of ocean systems.


Introduction
The first requirement of an ocean model is the definition of the system that the model is going to represent. As illustrated in Fig. 1, this usually amounts to defining an appropriate separation between the system (A) and the environment (B). For instance, in this study, we always use a stand-alone ocean model, which means that the atmosphere is not included in the system (A), but in the environment (B). A key property of any ocean model is also the separation between resolved scales (in A) and unresolved scales (in B), defining the spectral window that the model is going to represent. In a similar way, marine ecosystems are too complex to be entirely included in A. They can only be represented by a limited number of variables C i , i = 1, . . ., n, providing a synthetic pic-ture of the ecosystem, while the remaining biogeochemical diversity is included in B.
Even if the union of the two systems A and B could be assumed deterministic, this is in general not true for system A alone. The future evolution of A does not only depend on its own dynamics and initial condition, but also on the interactions between A and B. This means that the only two ways of obtaining a deterministic model for A are either to assume that the evolution of B is known (as is usually done for the atmosphere in stand-alone ocean models) or to assume that the effect of B can be parameterized as a function of what happens in A (as is usually done for unresolved scales and unresolved diversity). It is however important to recognize that this is always an approximation and that B is often an important source of uncertainty in the predictions made for A.
To obtain a reliable predictive model for A (in the sense given in Brier, 1950, andToth et al., 2003), a consistent description of this uncertainty should be embedded in the model itself. This transforms the deterministic model into a probabilistic model, which fully characterizes the quantity of information that the model contains about A. Two important advantages of this probabilistic approach are (i) to allow objective statistical comparison between model and observations (by providing sufficient conditions to invalidate the model; see for instance Candille and Talagrand, 2005), and (ii) to provide a coherent description of model uncertainty to data assimilation systems. The objective of the modeller also changes: instead of designing a deterministic model as close as possible to observations, a probabilistic model that Figure 1. Schematic of the separation between resolved and unresolved processes (systems A and B). Even if A ∪ B can be assumed deterministic, system A alone is not deterministic in general, because of the interactions with system B.
is both reliable (not yet invalidated by observations) and as informative as possible about A must be designed.
In practice, for a complex system, it is usually impossible to compute explicitly the probability distribution describing the forecast. In general, only a limited size sample of the distribution can be obtained through an ensemble of model simulations, as is routinely done in any ensemble data assimilation system (see Evensen, 1994). Ensemble simulations are produced by randomly sampling the various kinds of uncertainty (in the dynamical laws, in the forcing, in the parameters, in the initial conditions, etc.) in their respective probability distribution (Monte Carlo simulations). To allow objective comparison with observations or to correctly deal with model uncertainties in data assimilation problems, nondeterministic models are thus needed in many ocean applications. The most direct approach to introducing an appropriate level of randomness in ocean models is to use stochastic processes to mimic the effect of uncertainties. In the discussion above (summarized in Fig. 1), a specific focus was given to uncertainties resulting from the effect that unresolved processes (in B) produce on the system (A). However, there is a variety of other sources of uncertainty in ocean models (e.g. numerical schemes, machine accuracy, etc.) that do not enter this particular sketch, and that may also require a stochastic approach (Palmer et al., 2014).
Stochastic parameterizations explicitly simulating model uncertainty were first applied to ensemble weather forecasting by Buizza et al. (1999) about 15 years ago. Since then, stochastic parameterizations have emerged as a quickly developing area of research in meteorology (Palmer, 2001;Palmer et al., 2005). In oceanography, however, most stateof-the-art dynamical models are still deterministic. Up to now, the development of stochastic dynamical equations has been mainly focused on stochastic parameterization of Reynolds stresses in idealized ocean modelling systems (see Frederiksen et al., 2012a, andKitsios et al., 2013, for a review). Only a few exploratory studies have attempted to explicitly simulate uncertainties in realistic dynamical ocean models: this has been done for the ocean circulation (Brankart, 2013), for the ocean ecosystem (Arhonditsis et al., 2008), and for the sea ice dynamics (Juricke et al., 2013). These preliminary studies nonetheless already show that uncertainties can play a major role in dominant dynamical behaviours of marine systems.
In line with these studies, the objective of this paper is to propose a generic implementation of these stochastic parameterizations, and to investigate several applications in which the randomness of the ocean system may be an important issue. This is synthetically implemented in the ocean model (see Sect. 2) by adding one additional module providing appropriate random processes to any non-deterministic component of the system (circulation, ecosystem, sea ice). The method is designed to be simple enough to allow a quick check of the effect of uncertainties in the system, and flexible enough to apply to various sources of uncertainty (atmosphere, unresolved scales, unresolved diversity, etc.). Three applications are then illustrated in Sect. 3, showing that the explicit simulation of uncertainty can be important in a wide variety of ocean systems, by stimulating important nondeterministic dynamical behaviours. The first application (circulation model) is the same application as in Brankart (2013), but this previous paper only presented the average effect of the stochastic parameterization, whereas the focus is here on the randomness that is produced in the largescale ocean circulation. The second application (ecosystem model) is a first attempt to apply stochastic parameterizations and to explicitly simulate randomness in a basin-scale ocean ecosystem model. The third application (sea ice model) is an attempt to reproduce the parameterization developed in Juricke et al. (2013) in our ocean model using the generic implementation presented in Sect. 2, and to illustrate the randomness that is generated in the interannual variability of sea ice thickness.

Stochastic formulation of NEMO
The ocean model used in this study is NEMO (Nucleus for a European Model of the Ocean), as described in Madec et al. (2008). NEMO is the European modelling framework for oceanographic research, operational oceanography, seasonal forecast and climate studies. This model system embeds various model components (see http://www.nemo-ocean. eu/), including a circulation model (OPA, Océan PArallélisé), ecosystem models, with various levels of complexity (e.g. LOBSTER, LOCEAN Simulation Tool for Ecosystem and Resources), and a sea ice model (LIM, Louvain-la-Neuve Ice Model). The purpose of this section is to shortly describe the three kinds of stochastic parameterizations that have been implemented in NEMO, and to show that, from a technical point of view, they can be unified in one single new module in NEMO, feeding the various sources of

Order n autoregressive processes
The starting point of our implementation of stochastic parameterizations in NEMO is to observe that many existing parameterizations are based on autoregressive processes, which are used as a basic source of randomness to transform a deterministic model into a probabilistic model. A generic approach is thus to add one single new module in NEMO, generating processes with appropriate statistics to simulate each kind of uncertainty in the model (see examples in Sect. 3).
In practice, at every model grid point, independent Gaussian autoregressive processes ξ (i) , i = 1, . . ., m are first generated using the same basic equation: where k is the index of the model time step; and a (i) , b (i) , c (i) are parameters defining the mean (µ (i) ), SD (σ (i) ) and correlation timescale (τ (i) ) of each process: -for order 1 processes (AR(1)), w (i) is a Gaussian white noise, with zero mean and SD equal to 1, and the parameters a (i) , b (i) , and c (i) are given by -for order n > 1 processes (AR(n)), w (i) is an order n−1 autoregressive process, with zero mean, and SD equal to σ (i) , correlation timescale equal to τ (i) , and the parameters a (i) , b (i) , and c (i) are given by In this way, higher-order processes can be easily generated recursively using the same piece of code implementing Eq. (1), and using successively processes from order 0 to n − 1 as w (i) . The parameters in Eq. (3) are computed so that this recursive application of Eq. (1) leads to processes with the required SD and correlation timescale, with the additional condition that the n − 1 first derivatives of the autocorrelation function are equal to zero at t = 0, so that the resulting processes become smoother and smoother as n is increased. AR(2) processes (with other specifications) have already been applied in several studies (Berloff, 2005;Wilks, 2005), and will be used in this paper in the sea ice model application (see Sect. 3.3). Second, a spatial dependence between the processes can easily be introduced by applying a spatial filter to the ξ (i) . This can be done either by applying a simple filter window to the ξ (i) 2-D or 3-D matrices ξ (i) = F[ξ (i) ], or by solving an elliptic equation: L[ ξ (i) ] = ξ (i) . In both cases, the filtering operator could be made flow dependent, or more generally, the filter characteristics could be modified according to anything that is resolved by the ocean model (in system A in Fig. 1). Technically, this only requires that the description of the ocean model is made available to the filtering routines. This filtering option (using a simple Laplacian filter) is used in the sea ice application (see Sect. 3.3).
Third, the marginal distribution of the stochastic processes can also be easily modified by applying a nonlinear change of variable (anamorphosis transformation) to the ξ (i) before using them in the modelξ (i) = T [ξ (i) ]. This idea is similar to what is done in ensemble data assimilation methods to transform variables with non-Gaussian marginal distribution into Gaussian variables (Bertino et al., 2003;Béal et al., 2010;Brankart et al., 2012). For instance, this method can be very useful if the description of uncertainties in the model requires positive random numbers. In this case, anamorphosis transformation can be applied to transform the Gaussian ξ (i) into positiveξ (i) with lognormal or gamma distribution. This anamorphosis option (using a gamma distribution) is used in the sea ice application (see Sect. 3.3).
Overall, this method provides quite a simple and generic way of generating a wide class of stochastic processes. However, this also means that new model parameters are needed to specify each of these stochastic processes. As in any parameterization of lacking physics, a very important issue is then to tune these new parameters using either first principles, model simulations, or real-world observations. This key problem of assessing the parameters involved in Eq. (1) cannot be addressed in the present paper, and we can only provide a very brief overview of the nature of the problem. Many existing studies (e.g. Frederiksen et al., 2012b;Achatz et al., 2013;Grooms and Majda, 2013) already addressed the problem of choosing the coefficients of the AR(n) processes to simulate the Reynolds stresses in atmospheric and oceanic flows. Considerable progress has been made for this important problem, but not all unresolved processes have received so much attention, and it is often still difficult to figure out how to derive the parameters of the AR(n) processes.
Referring to the sketch presented in Fig. 1, the general idea to tune the parameters is to obtain reliable probabilistic information on what happens in system B, and to reduce this information to a simple statistical model (e.g. the autoregressive model described above). More precisely, the probability distribution simulating the effect of B should also be conditioned on what happens in system A. For instance, it can be very important that the probability distribution for the state of the atmosphere (e.g. surface winds) be conditioned on the state of the ocean model (e.g. mesoscale eddies), to simulate the interaction between A and B. Similarly, the probability distribution for unresolved scales or unresolved diversity usually depends on what happens in system A. This need to correctly simulate conditional probability distributions ex-plains why the tuning of the parameters is not easy, and why an extensive database to learn the statistical behaviour of the coupling between A and B is often necessary. In practice, this learning information can be obtained either from observations of the two systems or from other models explicitly simulating the coupling between A and B. For instance, high-resolution observations or high-resolution models can be used to tune a statistical model for unresolved scales; a model of the atmospheric boundary layer can be used to learn the statistical dependence of the state of the atmosphere on the ocean conditions; a generic biogeochemical model involving a large number of species can be used to understand the statistical effect that unresolved diversity can produce in a simple ecosystem model.
The identification of an appropriate statistical model is thus an important intermediate step that is far from straightforward, and for which it is difficult to provide very precise guidelines. Despite these difficulties, our point of view is that the tuning of the system is usually even more problematic with a deterministic parameterization of unresolved processes, since no deterministic simulation could exactly fit the real behaviour of the system. By explicitly simulating uncertainties, we can describe the actual random behaviour of the system (see Fig. 1); ensemble simulations can be objectively compared to observations (using probabilistic methods, see Brier, 1950;Toth et al., 2003;Candille and Talagrand, 2005); and the model (including the stochastic parameters) can be rejected as soon as the ensemble is not reliable. Unknown parameters could also be tuned by solving inverse problems, until ensemble reliability is achieved.

Stochastic perturbed parameterized tendency
A first way of explicitly simulating uncertainties in meteorological weather forecast was introduced about 15 years ago in the ECMWF ensemble forecasting system (Buizza et al., 1999). Their basic idea was to separate the model tendency (M) into non-parameterized (N P) and parameterized (P) tendencies (M = N P +P). The non-parameterized tendency (N P) contains all processes that are fully resolved by the model, and can be assumed free of uncertainties. The parameterized tendency (P) contains the parameterization of the effect of unresolved processes (system B in Fig. 1), which is essentially uncertain. The stochastic parameterization is then introduced by multiplying the parameterized tendency (P) by a random noise, explicitly simulating the uncertainties in P. The basic motivation was to produce ensemble forecasts with enhanced dispersion to improve their reliability (i.e. their consistency with available observations). This SPPT (for stochastic perturbed parameterized tendency) parameterization is still used today in the ECMWF ensemble forecasting system (Palmer et al., 2009).
This kind of stochastic parameterization is also meaningful in ocean models, and it can be directly applied in the model using the generic implementation described in Sect. 2.1. This can be done by using one or several of the ξ (i) given by Eq. (1) as multiplicative noise for the various terms of the parameterized tendency: where t is time; x, the model state vector; u, the model forcing; and p, the vector of model parameters. In this case, the mean of the ξ (i) must be set to 1, assuming that the model parameterized tendencies are unbiased, and the other statistical parameters (SD, time and space correlation structure, marginal distribution) are free to be adjusted to any reasonable assumption about the uncertainties. In ocean models, this stochastic parameterization can be applied to any parameterization of unresolved processes (see Fig. 1), as for instance the diffusion operators, simulating the effect of unresolved scales, the air-sea turbulent fluxes, the parameterization of the various functions of the ecosystem dynamics, usually describing the unresolved biologic diversity, etc. An example of this SPPT parameterization is given in the ecosystem application (see Sect. 3.2).

Stochastic parameterization of unresolved fluctuations
Another way of explicitly simulating uncertainties in ocean models is to directly represent the effect of unresolved scales in the model equations using stochastic processes. Unresolved scales can indeed produce a large-scale effect as a result of the nonlinearity of the model equations. Important nonlinear terms in ocean models are for instance the advection term, the seawater equation of state, the functions describing the behaviour of the ecosystem, etc. Concerning the advection term, the effect of unresolved scales is usually parameterized as an additional diffusion, while for the other terms it is most often ignored. However, in many cases, a direct way of simulating this effect would be to generate an ensemble of random fluctuations δx (i) with the same statistical properties as the unresolved scales, and to average the model operator over the ensemble: This corresponds to an averaging of the model equations over a set of fluctuations δx (i) representing the unresolved scales. The zero mean fluctuations δx (i) can produce an average effect (corresponding to an interaction between A and B in Fig. 1) as soon as the model M is nonlinear. In this parameterization, the number of independent fluctuations (m) and the statistics for each of them should be chosen to simulate the properties of the unresolved scales as accurately as possible.
Geosci. Model Dev., 8, 1285-1297, 2015 www.geosci-model-dev.net/8/1285/2015/ Obviously, the main difficulty with this method is to generate fluctuations δx (i) with the right statistics to faithfully correspond to the statistics of unresolved processes. As a first very simple approach, this can be done using one or several of the ξ (i) given by Eq. (1), either by assuming that the statistics of δx (i) can be directly approximated by the simple statistical structure of autoregressive processes ξ (i) , or by assuming that δx (i) can be computed as a joint function of the model state x and the autoregressive processes ξ (i) . For example, if the fluctuations can be assumed proportional to the large-scale gradient ∇x of the state vector, the fluctuations δx (i) could be computed as the scalar product of ∇x with random walks ξ (i) : This particular case corresponds to the stochastic parameterization proposed in Brankart (2013) to simulate the effect of unresolved scales in the computation of the horizontal density gradient because of the nonlinearity of the seawater equation of state. Examples of this parameterization are given in the circulation model application (Sect. 3.1) and in the ecosystem application (Sect. 3.2).
Before concluding this section, it is important to remember that the above discussion only provides one possible framework for simulating the effect of unresolved fluctuations, and that other approaches can be imagined. For instance, a specific stochastic parameterization is already routinely applied at ECMWF to simulate the backscatter of kinetic energy from unresolved scales to the smaller scales that are resolved by the model (Shutts, 2005). This scheme has been developed for atmospheric applications but might also be applicable to ocean models. On the other hand, the external forcing u (e.g. atmospheric data, river runoff, open-sea boundary conditions) can also be a major source of uncertainty in the model, which can be explicitly simulated using a formulation similar to Eq. (5): where the fluctuations δu (i) must be tuned to correctly reproduce the effect of uncertainties in the forcing. Introducing appropriate perturbations of the atmospheric data can for instance be useful to include them in the control vector of ocean data assimilation systems (Skandrani et al., 2009;Meinvielle et al., 2013).

Stochastic parameterization of unresolved diversity
Another general source of uncertainty in ocean models is the simplification of the system by aggregation of several system components using one single state variable and one single set of parameters. For instance, marine ecosystems always contain a wide diversity of species, which cannot be described separately by the model, and which must be aggregated in a limited number of state variables. In a similar way, sea ice can display a wide variety of dynamical behaviours, which cannot always be resolved by ocean models. As unresolved scales, unresolved diversity generates uncertainties in the evolution of the system, which can be explicitly simulated using a similar approach: where δp (i) are random parameter fluctuations representing the various possible dynamical behaviours that are simultaneously present in the system. The application of this method requires a statistical description of the uncertainties in the parameters; and again, as a first approach, this can be parameterized using one or several of the ξ (i) given by Eq. (1). As a particular case, this method includes the stochastic parameterization proposed in Juricke et al. (2013) to explicitly simulate uncertainties in ice strength in a finite element ocean model. It was thus very easy to apply the same scheme in the ice component of NEMO, as an example of this parameterization (see Sect. 3.3).

Impact on model simulations
The purpose of this section is now to illustrate the impact of the stochastic parameterizations presented in Sect. 2 in various components of NEMO: in the ocean circulation component in Sect. 3.1, in the ocean ecosystem in Sect. 3.2, and in the sea ice dynamics in Sect. 3.3. The focus of the discussion will be on the probabilistic behaviour of the system (A) as a result of the uncertainties (the interaction with B in Fig. 1). All applications have been performed using the same generic code implementing the stochastic formulation of NEMO described in Sect. 2.

Stochastic circulation model
As a result of the nonlinearity of the seawater equation of state, unresolved potential temperature (T ) and salinity (S) fluctuations (in system B) have a direct impact on the largescale density gradient (in system A), and thus on the horizontal pressure gradient through the thermal wind equation. As shown in Brankart (2013), this effect can be simulated using the scheme described in Sect. 2.3, by applying Eq. (5) to the equation of state: Table 1. Parameters of autoregressive processes for all applications described in this paper. The number of processes is the number of autoregressive processes used in each stochastic parameterization (sometimes multiplied by 3 to produce one process for each component of the random walks). The mean, SD and correlation timescale are the parameters µ (i) , σ (i) and τ (i) used in Eqs. (2) where δT (i) and δS (i) explicitly simulate the unresolved fluctuations of potential temperature and salinity. These fluctuations are generated using random walks following Eq. (6), with parameters for the ξ (i) given in Table 1 (i.e. the same parameterization as in Brankart, 2013). This stochastic parameterization simulates the exchange of potential energy between resolved and unresolved scales, which results from the nonlinearity of the equation of state (see Brankart, 2013, for more details). As for the Reynolds stresses, this should be strongly constrained by physical principles, but we will stick here to the parameters proposed in Brankart (2013), which were derived from a comparison with higher-resolution reanalysis data. It is interesting to note (as a complement to what is explained in Brankart, 2013) that there is a close similarity between this stochastic correction of the large-scale density and the semi-prognostic method proposed in Greatbatch et al. (2004) and Greatbatch and Zhai (2006). In both cases, indeed the only correction applied to the model occurs in the thermal wind equation through a direct correction of density, while the conservation equation driving the evolution of potential temperature, salinity and horizontal velocity are all kept unchanged. We can thus be certain that the stochastic parameterization displays the same nice conservation properties as the semi-prognostic method; in particular, there is no direct modification of the T and S properties of the water masses, no enhanced diapycnal mixing and thus no compromise with the fact that the ocean interior should primarily flow close to the neutral tangent plane. The modification of the thermohaline structure of the ocean is only produced indirectly through a modification of the main currents.
The first impact of the stochastic T and S fluctuations is indeed on the mean circulation simulated by the model. This mean effect in a low-resolution global configuration of NEMO (the ORCA2 configuration, see Madec and Imbard, 1996, for more detail) has been described in detail in Brankart (2013). In summary, the density correction is important (and quite systematically negative because of the convexity of the equation of state) along the main fronts separating the subtropical and subpolar gyres. The mean pathway of the mean current is thus modified, significantly reducing the biases of the deterministic model. In particular, the Gulf Stream pathway no longer overshoots and the structure of the northwestern corner becomes more realistic. The impact on the mean circulation is similar to what can be obtained with the semi-prognostic method (Greatbatch et al., 2004), in which the density correction is diagnosed from observations, whereas the stochastic model behaves as an autonomous dynamical system.
The second effect of the stochastic T and S fluctuations is to generate random variability in the system. Because of the nonlinearity of the equation of state, the small scales constantly modify the structure of the large-scale density, and thus the pathway of the large-scale circulation. There is a constant flux of information from system B (small scales) to system A (large scales), which is represented in the stochastic model by the random processes ξ (i) , and which is totally absent in the deterministic model. This effect is illustrated in Fig. 2, which displays the pattern of sea surface height (SSH) in several key regions of the Atlantic: the northwestern corner (top panels), the Brazil-Malvinas Confluence Zone (middle panels), and the Agulhas Current retroflection (bottom panels). In the non-stochastic simulation, in the absence of interannual variability of the atmospheric forcing (as in Brankart, 2013), the interannual variability is extremely Geosci. Model Dev., 8, 1285-1297, 2015 www.geosci-model-dev.net/8/1285/2015/ this is why only 1 typical year is shown, since all years would appear identical. In the stochastic simulation however, not only the mean SSH pattern is modified (as shown in Brankart, 2013), the interannual variability is also strongly enhanced, and thus becomes more compatible with the intrinsic largescale SSH variability that is obtained from higher-resolution models or from satellite altimetric measurements (as diagnosed in Penduff et al., 2011). This intrinsic variability (produced in the absence of any interannual variability in the atmospheric forcing) is a good proxy to the dispersion that would be observed in a truly probabilistic ensemble forecast. In a high-resolution model, this dispersion in the largescale behaviour can only result from the interaction with the mesoscale (as explained in Penduff et al., 2011). In the lowresolution ORCA2 configuration, this unpredictable and intrinsically variable behaviour of the large scales is here (at least partially) restored by a stochastic parameterization of the effect of the mesoscale (which is in system B) on the large-scale density. It must be mentioned however that such a small size sample is not sufficient to provide accurate quan-titative information on the magnitude of this effect. To give more precise quantitative results, further tuning and validation of the stochastic parameterization are required.
To further explore the effect of these uncertainties, we are currently applying the same stochastic parameterization to a 1/4 • resolution model configuration of the North Atlantic (NATL025). The results (obtained with the parameters listed in Table 1) indicate that the stochastic parameterization tends to produce a mean effect on the Gulf Stream pathway, and to decorrelate the mesoscale patterns produced in different members of the ensemble. The first questions that we would like to address with this kind of simulation are whether it is possible to better tune the stochastic parameterization using reference data, whether the ensemble dispersion can explain a substantial part of the misfit with altimetric observations, and thus whether this kind of ensemble can be used to assimilate SSH measurements in NATL025. And then, as a longerterm perspective, maybe the stochastic processes ξ (i) can be used as a control vector for data assimilation, which would therefore display the same nice conservation properties as the semi-prognostic method (Greatbatch et al., 2004).

Stochastic ecosystem model
There are many sources of uncertainty in marine ecosystem models. To simplify the discussion, only two classes of uncertainty will be considered here: uncertainties resulting from unresolved biologic diversity and uncertainties resulting from unresolved scales in biogeochemical tracers (see Fig. 1). On the one hand, the most common simplification in biogeochemical models (Le Quéré et al., 2005) is to aggregate the biogeochemical components of the ocean in a limited number of categories (defining system A in Fig. 1). This reduces the number of state variables and parameters, and introduces uncertainties in the model equations since the various components included in one single category (unresolved diversity, in system B) do not usually display the same dynamical behaviour. To simulate this first class of uncertainty, we will use the SPPT scheme described in Sect. 2.2 and multiply the "source minus sink" terms (SMS k ) of the ecosystem model by a multiplicative noise: where C l are the biogeochemical tracer concentrations, and ξ (k) are autoregressive processes obtained from Eq. (1), with parameters given in Table 1. To simulate unresolved diversity, the scheme described in Sect. 2.4 would probably have been more natural, but in view of the large number of parameters in the ecosystem model, the SPPT scheme is much easier to implement as a first approach. On the other hand, to simulate uncertainties resulting from unresolved scales, we will use the scheme described in Sect. 2.3, by applying Eq. (5) to the SMS terms: where δC (i) l explicitly simulate the unresolved fluctuations of biogeochemical tracer concentrations. These fluctuations are generated using random walks following Eq. (6), with parameters for the ξ (i) given in Table 1. (Since little is known about uncertainties in the ecosystem model, we just used here reasonable values for the parameters.) As a first approach, the impact of these two stochastic parameterizations has been studied in a low-resolution global ocean model, based on the ORCA2 configuration coupled to the LOBSTER ecosystem model (using exactly the same model settings as in the previous section). The ecosystem model (see Lévy et al., 2005 for more details) is a simple model including only six compartments (C k , k = 1, . . ., 6): phytoplankton, zooplankton, nitrate, ammonium, dissolved organic matter, and detritus. The behaviour of this model is here illustrated in Fig. 3 by the surface phytoplankton in the North Atlantic for 15 June (in the second year of simulation).
As compared to the deterministic simulation (top left panel), the stochastic simulation with the SPPT scheme (top right panel) does not modify very strongly the general behaviour of the system (despite the 50 % SD multiplicative noise), but substantially increases the patchiness of the surface phytoplankton concentration. This suggests the conjecture that uncertainties (in particular, unresolved diversity) may partly explain the patchiness of satellite ocean colour images. Conversely, the stochastic simulation of unresolved scales (bottom left panel) does not increase patchiness, but can significantly affect the local behaviour of the system, sometimes increasing or decreasing the production (whether the second derivative of the SMS term is positive or negative). At first sight, these two sources of uncertainty are thus insufficient to explain the considerable misfit between model simulation and ocean colour data.
As an additional experiment, the two stochastic parameterizations have then been used together (bottom right panel), by simply generating a sufficient number of autoregressive processes (corresponding to the two columns together in Table 1) to feed the two schemes. This result shows that there is a strong interaction between the two schemes, leading to a deep modification of the general behaviour of the system, and to enhanced patchiness as compared to the SPPT scheme alone. In our view, this directly leads to the idea that uncertainties may be an important ingredient to understand the dynamical behaviours of marine ecosystems, and to make the model distribution consistent with ocean colour observations (in magnitude and pattern). It must be noted however that these experiments only represent a first attempt to explicitly simulate uncertainties in the ecosystem component of NEMO, and that further studies are needed before any meaningful quantitative result can be obtained.

Stochastic sea ice model
One of the main difficulties of sea ice models is to correctly simulate the wide diversity of ice dynamical behaviours. Among ice characteristics, the most sensitive parameter is certainly the ice strength P * . In simple ocean models (as in LIM2 in NEMO), P * is assumed constant, whereas, in more complex models (as in LIM3 in NEMO), the variations of P * can be explicitly resolved as a function of the various types of ice simultaneously present at every model grid point. The impact of uncertainties in P * has already been studied in the work of Juricke et al. (2013) using a finite element ocean model (FESOM), coupled to a simple sea ice model similar to LIM2. The purpose of this section is to reproduce their parameterization in NEMO/LIM2 using the generic technical approach described in Sect. 2. This can be done very easily, almost without any additional implementation effort, using the scheme described by Eq. (8) with m = 1 and where ξ is one of the autoregressive processes given by Eq.
(1), with parameters given in Table 1. The parameters are chosen to be close to the stochastic parameterization in Juricke et al. (2013). Specificities are the use of order 2 instead of order 1 autoregressive processes, and the use of a gamma marginal distribution instead of another kind of positive distribution in Juricke et al. (2013). This stochastic parameterization has been applied to a lowresolution global ocean configuration of NEMO, again without interannual variability in the atmospheric forcing (using the same model settings as in Brankart, 2013). The behaviour of the model is here illustrated in Fig. 4 by the ice thickness in the Arctic at the end of March (when the ice extension is close to its maximum). As compared to the deterministic simulation (top left panel), the first impact of the stochastic parameterization is to systematically increase ice thicknesses, especially in the regions of old ice (north of Greenland and western Canada), and to slightly decrease the ice extension. This mean effect results from the nonlinearity of the model response to P * : during the periods of small P * , the ice thickness has the opportunity to increase, and this increase is not counterbalanced by a symmetric decrease of thickness during the periods of large P * . This behaviour is very similar to what is described in Juricke et al. (2013), and cannot be reproduced by a simple uniform modification of P * . On the other hand, the stochastic fluctuations of P * also generate random variability in the system. As for SSH in Sect. 3.1, the interannual variability of the ice thickness pattern is extremely weak in ORCA2 without interannual variability of the atmospheric forcing (which is why only 1 typical year is shown in Fig. 4). In the stochastic simulation however, not only is the mean ice thickness pattern modified (as for SSH in Fig. 2), but the interannual variability (which is again a good proxy to ensemble dispersion as explained in Sect. 3.1) is also strongly enhanced. What is expected from these results is thus that the explicit simulation of uncertainties can provide us with an adequate basis for probabilistic comparison with sea ice observations and help us in producing reliable ensemble forecasts for sea ice data assimilation problems. Consequently, it might also be that this stochastic approach represents a worthwhile alternative to explicit resolutions of sea ice diversity (as in LIM3).
In this paper, a simple and generic implementation approach has been presented, with the purpose of transforming a deterministic ocean model (like NEMO) into a probabilistic model. With this method, it is possible to easily implement various kinds of stochastic parameterizations mimicking the non-deterministic effect of unresolved processes, unresolved scales, unresolved diversity, etc. It has been shown indeed that ocean systems can often display a random behaviour, which needs to be explicitly represented in ocean models. Ensemble simulations are then required to sample all possible behaviours of the system. Getting a reliable overview of all dynamical possibilities is necessary to objectively compare models to observations, and to correctly apply the model constraint in ocean data assimilation problems.
Technically, what is proposed here is a very simple algorithmic solution that is easy to adapt to many kinds of models, and that is generic enough to deal with many different sources of uncertainty. This is obviously not intended to be the final theoretical or technical solution for simulating uncertainties. The algorithms and framework proposed in this study only provide a first-guess solution, which is simple enough to make a first quick evaluation of the effect of a given source of uncertainty, and flexible enough to easily evolve as a better understanding of the problem is progressively obtained. This technique has been applied to several applications, showing that randomness is ubiquitous in ocean systems: in the large-scale circulation (e.g. because of the effect of unresolved scales through the nonlinear equation of state), in the ecosystem model (e.g. because of the effect of unresolved scales and unresolved biogeochemical diversity), and in the sea ice dynamics (e.g. because of the unresolved diversity of sea ice characteristics). In each of these applications, uncertainty can be viewed as an essential dynamical characteristic of the system, which can modify our understanding of the ocean behaviour. As for any complex system, constructing ocean models using optimal (but imperfect) components can often be worse (less robust) than using unreliable components dealing explicitly with their respective inaccuracy. The ocean is like dice rolling on the table of a casino: we are unable to grasp all subtleties of their movements, and we can only sample from all possible outcomes of the game using probabilistic models.