An ensemble probabilisitic approach to reconstruct the biogeochemical state of the North Atlantic Ocean using ocean colour images

In this paper, we investigate the potential of using a probabilistic modelling approach in the prospect of ocean colour data assimilation. The main objective of the study is to assess the benefits of using error covariances based on an explicit simulation of model uncertainties. The relevance of this approach is evaluated by considering 3D observational updates of the ensemble (one update at one time step) performed every 5 10 days (over one year) using the statistics of a North Atlantic coupled NEMO/PISCES stochastic ensemble simulation involving 60 members, as previously described in Garnier et al. (2016). In this experiment, SeaWIFS ocean colour data are used to update the ensemble with a low rank ensemble Kalman Filter analysis scheme. The non-Gaussian behaviour of the model variables is taken into account using anamorphic transformations. Comparisons between the updated ensemble and the MERIS satellite 15 observations shows that the integration of high resolution SeaWIFS data significantly improves the representation and the ensemble statistics of chlorophyll concentrations. We also show that these improvements consistently cascade in the water column chlorophyll distributions and on non-observed variables closely linked with the primary production. In addition, we present first results illustrating the potential of our approach for biogeochemical forecasts. 20 The objective is to examine the model response to data assimilation in the perspective of future operational applications. For this purpose, we perform a 60 member simulation initiated from updated biogeochemical states. This forecast simulation shows that ocean colour data assimilation would be skillful considering integration cycles of the order of a day. Finally, the intend of this article is to point out the feasibility of operational biogeochemical data assimilation in the near future. 25

days (over one year) using the statistics of a North Atlantic coupled NEMO/PISCES stochastic ensemble simulation involving 60 members, as previously described in Garnier et al. (2016).
In this experiment, SeaWIFS ocean colour data are used to update the ensemble with a low rank ensemble Kalman Filter analysis scheme. The non-Gaussian behaviour of the model variables is taken into account using anamorphic transformations. Comparisons between the updated ensemble and the MERIS satellite 15 observations shows that the integration of high resolution SeaWIFS data significantly improves the representation and the ensemble statistics of chlorophyll concentrations. We also show that these improvements consistently cascade in the water column chlorophyll distributions and on non-observed variables closely linked with the primary production.
In addition, we present first results illustrating the potential of our approach for biogeochemical forecasts. 20 The objective is to examine the model response to data assimilation in the perspective of future operational applications. For this purpose, we perform a 60 member simulation initiated from updated biogeochemical states. This forecast simulation shows that ocean colour data assimilation would be skillful considering integration cycles of the order of a day. Finally, the intend of this article is to point out the feasibility of operational biogeochemical data assimilation in the near future.

