Vassilios D Vervatis
email: vervatis@oc.phys.uoa.gr

Pierre De Mey-Frémaux

Nadia Ayoub

John Karagiorgos

Malek Ghantous

Marios Kailas

Charles-Emmanuel Testut

Sarantis Sofianos

Assessment of a regional physical-biogeochemical stochastic ocean model. Part 1: Ensemble generation

Keywords: ensemble modelling, model uncertainties, stochastic physics-biogeochemistry, ocean colour, data assimilation, Bay of Biscay

In this article, Part 1 of a two-part series, we run and evaluate the skill of a regional physicalbiogeochemical stochastic ocean model based on NEMO. The domain covers the Bay of Biscay at 1/36° resolution, as a case study for open-ocean and coastal shelf dynamics. We generate ensembles of uncertainties from assumptions related to errors in the atmospheric forcing, the ocean model parameterizations and in the sources and sinks of the biogeochemical variables. The resulting errors are found to be mainly driven by the wind forcing uncertainties, with the rest of the perturbed forcing and parameters locally influencing the ensemble spread. Biogeochemical uncertainties arise from intrinsic ecosystem model errors and from errors in the physical state. Uncertainties in physical forcing and parameterization are found to have a larger impact on chlorophyll spread than uncertainties in ecosystem sources and sinks. The ensembles undergo quantitative verification with respect to observations, focusing on upper-ocean properties. Despite a tendency for ensembles to be generally under-dispersive, they appear to be reasonably consistent with respect to sea surface temperature data. The largest statistical sealevel biases are observed in coastal regions. These biases hint at the presence of high-frequency error sources currently unaccounted for, and suggest that the ensemble-based uncertainties are unfit to model error covariances for assimilation. Model ensembles for chlorophyll appear to be consistent with ocean colour data only at times. The stochastic model is qualitatively evaluated by analysing its ability at generating consistent multivariate incremental model corrections. Corrections to physical properties are associated with large-scale biases between model and data, with diverse characteristics in the open-ocean and the shelves. Mesoscale features imprint their signature on temperature and sea-level corrections, as well as on chlorophyll corrections due to the vertical velocities associated with vortices. Small scale local corrections are visible over the shelves. Chlorophyll information has measurable impact on physical variables.

Introduction