Introduction
Since the 1980's, observations from a sequence of ocean colour satellite missions (e.g. Eppley et al., 1985;Platt and Sathyendranath, 1988;Antoine et al., 1996;Behrenfeld and Falkowski, 1997) are available 30 to characterize the timing, spatial distribution and amplitude of ocean primary production and seasonal blooms. Modelling approaches have been developed in parallel, providing additional information to complement remote sensing data on the vertical (into the whole euphotic zone), the horizontal (between satellite swaths or to compensate cloud coverage) and temporal (between revisit observation times) dimensions (e.g., Doney et al., 1996;Aumont and Bopp, 2006;Gruber et al., 2006;Ilyina et al., 2013;Aumont et al., 2015). 35 Considering today's fundamental needs for estimating stocks of carbon and marine resources in a climate change context, a synthesis between these different sources of information including their respective accuracy is required to reduce uncertainties on the representation of the marine biogeochemistry (and in particular the primary production).
The generation of synthetic products describing the biogeochemical ocean state is one of the challenging 40 goals of the Copernicus Marine Environment Monitoring Service (CMEMS). One major issue is to develop a capacity to routinely assimilate ocean colour data into coupled physical-biogeochemical models, and thereof deliver real-time and reanalysis products of the global ocean biogeochemistry. Following the seminal modelling approach proposed by Garnier et al. (2016) to develop a stochastic version of the PISCES model, we examine in this paper the usefulness of ensemble biogeochemical simulations to make one step towards this 45 synthesis goal. So far, the data assimilated in biogeochemical models have varied from phytoplankton light absorption (Shulman et al., 2013), satellite diffusive attenuation coefficient (Ciavatta et al., 2014), ocean-colour plankton functional types (Skákala et al., 2018;Ciavatta et al., 2018) or colored dissolved organic carbon (Gregg and Rousseaux, 2016). However, most of the studies focused on the assimilation of chlorophyll a concentrations 50 derived from ocean colour data because other essential variables of the marine biogeochemistry are nearly unobserved. The assimilation techniques vary from optimal interpolation methods (Gregg, 2008), variational methods (Losa et al., 2003) or Kalman Filter methods (Ciavatta et al., 2011b(Ciavatta et al., , 2016Fontana et al., 2009;Simon and Bertino, 2009b;Simon et al., 2015).
Today, the performance of ocean colour data assimilation remains limited by the ability of biogeochem-55 ical models (BGCMs) to realistically represent the observed variability of biogeochemical properties. This limitation results, in part, from the many sources of uncertainties associated to the approximate description of the physical background, the reduced complexity of biogeochemical model, the unresolved scales or the empirical parameterizations of biogeochemical processes. In fact, the description of the model uncertainties (which are often referred to as "model errors") is a critical challenge of the data assimilation problem (Lahoz 60 et al., 2010). In this context, ensemble methods (e.g., Evensen, 1994Evensen, , 2003 provide a statistical means to describe uncertainties associated with a complex model system by simulating the evolution of the probability density function (pdf) of the system. The introduction of random processes into the model equations can be used to simulate the spread of the ensemble members introduced in a data assimilation scheme. In the numerical weather prediction context, an increasing number of operational systems rely on probabilistic ensemble simulations (e.g., Palmer, 2012;Berner et al., 2011;Houtekamer and Zhang, 2016), while ensemble simulations are still quite unusual in oceanography. For oceanic biogeochemical applications, Woods and Onken (1982) and Wolf and Woods (1988) addressed the question of biogeochemical modelling using an ensemble of Lagrangian particles to represent phytoplankton populations. More recently, several studies (e.g., Dowd, 2011;Weir et al., 70 2013) have used probabilistic configurations of coupled biogeochemical models (CPBMs) based on the introduction of stochastic parameterizations that simulate the effects of unresolved or poorly-resolve processes (see Leutbecher et al., 2017, for a review). With this approach, the effect of hypothesized uncertainties are directly implemented into the model equations by introducing random numbers to simulate uncertainties. In this context, the objective of the present study is to assess how the ensemble simulations, as produced by 75 Garnier et al. (2016), can be used in a data assimilation perspective. Compared to other studies, the originality of this approach is the use of ensemble statistics originating from the parameterization of modelling uncertainties to compute ocean colour assimilation updates. Finally, the overarching goal is to test whether the integration of ocean colour data into a probabilistic coupled simulation can improve the representation of primary production, as a first step to build a complete ensemble ocean data assimilation system.

80
To achieve this goal, we first use the probabilistic North Atlantic 1/4 o configuration of the NEMO-PISCES coupled model (hereinafter NATL025-PISCES) fully described in Garnier et al. (2016) to generate a 60-members prior ensemble simulation. Then, using high resolution SeaWiFS ocean colour data, we compute a sequence of 5-day ensemble updates to investigate if a low rank ensemble Kalman filter analysis scheme is able to reduce the model uncertainty. In particular we will verify if the analysis improves the 85 statistical consistency of the ensemble. I order to investigate the predictability of the system, some of the biogeochemical updates are thereafter used as initial conditions to produce ensemble forecasts. The objective is to assess whether the information brought by the observations is preserved with time during the forecast, which is a necessary condition prior to data assimilation.
The paper is structured as follows: section 2 introduces the model configuration and the prior ensemble 90 simulation. The section 3 presents the ensemble analysis method and a sequence of three-dimensional ensemble updates. The section 4 describes the potential of ocean colour observations in the perspective of biogeochemical forecasting. Finally, a summary including conclusions and future perspectives are proposed in Section 5.  Madec et al. (2012)) platform in its version 3.4 to realistically represent the ocean dynamics and its interactions with the biogeochemical components. The main part of the NEMO system 100 is the 3D hydrostatic OPA (Ocean PArallelised, Madec et al. (1998)) module which computes the oceanic circulation. For the purpose of this study, OPA is coupled with the sea ice LIM-v2 model (Louvain-La-Neuve sea Ice Model, Fichefet and Morales Maqueda, 1997) and the biogeochemical model PISCES-v2 (Aumont et al., 2015).
The coupled NEMO/PISCES model configuration is based on the NATL025 configuration, a sub-element 105 of 1/4 • global ORCA025 configuration developed within the framework of the DRAKKAR project (Barnier et al., 2006) to ensure a realistic representation of the ocean dynamic. Since the physical representation of the ocean is crucial for simulating primary production, an eddy-permitting 1/4 • horizontal resolution is used. This configuration includes 46 vertical levels and covers the North Atlantic basin from 20 • S to 80 • N and from 98 • W to 23 • E. In the context of biogeochemistry, the suitability of this eddy-permitting 110 framework has been demonstrated in previous studies using the LOBSTER model (e.g. Doron et al., 2011Doron et al., , 2013Fontana et al., 2013). The physical component is coupled with the biogeochemical PISCES-v2 model Aumont et al. (2015). The PISCES architecture consists on 24 biogeochemical variables from 4 major compartments including 2 classes of phytoplankton (diatoms and nanophytoplankton), 2 classes of zooplankton (mesozooplankton and 115 microzooplankton), 5 limiting nutrients (iron, silicate, phosphate, ammonium, and nitrate) and a detritus compartment composed of dissolved and particulate carbon. This model has been widely used in previous environmental studies such as Rodgers et al. (2008); Brasseur et al. (2009);Steinacher et al. (2010); Tagliabue et al. (2014). PISCES is computed to be able to simulate the large scale biogeochemical processes of the first trophic level of marine ecosystems  and is considered as a relevant solution to simulate 120 primary production at the North Atlantic basin scale.
To guarantee a realistic representation, the coupling between the physical and the biogeochemical models is performed on-line. The coupling frequency is 40 minutes. Note that the biogeochemistry has no feedbacks on the physics. The tracer advection is computed by the physical model using a MUSCL scheme (Monotone Upstream Scheme for Conservative laws, (LeVeque, 2002). Tracer vertical diffusivity fluxes are calculated 125 with the Turbulent Kinetic Energy (TKE, Blanke and Delecluse (1993)) closure scheme. The horizontal diffusivity is expressed with a laplacian operator.
All the physical parameters used in the coupling are issued from a NATL025 simulation previously validated in Candille et al. (2015). The biogeochemical part is initialized in January 2002 from the MERCATOR-OCEAN 1/4 • coupled NEMO/PISCES model BIOMER simulation. Between January 2002 and December 130 2004, an additional 3 years spin-up is performed to ensure a consistent biogeochemical initial state, prior to perform ensemble simulations used in this study.

The prior model ensemble chlorophyll estimation
In order to characterize model uncertainties, we introduce random processes (w tk ) into the model for-135 mulations. The evolution to the biogeochemical state x tk by the model M is then described by equation 1. The biogeochemical state temporal dynamic depends on the random perturbation sequence of w tk generated here from first order auto-regressive processes (AR1) including a spatial correlation of 5 grid points and a decorrelation time of 30 days. The w tk random perturbations are then introduced into the model formulation through the stochastic parameterizations proposed in Garnier et al. (2016) to simulate the uncertainties on 140 biogeochemical parameters and the uncertainties induced by unresolved scales in the presence of non-linear processes. Similar approaches (e.g., Ciavatta et al., 2016) use random perturbations over model forcing parameters in order to define the model error. With these methodologies, perturbations are only applied to the initial state. By contrast, using the present formulation, the co-variances statistics are time dependent and defined from the model dynamics which improves the prior model error definition.
Using this probabilistic approach, we perform a 1-year 60-member ensemble simulation during 2005 which is briefly shown here. For this purpose, figure 1 shows the ensemble quartiles of surface chlorophyll along 20 • W longitudinal sections together with MERIS observations. Quartiles are calculated at each grid point from the prior ensemble probability distributions. It therefore does not display a model solution but a composite between the 60-member of the ensemble.

150
Because of strong biogeochemical dynamics, the dispersion is higher during the spring bloom period. We also observe that the level of dispersion is very space and time dependent. It tends to be maximum during the spring bloom, in the northern part of the domain and in the equatorial current, where chlrophyll concentration are generally high. Note that around the 20 • latitude region, the very high dispersion level is due to the presence of the Mauritania's upwelling. In some particular geographical areas, the biogeochemical 155 activity can be very weak. Chlorophyll concentrations are close to zero and the ensemble is nearly nondispersive. In spite of this large spread, the median is nearly always very close to the deterministic standard simulation (detIni ) described in Garnier et al. (2016), proving that PISCES main large scale dynamics is well preserved. The main difference with respect to a deterministic simulation is that most of the observations (about 70% over the whole domain) are well included within the ensemble envelope, which means that the 160 stochastic parameterizations generate a level of uncertainty consistent with the variability of the observations. The ensemble quantiles also provides a rough vision of the shape of the probability distribution. The range of values of the 50% higher chlorophyll concentrations (above the median) is always larger than the 50% range of smaller concentrations (below the median) which underlines the non-gaussianity of the probability distribution. This point will be taken into account in the assimilation process using an anamorphosis 165 transformation (cf section 3.1) It is necessary to specify that envelopes presented in figure 1 depict the statistical information implemented as the model error in the analysis scheme presented in the next section. Because it is constructed from model solutions, it directly includes a spatial and temporal representation of its variability. Since these results indicate that the ensemble display a sufficient level of uncertainty to capture the information 170 content of the observations, the model error definition should also be consistent with the variability of the observations. In addition, such an evolutive characterization of model uncertainties present a robust technique which avoid to increase model co-variances when the dispersion tends to be very small (in particularly outside seasonal bloom periods). Red dots indicate SeaWIFS observations. the deterministic detIni standard PISCES simulation, performed without stochastic parameterizations (Garnier et al., 2016), is represented in cyan.

175
In this part, the prior probability distributions are used to update the model biogeochemical state. We present a sequence of analysis meaning that updated states are not propagated by the model. The objective is to assess the benefits of using prior covariances based on a explicit simulation of model uncertainties. The analysis presented in this article is based on the sequential Ensemble Transform Kalman Filter (ETKF, Bishop et al., 2001). The analysis scheme involves a transformation matrix to calculate the covariances in the space of anomalies δx f i (size s × N ) which are calculated from the ensemble mean state x f (equation 2) of the model forecasts {x f i : i = 1, 2..., N } here N = 60) following the equation 3.
These anomalies characterize the matrix S f of size i which defines the prior covariance matrix P f following equation 4. Note that with this formulation, error distributions are subject to the gaussian hypothesis.
The Kalman gain matrix (K) is calculated in the reduced space of anomalies from the product of the anomaly matrix S a . The change of space is made by the transformation matrix T expressed by equation 5, where H is the observation operator to pass through the observation space. In the observation space, the 190 anomalies are therefore expressed by Y f = HS f . R is the diagonal observation error covariance matrix.
We now define the transformation matrix T which calculates the updated anomalies and the ensemble mean.
The step of analysis as it will be done here is then only an updated of the ensemble members: In this scheme, the gaussianity of the prior probability distributions remains a strong hypothesis. To consider the non-gaussian shape of chlorophyll distributions, a solution is to apply a log-normal transfor-200 mation. However, using such a space in time invariant transformation is generally insufficient (Ciavatta et al., 2011a). In this study, the non-gaussianity of the probability distribution is taken into account using anamorphosis transformation ) depending on the shape of the probability distribution. The principle is to characterize a bijective transformation to redesign the percentiles of the distribution in such a way to produce a gaussian probability distribution. The ETKF analysis scheme is thereafter calculated 205 in this anamorphosis space. In the context of biogeochemical data assimilation, this method has already been successfully employed in several studies such as Simon and Bertino (2009a); Doron et al. (2011Doron et al. ( , 2013; Fontana et al. (2013).
Although a stochastic approach is a relevant solution to characterize model uncertainties, the sampling of the covariance matrix from 60 members can be a source of unrealistic large distance spatial correlations.

210
In an area like the North Atlantic where horizontal dimensions are much larger than the correlation lengths, this phenomenon must not be neglected. For this purpose, the spatial extent of corrections is reduced by introducing a domain localization matrix (Testut et al., 2003). Model covariances are multiplied by a correlation function γ(r) = exp(−r 2 /l 2 ) which exponentially attenuates spatial correlations as a function of the distance r. The cut-off radius l c defines the distance after which the covariances are imposed to zero.

215
The values l = 4 and l c = 16 grid points, inherited from the work of Candille et al. (2015), are used.

Ocean colour observations
In this study, chlorophyll concentrations assimilated in the coupled model are derived from the SeaW-iFS sensor launched on-board the SeaStar-Orbview2 satellite. More specifically, we use a high resolution 220 (1/24 • ) version of the SeaWiFS data set-up in the ESA (European Spatial Agency) project OC-CCI (Ocean Colour Climate Change Initiative). These data are accessible through the MyOcean-Copernicus web site (http://marine.copernicus.eu). In addition, 1/12 • MERIS-Envisat ocean colour observations are used to perform statistical assessments.
In the analysis scheme described in section 3.1.1, we suppose 1) a gaussian distribution of observa-225 tion errors and 2) that the observation are spatially uncorrelated. Consequently, R is a diagonal matrix characterized by a 30% observation error. In order to take into account the non-gaussianity of chlorophyll distributions, the observations y 0 and their associated standard deviations are also evaluated into the anamorphosis space defined by the ensemble simulation. Note that the anamorphosis transforms marginal distributions into gaussian distributions but it does not guarantee the gaussianity of joint distributions 230 (between variables). Using high resolution observations, a particular attention is given to the validity of the uncorrelated observation hypothesis. Indeed, in spite of its requirement in the assimilation scheme, errors between near observations points can not be entirely uncorrelated. Without considering these correlations, the model information is overwhelmed by the redundancy of the observations. To take into account correlation of 235 errors while maintaining a diagonal matrix, we inflate the diagonal terms of the observation covariances matrix proportionally to an estimated distance of the spatial correlations. A constant correlation distance of 2.3 • , determined experimentally, is set. The equation 10 presents these parameters where σ is the standard deviation and the α constant refers to the square root of the number of data observations contained within a 2.3 • radius circle around the grid point.

The Updated ensemble
Using the method described previously, the prior 60-member model simulation (cf. section 2.2) is corrected every 5 days during the entire year 2005. Hence, the results correspond to a 3D sequence of ocean colour data analysis (73 biogeochemical updated states).  May). Compared to the prior solution, the analysis increases chlorophyll concentrations over the subtropical gyre and decreases it in the Gulf Stream and the Northern regions. Consequently, the strong north to south gradient existing in the prior representation is smoothed, which results into an enlargement of the subtropical oligotrophic conditions. In general, the updates reduce deviations with observations but a more detailed comparison with observations is necesseray. It will be performed in next sections.

255
On average, the analysis do not induce abrupt variations of the large scale chlorophyll surface distributions and the solution always complies with the large scale biogeochemical dynamics defined by the prior model. In the context of data assimilation, this point is essential to maintain information brought by the analysis during propagation steps and does not prevent the model representation to include small scale variability. Thus, these results especially: 1) ensure the relevancy of the analysis method, and 2) confirm a major challenge in data assimilation. The first point to observe is that the analysis has an impact in the whole euphotic zone. In the equatorial zone and in the region located around 40 • N, the analysis 265 reduces the chlorophyll concentration in the first 50 meters. On the contrary, after 40 • N, it is reduced. A remarkable feature is that surface chlorophyll decreasings are associated with an enhancement the DCM (Deep Chlorophyll Maxima). These results highlight an important correlation between surface chlorophyll concentration and vertical chlorophyll distribution. Below the euphotic layer, the analysis does not produce adverse effects.

270
In spite of these impacts, the important point is that the analysis does not affect, on average, the vertical characteristics such as the depth of propagation and the localization of subsurface maxima. which shows that the method doesn't seem to induce unrealistic correlations on vertical fluxes (for example vertical organic matter fluxes). Such as for the surface, the analysis preserve main elements of the PISCES model vertical dynamics, to which the validity was previously verified in Garnier et al. (2016). The use of vertical 275 correlations based on our stochastic parameterizations is therefore able to describe relevant extrapolations of ocean colour data onto the vertical.