For a while now, the oceanographic community has been aware that regional and coastal ocean models have some quite specific requirements for both the methods and the observational data used for an effective data assimilation scheme (see extensive reviews e.g. in De Mey-Frémaux, 2000;[START_REF] Edwards | Regional Ocean Data Assimilation[END_REF][START_REF] Pinardi | From weather to ocean predictions: an historical viewpoint[END_REF][START_REF] Fujii | Observing System Evaluation Based on Ocean Data Assimilation and Prediction Systems: On-Going Challenges and a Future Vision for Designing and Supporting Ocean Observational Networks[END_REF]. But the success of data assimilation in a particular case will depend not only on the power of the methodology and data themselves, but on how realistic our estimates of observational errors and model errors are [START_REF] Oke | Representation Error of Oceanic Observations for Data Assimilation[END_REF]. Error processes in a model constitute a subspace of the complete space spanned by all of its variables as the model is integrated. Important sources of model errors in coastal domains are the presence of the coast and shallower water, strong bathymetry gradients, inputs from rivers, and are forced by (among other things) pressure and current fields from ocean-scale mass balances, circulation, tides and eddies, winds and air pressure variations, and non-uniform density (Kourafalou et al., 2015a;2015b).

In addition to the effects of local geography mentioned above, the dynamics of coastal (nested) models are largely controlled by the open boundary conditions [START_REF] Ghantous | Ensemble downscaling of a regional ocean model[END_REF]. [START_REF] Thacker | Propagating boundary uncertainties using polynomial expansions[END_REF] presented how model uncertainties in the open boundaries manifest as model sea-level uncertainties in the Gulf of Mexico. The Bay of Biscay's dynamics are influenced by the North East Atlantic circulation, especially along the southern slope, with the seasonal reversal of the Iberian poleward current. At depth, the entrance of Mediterranean water masses from the south has been shown to influence the Bay of Biscay's hydrology between 600m and 1500m depth [START_REF] Koutsikopoulos | Physical processes and hydrological structures related to the Bay of Biscay anchovy[END_REF], and J o u r n a l P r e -p r o o f potentially the circulation through interactions between eddies and deep, salty lenses (see for instance [START_REF] Carton | Interactions of surface and deep anticyclonic eddies in the Bay of Biscay[END_REF].

One common approach to realistically represent error processes and their time evolution is stochastic modelling (e.g. Adhler et al., 1996; and recent studies for regional configurations in [START_REF] Melsom | Forecasting search areas using ensemble ocean circulation modeling[END_REF][START_REF] Sakov | TOPAZ4: an ocean-sea ice data assimilation system for the North Atlantic and Arctic[END_REF][START_REF] Quattrocchi | Characterisation of errors of a regional model of the Bay of Biscay in response to wind uncertainties: a first step toward a data assimilation system suitable for coastal sea domains[END_REF]Vandenbulcke andBarth, 2015, 2019;[START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. [START_REF] Lucas | Stochastic study of the temperature response of the upper ocean to uncertainties in the atmospheric forcing in an Atlantic OGCM[END_REF] generated ensembles using NEMO (Nucleus for European Modelling of the Ocean; http://www.nemo-ocean.eu/; [START_REF] Madec | Nemo ocean engine[END_REF], based on a configuration for the North Atlantic at 1/4° resolution and multivariate EOF modes. The ensemble capabilities of NEMO have been discussed in the literature following the SANGOMA (http://www.data-assimilation.net/) and OCCIPUT projects [START_REF] Penduff | Ensembles of eddying ocean simulations for climate[END_REF], focusing on global and regional academic configurations spanning from 2° to 1/4° resolution and from seasonal to decadal time scales [START_REF] Brankart | Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modelling[END_REF][START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF][START_REF] Candille | Assessment of an ensemble system that assimilates Jason-1/Envisat altimeter data in a probabilistic model of the North Atlantic ocean circulation[END_REF][START_REF] Garnier | Stochastic parameterizations of biogeochemical uncertainties in a 1/4° NEMO/PISCES model for probabilistic comparisons with ocean color data[END_REF][START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: Application at eddying resolution[END_REF].

In this paper, we configure, run and assess a stochastic ocean model in the Bay of Biscay, consisting of physical-biogeochemical ensemble simulations with the NEMO-PISCESv2 model (Pelagic Interactions Scheme for Carbon and Ecosystem Studies volume 2; [START_REF] Aumont | PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies[END_REF], as a step towards coastal and regional data assimilation. The work is based on recent advances in NEMO explicitly simulating the effects of model uncertainties using an "ensemble generator" (e.g. the same as a scalar time-stepping numerical model can be seen as a "state generator"), modified here specifically for stochastic parameterizations in high-resolution configurations. Our "classic" stochastic modelling approach is based on introducing stochastic degrees of freedom of forcings and model parameterizations (the model response to those stochastic degrees of freedom being the so-called "system errors" whose second-order moments are usually modelled as the Q matrix in a data assimilation context). This stochastic model undergoes validation in Parts 1 and 2 of this study.

A generic perturbation approach based on first-order autoregressive processes -AR(1) -is proposed for the coupled physical and biogeochemical models. Autoregressive processes are based on statistical models operating under the hypothesis that the past state has an effect on the present state. First-order denotes that the current value is based on the immediately preceding value. A comprehensive analysis for the stochastic formulation of NEMO is given in [START_REF] Brankart | Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modelling[END_REF] and [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF]. Some theoretical background for probabilistic ocean modelling, with technical details on implementation strategies based on NEMO (e.g. online ensemble diagnostics, connection with observation operators and data assimilation systems) is provided by [START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: Application at eddying resolution[END_REF]. [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF] introduced two approaches for ensemble ocean modelling: the first is the stochastic perturbed parameterized tendencies (SPPT; [START_REF] Buizza | Stochastic representation of model uncertainties in the ECMWF ensemble prediction system[END_REF] and the second is the stochastic parameterization of unresolved fluctuations (SPUF). The SPPT implementation, which we use in this study, aims to generate perturbations on the models' parameterized tendencies (referred to here as the models' time derivative) and implements Monte-Carlo techniques to obtain a probability density function (pdf) of these tendencies. The stochastically derived parameterized tendencies are added to the models' non-parameterized tendencies (the latter assumed free of uncertainties). The SPUF implementation (not used in this study) is based on random walks sampling gradients from the state vector and adding them to the models' solution. In another study, [START_REF] Ollinaho | Towards process-level representation of model uncertainties: stochastically perturbed parametrizations in the ECMWF ensemble[END_REF] proposed a stochastic perturbed parameters (SPP) scheme, perturbing a number of model parameters whose values are presumed to contain errors. As in their work, model errors of physical parameters are introduced in this study, by applying spatiotemporally varying perturbations.

Marine biogeochemical data assimilation is increasingly being used in operational platforms as a tool to improve ocean forecasting systems. However, the subject is still immature with several challenges remaining, as for example using multivariate increments and validating non-assimilated variables in biogeochemical models. The Ensemble Kalman Filter (EnKF; [START_REF] Evensen | The Ensemble Kalman Filter: theoretical formulation and practical implementation[END_REF] was first used with a simple 1D ecosystem model by Eknes and [START_REF] J O U R N A L P R E -P R O O F Eknes | An ensemble Kalman filter with a 1-D marine ecosystem model[END_REF], and later, by [START_REF] Allen | An Ensemble Kalman Filter with a complex marine ecosystem model: hindcasting phytoplankton in the Cretan Sea[END_REF], to control the evolution of zooplankton and nutrients by assimilating chlorophyll. [START_REF] Simon | Application of the Gaussian anamorphosis to assimilation in a 3-D coupled physical-ecosystem model of the North Atlantic with the EnKF: a twin experiment[END_REF] extended the EnKF to include a Gaussian anamorphosis transformation, accounting for non-Gaussian biogeochemical distributions. [START_REF] Ciavatta | Assimilation of remotely-sensed optical properties to improve marine biogeochemistry modelling[END_REF], instead of assimilating chlorophyll, adapted the J o u r n a l P r e -p r o o f

EnKF to assimilate the SeaWiFS (https://oceancolor.gsfc.nasa.gov/SeaWiFS/) light attenuation coefficient incorporating a bio-optical model. Biogeochemical data assimilation mostly relies on satellite chlorophyll a. Recent studies have shed light on the assimilation of ocean colour plankton functional types, further improving marine ecosystem simulations [START_REF] Ciavatta | Decadal reanalysis of biogeochemical indicators and fluxes in the North West European shelf-sea ecosystem[END_REF][START_REF] Ciavatta | Assimilation of oceancolor plankton functional types to improve marine ecosystem simulations[END_REF].

In light of these advances, the ocean modelling and data assimilation communities have put significant effort in improving the predictive skill of biogeochemical models [START_REF] Song | Data assimilation in a coupled physicalbiogeochemical model of the California current system using an incremental lognormal 4-dimensional variational approach: Part 3-Assimilation in a realistic context using satellite and in situ observations[END_REF][START_REF] Mattern | Data assimilation of physical and chlorophyll a observations in the California Current System using two biogeochemical models[END_REF][START_REF] Mattern | Improving Variational Data Assimilation through Background and Observation Error Adjustments[END_REF][START_REF] Kaufman | Assimilating biooptical glider data during a phytoplankton bloom in the southern Ross Sea[END_REF][START_REF] Yu | Insights on multivariate updates of physical and biogeochemical ocean variables using an Ensemble Kalman Filter and an idealized model of upwelling[END_REF][START_REF] Ford | Assessing the role and consistency of satellite observation products in global physicalbiogeochemical ocean reanalysis[END_REF][START_REF] Goodliff | Temperature assimilation into a coastal ocean-biogeochemical model: assessment of weakly and strongly coupled data assimilation[END_REF]. Meanwhile, there is an increasing interest in biogeochemical model optimization [START_REF] Wang | Assessing the value of biogeochemical Argo profiles versus ocean color observations for biogeochemical model optimization in the Gulf of Mexico[END_REF], as well as for coupled probabilistic systems and methods for their evaluation [START_REF] Candille | Assessment of an ensemble system that assimilates Jason-1/Envisat altimeter data in a probabilistic model of the North Atlantic ocean circulation[END_REF][START_REF] Garnier | Stochastic parameterizations of biogeochemical uncertainties in a 1/4° NEMO/PISCES model for probabilistic comparisons with ocean color data[END_REF][START_REF] Santana-Falcón | Assimilation of chlorophyll data into a stochastic ensemble simulation for the North Atlantic ocean[END_REF]. Errors in biogeochemical models stem from parameterizations and unresolved scales, and can be investigated by uncertainties in initial and atmospheric conditions [START_REF] Verdy | A data assimilating model for estimating Southern Ocean biogeochemistry[END_REF][START_REF] Fransner | Ocean Biogeochemical Predictions-Initialization and Limits of Predictability[END_REF], in biogeochemical tracers and their sources and sinks [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF], and model parameters (e.g. nutrient limitations, growth and mortality rates, grazing etc.; [START_REF] Garnier | Stochastic parameterizations of biogeochemical uncertainties in a 1/4° NEMO/PISCES model for probabilistic comparisons with ocean color data[END_REF][START_REF] Gharamti | Ensemble data assimilation for ocean biogeochemical state and parameter estimation at different sites[END_REF]. Moreover, biogeochemical model performance is strongly dependent on ocean dynamics and on choices made in the assimilation scheme. The use of (pre-) operational biogeochemistry data assimilation is of vital importance to assess these systems and their products, and advance on the management of marine ecosystems [START_REF] Aumont | PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies[END_REF].

The most useful statistical properties of ensembles include their mean, spread, and -especially for data assimilation -their covariances. Higher order moments can also sometimes prove useful (e.g. [START_REF] Quattrocchi | Characterisation of errors of a regional model of the Bay of Biscay in response to wind uncertainties: a first step toward a data assimilation system suitable for coastal sea domains[END_REF]. Where observations are available the ensemble mean can reveal biases, while the ensemble spread can measure model sensitivity to perturbed quantities used to generate the ensemble. When data are available and their uncertainties reasonably well known, consistency analyses can be performed between ensemble-based uncertainties and observational errors [START_REF] Edwards | Regional Ocean Data Assimilation[END_REF], using differences between modelled and observed variable values (the innovation vectors).

We use the above techniques to validate our ensembles, and hence the underlying stochastic model. We describe and implement consistency analysis methods to produce a quantitative assessment of the stochastic modelling of coupled physical and biogeochemical processes. Two other methods to assess the empirical consistency of ensembles will be discussed in a companion article [START_REF] Vervatis | Assessment of a regional physical-biogeochemical stochastic ocean model. Part 2: empirical consistency[END_REF]. We complement this quantitative assessment with a qualitative approach, in which the stochastic model is evaluated by exploring its ability to generate consistent multivariate incremental model corrections. We illustrate this by means of multivariate representers and "stochastic" EnKF incremental analyses for a few members.

The scientific objectives of this paper focus on the generation of ensembles in high-resolution coastal and regional models. The study aims at guiding future ensemble-based modelling strategies, in support of data assimilation and probabilistic forecasting approaches. On the basis of prior knowledge in the literature for stochastic approaches and on what is feasible in terms of computational resources, we investigate physical-biogeochemical model uncertainties and what the impacts of the choice of different sets of perturbations and observational networks are on the analysed ocean state. In Section 2, we begin by describing the coupled physical-biogeochemical model and the methods to evaluate the stochastic model with respect to satellite observations. The ensembles and the model-data consistency assessment are discussed in the results section (Section 3). A discussion and concluding remarks are presented in Section 4.

Methods and data

The coupled physical-biogeochemical deterministic ocean model

We use the release 3.6 of the NEMO ocean model in a regional configuration encompassing the Bay of Biscay and the western part of the English Channel (Fig. 1). The configuration, named BISCAY36, is a subgrid of the operational Copernicus Marine Environment Monitoring Service (CMEMS) Iberia-Biscay-Ireland (IBI) configuration [START_REF] Sotillo | The MyOcean IBI Ocean Forecast and Reanalysis Systems: operational products and roadmap to the future Copernicus Service[END_REF] previously applied and validated [START_REF] Quattrocchi | Characterisation of errors of a regional model of the Bay of Biscay in response to wind uncertainties: a first step toward a data assimilation system suitable for coastal sea domains[END_REF][START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. The horizontal resolution is 1/36° (about 3 km and 2 km in the meridional and zonal directions respectively). For a complete description of the original BISCAY36 configuration, and for validation details, the reader is referred to [START_REF] Maraldi | NEMO on the shelf: Assessment of the Iberia-Biscay-Ireland configuration[END_REF].

J o u r n a l P r e -p r o o f

In our BISCAY36 implementation, the NEMO ocean engine OPA (Océan Parallélisé) is coupled on-line with the passive tracer package TOP2 and the biogeochemical model PISCESv2. The model solves for 24 prognostic variables simulating the biogeochemical cycles of oxygen, carbon and the main nutrients controlling phytoplankton growth: nitrate, ammonium, phosphate, silicic acid and iron. The model distinguishes four plankton functional types based on size, including two phytoplankton compartments (nanophytoplankton and diatoms) and two zooplankton classes (microzooplankton and mesozooplankton). The distinction of the two phytoplankton size classes, along with the description of multiple nutrient co-limitations, allows the model to represent ocean productivity across different biogeographic ocean provinces [START_REF] Longhurst | Ecological Geography of the Sea, 398[END_REF].

The biogeochemical model tracks phytoplankton biomass in units of carbon, iron, silicon (the latter only for diatoms) and chlorophyll. The ratios between carbon and the latter three variables can change, which allows for a more accurate conversion from phytoplankton to chlorophyll concentrations, which is of great importance for comparisons with ocean colour satellite data. The ratio between carbon, nitrogen and phosphorus is kept constant at the Redfield values of 122:16:1 following [START_REF] Takahashi | Redfield ratio based on chemical data from isopycnal surfaces[END_REF].

PISCESv2 also distinguishes three non-living pools for organic carbon: small Particulate Organic Matter (sPOM), big Particulate Organic Matter (bPOM; different settling velocities from sPOM) and semi-labile Dissolved Organic Carbon (DOC).

The physical and biogeochemical models are coupled on-line, using a high coupling frequency for the conservation of tracers of once every two time-steps, i.e. 150 s for physics and 300 s for biogeochemistry.

The coupling is one-way, so that the ocean forcing is applied to the biogeochemical model but there is no feedback. The primitive equations and the tracer transport model are discretized on the same 1/36° curvilinear Arakawa C-grid. The TOP2 package controls the advection-diffusion equations of the passive and biogeochemical tracers. The numerical scheme for the biogeochemical processes is forward in time (Euler) and differs from the classical leap-frog scheme used for the physics. The advection scheme is the same as for the physics, i.e. QUICKEST [START_REF] Leonard | A stable and accurate convective modelling procedure based on quadratic upstream interpolation[END_REF], but using the limiter of [START_REF] Zalesak | Fully multidimensional flux corrected transport algorithms for fluids[END_REF]. These options have been tested by [START_REF] Gutknecht | Evaluation of an operational ocean model configuration at 1/12 spatial resolution for the Indonesian seas (NEMO2.3/INDO12)-Part 2: Biogeochemistry[END_REF] and are now used in the IBI-MFC operational system (http://marine.copernicus.eu/).

The meteorological fields are provided by the ECMWF (European Center for Medium-Range Weather Forecasts). The initial state and open boundary conditions are acquired from the daily archives of the CMEMS infrastructure for physics, and the weekly archives for biogeochemistry. Details of the physical model set-up are described in [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. The open boundaries of the biogeochemical model are forced by the global system BIOMER4V1R1 (resolution: 1/2°; http://www.mercator-ocean.fr/), which provides (via the CMEMS archives) a 3D global weekly mean analysis for dissolved iron, nitrate, phosphate, silicate, oxygen, chlorophyll and phytoplankton concentrations expressed as their carbon contents.

The stochastic model

The main issues to be faced when building a stochastic model is first to identify and select the relevant quantities to perturb (depending on the dynamics and on the objectives of the stochastic model) and secondly to generate the perturbations. This is what we propose to investigate in this study, the general context being an operational system assimilating satellite data in coastal and regional domains.

The ensembles of simulations are generated by perturbing some quantities that are considered as major sources of errors for the model. These quantities consist of forcing fields at the boundaries, parameters and state variables; these are referred to as "parameterized tendencies" (see [START_REF] Buizza | Stochastic representation of model uncertainties in the ECMWF ensemble prediction system[END_REF], in the sense that the time evolution of the model state variables can be written as follows: 𝜕 𝑀 𝒙, 𝒖, 𝒑, 𝑡 𝐷 𝒙, 𝑡 𝑃 𝒙, 𝒖, 𝒑, 𝑡 , where 𝜕 𝑀 𝒙, 𝒖, 𝒑, 𝑡 is the time derivative of the model 𝑀 𝒙, 𝒖, 𝒑, 𝑡 , as a function of the model state vector 𝒙 (e.g. temperature, salinity or currents), the forcing 𝒖 and the vector of the model parameters 𝒑. 𝐷 𝒙, 𝑡 is the dynamics tendency of the ocean state for well-resolved, nonparameterized processes, equal to the advection term 𝛻 ⋅ 𝒙𝑼 , where 𝑼 the vector of the Eulerian velocity, and 𝑃 𝒙, 𝒖, 𝒑, 𝑡 is the tendency of parameterized processes. The latter is perturbed based on assumptions on the amplitude and space-time structure of atmospheric forcing uncertainties, model errors in physical parameterizations and sources-minus-sinks 𝑆𝑀𝑆 biogeochemical model errors (state J o u r n a l P r e -p r o o f variables are not explicitly perturbed in our experiments); the scientific basis of these assumptions is discussed in the next sections.

In this study, the perturbations, hereafter noted 𝜉, are computed using the stochastic perturbed parameterized tendencies (SPPT; [START_REF] Buizza | Stochastic representation of model uncertainties in the ECMWF ensemble prediction system[END_REF] and the stochastic perturbed parameters (SPP; [START_REF] Ollinaho | Towards process-level representation of model uncertainties: stochastically perturbed parametrizations in the ECMWF ensemble[END_REF] schemes. Assuming that the past value of the uncertainties affects the present one, the 𝜉 are modelled as first-order auto-regressive processes. In practice, at every time-step, Gaussian AR (1) processes are generated following the expression defining:

𝜉 𝑒 • 𝜉 𝜎 • 1 𝑒 • 𝑤 𝜇 • 1 𝑒 (1)
where 𝑘 is the model time-step, 𝑤 is white Gaussian noise, 𝜇, 𝜎 and 𝜏 the mean, standard deviation (uncertainty amplitude) and correlation timescale, respectively. Here, we assume the errors are unbiased, so 𝜇 is set to zero. The perturbation fields are different for each perturbed tendency or parameter, and vary under the assumption of spatial and temporal correlated scales. Also, the fields are bounded 𝜉 ∈ 1 1 to retain the sign of the tendency [START_REF] Palmer | Stochastic Parameterization and Model Uncertainty[END_REF]. [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF] generated independent Gaussian autoregressive processes at every model grid point and introduced spatial dependence in the 𝜉 2D or 3D fields, by applying a Laplacian filtering operator ℱ to the white noise 𝑤 in Eq. ( 1) as 𝜉 ∝ ℱ 𝑤 . The Laplacian filter implemented in NEMO is equivalent to the five-point smoothing operator described by [START_REF] Shapiro | Smoothing, filtering, and boundary effects[END_REF]. This is a low-pass, first-order Shapiro filter used to remove small scale grid noise. When applied to white noise on a coarse-resolution domain (typically of a few tens of kilometres), a few Laplacian passes are sufficient to introduce correlation at scales of a few hundreds of kilometres. In this context, the Laplacian-Shapiro filter acts to "blur" the white noise at every model grid point and is computationally efficient [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF][START_REF] Garnier | Stochastic parameterizations of biogeochemical uncertainties in a 1/4° NEMO/PISCES model for probabilistic comparisons with ocean color data[END_REF]. By contrast, in high-resolution configurations of a few kilometres (e.g. BISCAY36), the Laplacian-Shapiro filter is not an optimal approach to attenuate noise and produce long-range correlations (cf. Appendix A; Fig. A1); therefore, another approach is required.

Here, long-range spatial correlations (on the order of 10 to 100 km) are explicitly calculated by applying a 2D Gaussian function to estimate iso-correlation contours with a pdf given by the mathematical expression:

𝑝 𝑿|𝑿 , 𝜮 ∝ 𝑒𝑥𝑝 𝑿 𝑿 𝜮 𝑿 𝑿 (2)
where 𝑿 𝑥, 𝑦 represents the whole model domain, 𝑿 𝑥̅ , 𝑦 is the central grid point of the distribution (the mode) and 𝜮 𝑑𝑖𝑎𝑔 𝜎 , 𝜎 is the bivariate diagonal covariance kernel controlling the length scales of the perturbation. The spatial covariance kernel has anisotropic (important in coastal domains) e-folding length scale variances in the meridional 𝜎 and zonal 𝜎 directions, varying randomly around a typical correlation length per variable and across the ensemble members. In order to generate a multimodal pattern (with two or more local maxima; 𝓂 2), we recursively implement the Gaussian function at 𝓂 random central grid points 𝑿 in the model domain, resulting in a finite sum of Gaussian distributions 𝑝 𝑿|𝑿 𝓂 , 𝜮 𝓂 ∝ 𝑝 𝑿 𝑿 𝓂 , 𝜮 𝓂 ∑ 𝑝 𝑿 𝑿 , 𝜮 𝒊 𝓂 . The recursive implementation of the function further increases the perturbation's anisotropy. The number of recursive iterations 𝑖 1, … , 𝑚 (and therefore the number of modes) depends on the correlation length within a given model domain. For instance, a correlation length half the size of the model domain yields approximately a bimodal distribution. The Gaussian function in Eq. ( 2), is used as a squared exponential (apodization) covariance function, so instead of the normalization factor 1 2𝜋|𝛴| / ⁄ we used the perturbation uncertainty amplitude defined by 𝜎 in Eq. (1). Stochastic perturbations are generated by substituting the Laplacian convolutions in the white noise 𝑤 in Eq. ( 1), with a multimodal Gaussian pattern spanning the whole model domain i.e. ℱ 𝑤 → 𝑝 𝑿|𝑿 𝓂 , 𝜮 𝓂 . AR (1) processes estimated by Eq. (1), are performed at every model grid point shifting and altering the Gaussian modes in space and time. The implementation is compatible with the NEMO MPI double parallelization environment in the J o u r n a l P r e -p r o o f spatial domain and across an ensemble [START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: Application at eddying resolution[END_REF]. The algorithm is available from the opensource software platform Zenodo with doi:10.5281/zenodo.2556530.

Atmospheric forcing uncertainties

The atmospheric forcing in coastal and regional applications constitutes a major source of ocean model uncertainty and is likely to lead to large-scale biases [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. We investigate here uncertainties in the wind velocities 𝑈 𝑢 , 𝑣 , the sea level pressure 𝑆𝐿𝑃) and the air temperature 𝑇 , i.e. 𝑢 ∈ 𝑈 , 𝑆𝐿𝑃, 𝑇 . The ECMWF atmospheric fields 𝑢 are multiplied by AR(1) stochastic processes 𝜉 (Figs. 2a-f) following the SPPT scheme:

𝑢 → 1 𝛼 • 𝜉 • 𝑢 (3)
at every time-step 𝑘, where 𝛼 is an optional tapering value in the interval 0 1 [START_REF] Buizza | Stochastic representation of model uncertainties in the ECMWF ensemble prediction system[END_REF]. To calculate the AR(1) stochastic parameterization 𝜉 using Eq. ( 1), we select representative values for the uncertainty amplitude 𝜎, the average spatial correlation length 𝜎 , used in Eq. ( 2) and the correlation timescale 𝜏 for each tendency.

Output diagnostics from a one-year control run are analysed to estimate the distribution of pdfs (not shown) and tune the stochastic parameterizations for all atmospheric variables of interest. Wind velocity 𝑈 and 𝑆𝐿𝑃 are distributed normally, whereas the distribution 𝑇 is bimodal due to the seasonal cycle.

We assume that uncertainties are related to synoptic timescales, e.g. atmospheric phenomena such as storms, and we set a temporal correlation length of a few days for all atmospheric variables. The synoptic timescales are also verified by time-lagged autocorrelation methods applied to the control run. The spatial scales of the atmospheric fluctuations are determined over synoptic timescales. The signal-to-noise ratio is assigned according to the statistical laws of the pdfs. These stochastic parameterizations for Eq. ( 1) are in agreement with other studies in the literature (e.g. [START_REF] Palmer | Stochastic Parameterization and Model Uncertainty[END_REF].

Model uncertainties in physical parameterizations

Fluxes of momentum, heat and mass are the key quantities linking the air and sea. The physical processes related to them are parameterized in terms of bulk coefficients, which are deduced from empirical laws, incorporating wind speed dependent coefficients and feedback from the sea state. In this study, we assume that model errors are in part due to limitations of these parameterizations of air-sea interaction. Stochastic perturbations varying spatiotemporally are imposed on the models' momentum drag 𝑐 , latent heat 𝑐 and sensible heat 𝑐 coefficients. The AR(1) distribution and temporal scales are the same as those we prescribed for the wind; the spatial scales on the other hand, are assumed to be that of the ocean state and are set to a few Rossby radii of deformation. The same stochastic pattern is applied to all coefficients in the CORE bulk formulae [START_REF] Large | Diurnal to decadal global forcing for ocean and sea-ice models: the data sets and flux climatologies[END_REF]. The positiveness of the coefficients is verified by tapering methods and for different stability conditions and wind speed regimes in the bulk formulae. Like Eq. (3), the SPP scheme introduces spatiotemporal perturbations following a normal distribution of AR (1) processes (other options may include simpler perturbations of white noise uncorrelated in space and time), expressed as:

𝑝 → 1 𝛼 • 𝜉 • 𝑝 (4)
In order to represent model uncertainties for tidal mixing over the shelves, stochastic fluctuations based on flux boundary conditions are also imposed at the bottom layer. The bottom drag parameterizations 𝑐 are based on model assumptions for the vertical shear, the mixing scheme and the nature of the seabed (rocky, sandy or muddy), when known, which modify the bottom boundary layer. Because in many ocean models the bottom drag is approximated as a permanent feature (e.g. constant minimum values in the abyssal plain as in [START_REF] Maraldi | NEMO on the shelf: Assessment of the Iberia-Biscay-Ireland configuration[END_REF], large temporal scales up to one month are imposed in Eq. ( 1). The true nature of the seabed is unknown and consequently so are the dominant scales in the bottom layer, so we simply apply white noise and Laplacian filtering to introduce AR(1) spatial scales. The formulation of the bottom drag follows a quadratic log law, with minimum positive values clamped at 2.5 10 in the abyssal plain and maximum values observed in the shallow areas of the English Channel.

Biogeochemical sources-minus-sinks model uncertainties

J o u r n a l P r e -p r o o f

Marine biogeochemical models encompass many sources of uncertainty stemming from unresolved model processes and poorly constrained parameterizations. These uncertainties fall into two broad categories: unresolved biodiversity and unresolved scales. The first category refers to the biodiversity restriction of only a few tens of species resolved by the model in an effort to reduce state variables [START_REF] Quéré | Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models[END_REF]. This category also includes errors in the parameterization or missing biogeochemical processes, those that control the feedback between sub-systems (ecosystem, chemistry, oxygen and carbonate models; [START_REF] Gehlen | Biogeochemical consequences of ocean acidification and feedbacks to the earth system[END_REF]. These errors emerge from having a limited number of compartments, which often leads to a crude parameterization of their processes. The second categoryerrors due to unresolved scales -includes errors related to the model resolution and the species included in the model (e.g. lack of diurnal migration in the model), and the method for handling them is similar to that for small scale unresolved physical processes.

The SPUF scheme appears to be the most natural method to simulate uncertainties for both categories. However, one must consider that performing random walks with a large state vector could be computationally expensive. Hence, as with the other sources of error, unresolved biodiversity is here explored via the SPPT AR( 1) scheme [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF][START_REF] Garnier | Stochastic parameterizations of biogeochemical uncertainties in a 1/4° NEMO/PISCES model for probabilistic comparisons with ocean color data[END_REF], introducing uncertainties in the sources-minus-sinks terms for each of the 24 biogeochemical tracers 𝐶, using the stochastic field 𝜉 described in Eq. ( 1) at every time step 𝑘, as

𝑆𝑀𝑆 𝐶 → 𝑆𝑀𝑆 𝐶 • 𝑒 . ( 5
)
As an example of this term, phytoplankton sources and sinks are photosynthesis, respiration, death and grazing. The biogeochemical model parameters themselves are not perturbed. 2016) used SeaWIFS at 1/12° resolution to assess the unresolved scales of a coarser 1/4° model and tune the AR(1) stochastic parameterizations. In this study, we are unable to assess the unresolved scales of the high-resolution BISCAY36 model, with respect to a dense observational network. Therefore, stochastic parameterizations are deduced from a one-year control run, mainly focusing on the statistical properties of chlorophyll. The sources-minus-sinks perturbations are assumed to have a lognormal distribution and large uncertainty amplitude 𝜎 of up to 60%. The bias correction term 𝜎 2 ⁄ [START_REF] Simon | Application of the Gaussian anamorphosis to assimilation in a 3-D coupled physical-ecosystem model of the North Atlantic with the EnKF: a twin experiment[END_REF] is part of the model tuning, implying that the ensemble mean is equal to the unperturbed 𝑆𝑀𝑆 term. The spatial scales are representative of a few Rossby radii of deformation and the correlation time-scale is set to10 days, typical for the underlying mesoscale dynamics and associated upwelling patterns. We perturb the state of the sources-minus-sinks terms in all three dimensions and across all vertical levels with the same stochastic pattern. Sensitivity experiments showed that this approach yielded uncertainty regimes with an ensemble spread increasing in time. We did test alternatives by applying different stochastic patterns per level or per tracer but these degraded the impact of the method.

Garnier et al. (

Experimental design

In order to test which error sources are most significant for the biogeochemical uncertainties, and to better understand the couplings within the stochastic model, we designed three ensemble experiments: perturbing only the physics (EnsP), only the biogeochemistry (EnsB), and perturbing both simultaneously (EnsPB). In the coupled simulation, the evolution of the biogeochemical tracers is described by the advection-diffusion equation:

𝜕 𝐶 𝛻 𝑢 • 𝐶 𝐾 𝛻 𝐶 𝜕 𝐾 𝜕 𝐶 𝑆𝑀𝑆 𝐶 (6)
where on the right hand of Eq. ( 6) the first term represents the advective transport of tracers along isopycnals, and the second and third terms the 3D parameterized diffusion processes. The last term denotes all biological processes affecting the concentration of tracers 𝐶, due to the sources-minus-sinks terms 𝑆𝑀𝑆 𝐶). The three experiments defined above are indicated schematically in Eq. ( 6).

J o u r n a l P r e -p r o o f

In Table 1, we summarize our settings for the stochastic parameterizations in the sensitivity experiments perturbing individual sources of model uncertainty and the seasonal-range ensembles EnsP, EnsB, EnsPB. We also include a reference ensemble (Ens0) and a control run (CR). A detailed presentation of the Ens0 ensemble is given by [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF], but it suffices to say here that it differs fundamentally from the EnsP, EnsB, EnsPB ensembles in that it was generated by performing stochastic modelling of the wind forcing based on empirical orthogonal functions (EOF); there is no other variable perturbed in Ens0.

A major challenge in ensemble forecasting is to identify the most important sources of model error and properly configure the stochastic parameterizations. In order to tackle this, we adopted the following twostep procedure: firstly, we identified which error sources the model was most sensitive to in the Bay of Biscay, by running medium-range ensembles of 20 members for one month during the spring bloom, namely April, 2012 (S1-5, SP, SB, SPB experiments in Table 1). Several stochastic parameterizations, deduced from the output diagnostics of a one-year control run, were tested on a case-by-case basis for each sensitivity experiment. A final selection of perturbations was established, based on their impact on the spread of upper-ocean variables as well as on their computational efficiency (e.g. the requirement that no members blow up during the ensemble run). Based on these experiments, we decided not to perturb variables and parameters that the model in this configuration is insensitive to, such as, for example the photosynthetically active radiation coefficient 𝑘 for the penetrative solar radiation (sensitivity experiments not shown). Results from the sensitivity experiments are discussed in Section 3.1. For the second step, we used these perturbations to run seasonal-range ensembles of 40 members each (EnsP, EnsB, EnsPB experiments in Table 1). The stochastic model skill of the seasonal-range ensembles is evaluated with respect to observations in Sections 3.2 and 3.3.

Figure 2 shows examples of AR( 1) stochastic patterns drawn for several variables of the coupled system, in order to introduce model errors. We generate wind uncertainties with smaller spatial scales, varying around an average of 𝜎 , ~1° (Figs. 2a-d), compared to the work presented by [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF], who used in their study broader patterns based on dominant EOF modes (cf. Appendix A; Fig. A2). Figures 2a-c depict the spatial variations of the wind stochastic pattern for member-001, over a temporal correlation period of three days. In the same simulation stochastic patterns with larger spatial scales can be applied to other atmospheric variables: for example, scales of 𝜎 , ~2° and 𝜎 , ~3° have been adopted for the air temperature and sea level pressure respectively (Figs. 2e-f). Figure 2h illustrates the effect of the Laplacian filtering after 100 iterations in BISCAY36. It shows noisy spatial patterns, not representative of most oceanic processes.

In the same vein as [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF], a free run (i.e. the CR in Table 1) is carried out with no stochastic parameterizations. Starting from the PSY2V4R2 analysis (1/12° resolution; [START_REF] Lellouche | Evaluation of global monitoring and forecasting systems at Mercator Océan[END_REF] http://www.mercator-ocean.fr/) the free run is used to generate a five-month spin-up, from July to November 2011, to allow the ocean model to develop coherent structures and represent the main physical processes in the Bay of Biscay. The same period is used for the spin-up of the on-line coupled biogeochemical model. The post spin-up solution is used as the initial conditions for all other experiments. The free run is then extended from December 2011 to June 2012 and serves as a reference for the abovementioned ensemble experiments; hence, we refer to it as the one-year control run. In Fig. 3, we present the ocean state of the control run on April 30, 2012, for a number of surface variables: sea surface height (SSH), sea surface temperature (SST), sea surface salinity (SSS), total chlorophyll, and nanophytoplankton and diatom concentrations. In supplementary material, we verify the model's skill to simulate subsurface nutrient conditions compared to world ocean atlas climatology [START_REF] Garcia | World Ocean Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate and nitrate+nitrite, silicate). A. Mishonov Technical Ed[END_REF].

The seasonal-range ensembles were run from December 2011 to June 2012. The computational resources for the ensemble experiments are presented in Appendix B. All sensitivity (S1-5, SP, SB, SPB) and seasonal-range (EnsP, EnsB, EnsPB) ensembles were initialized by using the ocean and the biogeochemical states of the control run, with the exact same initial conditions. No data was assimilated and the perturbation mechanism remains at work throughout the whole simulation period. A total number of 40 ensemble restarts was archived, where pseudo-random seed numbers (with different random sequences) were saved for all stochastic parameterizations. These ensemble restarts were used to ensure the reproducibility of the stochastic patterns upon initialization of the ensembles. This capability allowed J o u r n a l P r e -p r o o f us to investigate the growth rate of model uncertainties, with the same forecast lead time, under different atmospheric and ocean states.

Datasets

We focus on upper-ocean properties for sea level anomaly (SLA), SST and surface total chlorophyll. The observational networks were accessed via the CMEMS infrastructure (http://marine.copernicus.eu/) and are summarized in Table 2. The observational errors include both measurement and representativity errors, and are usually unknown. Here, the observational errors were taken from those reported at CMEMS and averaged to create a single representative value for each product (Table 2).

High-resolution SST observations are necessary when it comes to validating ensembles at eddy-resolving scales, such as in BISCAY36 at 1/36° resolution. We use the OSTIA SST L4 gap-free gridded dataset which is a near-real-time daily-mean product of foundation SST free of diurnal variability. The data set is provided with an estimate of their uncertainty which is 0.5°C on average. The model proxy for the foundation SST has been chosen to be the temperature interpolated at 10 metres depth.

For the sea-level, we use the CMEMS L3 along-track product based on several altimetry missions and which provides SLA at a resolution of 14 km with uncertainties of 0.05 m. Both data and model include tides and both observed and simulated SLA are detided. The model includes pressure forcing, and therefore an inverted barometer correction is applied to the model and observations (a dynamic atmospheric correction is applied to the sea-level observations). In order to calculate the SLA model equivalent, we use the mean dynamic topography of the parent model IBI36 (Benkiran, pers. comm.).

Of central importance to this work is the consistency analysis of biogeochemical ensembles, where physics and biogeochemistry are both perturbed, against ocean colour products. The ocean colour method exploits different radiation wavelengths and reflectances emitted from the sea surface, affected by phytoplankton and corresponding to different water types. In this study, we use surface total chlorophyll produced for the Global Ocean in the framework of the ESA ocean colour climate change initiative (OC-CCI), made available through CMEMS. This is a merged data records product collected from multiple sensors and ocean satellite passages provided in gridded format at 4 km resolution. The ocean colour L4 product is reconstructed from L3 reprocessed daily composites applying 8-day temporal averaging to fill in missing data, though some gaps remain.

The satellite chlorophyll a data are provided with a scaled (%) observational error in comparison to the signal and therefore, the error has spatial distribution. We decided not to use this information because it would be harder to interpret pattern-consistency and because here we mostly focus on the investigation of model errors. Instead, we estimated a static error at 0.3 mg/m 3 , which is moderately larger (smaller) in the abyssal plain (shelf areas) compared with the scaled error. The ocean colour chlorophyll proxy in models is often taken as an average of chlorophyll over the top 10% of the euphotic layer [START_REF] Prunet | Assimilation of surface data in a onedimensional physical-biogeochemical model of the surface ocean: 1. Method and preliminary results[END_REF]. In our case, estimating a euphotic layer of approximately 50 m, we decided to tune the observation operator to return a model proxy as the mean value of the first 5 m of the water column. This is in agreement also with calculations based on 1/𝑘 .

2.5

Consistency analysis framework

Ensemble-based innovation metrics

In Section 3.2, we evaluate the stochastic model, by quantifying the ability of the ensembles at representing the model's errors with respect to observations. Our method assesses the statistical consistency of the following distributions: (a) ensemble-based uncertainties, and (b) empirically-based uncertainties for upper-ocean properties. As the system is dynamically heterogeneous, the consistency is checked locally.

The empirically-based uncertainties are represented by an ensemble of observations for each dataset. The observations are perturbed with a Gaussian random number to generate data distributions with the proper error standard deviation for each network. All observational errors are considered independent (no crosscorrelations) for each dataset so that their error covariance matrix is diagonal (a useful, but not entirely true assumption). For the consistency analysis of EnsP perturbing only the physics, we used the SST and J o u r n a l P r e -p r o o f SLA observational networks. The chlorophyll a product is used for the consistency check of EnsB and EnsPB biogeochemical model uncertainties. The statistical consistency between the simulated and empirical ensembles is evaluated by comparing the shape of the ensemble envelopes locally and as a function of time. For each dataset, we also compute the innovation vector, that is the model-data misfits, and compare the second-order moments: the innovation variance across the ensemble is the result of uncertainties in both the model and observations. Under the assumption that errors in the model and in observations are uncorrelated, the innovation variance should be close to the sum of the model uncertainty variance (estimated from the ensemble variance) and the observational error variance. The diagnosis proposed by [START_REF] Desroziers | Diagnosis of observation, background and analysiserror statistics in observation space[END_REF] for data assimilation is based on a similar assumption.

Denoting 𝑦 , the measured value of observation 𝑛 for ensemble member 𝑗, and 𝑦 , the model forecast equivalent for observation 𝑛 for member 𝑗, the following ensemble-based consistency statistics in data space may be calculated:

𝑂𝑚𝐸 𝑦 , 𝑦 , , (7)
𝑂𝑚𝐸 𝑦 , 𝑦 , , (8)
𝐷𝑂𝑚𝐸 𝑦 , 𝑦 , 𝑦 , 𝑦 , (9)
where the 𝑗 overbar averages are taken over the ensemble members 𝑗 1, … , 𝑁 and the 𝑛 overbar averages are taken over all observations 𝑛 1, … , 𝑝 of a given type within a given region, and the combined 𝑗, 𝑛 overbar denotes a double average taken first over the ensemble and second over all observations of a given type within a given region. 𝑂𝑚𝐸 stands for Observation minus Ensemble metric (it is an ensemble-wide extension of the so-called Norm of the innovation vector used in data assimilation). 𝐷𝑂𝑚𝐸 stands for Debiased 𝑂𝑚𝐸 metric defining the mean spread of innovation, which should be larger than the ensemble spread, and can be used even if the 𝑂𝑚𝐸 is not small. For an unbiased system, i.e. no biases between observations 𝑦 , and model 𝑦 , , ∀𝑛, we get 𝑂𝑚𝐸 0 and 𝐷𝑂𝑚𝐸 𝑂𝑚𝐸. In Appendix C, we provide a detailed description of the notation for the ensemble-based innovation metrics.

In addition, we calculate the ensemble spread and the quantiles in data space to assess the ensemble median Q2(50%), the mid-spread Q1(25%)-Q3(75%) and the ensemble outermost quantiles Q0(1%)-Q4(99%).

Let us also define the support of a probability density function (pdf) as the smallest closed set outside of which the pdf vanishes. For a pdf defined in ℝ, the pdf envelope is defined as the range between the minimum and maximum values of the support.

In order to illustrate the dependency of assessment results on geographic region, we focus on two distinct areas in the Bay of Biscay, namely the abyssal plain and the Armorican shelf. Both areas exhibit strong spatiotemporal variability, governed by diverse physical-biogeochemical open-ocean and coastal shelf processes, respectively.

Multivariate incremental analysis

The stochastic model skill is qualitatively evaluated by means of multivariate representers (defined below) and EnKF-type incremental analyses [START_REF] Evensen | The Ensemble Kalman Filter: theoretical formulation and practical implementation[END_REF]. The incremental analysis corresponds to the correction step in the sequential data assimilation scheme; in other words, at a given date, we compute the correction to the controls using an estimate of the model and data errors and assuming a given value for the innovations, but then we do not restart the model from the corrected state. This is what makes it different from a data assimilation experiment. A by-product will be to illustrate the potential impact of J o u r n a l P r e -p r o o f our ensemble-modelled uncertainties on data assimilation at analysis time. In all cases, we investigate the impact of observations on unobserved variables, such as other data types or subsurface variables. The novelty in this exercise is to showcase how stochastic modelling of different components in the coupled physical-biogeochemical system contributes to incremental analyses. We emphasize the fact that biogeochemical model errors in Eq. ( 6) have been generated in simulation experiments of increasing complexity, perturbing separately the physics and biogeochemistry. The calculations are conducted with SDAP (Sequoia Data Assimilation Platform, https://sourceforge.net/projects/sequoia-dap/), whose functions were expanded to interface with the NEMO platform and its biogeochemical component PISCESv2.

Let us define 𝑯 to be the observational operator and 𝑷 the error covariance matrix of the prior controls. In this study, the controls are the following state variables: 𝑆𝑆𝐻, 𝑇, 𝑆, 𝐶ℎ𝑙. In an ensemble-based context, the prior error covariance matrix is approximated by 𝑁 samples (i.e. members) in a decomposed form as 𝑷 𝑺 𝑺 , where

𝑺 √ ∏ 𝒙 𝒙 (10)
is the square root matrix of the error-subspace, with ensemble mean

𝒙 ∑ 𝒙 , ( 11
)
𝒙 being the control vector.

The representers are the influence functions of single observations: the correction added to the forecast during the analysis step can be written as a linear combination of representers weighted by the innovation [START_REF] Bennett | Generalized inversion of a global numerical weather prediction model[END_REF]. The matrix of representers is 𝑯𝑷 𝑯 . In Section 3.3, we analyze zero-lag representers, i.e. where observations and controls are at the same time (operator 𝑯 has no time component).

The corrections are computed from the Kalman gain matrix

𝑲 𝜚 ° 𝑷 𝑯 𝑯 𝜚 ° 𝑷 𝑯 𝑹 (12)
multiplied by the member-dependent innovation vector 𝒅 𝒚 𝐻 𝒙 , where 𝒚 denotes the (memberindependent) observation vector. The symbol ° denotes the Schur (element-wise) product of two matrices and 𝜚 is the localization function [START_REF] Gaspari | Construction of correlation functions in two and three dimensions[END_REF]. Here, unlike the analysis step of a "stochastic" EnKF (and from the previous Section 2.5.1), we chose not to perturb the observations (i.e. as in a "deterministic" EnKF, cf. [START_REF] Sakov | A deterministic formulation of the ensemble Kalman Filter: an alternative to ensemble square root filters[END_REF][START_REF] Yu | Insights on multivariate updates of physical and biogeochemical ocean variables using an Ensemble Kalman Filter and an idealized model of upwelling[END_REF]. There is a trade-off -using unperturbed observations does lead to the loss of statistical consistency for the second-order moments and forbids calculating innovation-based diagnostics, but properly calculated EnKF increments for individual members would be harder to interpret physically because of the observation noise. Convolution with a localization function in Eq. ( 12) constrains spurious long-distance correlations resulting from the small ensemble size. We used a localization Gaussian function with cut-off radius 3° and e-folding scale 0.2°. Vertical localization is not applied to make results as "raw" and easy to interpret as possible.

In Section 3.3, we compute and analyse the multivariate EnKF-type incremental analyses 𝑲 • 𝒅 for the first two members (001/002) of the three long-range ensembles (EnsP, EnsB and EnsPB), as well as for different subsets and seasons, using multiple observations together. Thinning techniques are applied to high density SST and ocean colour L4 observational networks [START_REF] Zuo | The ECMWF operational ensemble reanalysis-analysis system for ocean and sea ice: a description of the system and assessment[END_REF], in order to reduce the computational time and optimize data assimilation performance.

We define specific criteria to assess the skill of the multivariate representers and the incremental analyses, in the context of different perturbation methods (i.e. Ens0 and EnsP, EnsB, EnsPB). First, we examine whether the increment correction is large constraining the model ensemble. A subsequent criterion is whether the increment (regardless of its magnitude) has a physical-biogeochemical meaning attributed to specific dynamical processes in the area, showing also potential to constrain the ensemble. Additional criteria are used to investigate the consistency of the increments (a) when calculated from a smaller number of members (testing whether we have convergence of covariances), (b) for different periods (e.g.

J o u r n a l P r e -p r o o f flow dependent errors relevant to different seasons), (c) when assimilating different datasets (e.g. SST and chlorophyll a) and (d) their impact on unobserved variables (e.g. SSS).

Results

Ensemble sensitivity experiments

To assess error patterns prior to the seasonal-range ensemble simulations, we ran a series of mediumrange ensemble simulations. These small ensembles are useful for tuning the stochastic parameterizations and to verify that their spreads are moderate compared to the seasonal-range ensembles with longer spinup period.

In Fig. 4, we present maps for April 30, 2012, of model uncertainties for all medium-range ensembles, S1-5, SP, SB and SPB, after a one-month spin-up. If we compare the effect on SST of wind perturbations alone (Fig. 4a) to all perturbations together (Fig. 4f) we see little difference; likewise, for the SSH (Figs. 4g,l). This implies that error regimes for SST and SSH are mainly driven by the wind forcing. Wind uncertainties have a large impact on upper-ocean uncertainties for both the geostrophic and the Ekman components, defined by Sverdrup dynamics. Imposing the same perturbation field for both u and v wind velocities, i.e. not changing the wind azimuth, results in similar uncertainties for the vorticity and Ekman pumping, further enhancing model errors.

The rest of the perturbed variables locally augment the ensemble spread in filament-like patterns in the periphery of eddies, near river plumes and over the shelf slope due to energy trapping (cf. Fig. 4 with respect to Fig. 3). After the wind, air temperature uncertainties have the next greatest impact on SST (Fig. 4b) and chlorophyll (Fig. 5b). Uncertainties in the wind drag coefficient have moderate impact on the wind stress and consequently on the spread of surface variables. Compared to the spread associated with uncertainties in the wind forcing they are approximately an order of magnitude smaller (Figs. 4d,j). This is because the expression for the wind drag coefficient 𝑐 has a correction dependency based on the different wind speed regimes (i.e. discrete 𝑐 values that the perturbation may be ineffective to change the class of these values). The ocean response to 𝑆𝐿𝑃 forcing has two components: the first-order Inverse Barometer (IB), which is isostatic and dominant at large scales, and the second-order non-IB, which depends mostly on the geographic region. The IB response to 𝑆𝐿𝑃 perturbations has spatial scales that are equal to or larger than the external Rossby radius and the IB pumping on the abyssal plain [START_REF] Tai | On the Quasigeostrophic oceanic response to atmospheric pressure forcing: the Inverted Barometer pumping[END_REF]. Subtracting the IB response from the SSH leaves us with the non-IB, which in our case has a limited impact on SST and SSH (Figs. 4c,i). Uncertainties in the bottom drag coefficient amplify error regimes, mostly for SST and less for SSH, along the shelf break and on the shelves, and especially on the macrotidal area of the English Channel dynamically controlled by strong tidal currents and fronts (Figs. 4e,k).

Biogeochemical uncertainties arise from the ecosystem model's intrinsic errors and from errors in the physical state variables (Fig. 5). All sensitivity experiments perturbing the physics leave an imprint on chlorophyll uncertainties, which are sometimes significant compared to the values of the control run and to other model variables, like for instance the SSH and SST. When perturbing only the biogeochemical model, we implement an identical stochastic pattern across all variables and vertical levels in order to increase the impact of the method. Since the physics are not perturbed, the biogeochemical uncertainties are passively advected via the ocean circulation. In all cases, physical model errors are found to have a larger impact on chlorophyll spread than those of the biogeochemical perturbations. It is straightforward to see why: if the ocean physics are identical for all members (as is the case for the SB ensemble) it behaves as a strong dynamical attractor for the system; vertical nutrient fluxes are of course sensitive to the velocity field and in particular its divergence, the latter scaling over the continental shelf break [START_REF] Karagiorgos | The Impact of Tides on the Bay of Biscay Dynamics[END_REF], a process which might dominate other sources of uncertainty, such as those of the 𝑆𝑀𝑆 of the biogeochemical model. In general, the spread is largest in SPB where both physics and biogeochemical perturbations are applied simultaneously. Biogeochemical uncertainties resulting from errors in the coupled system depend on the geographic location as well. An example is the chlorophyll uncertainties in the English Channel compared to the mesoscale field in the abyssal plain.

There is a stark geographic difference in spread between nanophytoplankton chlorophyll and diatom chlorophyll (Figs. 5i-n). This is expected from the different model parameters of the phytoplankton J o u r n a l P r e -p r o o f

classes and the additional requirements in nutrient supply of Si for the primary production of diatoms. Larger uncertainties are observed for the nanophytoplankton, especially in the open-ocean, whilst diatoms have larger uncertainties over the shelves and in the English Channel, where errors are dynamically controlled by tides. It is evident that the uncertainty regimes for each class follow the different size class chlorophyll abundance patterns in the region (cf. Fig. 5 with respect to Fig. 3), especially during the spring bloom.

Interestingly, even though model uncertainties are in general larger when both physical and biogeochemical models are perturbed, a small decrease in chlorophyll spread may be observed over the abyssal plain in the presence of coherent eddies. This effect is attributed to two facts: (i) the ensembles are not constrained by data assimilation and therefore model errors with long forecast lead time from different periods, may partially overlap or even cancel each other with "flow-dependent" errors (also called "errors of the day"); (ii) sub-mesoscale physical and biogeochemical processes with length scales smaller than 10 km can be resolved in the abyssal plain (Levy et al., 2012a;2012b;[START_REF] Charria | Interannual evolution of (sub)mesoscale dynamics in the Bay of Biscay[END_REF], therefore not contributing to the mesoscale error patterns being different in EnsB and EnsPB. This does not change the fact that the dominant spatial scales in the abyssal plain are characteristic of the underlying quasigeostrophic mesoscale features, due to the mesoscale vertical velocity field.

3.2

Quantitative evaluation of stochastic model skill

SST L4 gap-free gridded observations

Figure 6 shows examples of consistency metrics, as defined in Section 2.5.1, based on daily innovation samples. An ensemble of OSTIA SST is compared with EnsP (Figs. 6a,c). The EnsPB ensemble is statistically identical, since there is no feedback from the biogeochemical model into the circulation model. The calculations are carried out through the whole period of the seasonal-range ensemble. For the SP ensemble in Figs. 6b,d, calculations are performed for April 2012.

The model and data ensembles in both regions are consistent with each other since the supports of the pdfs overlap (Figs. 6a,c). Assimilating these observations, with these error estimates, with an ensemble or Bayesian filter would probably be well-posed, at least in data space (additional criteria to the pdf supports may need to be verified). In both regions, the model error estimate is lower than the innovation spread 𝐷𝑂𝑚𝐸, and the ensemble mean innovation vector 𝑂𝑚𝐸 is mostly contained within the observational error interval 0.5 °C (Figs. 6b,d), showing that the ensemble is consistent. The spread indicator slowly evolves in time, mainly increasing, consistently with the fact that the perturbation mechanism remains active throughout the period. However, there are occasions where the EnsP spread is reduced and the ensemble envelope is under-dispersive and sometimes biased with respect to the observational pdf. This occurs over a short period during the spring shoaling of the thermocline and mainly over the Armorican continental margin. This model overconfidence could be associated with missing re-stratification processes in our stochastic protocol, such as the vertical subgrid scale physics like the vertical eddy viscosity and diffusivity coefficients calculated in the turbulent closure scheme, which are not perturbed.

An interesting feature when comparing the seasonal-range EnsP and medium-range SP ensembles, is that the SST model errors over the abyssal plain increase with similar rates given a forecast lead time of a few days (Fig. 6d), despite the fact that they have different numbers of members and initial conditions (i.e. December 2011 for EnsP and April 2012 for SP). The stochastic restarts in these two experiments are identical and therefore the stochastic fields applied in the perturbed variables and parameters are alike. However, after this initial period of a few days' forecast lead time the model errors appear to increase with different rates, pertaining to the different ocean-atmosphere states during December for EnsP and April for SP. The latter is more apparent in the Armorican margins area (Fig. 6b).

In Figs. 6b,d, we compare the EnsP SST spread with the reference ensemble Ens0 [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. Both stochastic approaches in Ens0 and EnsP ensembles can be considered as variants of the perturbed tendency scheme. In the case of the AR(1) processes we have the option to perturb several variables in EnsP, compared to the stochastic modelling of only the wind forcing in Ens0. Therefore, the former yields a larger SST spread during winter and at the end of the run, with their differences ranging between about 0.05 °C to 0.1 °C, and with larger differences observed in the open-ocean. At the end of the run, both J o u r n a l P r e -p r o o f ensembles become biased with respect to the observations (more so over the shelves; cf. also Vervatis et al., 2016 for Ens0 and Figs. 6b, d for EnsP), since they are not constrained by data assimilation. During spring, in the presence of a strong thermocline, the spread is reduced in both ensembles and regions (though more over the shelf), suggesting that there are missing error processes in both stochastic protocols (i.e. in addition to the wind for Ens0 and all variables perturbed in EnsP).

SLA L3 along-track observations

The data space consistency analysis for EnsP SSH is carried out with the sea-level L3 along-track product in Table 2. Figure 7a shows the distributions corresponding to averages of all SLA along track observations crossing the BISCAY36 domain, as well as the EnsP and Ens0 supports of pdfs in data space. Both model ensembles appear to have more energy at the weekly timescale, occasionally at shorter timescales, and overall at the seasonal timescale compared to the observations. The minimum levels of spread in both model and data are reached in mid-February, 2012; for the rest of the simulation the SLA model equivalent shows larger variability than the data. There is a notable difference in spread between the EnsP and Ens0 ensembles, with the former's being greater possibly because of the SLA response to 𝑆𝐿𝑃 perturbations, as well as the different stochastic approaches in the wind perturbations. EnsP model uncertainties are comparable with those of the observational error at 0.05 m, reaching a few cm in magnitude.

We also find that in about one-third of cases the 𝑂𝑚𝐸 stays within the observational error interval (Fig. 7b). The most consistent configuration is over the abyssal plain, where the spread is increased due to mesoscale decorrelation of eddies after spin-up. In contrast, the largest statistical biases are observed in coastal regions and in the English Channel. This is also confirmed by the box-and-whisker plots that are plotted for several regions (Figs. 7c-e). These serve as a means to visualize both the distributions and their consistency. There are many cases where the joint uncertainty associated with both sources of information is clearly nonzero. However, several instances of strong bias are evident in two of the regions, sometimes to the point that the supports of the pdfs become disjoint, e.g. during the week between 06 to 12 March, 2012 (Figs. 7c-e). If those pdfs are used in an Ensemble filter, disjoint supports will lead to illposed assimilation problems, meaning that a solution will be obtained since all analysis schemes are convex, but the result will be meaningless.

When inconsistencies occur, it is generally because other error processes are active in the model in addition to the ones generated by the range of EnsP perturbations. This is more apparent in Ens0, since in this particular ensemble there are no perturbations in 𝑆𝐿𝑃. Both ensembles appear to be occasionally under-dispersive. In the English Channel, there are missing error processes such as the residual tidal error, which is enhanced by the presence of local tidal fronts and occasional Kelvin waves propagating along the coasts, and the non-isostatic response to atmospheric pressure (e.g. the non-IB response). The missing error sources are also particularly evident in the English Channel, hinting at the presence of highfrequency errors currently unaccounted for, although we cannot currently identify them with the tools at hand. Another possible reason for the observed statistical inconsistency is the low coverage of altimeter observations. This is exacerbated by the larger errors which altimetric data are subject to in coastal areas, which leads to less available data there. Future wide-swath altimetry products (e.g. SWOT, https://swot.jpl.nasa.gov/) are likely to provide better coverage in coastal regions and stronger constraints on models overall.

Ocean colour L4 gridded observations

Figure 8 shows results of biogeochemical ensembles and innovation statistics in the data space assuming an observational error of 0.3 mg/m 3 . Innovation statistics are calculated in log space applying an anamorphosis function to transform chlorophyll lognormal distributions into Gaussian distributions. The lognormal assumption for chlorophyll a has shown advantages in data assimilation systems [START_REF] Song | Data assimilation in a coupled physicalbiogeochemical model of the California current system using an incremental lognormal 4-dimensional variational approach: Part 3-Assimilation in a realistic context using satellite and in situ observations[END_REF]. Among the seasonal-range ensembles the EnsPB exhibits the largest chlorophyll spread, with EnsB being the least dispersive, as seen in the medium-range ensembles (Fig. 5). The model ensemble in the Armorican shelf shows a strong bias with respect to the ocean colour data, with disjoint pdfs supports that may lead to an ill-posed assimilation problem if those error estimates were used in a data assimilation system. On the other hand, in Fig. 8c, the EnsPB model-data samples appear to be marginally compatible J o u r n a l P r e -p r o o f with each other over the abyssal plain (also for EnsP and EnsB, not shown). There are also periods with disjoint pdfs supports with the most prominent in late March. During this period, observations show a strong phytoplankton bloom in the open-ocean, weakly present in all ensemble members (but with the correct phase), resulting in disjoint pdfs. Subsequent, less intense blooms in the open-ocean just before and after this event are also found in both model and data samples, and in those cases the pdfs overlap partially. Just as with the previously discussed SST and SLA networks, assimilating these ocean colour observations with these error estimates in the open-ocean would probably be marginally well-posed, in the sense that the joint probability associated with both sources of information appears to be nonzero only at times. The most consistent configuration for EnsPB after the spin-up period appears to be May-June and only for the abyssal plain.

In Figs. 8b,d, we illustrate the innovation metrics defined in Section 2.5.1. All biogeochemical model ensembles are under-dispersive during the first three months of the simulation, associated also with the fact that primary production is low in winter. After this period dispersion slowly increases over time. Towards the end of the run chlorophyll uncertainties in the Armorican shelf exceed 0.15 mg/m 3 (Fig. 8b). These model errors are small compared with the chosen observational error and with ocean colour data measured by satellite, often yielding a model-data misfit at about ~1 mg/m 3 . However, considering that chlorophyll has a lognormal distribution with near-zero positive values, they are not small with respect to the model's overall chlorophyll concentration of about 0.3 mg/m 3 . In the abyssal plain, the model-data ensemble mean misfit expressed by the 𝑂𝑚𝐸 metric, is well defined within the observational error interval for most of the period, except during the phytoplankton bloom observed in spring (Fig. 8d). In general, the control run underestimates chlorophyll abundance compared with ocean colour data and subsequently leads to under-dispersive biogeochemical model ensembles. Inconsistency is mostly attributed to the low bias. An examination of the temporal evolution of the spread resulting from different stochastic protocols, reveals that chlorophyll uncertainty variations are mainly controlled by physical processes and their errors. Biogeochemical processes and their uncertainties have only a moderate impact on model errors, one reason could be the unaccounted errors in nutrient initial conditions (cf. supplementary material using the world ocean atlas nutrient climatology; [START_REF] Garcia | World Ocean Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate and nitrate+nitrite, silicate). A. Mishonov Technical Ed[END_REF], except during periods of phytoplankton blooms where both components of the coupled system are important. In line with this, biogeochemical model errors are higher during periods of high biological productivity, than in periods of low biological productivity (e.g. comparison of EnsPB and SPB in Figs. 8b,d).

Examining the mean model spread for all ensembles against the mean spread of the innovation 𝐷𝑂𝑚𝐸, the former is always contained within the latter (Figs. 8b,d). However, this is not sufficient to verify consistency. This is because there is a constant bias of greater amplitude than the model dispersion, especially in coastal regions.

3.3

Qualitative assessment of stochastic model skill through its impact on multivariate incremental analysis

Ensemble-based single observation representers

In Figs. 9a-f, we illustrate examples of zero-lag representers for three individual observations in different locations on May 07, 2012. In order to address the qualitative criteria in Section 2.5.2, the representers are calculated from 40 members of EnsP (Figs. 9a-d), also from a 20-member subset of EnsP (Fig. 9e) and finally from 40 members of Ens0 (Fig. 9f), and are shown as correlations between SST and other surface variables of the control vector. The correlation structures reveal differences between the abyssal plain and coastal areas, as well as between variables. In the English Channel and over the south Armorican shelf, the filament-shaped structures for SST, SSS and surface chlorophyll are linked to near-shore features, such as river discharges (e.g. Loire river plume), mid-shelf thermal fronts and tidal fronts. The SSS pattern is dipolar, which could be explained by meridional plume displacement (Fig. 9b). Spring bloom is suggested by the negative correlation between SST and chlorophyll over the shelf (Fig. 9d). As the surface layers heat up during spring plankton depletion follows the bloom.

SST and SSH are largely decorrelated in the domain, probably because of the large-scale atmospheric forcings directly influencing SST in the spring, as well as low-frequency mesoscale variability (Fig. 9c). Due to the mixed conditions over the inner shelf at that time of year, the SSH response is relatively large-J o u r n a l P r e -p r o o f scale as it is associated with barotropic processes at the scale of the external Rossby radius. In the case of smaller ensembles, i.e. 20 members (Fig. 9e) or under-dispersive ensembles, i.e. SSH in Ens0 (Fig. 9f), correlations due to partially-converged statistics increasingly contaminate the pattern of representers. These patterns are however not very different from those of larger or more dispersive ensembles, at least in their general shapes and signs, and mostly the correlations or negative correlations are amplified. Similar results of rather broad and symmetrical structures are also found in the case of SSH single observation representers (not shown). Model errors associated with SST and SSH, at that time of the year and for those ensembles, may have a limited impact on the ensemble model skill.

Zero-lag representers of single-observation ocean colour data are calculated from EnsPB. The most important finding is that correlations behave differently in the open-ocean and over the shelves (Figs. 9gi). Chlorophyll autocorrelation structures appear broad and symmetrical (with respect to the singleobservation location) in the open-ocean, with scales dictated by the vertical velocity field of the underlying mesoscale features. This is in agreement with the biogeochemical model error regimes for the medium-range ensembles depicted in Fig. 5. Over the shelves, chlorophyll correlations are more dipolar in nature, when calculated for instance with respect to the SST and SSS fields, representing smaller-scale local conditions. Filament-shaped structures of negative correlations between chlorophyll and SST are seen again with the EnsPB stochastic protocol, indicating that model errors in primary production are mainly controlled by model uncertainties in physical processes, implying once again the possibility to improve ecosystem model skill performance by perturbing the physics alone.

Incremental analysis using SST L4

In this section, we qualitatively assess the data assimilation impact of the stochastic model in correcting the surface variables at analysis time. In Figs. 10a,b, we illustrate the SST correction of the first two members 001 and 002 on May 07, 2012, using the prior error covariance matrix 𝑷 (in a decomposed form of the square root matrix 𝑺 as in Eq.( 10)) from all 40 members of EnsP and assimilating the OSTIA SST dataset. SST corrections are consistent between members, suggesting that they are meaningful and that they could enhance model skill. The increments reveal a large-scale north-south pattern, which is consistent with the dynamics of the region and the model's cold bias with respect to the observations for this period over the Irish shelf and the English Channel. They also correct mesoscale processes over the abyssal plain and shelf processes near the river plumes. An interesting feature is the correction of the Bay of Biscay sub-gyre located between 4°-6°W and 44°-46°N, confirmed by consistent increments in SSH, SSS and especially in surface total chlorophyll, this last hinting at sub-gyre scale changes in the vertical velocity (Figs. 10c-e). Increments of opposite signs for SST and chlorophyll indicate that physical processes in the Bay of Biscay, such as tidal mixing, slope currents, river plumes and open-ocean mesoscale activity, play an important role in the biological productivity of the area. Nonzero chlorophyll corrections, calculated by cross-covariances of the coupled system, offer potential for improving ecosystem model skill, even if only physical data are assimilated.

In order to qualitatively evaluate the impact of EnsP on the 3D temperature, salinity and chlorophyll model update, we examine the increment profiles for the first member 001 at two specific locations over the abyssal plain and the Armorican shelf (Fig. 11). The increments are able to capture the region's spatiotemporal variability over two different seasons and in areas with different characteristics, suggesting that model skill performance depends on flow-dependent errors in the physics and biogeochemistry. The correction profiles reflect a fairly deep mixed layer in winter and a shallower mixed layer in spring, with chlorophyll changes at sub-surface layers hinting at ongoing bloom-related changes. The increments are associated with the shallow thermocline depth during spring (~10-30 m) and with the depth of the euphotic layer (~40-50 m) influencing the sub-surface vertical corrections of chlorophyll. Over the abyssal plain, at depths greater than 1000 m the vertical T/S corrections are possibly linked to the upper-ocean low-frequency mesoscale circulation, which affects the deep vortex dynamics. This is an intriguing result showing that deep model errors can be controlled by ensemble-based data assimilation methods, in which model ensembles are generated by perturbing surface variables and assimilating data of near-surface ocean properties.

Incremental analysis using ocean colour L4

J o u r n a l P r e -p r o o f

In Fig. 12, we see that most of the total chlorophyll correction on May 07, 2012, as seen from the increments of the first member, arises from uncertainties in the physics (i.e. EnsP and EnsPB), though biogeochemical model uncertainties in most areas in the domain are still significant. We focus on that specific date, because of the good spatial data coverage (ensembles are projected in data space in Figs. 12d-f), in order the analysis to be meaningful. When the physics are not perturbed (i.e. EnsB) one can see coherent mesoscale features in the correction (e.g. over the abyssal plain), since the physics and hence the dynamics are the same for all members. At 15 m depth, a signature of the ocean bottom relief can be seen in the correction fields, perpendicularly to the seabed features (resulting in a parallel crest signal) over the Celtic shelf near 7°-8°W and 48°-49°N. Biogeochemical processes and consequently ecosystem model performance are sensitive to bottom Ekman flow transport due to nutrient upwelling, especially in shallow shelf regions, and possibly to internal tides and waves induced by barotropic tides and winds.

High positive increments over the shelves for all three ensembles suggest that chlorophyll abundance is underestimated in those areas during subsequent spring blooms (Figs. 12a-c). This holds true over the continental shelf break near 3°-4°W and 46°-47°N, an area dynamically controlled by both barotropic and baroclinic tidal processes, which contribute to vertical mixing and enhance primary production. It is worth noting that in some coastal areas where the coupled system appears to underestimate chlorophyll abundance, a very small correction is applied instead, due to model overconfidence (i.e. under-dispersive ensembles lacking model skill) with respect to observational errors. In order to qualitatively assess the ensemble model skill, we may compare the chlorophyll increments contributed by different observational networks. Figure 10e and Fig. 12a present different correction patterns when different observations are assimilated, hinting at dissimilar processes captured by the prior error covariance matrices in combination with each network, in these cases by SST or ocean colour respectively. Here the chlorophyll correction is rather small when ocean colour is not assimilated, possibly because of weak cross-covariances between ocean physics and biogeochemical properties. Interestingly, we observe different signs in the increments, e.g. over the Celtic shelf, as calculated from the different model ensembles of the coupled system (Figs. 12a-c), and in particular between EnsP perturbing only physics and EnsB perturbing only biogeochemistry sources and sinks. In the present configuration, we conclude that the assimilation of ocean colour is important, and the specific stochastic method less so, when we seek to enhance model performance in marine ecosystems.

Incremental analysis using both SST L4 and ocean colour L4

In this section, we qualitatively assess the multivariate impact of the stochastic model in both temperature and total chlorophyll, from two different observational networks assimilated simultaneously, namely the OSTIA SST L4 and the ocean colour L4 (cf. Table 2). The changes that the extra observation network brings to the analyses vary from moderate to locally large for all variables. In Figs. 13a-d, we depict the correction fields for May 07, 2012 of the state vector surface variables for the first member, based on ensemble covariances of EnsPB with 40 members.

The north-south SST correction pattern is not as distinct as that when assimilating only OSTIA SST, especially in the English Channel area (Fig. 13a vs. Fig. 10a). In addition, the increments are greater, especially over the Irish shelf. This last aspect is also true for other surface variables, such as the SSH and SSS . We attribute this effect to the presence of ocean colour data in a situation where error pdfs might be disjoint, even when the data assimilation problem is convex, which suggests weak ensemble model skill. In these situations, the analysis scheme possibly impacts biogeochemical properties more than the ocean physics do. Similarly, chlorophyll correction values are moderately decreased when both observational networks are assimilated, compared to the analysis when assimilating only ocean colour (Fig. 13d vs. Fig. 12c). Compared to the correction of physical properties, the chlorophyll increments seem to be less sensitive to changes when multivariate observing networks are brought together. Thus, the ensemble ecosystem model skill appears to be less sensitive compared to the ocean model skill when ocean colour is assimilated, along with other more traditional networks like SST.

We investigated the convergence of covariances and its impact on the incremental analysis to assess model performance, incorporating different ensembles and ensemble sizes (i.e. how effective are the J o u r n a l P r e -p r o o f ensembles EnsP vs. EnsPB and 10 vs. 40 members at extracting reliable covariances; Figs. 13a,d,). Differences in the SST analyses between EnsPB and EnsP with 40 members in a multivariate context are minor (Fig. 13a vs. Fig. 13e). If we use fewer ensemble members, the analyses for both ocean physics and biogeochemistry properties resemble the correction patterns of larger ensembles. However, the increments are noticeably less smooth because covariances are calculated from partially converged statistics . As a final remark, we note that when chlorophyll is assimilated in conjunction with OSTIA SST, the analysis scheme moderately affects small scale processes for all variables, augmenting the increment values around coherent ocean dynamical features; this is suggestive of good stochastic model skill.

4

Discussion and conclusions

In this paper, Part 1 of a two-part series, we evaluated the skill of a NEMO-based physicalbiogeochemical stochastic model. To build this stochastic model, existing NEMO stochastic modules were complemented by integrating a 2D Gaussian function to introduce spatial correlations in the stochastic perturbations. [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF] generated correlations by applying a spatial filter in global and regional coarse domains. Here, the correlations have variable and anisotropic spatial scales, which is of particular importance in high-resolution coastal configurations. Another option would have been to use a tensor form that would make those patterns coast-aware [START_REF] Barth | Dynamically constrained ensemble perturbations-application to tides on the West Florida Shelf[END_REF]. We decided not to follow that route here, since this would have added significant complexity while, we felt, not being critical for the specific BISCAY36 configuration (e.g. there are not many islands and islets) and for the variables we chose to perturb (e.g. atmospheric forcing). Recursive filters and Gaussian functions have been used in the literature to calculate correlation lengths, but most of them focused on large spatial scales in coarse global domains and few of them have been used recently for biogeochemistry [START_REF] Storto | Estimation and impact of nonuniform horizontal correlation length scales for global ocean physical analyses[END_REF][START_REF] Kuhn | Temporal and spatial scales of correlation in marine phytoplankton communities[END_REF].

Our stochastic implementation was based on the SPPT and SPP schemes, and first-order autoregressive processes [START_REF] Brankart | A generic approach to explicit simulation of uncertainty in the NEMO ocean model[END_REF][START_REF] Ollinaho | Towards process-level representation of model uncertainties: stochastically perturbed parametrizations in the ECMWF ensemble[END_REF], applied to several sources of model uncertainties in the coupled system. These uncertainties emerge mainly from atmospheric forcing uncertainties, model errors in physical parameterizations and biogeochemical model sources and sink term uncertainties. Wind uncertainties were found to dominate all other atmosphere-ocean sources of model errors. Biogeochemical model uncertainties resulting from perturbations in physics were greater than those perturbing the 𝑆𝑀𝑆 concentration of the biogeochemical tracers.

Sensitivity studies of the stochastic model, e.g. with respect to (1) initial and boundary conditions (biases in this stochastic context), ( 2) parameter lognormal perturbations [START_REF] Ollinaho | Towards process-level representation of model uncertainties: stochastically perturbed parametrizations in the ECMWF ensemble[END_REF], and (3) perturbations of the unresolved scales, would be an ulterior step and are not included in this study.

The quantitative assessment of ensembles with respect to the gridded gap-free OSTIA SST L4, assuming an observational Gaussian error of 0.5 °C standard deviation, suggests that the seasonal-range ensembles where only the physics were perturbed are fairly consistent with the data distribution. Analysing the consistency of the ensembles with respect to the along-track SLA L3 CMEMS product (observational Gaussian error of 0.05 m standard deviation), we could see the presence of strong biases between the model and along-track data distributions. All ensembles were under-dispersive for sea-level, especially in coastal regions. The sea-level model-data misfits were found to be associated with strong SSH spatial gradients, in particular in the shelf regions such as the English Channel and Celtic Sea. High-frequency error processes are currently unaccounted for in the English Channel (e.g. open boundary conditions in the English Channel, which do not deal properly with high-frequency processes). Probabilistic attribution approaches (see e.g. [START_REF] Hannart | DADA: data assimilation for the detection and attribution of weather and climate-related events[END_REF], in addition to the quantitative assessment presented in this study, may provide more insight.

Our consistency analysis in log-space showed a statistical spin-up time for chlorophyll on the order of 3 months. During this spin-up period the biogeochemical model ensembles EnsP, EnsB and EnsPB were under-dispersive and biased with respect to gridded ocean colour L4 data to which a log-transformed, observational Gaussian error of 0.3 mg/m 3 standard deviation was applied. Statistical consistency was not always maintained for chlorophyll as it was for SST and, to a lesser extent, for SSH. It is difficult to attribute those error patterns to specific physical or biogeochemical processes, without a further J o u r n a l P r e -p r o o f probabilistic attribution analysis. The investigation of uncertainty regimes for other state variables of PISCESv2 is beyond the scope of this study, however we mention in passing that the observed error patterns in upper-ocean ecosystem properties, for instance for zooplankton and nutrients, are similar to those of chlorophyll (not shown).

The low chlorophyll bias may be a sign that the biogeochemical model's initial conditions for our ensembles (after the 5-month control run spin-up) are incorrect (cf. supplementary material using the world ocean atlas nutrient climatology; [START_REF] Garcia | World Ocean Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate and nitrate+nitrite, silicate). A. Mishonov Technical Ed[END_REF]. Perturbing the physical forcing and the biogeochemical 𝑆𝑀𝑆 terms cannot properly account for uncertainties in the initial biogeochemical state. These initial uncertainties could lead to unrealistic adjustments in the properties that are being perturbed.

Finally, we used the stochastic model output to generate multivariate increments, by means of representers and EnKF-type analyses with the aim to provide a qualitative assessment step towards developing an assimilation scheme. One objective was to assess the impact of observations on unobserved variables, such as other data types or subsurface variables. Specific qualitative criteria were used to assess the stochastic model skill at analysis time.

Based on these qualitative criteria, we summarize here the most important findings in the incremental analysis, arranging them into three broad categories. (i) Concerning the physics alone: corrections to physical properties are associated with large-scale biases between the open-ocean and the shelves; small scale local corrections are mainly visible over the shelves in near-shore coastal areas, which may be explained by meridional river plume migration, mid-shelf thermal fronts, barotropic tides and possibly baroclinic tides. (ii) Concerning the vertical profiles: incremental analysis on the water column structure reveals vertical changes linked to seasonal variability of the thermocline, such as for example the extended mixed layer depth in winter and the spring shoaling of the thermocline; increments are also observed at depths near the euphotic layer (~40-50 m) controlling the sub-surface vertical corrections of chlorophyll.

(iii) Concerning the interplay between the physics and the biogeochemical model: the scales of the correction patterns in the abyssal plain are often characteristic of the underlying quasigeostrophic mesoscale features associated with vortices; most of the chlorophyll correction arises from uncertainties in the physics, but biogeochemical model errors tend to enhance this effect; as expected, assimilating chlorophyll has a measurable impact on physical variables, e.g. when assimilating both ocean colour and SST datasets, the increments are larger than when assimilating SST only.

Results from the two assessment approaches in Sections 3.2 and 3.3 are analysed together in an attempt to evaluate the model skill in physics and biogeochemistry. The skill of the perturbation method at properly estimating model errors is in general improved for physics, though data assimilation performance to correct these errors depends strongly on the multivariate analysis. The biogeochemical model is not as sensitive as the physical model to the multivariate analysis and its skill is mainly defined by the assimilation of chlorophyll. The biogeochemical model errors appear to be complementary to those of the dynamical model alone and they contribute moderately to the incremental analysis. However, both increments derived from stochastic physics and biogeochemistry have marked differences in their spatial patterns. In some specific areas, the stochastic perturbations in physics and biogeochemistry lead to an under-dispersive spread. This is a weakness of the stochastic model which results in a limited impact of data assimilation in these areas, especially on the biogeochemical variables.

As our ensembles were sometimes found to be under-dispersive, for future efforts additional approaches should be envisaged to augment ensemble spread, in particular for biogeochemical variables; to list just a few: (a) inflation techniques [START_REF] Hamill | Distance-Dependent Filtering of Background Error Covariance Estimates in an Ensemble Kalman Filter[END_REF][START_REF] Anderson | Spatially and temporally varying adaptive covariance inflation for ensemble filters[END_REF] in the initial conditions may open some degrees of freedom in the first few time-steps of the model run, (b) perturbing the biogeochemical parameters (see also [START_REF] Garnier | Stochastic parameterizations of biogeochemical uncertainties in a 1/4° NEMO/PISCES model for probabilistic comparisons with ocean color data[END_REF] in addition to 𝑆𝑀𝑆 concentrations, (c) activation of the feedback of biology onto the physics in the NEMO-PISCESv2 coupled system, and (d) incorporating atmospheric ensembles, such as the ECMWF Ensemble Prediction System (https://apps.ecmwf.int/archivecatalogue/).

Uncertainties in the open boundary conditions, either in the numerical scheme or in the prescribed values for active boundaries, are expected to contribute significantly to the model error budget (see also [START_REF] Kim | Toward an uncertainty budget for a coastal ocean model[END_REF]. So, it is natural to consider perturbing the open boundary conditions. This would very likely J o u r n a l P r e -p r o o f lead to an increase in the ensemble spread from the surface to at least 1500 m. However, perturbing the open boundaries is not straightforward: first because of the need to ensure physical consistency between the perturbed variables, and second because the errors in the prescribed fields at the open boundaries are usually unknown. A favorable solution is when an ensemble of nesting (or parent) solutions is available and provides an estimate of open boundary uncertainties to the child model [START_REF] Ghantous | Ensemble downscaling of a regional ocean model[END_REF]. However, in the absence of a larger-scale ensemble to force the boundaries, a boundary error would appear to the stochastic model as a systematic error (bias) across all members, which was beyond the scope of this paper, but will be important to address in the future.

The success of physical-biogeochemical data assimilation will also depend on different observational networks brought together, e.g. remote sensing in synergy with in-situ observations, and their multivariate nature [START_REF] Song | Data assimilation in a coupled physicalbiogeochemical model of the California current system using an incremental lognormal 4-dimensional variational approach: Part 3-Assimilation in a realistic context using satellite and in situ observations[END_REF][START_REF] Verdy | A data assimilating model for estimating Southern Ocean biogeochemistry[END_REF][START_REF] Mattern | Data assimilation of physical and chlorophyll a observations in the California Current System using two biogeochemical models[END_REF][START_REF] Mattern | Improving Variational Data Assimilation through Background and Observation Error Adjustments[END_REF][START_REF] Kaufman | Assimilating biooptical glider data during a phytoplankton bloom in the southern Ross Sea[END_REF][START_REF] Ford | Assessing the role and consistency of satellite observation products in global physicalbiogeochemical ocean reanalysis[END_REF][START_REF] Goodliff | Temperature assimilation into a coastal ocean-biogeochemical model: assessment of weakly and strongly coupled data assimilation[END_REF]. Regarding in-situ observations, the increasing availability of high-resolution profiles provided by Biogeochemical-Argo floats (http://biogeochemical-argo.org), enables us to investigate the ocean interior across a wide range of spatiotemporal scales. The development of the Biogeochemical-Argo network is expected to substantially improve the quality of biogeochemical models, by incorporating vigorous validation, improved parameterizations [START_REF] Wang | Assessing the value of biogeochemical Argo profiles versus ocean color observations for biogeochemical model optimization in the Gulf of Mexico[END_REF], and novel data assimilation techniques, including ensemble-based methods [START_REF] Yu | Insights on multivariate updates of physical and biogeochemical ocean variables using an Ensemble Kalman Filter and an idealized model of upwelling[END_REF]. 1) using a 2D Gaussian function with length scales 𝜎 , ~1° and in Ens0 by using EOF modes [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. Colorbar units in m/s.

The model output consists of daily files of the ocean state vector and the two classes of chlorophyll, as well as three-day averages of 14 3D-biogeochemical variables.

We made used of NEMO's enhanced MPI strategy whose features allow for parallelization in both the spatial domain and across ensemble members [START_REF] Bessières | Development of a probabilistic ocean modelling system based on NEMO 3.5: Application at eddying resolution[END_REF]. Ensemble simulations were then carried out by just one call of the executable of the coupled NEMO-PISCESv2 system. BISCAY36 scales out using 96 processor cores of domain decomposition per ensemble member, excluding land processors. The configuration used the NEMO I/O and was connected to an external server (i.e. XIOS controlled by an XML file), thus increasing the total number of processors for the ensemble simulations, including those handling the I/O specifications (i.e. model variables, domains, grid, output frequencies etc.). Our ensemble experiments were designed to fit the scalability limits encountered in current operational systems. [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF] showed that with increasing ensemble size, results for BISCAY36 converge in the range of 20 to 40 members. We used the same range of members resulting in a scalability problem of the order of 10 cores.

Taking into account the ECMWF's hardware and software specifications, we tested the following resource geometries: (a) for 10 members, we used 960 NEMO processors and 48 XIOS servers filling a total of 28 nodes, and (b) for 20 members, we used 1920 NEMO processors and 24 XIOS servers filling a total of 54 nodes. The ensemble simulations were submitted as batch jobs for a 30-day run. For these examples, the ECMWF's job epilogue during production indicated a runtime average of about 489 minutes, with runtime standard deviation of approximately 29 minutes, including the first and last reading and writing time-steps.

Appendix C: Notation for ensemble-based innovation metrics

Let us note 𝑌 𝑦 , 𝑦 , the ensemble mean of the innovation for observation 𝑛 as:

𝑌 ∑ 𝑦 , 𝑦 , (C1)
where 𝑗 1, … , 𝑁 the size of the ensemble and 𝑛 1, … , 𝑝 the number of observations. With this notation, 𝑂𝑚𝐸 in Eq. ( 7) writes:

𝑂𝑚𝐸 ∑ 𝑌 (C2)
as the sum over 𝑝 observations of the ensemble mean of the innovation. From Eq. ( 8) the square of 𝑂𝑚𝐸 writes:

𝑂𝑚𝐸 ∑ ∑ 𝑦 , 𝑦 , (C3)
We now take the square of 𝐷𝑂𝑚𝐸 from Eq. ( 9) using the 𝑌 notation in (C1) as:

𝐷𝑂𝑚𝐸 ∑ ∑ 𝑦 , 𝑦 , 𝑌 (C4)
This can be written:

𝐷𝑂𝑚𝐸 ∑ ∑ 𝑦 , 𝑦 , ∑ 𝑌 (C5)
and using Eq. (C3) we get:

𝐷𝑂𝑚𝐸 𝑂𝑚𝐸 ∑ 𝑌 (C6)
If 𝑂𝑚𝐸 0 then from Eq. (C2) we get ∑ 𝑌 0. In general, this yields ∑ 𝑌 0 and subsequently 𝐷𝑂𝑚𝐸 𝑂𝑚𝐸 as seen from Eq. (C6), unless we have an unbiased system, i.e. no biases between model and observations with 𝑌 0, ∀𝑛. J o u r n a l P r e -p r o o f amplitude (1std: 𝜎 no units) of the 2D Gaussian distributions (prior to the anamorphosis function if applied). The spatial correlation value for the bottom drag coefficient 𝑐 is an approximation after 100 passes of the Laplacian operator. Ens0 is a seasonal-range ensemble performing stochastic modelling of the wind forcing based on EOF modes [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. J o u r n a l P r e -p r o o f Figure 5 (a-f) S1-5 and SP model uncertainties (i.e. 1std ensemble spread) of total surface chlorophyll (mg/m 3 ) perturbing physics (PHY), (g) SB chlorophyll spread perturbing 𝑆𝑀𝑆 𝐶 , (h) SPB chlorophyll spread perturbing all physical-biogeochemical variables. (i-k) and (l-n) are the same as (f-h) but for nanophytoplankton and diatom spread.

J o u r n a l P r e -p r o o f Figure 6 (a, c) OSTIA SST L4 (°C) observation distribution and EnsP ensemble inter-quantile ranges in data space averaged over the Armorican shelf and the abyssal plain, (b, d) innovation statistics and spread for the two regions, respectively; the thin horizontal line denotes the observational error 0.5 °C (the legend is split across two subplots).

Figure S4

Hovmöller plot of rank histograms (same as in Figs. 5b-c) between EnsPB and PFT (a) nano not combined with pico, and (b) diatoms not combined with dino.

  Figure A1Response function 𝑅 (no units) and noise attenuation (%) of the Laplacian-Shapiro filter for different number of iterations 𝑚, as a function of the correlation length scale L (km) and for different model resolutions: (a) 𝑑𝑥~25 km, (b) 𝑑𝑥~2.5 km.

  Figure A2(a, b) EnsP and Ens0 wind ensemble spread, and (c, d) wind ensemble mean on May 07, 2012. Stochastic modelling in EnsP is by means of SPPT AR(1) using a 2D Gaussian function with length scales 𝜎 , ~1° and in Ens0 by using EOF modes[START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. Colorbar units in m/s.

Table 2 . 1

 21 abbreviations: flt-filter; ens-ensemble; mem-members; 𝑈 -𝑢, 𝑣 wind velocities; 𝑇 -air temperature; 𝑆𝐿𝑃-sea level pressure; 𝑐 , 𝑐 , 𝑐 -wind drag and turbulent coefficients; 𝑐 -bottom drag; 𝑆𝑀𝑆 𝐶 -sources minus sinks of biogeochemical tracers 𝐶. *100 passes of the Laplacian filter in a Gaussian distribution per model grid point (spatial scales estimated by the response function of the first-order Shapiro filter). **A lognormal anamorphosis function is applied to the 𝑆𝑀𝑆 𝐶 of the 24 PISCESv2 prognostic variables.J o u r n a l P r e -p r o o f

J

  Figure 4Model uncertainties (i.e. 1std ensemble spread) of medium-range ensembles S1-5 and SP on April 30, 2012: (a-f) SST spread (°C) and (g-l) SSH spread (m) perturbing the (a, g) wind, (b, h) air temperature, (c, i) sea level pressure, (d, j) wind drag and turbulent coefficients, (e, k) bottom drag coefficient, (f, l) all variables together. Note the different colorbars with units varying up to 0.2 °C in SST and up to 0.02 m in SSH.

J

  Figure 7 (a) SLA along track L3 observation distribution (m) and EnsP and Ens0 ensembles in dataspace, (b) 𝑂𝑚𝐸 map using EnsP for the period starting on February 25, 2012 and for three consecutive weeks, (c) EnsP box-whisker plots and observation error bars averaged over the abyssal plain, (d) the Armorican shelf and (e) the English Channel. J o u r n a l P r e -p r o o f

  Figure 11Incremental analysis using OSTIA SST L4 on February 01, 2012 (grey lines) and on May 07, 2012 (black lines), using EnsP 40 members; (a-c) vertical correction of member 001 on, respectively, T (°C), salinity and total chlorophyll (mg/m 3 ) over the Armorican shelf, and (d-f) the same for the abyssal plain (cf. Fig.10a).

  Figure S1Ocean colour L4 (8-day frequency) and PFT (daily) chlorophyll concentration (𝑚𝑔/𝑚 3 ) in the Bay of Biscay from December, 2011 to June, 2012.

  Figure S2Spatial distribution of (a) PFT total chlorophyll and (b) pico, (c) nano, (d) diatoms, (e) dino concentrations (𝑚𝑔/𝑚 3 ) on March 28, 2012.

  Figure S3Scatter plots of chlorophyll concentrations in 𝑚𝑔/𝑚 3 on March 28, 2012: (a) PFT (nano and pico) vs. only nano, (b) PFT (diatoms and dino) vs. only diatoms. r is the correlation coefficient, RMSE the root mean square error (𝑚𝑔/𝑚 3 ) and with red the 1:1 line.

Table 1 .

 1 Table of simulations: control run (CR), ensemble medium-range sensitivity experiments (S1-5, SP, SB, SPB) and seasonal-range ensembles (EnsP, EnsB, EnsPB). The table shows the stochastic parameterizations of first-order autoregressive processes AR(1), deduced from the statistical properties of the one-year CR. Spatial correlations are explicitly calculated by applying a 2D Gaussian function with variable anisotropic length scales (1std: 𝜎 , in degrees), temporal correlations (𝜏 in days) and uncertainty

Acknowledgments. This work was carried out as part of the Copernicus Marine Environment Monitoring Service (CMEMS) "Stochastic Coastal/Regional Uncertainty Modelling (SCRUM)" Service Evolution project. CMEMS is implemented by Mercator Ocean International in the framework of a delegation agreement with the European Union. Part of this research was also made possible through the IKY Scholarships Programme and co-financed by the European Union (European Social Fund-ESF) and Greek national funds through the action entitled "Reinforcement of Postdoctoral Researchers", in the framework of the Operational Programme "Human Resources Development Programme, Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) 2014-2020. The contribution of P. De Mey-Frémaux and N. Ayoub is supported by the Centre National de la Recherche Scientifique (CNRS). We acknowledge the use of the ECMWF's computing and archive facilities in this research. This work was also supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facility -ARIS -under project ID PA002007. We thank four anonymous reviewers for their constructive comments.

Appendix A: Spatial correlation length scales

The response function of the first-order Shapiro filter [START_REF] Shapiro | Smoothing, filtering, and boundary effects[END_REF] measures the skill of the filter to attenuate noise and is given by the equation 𝑅 𝑐𝑜𝑠 • 𝜋 • 𝑑𝑥 𝐿 ⁄ , where 𝑑𝑥 the mesh grid resolution and 𝐿 the correlation length scale produced after repeated iterations 𝑚. Figure A1 shows the performance of the Laplacian-Shapiro filter to attenuate noise for two different mesh grid resolutions. Figure A1a refers to a coarse-resolution domain at 1 4°⁄ (e.g. ORCA025, 𝑑𝑥~25 km) and Fig. A1b refers to a highresolution domain at 1 36°⁄ (e.g. BISCAY36, 𝑑𝑥~2.5 km). For the coarse-resolution domain the Shapiro filter is able to attenuate noise efficiently (e.g. in most cases more than 50%) and produce long-range correlations of a few hundreds of kilometres. For the high-resolution domain the Shapiro filter may introduce length scales of a few tens of kilometres, but the noise is not attenuated for the large scales (cf. also Fig. 2h).

Figure A2 shows the first and second order statistical moments of two different methods for performing stochastic modelling of the wind forcing. In this study, we introduced wind uncertainties by means of SPPT AR (1) processes and a 2D Gaussian function with length scales 𝜎 , ~1° (Fig. A2a). In a previous study, we used EOF modes to generate wind forcing uncertainties [START_REF] Vervatis | Data assimilative twin-experiment in a high-resolution Bay of Biscay configuration: 4D EnOI based on stochastic modelling of the wind forcing[END_REF]. The most important attribute of both methods is that they have similar ensemble means (and their ensemble means are also similar to the unperturbed wind; not shown), whilst there are marked differences in their ensemble spreads.

J o u r n a l P r e -p r o o f Authors' contribution statement using CRediT 1) Vassilios D. Vervatis: conceptualization, methodology, software, visualization, formal analysis, investigation, writing -original draft, writing -review and editing. 2) Pierre De Mey-Frémaux: conceptualization, methodology, software, formal analysis, writing -original draft, writing -review and editing. 3) Nadia Ayoub: conceptualization, formal analysis, writing -original draft, writingreview and editing. 4) John Karagiorgos: software, visualization, formal analysis. 5) Malek Ghantous: writing -original draft, writing -review and editing. 6) Marios Kailas: software, formal analysis. 7) Charles-Emmanuel Testut: conceptualization, methodology. 8) Sarantis Sofianos: conceptualization, resources.

The corresponding author, Vassilios Vervatis

J o u r n a l P r e -p r o o f

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

J o u r n a l P r e -p r o o f

Supplementary material

Fig. S1 presents the temporal variability of the PFT data in the Bay of Biscay, during the period December, 2011 to June, 2012. We also show the ocean colour L4 total chlorophyll product used in the Part 1 article [START_REF] Vervatis | Assessment of a regional physical-biogeochemical stochastic ocean model. Part 2: empirical consistency[END_REF] and we verify that concentrations are of the same order to those of the PFT total chlorophyll. According to this, we categorised the four satellite PFT (pico, nano, diatoms and dino) into the two broad size groups in PISCES (nano or diatoms), in a manner most representative, ensuring that the total biomass (chlorophyll) from the model and satellite data can be compared like-for-like.

PFT diatoms and nanophytoplankton contribute together approximately more than 80% in total chlorophyll, whereas picoplankton contributes at about 10% and dinoflagellates less than 10% (Fig. S1). PFT diatom chlorophyll concentration is an order of magnitude larger compared with dinoflagellates and nanophytoplankton is about three times larger than picoplankton chlorophyll concentration (Fig. S1). In Fig. S2, we show the spatial distribution of the four satellite PFT and the total chlorophyll during the peak of the spring bloom on March 28, 2012. We confirm that the satellite micro class (i.e. diatoms and dino) is driven primarily by diatoms, far more abundant in the satellite data, with the two functional types being highly correlated in spatial. Fig. S3 presents scatter plots of combined vs. non-combined PFT chlorophyll, verifying the close relationship between functional types in a size class-based approach.

We also present results from one-on-one comparisons between model and PFT data, as opposed to the size class-based categorization merging different functional types. Figure S4 shows Hovmöller plot of rank histograms between EnsPB and PFT, in the same way as Figs. 5b-c, with one main difference: in Figs. S4a-b we do not combine the nano functional type with pico, nor we combine diatoms with dino. Rank histogram results for the nano class are degraded when pico and nano PFT data are not combined together in late-winter and early-spring when a primary bloom occurs (Fig. 5b vs. Fig. S4a). The latter may suggest that PISCES nano can be representative of a broader phytoplankton community, accounting also for smaller size classes. Rank histogram results are almost identical throughout the whole period for the micro class, regardless of whether dino and diatoms are combined together or not (Fig. 5c vs. Fig. S4b). Model-data one-on-one array mode consistency results (not shown) are in practice indistinguishable by visual inspection with the results presented in Figs. 10 and 11, confirming the validity of the size class-based approach.