Uncertainty of the analysis
Using a probabilistic approach, a measure of the uncertainty must take into account the various shapes tainty. This point is fundamental as it will characterize the prior error model in an assimilation scheme. Rationally, the analysis mainly impacts high chlorophyll concentration patterns for which prior uncertainties are higher. Therefore, the system acts mainly during the high biogeochemical activities of the spring bloom period. Observations remain most of the time close to the median solution and are nearly always captured by updated ensemble spread. Note that updated members are also included within the prior ensemble en-290 velope. It means that chlrophyll is only consistently corrected from a prior distributions in good agreement with observations. In this sense, a good definition of the stochastic parameterization appears is essential.
As regards to the vertical dispersion, figure 5 presents vertical chlorophyll profiles of the 60-member of the prior and the analyzed ensembles at one grid point (STAT C ) located in the Gulf Stream region. As for the 295 surface, the analysis reduces the dispersion while keeping a significant degree of uncertainty throughout the euphotic layer. It is relevant to observe that the 2 different biogeochemical regimes (in spring and automn) are kept during the analysis. The dispersion seems to be higher for the depths of subsurface maxima and during the spring bloom, when model uncertainties are higher. At a global scale, the level of uncertainty is only slightly reduced which is related to the strong level of biogeochemical uncertainties rather than a 300 failure of the system. In this case, this is quite positive since a small dispersion would have been detrimental to the assimilation process. Finally, assuming relevant prior ensemble definition, i.e relevant stochastic parameterization, surface chlorophyll corrections should be consistentely reflected into vertical chlorophyll distributions.

305
Using the proposed probabilistic method in assimilation processes, the introduction of ocean colour data should be able to consistently correct chlorophyll fields. Since stochastic parameterizations characterize an prior ensemble with relevant statistics, co-variances computed from the model dynamics avoid to drive the model toward incoherent biogeochemical states. Obviously, we also need to ensure that this method can sufficiently constrain the biogeochemical system to be considered for data assimilation.

Probabilistic assessment
In this section, impacts of the analysis on the ensemble distribution are investigated from a statistical comparison with MERIS ocean colour observations. Two statistical properties are used: the reliability, which evaluates the level of consistency between the ensemble probability distributions and the statistical distri-315 bution of observations and the resolution which measures the distance between the cumulative probability distribution of the ensemble and the observations expressed in terms of probability distributions. Further information is available in Candille and Talagrand (2005) and Candille et al. (2015). Here, these properties are deduced from the rank of observations (rank histograms, Anderson (1996)) and from the Continuous Rank Probability Score (CRPS, Brown (1974)).

Ranks of MERIS observations
The reliability is evaluated from the repartition of MERIS observations within the ensemble range of chlorophyll concentrations (the rank of observations), displayed in the form of histograms. The flatness indicates the level of reliability: a flat histogram depicts a perfect reliability. As decribed in Saetra et al. 325 (2004), ranks are calculated taking into account a 30% observation error by increasing the ensemble dispersion. The reliability investigated here is a spatial reliability in the sense that the distribution of observations is deduced from their spatial variability at a given date.
In figure 6, the shape of the prior ensemble rank histogram is compared to that of the updated ensemble during the spring bloom. All available observation over the domain are considered. The rank histogram of the prior ensemble display an underdispersive behaviour (about 80% of the observations are included) with a distinctive ∪ shape reflecting the difficulties of the ensemble to encompass the observations. The analysis clearly reduces this underdispersion. The updated ensemble capture approximatively 10% more observation and its histogram is strongly flattened, which indicates that the introduction of SeaWiFS information improves the reliability.

335
In addition, the spatial distribution of the ranks is much more heterogeneous for the updated ensemble. Depending on the geographical zone, this highlights local improvements on the reliability and a better midrange representation of chlorophyll concentrations. The updated ensemble ranks also present higher small scale variability that seems to be in better agreement with the variety of biogeochemical behaviours existing at a basin scale. Nevertheless, as for the prior simulation, many observations (found around the Gulf Stream 340 and the subtropical gyre regions) still remain outside the envelope of the ensemble. Top panels present the ranks from histograms while the bottom panels present spatial maps of the ranks. The first column displays the prior ensemble ranks (stoUpd-xf ) and the second column shows ranks of the updated ensemble (stoUpd-xa). The frequency of occurence is given in percent (over the total number of observations). The blue dotted lines in rank histogram figures refer to the rank repartition of a perfectly reliable ensemble.

The continuous ranked probabilistic score
The continuous ranked probability score (CRPS) is a generalization to continuous events of the Brier score (Brier, 1950). As expressed by equation 11, it corresponds to the distance between the ensemble 345 cumulative probability density function F p (x) of a x variable (here the chlorophyll) and the probability distribution F o (x) of its observation y o (here the ocean colour data).
In order to differentiate the statistical properties it is possible to decompose the CRPS as the sum of 2 terms: CRPS=RELI+RESO. The CRPS decomposition is fully documented in Hersbach (2000). RELI refers to the ensemble reliability and RESO refers to the potential resolution, that is the resolution of the ensemble in the case of a perfect reliability. The closer to 0 is the value the better is the score.  Using all MERIS observations available every 5 days during 2005, we calculate the RELI and the RESO terms of the CRPS decomposition. Results are given in time series in figure 7.
As already observed with rank histograms, it shows that the analysis systematically improves the global reliability of the ensemble. This result is more moderate during the first months of the year because the dispersion coming from the stochastic parameterizations has not yet been entirely generated. In average 360 the reliability is improved of about 70% in summer and of about 30% in winter which confirms the results observed with the rank histogram. The higher biogeochemical activity (and so the ensemble dispersion) during the spring and the summer period probably mainly explain this difference. Indeed, deviations between the model and the observations are higher during this period and the analysis has also higher impacts.
Concerning the resolution, improvements due to the analysis are smaller and remain relatively constant.

365
The analysis improves the resolution for about 15-20% This is the result of a slight narrowing of the probability distribution extent which mainly causes a decreasing of extreme events probability of occurrence. Finally, the analysis never induces negative impacts on the ensemble statistics, which underline the consistency of our method.

correction on non-observed variables 370
So far, we have focused on chlorophyll concentrations. Indeed, the stochastic parameterizations used to generate the ensemble aimed to simulate some source of uncertainties centered on primary production. One of the major issue in data assimilation is,nevertheless, the correction of non-observed variables. This means that we need to verify that the integration of chlorophyll data could cascades onto these non-observed 375 variables To address this issue, the figure 8 presents 60-member surface time series for 3 variables closely linked with primary production, the nitrates, the phosphates and the microzooplankton. The aim is only to ensure that the co-variances inherited from the stochastic parameterizations produce updated biogeochemical states in accordance with the prior information. This figure demonstrates that the stochastic parameterizations also induce dispersion onto the non-380 observed variables (of the prior simulation). In return, the method of analysis is able to affect their statistical distributions and to reduce the ensemble dispersion. As for the chlorophyll, all updated members are kept inside the prior envelope, which means that prior characteristics are conserved after the analysis and highlight the importance of an appropriate definition of prior model statistics. Finally it is relevant to underline that the correlations established from the stochastic parameterizations are able to correct some 385 non-observed variables closely related to primary production. This point is crucial in order to constrain the biogeochemical state in the North Atlantic basin scale using only chlorophyll satellite data.
Nonetheless, a comparison with observation, as it has been done for the chlorophyll, would be necessary to assess the validity of the prior statistics and the following corrections. As long as the stochastic 390 parameterizations would be defined in such a way they can generate consistent dispersion for the concerned non-observed variables, the corrections should be relevant. The question which still held open: would the analysis performs with only one observed variable be sufficient to correctly constrain a biogeochemical state characterized by many variables ? In order to give first elements of response, we perform forecast expriments presented in the next section. The approach presented here consists in producing ensemble forecasts initialized from updated biogeochemical states. The objective is to investigate the model response to the integration of ocean colour information in order to define data assimilation approaches adapted to the context of marine biogeochemistry. Furthermore, another aspect addressed here is the ability of the system to accomplish relevant short 400 and/or long term biogeochemical forecasting.

Impact of the analysis on the prior ensemble distributions
Updates of the ensemble were previously computed every 5 days. In order to make consistent forecast experiments, a finer temporal resolution is essential. Indeed, improvements due to the analysis could be 405 rapidly deteriorated because of the strong level of uncertainty assumed for the PISCES model. We therefore compute daily outputs during a 1-month period of the spring bloom (5 April 2005 to 5 May), known to display strong discrepancies between model and observations. For this experiment, two 60-member ensemble simulations have been carried out. The first simulation is initialized by the multivariate updated biogeochemical state presented in the previous section, while the second is set to be identical to the prior 410 simulation (cf. section 2.2), except that the outputs are daily.  Figure 9 presents 1-month daily time-series of the ensemble quartiles for the prior simulation (dotted lines) and the forecast simulation (full line) in two grid points highlighting typical behaviours. Depending on the geographical area (more precisely on the biogeochemical dynamics), the ocean colour information is generally sustained during several days by the ensemble (e.g at STAT G). A maximum of about 1 month 415 (e.g at STAT D) was detected. This seems to reveal a short term predictability of the system, probably due to the strong biogeochemical model uncertainties. The tendency to rapidely converge towards the prior solution illustrates the stability of the PISCES model. Indeed, PISCES is strongly constrained by the value 15 Ocean Sci. Discuss., https://doi.org/10.5194/os-2018-153 Manuscript under review for journal Ocean Sci. Discussion started: 24 January 2019 c Author(s) 2019. CC BY 4.0 License. of many parameters and the model mostly tends to converge toward the solution for which it was initially calibrated. Note that forecast members tend toward their free corresponding member because the exact 420 same random numbers are used in the free simulation and the forecast.

Comparison with observations
Taking into account our statistical approach, it is essential to assess the consistency between the forecast ensemble and ocean colour observations. Two issues have to be investigated: (1) for how long an update is 425 able to reduce the prior level of uncertainty, and (2) how does it affect the statistical characteristics of the prior probability distributions. To investigate the first issue, figure 10 shows a 1-month time serie of the root mean square error (RMSE) for the forecast and the prior simulation. The RMSE is calculated every day using the mean of the two ensembles and the mean of all the available MERIS observations. Note that observations used for the calculation are spatially varying with time.

430
As expected, the figure shows that the integration of ocean colour data reduces the system level of uncertainty as compared with the prior simulation. It also confirms the short-term predictability of the system. According these results, the decreasing of the RMSE due to the analysis vanishes after 5 days but a significant reduction of the uncertainty is only kept for 2 days. In a context of data assimilation, this indicates that it is necessary to integrate ocean colour data with high frequency to consistently constrain 435 the model. In section 3.4, we have shown that the analysis improves the statistical coherence between the ensemble and satellite observations. The second question to address is whether this feature can be preserved during the forecast. As before, we compare the ensemble simulations with observations using rank histograms of MERIS observations which are presented in figure 11. Note that calculations of the ranks are identical to 440 those presented in figure 6.
After 1 day of forecast (6 April 2005), it is remarkable to see that the shape of the histogram remains nearly unchanged. Statistical improvements due to the analysis are then preserved during this period. Beyond this time, histograms become comparable again with those of the prior ensemble. After 5 days of prediction, the rank histograms of the forecast and the prior ensemble are almost identical. The gain of 445 reliability coming from the analysis has been lost. Because of a high level of uncertainty, the biogeochemical model very rapidly generates some dispersion which make ensemble forecasts no more reliable than the free run after 1 day.
In summary, the analysis improves the predictions of the mean representation of surface chlorophyll for a 450 few days but the ensemble forecast remains reliable only after one day of forecasting. Using our framework, this means that data assimilation could only be efficient performing daily cycles. Next step would be to set-up a complete data assimilation system that integrates ocean colour information to daily update the forecast. In this manner we could evaluate whether statistical improvements seen in the present study are preserved in time or not, even though nowadays satellite products lack of a complete spatial coverage.

Summary and perspectives
The ensemble approach developed in this paper is a first step towards ocean colour data assimilation for ocean biogeochemistry, with a focus on the representation of primary production. Using the probabilistic version of the PISCES model assessed in Garnier et al. (2016), a time-dependent integration of ocean colour images based on stochastic parameterizations has been documented using a 60-member ensemble 460 simulation. The main objective is to demonstrate the potential of this approach to assimilate ocean colour satellite observations into realistic coupled physical/biogeochemical model configurations.
Our results demonstrate that the proposed model error parameterization is appropriate for biogeochemical data assimilation context. In particular, the method is shown to significantly improve the statistical consistency (indicated by reliability and resolution) of the ensemble. Furthermore, the uncertainty of the 465 updated biogeochemical states is reduced while keeping the dynamical characteristics defined by the prior ensemble simulation. These outcomes indicate that the methodology is promising for future operational applications. Nonetheless, some issues need to be further explored.
According to our forecast experiments, the implementation of ocean colour data assimilation is skillful 470 over short time cycles of about one or two days. Beyond this time, improvements due to the integration of observations into the biogeochemical model become negligible with the currently used PISCES formulation. The dynamics of chlorophyll so appears to be strongly constrained by the prescribed biogeochemical PISCES model parameters and forecasts tend to quickly converge towards the prior ensemble trajectories. An important feature highlighted is the importance of the initial definition of stochastic parameterizations, that 475 Ocean Sci. Discuss., https://doi.org/10.5194/os-2018-153 Manuscript under review for journal Ocean Sci. Discussion started: 24 January 2019 c Author(s) 2019. CC BY 4.0 License. define the model covariances. The efficiency of the analysis strongly depends both on a relevant choice of the biogeochemical sources of uncertainties, and of an appropriate level of uncertainties. In order to improve the assimilation of biogeochemical data, it is therefore suggested to refine stochastic paramaterizations, relatively to the integrated observations. Further experiments should allow to investigate whether using parameter estimation techniques could sustanaibly constrain the primary production.

480
In practice, the proper duration of the ocean colour data assimilation window in an operational system should in the order of a day, while it is in the order of the week for altimetric data assimilation in a mesoscale circulation model. A question that arises from the present study is therefore the controllability of biogeochemical systems. First, it is important to remind that the surface chlorophyll is the only observed variable over 24 that make up the PISCES model. Despite the use of a multivariate scheme with a relevant 485 definition of correlations, fundamental processes such as the vertical nutrient fluxes are not be sufficiently monitored to constrain the phytoplankton growth in the whole euphotic zone. An issue that need to be investigated is therefore the introduction of vertical profiles in addition to surface information. The integration of in situ physical variables data into ocean circulation models have already been investigated in earlier studies (e.g., Griffa et al., 2006). Similarly, other biogeochemical variables should be incorporated to 490 coupled systems. For this purpose, data provided by BIO-ARGO floats (Johnson et al., 2009;Johnson and Claustre, 2016) are essential to consolidate the observing system. Meanwhile, this kind of requirements for data assimilation is an important feature to promote the development of observation networks. The main obstacle for using verticale profiles is the need for sufficient space and time coverage to generate large scale impacts From an algorithmic point of view, the localization scheme used in our analysis will be an issue to 495 assimilate vertical profiles as it could only impact a small radius around the observation points whereas an effective large scale constrain would be required to extract the information of various scales contained in these observations. However, the correction of large scale nutrient fluxes and surface chlorophyll concentrations would presumably improve the model representation of primary production.
On the other hand, uncertainties coming from the physical part of the system are not taken into account 500 in the present configuration despite the critical role of the physics on biogeochemical fluxes. In our system, physics is considered as deterministic and certain. However, the vertical upcoming of nutrients in the euphotic zone depends on the vertical oceanic currents and vertical mixing induced by various physical sources. Thus, it is necessary to correct the physics in order to constrain biogeochemical vertical fluxes. In this sense, uncertainties coming from both the physical and the biogeochemical parts of the system have to 505 be considered. Note that ensemble twin experiments are being develop to address this question (e.g., Yu et al., 2018), and future attempts with the present methodology would help in that task.