Building the capacity for forecasting marine biogeochemistry and ecosystems: recent advances and future developments

Building the capacity for monitoring and forecasting marine biogeochemistry and ecosystem dynamics is a scientific challenge of strategic importance in the context of rapid environmental change and growing public awareness of its potential impacts on marine ecosystems and resources. National Operational Oceanography centres have started to take up this challenge by integrating biogeochemistry in operational systems. Ongoing activities are illustrated in this paper by presenting examples of (pre-)operational biogeochemical systems active in Europe and North America for global to regional applications. First-order principles underlying biogeochemical modelling are briefly introduced along with the description of biogeochemical components implemented in these systems. Applications are illustrated with examples from the fields of hindcasting and monitoring ocean primary production, the assessment of the ocean carbon cycle and the management of living resources. Despite significant progress over the past 5 years in integrating biogeochemistry into (pre-)operational data-assimilation systems, a sustained research effort is still needed to assess these systems and their products with respect to their usefulness to the management of marine systems.


Introduction
Global warming, in synergy with ocean acidification, eutrophication, deoxygenation, as well as the ongoing exploitation of living marine resources drive major changes in marine biogeochemistry and ecosystems (Bopp et al. , 2013Orr et al. 2005;Keeling et al. 2010;Gehlen et al. 2011;Stock et al. 2011). The awareness of ongoing and future changes prompted a considerable community effort over the past decade, resulting in the development of increasingly complex biogeochemical ocean general circulation models (BOGCM) (Moore et al. 2002;Le Quéré et al. 2005;Aumont & Bopp, 2006;Follows et al. 2007;Vichi et al. 2013). Today's BOGCMs capture the major features of large-scale distributions of dissolved inorganic carbon, alkalinity, nutrients, primary and export production (Schneider et al. 2008;Lazzari et al. 2010;Séférian et al. 2012). These models gradually developed into lower trophic level ecosystem (LTE) models by explicitly representing the first levels of the marine food-web from phytoplankton to zooplankton. They are increasingly used in a variety of studies encompassing the monitoring and short-term forecasting of ecosystem health up to the assessment of potential impacts of climate change on higher trophic levels (Stock et al. 2011). Biogeochemical-LTE models are becoming essential tools for: (1) producing decadal reanalyses (Ford et al. 2012;Fontana et al. 2013) as well as projecting trends in ocean biogeochemistry and ecosystems in a changing global environment (Orr et al. 2005;Schneider et al. 2008;Bopp et al. 2013); (2) designing the observation system needed to monitor the current marine biogeochemical status (Lin et al. 2010;Xue et al. 2012); (3) quantifying anthropogenic impacts on biogeochemical cycles and marine ecosystems Orr et al. 2005); (4) assessing and forecasting changes on time-scales relevant to management purposes (Stumpf et al. 2009; While there is an increasing demand for the last item, the development of a predictive capability for ocean biogeochemistry and LTEs on time-scales ranging from weeks to seasons remained, until recently, impeded by the lack of appropriate modelling tools and observational products. Past decades saw the rapid expansion of global ocean observing capabilities building on remote sensing and in situ technologies [e.g. Argo floats, moored arrays such as Tropical Ocean-Global Atmosphere (TOGA)/Tropical Atmosphere Ocean (TAO) array, Prediction and Research Moored Array in the Tropical Atlantic (PIRATA) or Research Moored Array for African-Asian-Australian Monsoon Analysis and Prediction (RAMA)]. The improved data coverage, both in time and in space, allowed the successful implementation of ocean data-assimilation platforms in parallel with the development of operational oceanography applications. The development, implementation and demonstration of feasibility were at the centre of the 'Global Ocean Data Assimilation Experiment (GODAE)'. This important 10-year international effort proved the feasibility and usefulness of integrated observation-model systems for physical variables. Over the recent years, biogeochemical variables were progressively added to global observing systems (e.g. Bio-Argo, ocean colour sensors and autonomous nutrient analysers). It was paralleled by the first 'green' applications of operational oceanography aiming at the near-real-time assessment of the ocean biogeochemical state ).
Here, an overview of several (pre-)operational biogeochemical systems presently active in Europe and North America for global to regional and sub-regional applications is presented. Systems were selected such as to illustrate regional versus global scales of applications, complexity of biogeochemical models, various approaches to data assimilation, as well as status of implementation (e.g. hindcast simulation without biogeochemical data assimilation to fully coupled operational systems). Emphasis is put on developments over the past decade, and a status report highlighting new opportunities for applications is presented. A brief overview of physical-biogeochemical systems for operational applications is first discussed. The overview distinguishes between biogeochemical models and the forecasting systems in which these are embedded. Next, (pre-)operational applications for topics covering biogeochemistry to ecosystems are illustrated, and to conclude, major blocking points and prospects for future evolution are identified.

Physical-biogeochemical systems for operational applications Biogeochemical components
Biogeochemical components embedded in large-scale physical operational systems are simplified representations of the complex and poorly constrained network of biological interactions governing the cycling of matter between reservoirs of the Earth system. These models do not aim at a detailed representation of marine ecosystems, but rather focus on processes and stocks relevant to biogeochemical cycles up to the dynamics of the first trophic levels. Organisms are often grouped into 'Plankton Functional Types' (PFTs) according to their specific function in carbon and nutrient cycles (Le Quéré et al. 2005;Hood et al. 2006). Next to function, size is an important structuring element in these models. Size participates in modulating nutrient acquisition, as well as grazing pressure and hence both bottom-up and top-down control acting on LTEs. The recognition that fluxes within the marine foodweb are characterized by near-constant elemental ratios (C:N:P = 106:16:1) (Redfield et al. 1963) when averaged over space and time (e.g. seasonal cycle) allows models to use a single basic currency (e.g. carbon, nitrogen or phosphorus). Fluxes of the remaining elements are computed from fixed stoichiometric relationships. While this approach limits the number of state variables, it does not allow representation of the decoupling of elemental cycles observed during blooms and subsequent export events (Engel et al. 2002). Some models are thus based on plastic stoichiometry (Baretta et al. 1995;Vichi et al. 2013) which allows for flexible ratios between the major elements in response to environmental conditions and physiological requirements. The stoichiometric plasticity is generally implemented through a modified internal cell-quota model (Droop, 1973), which relates the specific growth rate and nutrient uptake to the intracellular nutrient to carbon ratio (Baretta-Bekker et al. 1997), by adding a Michaelis-Menten limitation term for the external nutrient availability. Flexible stoichiometry improves the representation of co-limited oligotrophic regimes in biogeochemical models.
The biogeochemical models (Table 1, supplementary material for a detailed presentation of individual models) considered here are of the NPZD type, where NPZD denotes nutrient, phytoplankton, zooplankton and detritus. While HadOCC, PISCES and BFM have been developed for open ocean applications, NORWECOM and GSBM are regional in scope. The European Regional Seas Ecosystem Model (ERSEM) (www.shelfseasmodelling.org) was originally developed and routinely applied in the shelf-seas context but has since been extended for the global ocean and basin scale applications. It is consequently applicable for shelf seas and open ocean. The contrast in complexity between these models (Table 1) reflects at first order the level of complexity required by the precise context of the application, e.g. biogeochemical cycling versus ecological application. The trade-off between complexity and computational cost is an additional important consideration. With the exception of HadOCC, all models distinguish between large and small plankton or, in other terms, between fast-growing small cells and slower-growing large cells. The former are favoured under low-nutrient conditions, while the latter tend to dominate under high-nutrient conditions, yielding high biomasses and efficiently driving carbon export. All models except NORWECOM include a formulation of the carbon cycle. In GSBM, alkalinity is not a prognostic tracer but is derived from the relationship between observed alkalinity and modelled salinity. The number of nutrients is another important between-model difference, in particular the inclusion or omission of iron (Fe). Fe limits primary production across large oceanic regions, i.e. the high-nutrient low-chlorophyll (HNLC) regions, such as the Southern Ocean, the Equatorial Pacific. The availability of Fe (along with phosphate) also modulates atmospheric N 2 fixation. Both NORWECOM and GSBM do not include Fe cycling, since Fe is not a limiting micro-nutrient in their target regions. BFM, as well as PISCES, consider Fe biogeochemistry along with corresponding sources of external input.
Despite all models being of the NPZD type and hence sharing similarities in their underlying conceptual framework, there is a progression in complexity from HadOCC to NORWECOM, respectively GSBM, and PISCES to BFM. With the widespread use of biogeochemical models and the rapid increase in state variables, the identification of optimal complexity has become an active area of research (Friedrichs et al. 2007;Ward et al. 2013;Kwiatkowski et al. 2014;Xiao & Friedrichs, 2014a. In the context of this paper, we are more specifically concerned with how model complexity plays out in model forecast skill. The models presented here are coupled to different physical models, and any coordinated model intercomparison exercise would highlight differences in skill resulting from differences in both the biogeochemical and physical components. It would fall short in identifying the link between the structure of the biogeochemical model and forecast skill of the system. Addressing the question of optimal complexity of biogeochemical models calls for an approach in which the complexity of the biogeochemical component is increased in a stepwise fashion within an unchanged physical environment (Kwiatkowski et al. 2014).
Published studies suggest that models with two phytoplankton and two zooplankton size classes would have the optimal complexity for representing biogeochemical dynamics across a variety of regimes together with a good predictive ability (Friedrichs et al. 2007;Ward et al. 2013;Xiao & Friedrichs, 2014a. Increasing the number of state variables might increase model skillin particular when compared with local 1D or regional data setsbut goes along with a decrease in portability and, possibly, reduced predictive skill. Finally, model complexity is a function not only of the number of PFTs, but also of the degree of sophistication of process formulation. Future studies are required to address the benefits of the  inclusion of multiple nutrients and chemical speciation, of flexible stoichiometry and hence internal cell quota, etc. to model forecasting. For operational applications, trade-offs between model complexity, computational costs and forecast skill have to be carefully evaluated, keeping in mind that models should be as simple as possible, but as complex as needed. Ultimately, the required level of complexity will depend on the scientific questions the models have been/ are being designed to answer. Operational systems are deployed over a variety of spatial scales from regional to basin-wide to global applications. Models designed to address coastal processes often require a higher model complexity because key processes that govern plankton dynamics in the shallow waters (such as tidal processes, resuspension, sediment dynamics, benthic dynamics, etc.) are not of critical importance in the open ocean. Complexity in the number of plankton variables in coastal waters is also justified in order to capture, for example, occurrences of harmful algal blooms that affect fisheries, aquaculture, human health and local economies, whereas this is not an issue in the open ocean. There is thus no unique answer to the question of optimal model structure, and the preceding rather pleads in favour of a diversity of biogeochemical models.

Forecasting systems
The major characteristics of the forecasting systems in which the biogeochemical components are embedded are summarized in Table 2.

FOAM-HadOCC
HadOCC has been coupled online to the global configuration of the Met Office's operational Forecasting Ocean Assimilation Model (FOAM) (Storkey et al. 2010;Blockley et al. 2014), which is based on the Nucleus for European Modelling of the Ocean (NEMO) hydrodynamic model (Madec, 2012) and the Los Alamos sea ice model (CICE) (Hunke & Lipscomb, 2010). A key feature of FOAM is the ability to assimilate remotely sensed and in situ observations of temperature, salinity, sea-level anomaly (SLA) and sea ice concentration. The data assimilation has recently been upgraded from the analysis correction scheme of Martin et al. (2007) to the 3D-Var NEMOVAR scheme .
The data-assimilation capability has been extended to also assimilate chlorophyll data derived from remotely sensed ocean colour and in situ observations of sea surface partial pressure of carbon dioxide (pCO 2 ). This has been used for both reanalysis and near-real-time forecasting applications, with the aim of improving estimates of air-sea CO 2 flux and other carbon-cycle variables, as well as biological variables. Both the chlorophyll (Ford et al. 2012) and pCO 2 ) assimilation schemes follow a two-stage process. The first stage is to generate a set of sea-surface chlorophyll or pCO 2 increments. This is currently done using the analysis correction scheme of Martin et al. (2007) but will be updated to use the same NEMOVAR scheme as the operational physics system. Both employ an observation operator that uses a first-guess-at-appropriate-time (FGAT) technique. The second stage is to generate a set of multivariate increments, which are then applied to the model. Surface chlorophyll increments are converted to a set of increments for all biogeochemical model state variables, at all depths, using the nitrogen-balancing scheme of Hemmings et al. (2008). This uses a principle of conservation of nitrogen and carbon to propagate the information from the chlorophyll observations as fully and realistically as possible to the rest of the model. The scheme is designed to be computationally efficient, making it suitable for operational applications. In the pCO 2 assimilation scheme, the surface pCO 2 increments are converted to dissolved inorganic carbon and alkalinity increments using the method described in While et al. (2012), and applied throughout the mixed layer.
The coupled FOAM-HadOCC system is currently preoperational, and different configurations are used for different experiments. The open ocean physical FOAM system is fully operational, and details of its running are given in Blockley et al. (2014). There is also an operational FOAM configuration covering the northwest European shelf seas O'Dea et al. 2012). This is coupled online to ERSEM, and both physical and biogeochemical analysis and forecast products are provided operationally via the MyOcean project (http://www. myocean.eu/).

TOPAZ-NORWECOM
TOPAZ-NORWECOM is used in the Arctic Forecasting Center in the MyOcean project. The modelling system is the result of a collaboration between the IMR and NERSC for the developments and MET Norway, which exploits the system operationally. The coupled model is HYCOM-NORWECOM (Pätsch et al. 2009;Samuelsen et al. 2009). It uses the HYbrid Coordinate Ocean Model [HYCOM (Bleck, 2002)] as the physical model and NOR-WECOM (Oki & Sud, 1998;Skogen & Søiland, 1998) as the ecosystem model. The physical and biological models are coupled online. HYCOM uses a combination of isopycnal and z-coordinates, which allows for both good conservational properties in the deep ocean and high vertical resolution in the upper mixed layer. The present model configuration has 28 vertical layers, of which the five upper layers are in z-coordinates, and the lower 23 layers are hybrid layers.
The river forcing is generated using a hydrological modeltotal runoff integrating pathways (TRIP) (Oki & Sud, 1998). The river outflow calculated by TRIP was combined with data from Global Nutrient Export from Watersheds [GlobalNEWS (Beusen et al., 2009;Seitzinger et al. 2005)] and used as nutrient river input to the model. Nutrients and oxygen are relaxed to climatology at the lateral boundaries. For initialization of NORWECOM, climatological values of nutrients and oxygen were used; all other variables were set to a constant low value. The capacity to assimilate ocean colour data has been developed using a version of the Ensemble Kalman Filter that includes both a Gaussian anamorphosis in order to handle the non-Gaussian model variables (Simon & Bertino, 2009) and a state augmentation in order to estimate biological model parameters together with the model variables (Simon & Bertino, 2012). These assimilation procedures have been applied with the HYCOM-NORWE-COM system coupled online at a coarser resolution than that of the real-time forecast system -50 km instead of 12 kmfor the period 2007-2010. An ensemble of 100 members has been integrated. The data are available as part of the MyOcean service as the first Arctic biological reanalysis that assimilates jointly physical (altimeter, SST, sea ice) and biological observations.

MERCATOR-OCEAN/BIOMER
Since 2012, the French centre for operational oceanography Mercator Océan has added near-real-time assessment of ocean biogeochemistry called BIOMER to its suite of products. BIOMER is based on the ecosystem model PISCES forced offline with ocean physical fields provided by the global operational system PSY3V3 (Lellouche et al. 2013) at a horizontal resolution of 1/4°cos(lat) (NEMO 3.1, 50 vertical levels, atmospheric forcings from ECMWF operational analysis at 3 h, CORE bulk formulation), with assimilation of temperature, salinity and sealevel data via the SAM2 assimilation scheme based on the SEEK method , an Incremental Analysis Update (Bloom et al. 1996;Ourmières et al. 2006) and a large-scale bias correction for temperature and salinity based on a 3D-var approach. To decrease the computational burden, the spatial resolution of physical fields was coarsened to 1°, and temporal output was aver-  Figure 1(b)] and nutrient distributions (not shown) in the equatorial band compared with a simulation without assimilation of physical data. The overestimation of chlorophyll and nutrients was attributed to spurious vertical fluxes in response to unbalanced physical data assimilation ) associated with biases in assimilated MSSH (mean sea surface height).
The recent version of the system, BIOMER4, uses the same spatial resolution for the physical and biogeochemical systems [1/4°cos(lat)]. Besides its higher spatial resolution, BIOMER4 [ Figure 1(a)] relies on a daily forcing of the offline biogeochemical model with ocean physical and atmospheric fields, and assimilation of MSSH is changed to hybrid MSSH. Modifications of the biogeochemical component include the transport of passive tracers (coefficient of diffusion) and the tuning of biogeochemical parameters (specifically the photosynthetic response to changes in light intensity and the sinking speed of large particles).  Figure 2, for North Atlantic chlorophyll concentrations. For a complete evaluation, please refer to MyOcean quality information document MYO2-GLO-QUID-001014 (http://www.myocean.eu/web/69-myocean-interactivecatalogue.php). Improvements between BIOMER1 and BIOMER4 result from the combination of changes made to the physical-biogeochemical system: (1) increasing the spatial resolution from 1°to 1/4°improved the representation of productive shelf regions and Eastern Boundary Upwelling Systems; (2) shifting from assimilation of MSSH to hybrid-MSSH reduced spurious vertical velocities; and (3) adjusting selected parameters of PISCES reduced its productivity.
Model output will be made available in near-real-time with a lag of one day for nitrate, phosphate, dissolved oxygen, iron, dissolved silica, phytoplankton biomass, net primary production (NPP) and chlorophyll on the MyOcean portal in September 2014. The near-real-time assessment of ocean biogeochemistry is completed by an interannual hindcast simulation from 2007 to the present.

MFS-OPATM-BFM
The OPATM-BFM model system (Lazzari et al. 2010) has been developed at National Institute for Oceanography and Experimental Geophysics (OGS) and is the biogeochemical component of the MyOcean Mediterranean Monitoring and Forecasting System (MedMFC). Since the start of the preoperational implementation in 2007 within the MERSEA project, substantial improvements have been made. In the framework of MyOcean, the system was upgraded by adding a 3D variational assimilation scheme (Teruzzi et al. 2014a) for satellite (MODIS)-derived chlorophyll and apparent optical properties (AOP). While satellitederived chlorophyll constrains the C:chlorophyll ratio of phytoplankton, AOP allow the estimation of the optical model. The transport model is based on OPA 8.1 (Foujols et al. 2000) including additional enhancements such as a flux-corrected advection scheme. The horizontal spatial resolution is 1/16°. The physical forcing for biogeochemical processes is pre-computed at the same resolution by the ocean general circulation model MFS based on NEMO3.4. MFS supplies the daily averaged fields of horizontal and vertical current velocities, vertical eddy diffusivity, potential temperature and salinity, in addition to surface data for solar shortwave irradiance and wind stress.
The Med-MFC operational chain runs unattended twice a week in an offline mode. It is composed of four steps: (1) downloading of forcing and daily MODIS chlorophyll data and pre-processing, (2) 3D-var satellite chlorophyll assimilation, (3) 7 day forecast of the biogeochemical state of the sea and (4) MyOcean catalogue update. The biogeochemical run is composed of (1) 7 days of hindcast (the first day: analysis with data assimilation on Wednesday run only) and (2) 10 days of forecast for the Wednesday run. Operational output includes daily 3D fields of chlorophyll, oxygen, nitrate, phosphate, aggregated phytoplankton biomass and NPP.

CANOPA-GSBM
CANOPA-GSBM is being developed in collaboration with three Fisheries and Ocean Canada Institutes. The ocean circulation model, CANOPA, is described in Brickman & Drozdowski (2012). The modelling system is based on the ocean code OPA version 9.0 (Madec, 2012) coupled to the sea-ice model LIM2 (Madec et al. 1998). The physical and biochemical (GSBM) models are coupled online. The grid of the model covers the Gulf of St. Lawrence, the Scotian Shelf and the Gulf of Maine with a spatial resolution of 1/12°and 46 layers of variable thickness in the vertical. The model includes tidal and river forcing. It is a prognostic model, meaning that the temperature and Lawrence river is estimated from water-level measurements, while runoff from the other rivers is obtained from a simple hydrological model (Lambert et al. 2013).
Monthly climatologies for temperature, salinity, nitrate, oxygen and dissolved inorganic carbon are used to initialize the model and to provide lateral boundary conditions. Monthly climatologies of nitrate are also used for river input. Surface forcing for CANOPA is obtained from the Canadian Meteorological Centre Global Environmental Multiscale (CMC-GEM) atmospheric model (Pellerin et al. 2003). CANOPA-GSBM is used to produce an annual state of the ecosystem for part of the Gulf of St. Lawrence. Within the Canadian National Approach under GODAE Ocean-View, CONCEPTS, EC and DFO plan over the next 4 years to couple GSBM to a subset of their Global Ice Ocean Prediction Systems (GIOPS) currently operational at Environment Canada as well as to a planned Regional Ocean Ice Prediction System (RIOPS) to be implemented in 2015 (Ryan et al. 2015;Tonani et al. 2015).

Synthesis of present status of coupled physicalbiogeochemical systems
The review of (pre-)operational physical-biogeochemical systems presented here highlights the diversity of biogeochemical models and data-assimilation methods used across the community (Tables 1 and 2). CANOPA-GSBM, in its present set-up, produces a yearly biogeochemical hindcast assessment at a regional scale. Two systems (MERCATOR-OCEAN/BIOMER and OPATM/ BFM) take advantage of physical data-assimilation products to drive ocean biogeochemistry offline. This approach might be seen as an intermediate step towards online coupled physical-biogeochemical systems as exemplified by the FOAM-HadOCC and TOPAZ-NORWECOM model suites. Since the review by Brasseur et al. (2009), several GODAE systems are in the process of moving, or have moved, to operational mode (Table 2). Satellite ocean colour data are still the only biogeochemical data set meeting the requirements of spatial and temporal coverage for assimilation into operational systems. However, other sources of observations can be useful depending on the aims of the experimentfor instance, a pCO 2 assimilation scheme is described in the following section of this paper.
A number of different techniques have been used to assimilate surface ocean chlorophyll, and these consistently improve modelled distributions compared with both the assimilated satellite data and independent in situ observations. However, very much as for the issue of model complexity discussed above, no unique recommendation emerges at this stage. The diversity of approaches guarantees the necessary degree of scientific diversity for the development of targeted applications and variable spatial scales. The increasing number of systems producing biogeochemical forecasts and reanalyses will enable coordinated evaluation and model intercomparison exercises in the near future.

Physical-biogeochemical model applications in operational oceanography
Marine primary production NPP is at the base of the marine food web. The reliable simulation of its spatial distribution and temporal variability is of importance for understanding and predicting ecosystem dynamics, as well as for the development of downstream applications such as the operational management of marine resources. The evaluation of modelled NPP is not straightforward. Estimates derived from remote sensing rely themselves on empirical models that derive NPP as a function of chlorophyll, temperature and light. In situ measurements are scarce and fail to resolve spatial and temporal variability. This explains why biogeochemical models are evaluated against chlorophyll and nutrient distributions. Three-dimensional nutrient distributions reflect the combination of physical transport and biological activity in surface waters and across the water column (Séférian et al. 2012). In addition, data from Eulerian time-series stations and cruises provide complementary information on plankton biomasses and biological rate estimates, including NPP. In order to illustrate the skill of (pre-)operational physical-biogeochemical model systems, we present hereafter two biogeochemical reanalyses: one at the global scale and one for the Mediterranean Basin.
FOAM-HadOCC has recently been used to produce a 1°global daily biogeochemical reanalysis covering the period September 1997 to July 2012, assimilating daily merged chlorophyll products from GlobColour (http:// www.globcolour.info). For this particular work, only chlorophyll data were assimilated, not physics or pCO 2 . The data assimilation has been shown to improve the fit of modelled chlorophyll not only to the assimilated satellite data but also to independent in situ chlorophyll observations, both at and beneath the surface. The changes to model chlorophyll result in corresponding changes to the primary production fields. Comparisons with in situ observations of nitrate show no degradation of the model nutrient fields as a result of the assimilation. These results are similar to those detailed in Ford et al. (2012), to which the reader is referred for a more detailed validation of the chlorophyllassimilation scheme. The reanalysis is being used to investigate seasonal and interannual variability in the model results, such as relationships between biogeochemical variables and known climate drivers such as the El Niño Southern Oscillation and North Atlantic Oscillation, and how the data assimilation impacts this. Work is also under way to assess reanalysis results using initial products from the European Space Agency's Climate Change Initiative project (http://www.esa-oceancolour-cci.org).
A 12-year (1999-2010) reanalysis of Mediterranean Sea biogeochemical tracer fields at 1/8°has been made available as part of the recent MyOcean release V3.1. It provides monthly 3D mean fields of chlorophyll, phosphate, nitrate, dissolved oxygen, NPP and phytoplankton biomass. A database of satellite derived (SeaWIFS) surface chlorophyll concentrations delivered by the MyOcean Ocean Colour Thematic Assembly Center was used for data assimilation. Data assimilation of surface chlorophyll concentration is performed weekly through a 3D variational scheme (3D-var) following Dobricic and Pinardi (2008) and modified as in Teruzzi et al. (2014a). To statistically project surface information to a fully 3D chlorophyll field, EOF analysis has been applied to a multi-annual (1998)(1999)(2000)(2001)(2002)(2003)(2004) OPATM-BFM run (Lazzari et al. 2012) identifying 12 months and nine sub-regions. The efficiency of the procedure was assessed by visual inspection (see Figure 3), as well as by comparing predicted surface chlorophyll fields to satellite-derived observations and in situ data (Teruzzi et al. 2014a).
For the purpose of model validation and assessment of the impact of monovariate chlorophyll assimilation on the overall model prediction skill, an extensive dataset was collated for dissolved oxygen, nitrate, phosphate, dissolved silica and NPP (Lazzari et al. 2012). It combines biogeochemical measurements available at the OGS National Oceanographic Data Centre with additional recent cruises. The metrics introduced are based on nearestpoint statistics. Model bias and correlation were calculated according to Allen et al. (2007) and summarized in Table 3 for chlorophyll, dissolved oxygen and nutrient data for the reanalysis and the corresponding free simulation (hindcast). The mostly lower bias at depth reflects the realistic initial conditions. Relative model output/data mismatch is amplified when expressed as relative error for low to very low nutrient levels. This explains in part the larger bias obtained for nitrate and phosphate for oligotropic surface waters compared with the ocean interior. Dissolved silica distributions are poorly represented. Basin-wide correlations between observed and modelled tracer distributions are good at depth (again reflecting initial conditions) and for surface phosphate and dissolved oxygen fields. It follows from the comparison between the two simulations that there is a substantial improvement in surface chlorophyll. The bias on nutrient estimates is not adversely affected by data assimilation with improvements for four out of nine fields, while reanalysis correlation outperforms hindcast in seven out of nine cases. For this reanalysis, chlorophyll was assimilated and projected on phytoplankton biomass only, that is, without correcting nutrient concentrations. While this approach improved surface chlorophyll distribution, 3D nutrient fields and oxygen distribution only partially benefit from univariate assimilation of surface chlorophyll. This example argues in favour of the development of multivariate approaches to the assimilation of biogeochemical data in coupled physical-biogeochemical models. s176 Figure 3. Comparison of surface chlorophyll as (1) observed by SeaWiFS (first column), (2) predicted by a hindcast (control run) (3) as in (2) but with surface chlorophyll-a data assimilation (Assimilation Run). Three model daily means (3, 8 and 13 March 2008) exhibit better predicting capabilities of the DA run than with the control run. The SeawiFS chlorophyll-a field is assimilated to produce the prediction of the next time (blue arrows), (modified after . Table 3. Uncertainty of Mediterranean biogeochemical reanalysis. Basin-wide bias (percentage) and correlation have been produced for the upper ocean (0-200 m, italics) and for the ocean interior (>200 m). The statistics are computed against in-situ profiles for dissolved oxygen, nitrate, phosphate, dissolved silica while surface Chl-a is compared with satellite observations. The best performances are indicated in bold [modified after Teruzzi et al (2014a)].

Ocean carbon cycle
Accurate simulation of the marine carbon cycle, particularly air-sea CO 2 fluxes, is desirable for ocean carbon cycle monitoring and reanalysis purposes. As such, the aim of producing a reanalysis with FOAM-HadOCC based on the assimilation of ocean colour data is not just to simulate biological variables such as chlorophyll and primary production, but to assess the impact that assimilating satellite chlorophyll data have on carbon-cycle variables. When comparison is made with in situ observations of sea surface pCO 2 in the open ocean, the overall impact of the chlorophyll data assimilation on the model carbon cycle is small, with similar error statistics seen for both a free run and an assimilation run. This is due in part to pCO 2 being largely physically controlled. However, in regions where there is a strong biological component to the carbon cycle, for instance during the North Atlantic spring bloom, the chlorophyll assimilation can have a beneficial impact. An example is shown in Figure 4, which compares the June mean North Atlantic air-sea CO 2 flux from equivalent free and chlorophyll assimilation runs of FOAM-HadOCC with the climatology of Takahashi et al. (2009) for June. In the free run, the airsea CO 2 flux is too weak across the northern part of the basin, with an area of spurious outgassing seen in the central North Atlantic. This error is much reduced in the chlorophyll assimilation run, which also shows increased CO 2 uptake in the northeast Atlantic, better matching the climatology. This demonstrates the potential of ocean colour assimilation to improve estimates of non-assimilated variables.
Next, a scheme to assimilate in situ pCO 2 observations has been developed by While et al. (2012) and implemented in FOAM-HadOCC. The assimilation was shown to reduce model pCO 2 errors, which has a direct impact on air-sea CO 2 fluxes. Impacts were also seen in the dissolved inorganic carbon and alkalinity fields. Where systematic pCO 2 errors were not reduced, it was concluded that this was due to sparse observational coverage in those regions. However, the assimilation had a long memory, with information from the increments retained by the model for several months. This demonstrates that even a limited number of observations can have a substantial and lasting impact on the model results. In future, this method can be applied to larger pCO 2 observation databases such as the Surface Ocean Carbon Atlas (SOCAT) (Bakker et al. 2014) in order to produce long-term reanalyses. A development could be to assimilate physical, chlorophyll and pCO 2 data simultaneously, with cross-covariances or balance relationships between the physical and biogeochemical variables, thereby making maximum use of available observations.
In order to achieve successful pCO 2 forecasts, the estimation of poorly known static and globally constant biological model parameters is a necessary way ahead. However, affordable parameter estimation methods still need to be proven for non-linear models with positivevalued variables. An attempt was therefore made at NERSC to estimate two parameters of the MICOM-HAMOCC model (Assmann et al. 2010) using a local Ensemble Kalman Filter with state augmentation and Gaussian anamorphosis, similarly to the aforementioned HYCOM-NORWECOM system. The two model parameter values (namely phytoplankton growth rate and air-sea gas transfer) were also estimated locally. The surface pCO 2 observations from SOCAT were assimilated for 4 years 1999-2003 and did improve marginally, but systematically, the forecasts of pCO 2 (Simon & Bertino, 2013) (Figure 5). The proof of concept was only a partial success: a perfect twin experiment showed that the air-sea gas transfer coefficient could be retrieved everywhere, but the estimates of phytoplankton growth rate diverged in extra-tropical areas such as the North Pacific, where the nonlinearities seemed more severe than elsewhere, and denser observations would be necessary. For real-time forecasting of pCO 2 , the ultimate system should combine assimilation of physical data, chlorophyll data and surface pCO 2 as well as online estimation of a few model parameters, for which the scientific basis is being progressively established.  M. Gehlen et al.

Operational management of fisheries
As illustrated throughout this paper, coupled physical-biogeochemical model systems are becoming increasingly skilful and starting to provide useful information to realworld applications, such as the operational management of fisheries. Preliminary steps towards this goal have been initiated, for example, in the case of the high-seas pelagic tuna fisheries. Operational physical model data are used to drive the spatial ecosystem and populations dynamics model SEAPODYM, to simulate functional groups of organisms at the intermediate trophic levels (Lehodey et al. 2010) and the dynamics of their predators, i.e. tuna . The model also includes a quantitative estimation approach ) allowing the use of this model for fish stock assessment and testing management scenarios (Sibert et al. 2012). The INDESO (Infrastructure Development of Space Oceanography) project, which includes the development of an operational ecosystem and tuna model (SEAPO-DYM) in a regional Indonesian configuration, provides a proof of concept. A global operational configuration driven by Mercator-Ocean PSY3V3 model outputs (temperature and currents), and global primary production derived from satellite ocean colour data (Behrenfeld & Falkowski, 1997) provide the boundary conditions for the regional configuration. The first version of the model chain uses dynamical output from a regional configuration of the ocean general circulation model NEMO/OPA at 1/ 12°resolution and satellite derived NPP. In the near future, the regional coupled physical-biogeochemical model based on NEMO/PISCES will provide physical and biogeochemical forcings for SEAPODYM (e.g. temperature, salinity, euphotic depth, oxygen and NPP). Environmental constraints on tuna population dynamics are complemented by best estimates of fishing mortality. The chain of production runs every week and delivers one week of hindcast, one week of nowcast and 10 days of forecast. These outputs will be used by the Fishing Authority to implement a better collection and verification of fishing data. The real-time and high-resolution configuration should help to quickly improve the model calibration and consequently abundance estimates needed to establish the optimal level of exploitation (total allowed catch) and the conservation measures (e.g. identification and protection of spawning grounds and nurseries) required for the sustainable exploitation of this essential resource.
Another example is given by Christensen et al. (2013), who, within the first MyOcean project, developed a shortterm (seasonal up to a couple of years) forecasting system for sand-eel population based on larval dispersal and population analysis covering several life cycles, in order to support fisheries management providing information on the maximum sustainable yield and the effects of alternative fishing scenarios.

Future challenges, perspectives and blocking points Biogeochemical data assimilation
Today, ocean colour observations and derived products (e.g. surface estimates of chlorophyll concentration, NPP, particulate organic carbon) are the only source of data collected routinely at a global scale that provide a quantitative description of the surface signature of marine primary production at relevant space and time-scales. However, when the models reviewed in the previous section are used as deterministic tools to hindcast/forecast the biogeochemical properties of the ocean without data to constrain the simulations, they quite often display systematic biases as compared with ocean colour patterns observed from satellites. Misfits between ocean colour observational products and simulated field properties may have many different origins, e.g. inaccurate forcing fields (heat/momentum surface fluxes, photosynthetically active radiation fluxes, precipitations) used to drive the simulations, unrealistic physical model response in terms of advection/diffusion by ocean currents or unresolved spatial scales, oversimplified complexity of biogeochemical model structure and process formulation (e.g. no day-night cycle for photosynthesis) and uncertain parameterizations of interactions between biological species. Satellite-derived chlorophyll fields in turn suffer from a low signal-to-noise ratio, sensitivity to cloud cover and aerosols, as well as water transparency. They depend on (regional) empirical algorithms that transform the water-leaving radiance into average chlorophyll concentration of surface waters. As both models and observation products are uncertain, the assimilation of data into coupled physical-biological models provides a useful means to reduce uncertainty in modelled properties and help identify the source of errors in modelling and data processing chains. Data assimilation into coupled models can be used to (1) refine the forcing functions of physical models (Skandrani et al. 2009) that drive biology, (2) improve the representation of mesoscale features into eddy-permitting or eddy-resolving models (e.g. Martin et al. 2015), (3) reduce the range of uncertainty of model parameters in biological interaction terms (Doron et al. 2013) or (4) produce statistical estimations of the biological fields (Fontana et al. 2013). Data assimilation, however, cannot be a substitute for model improvement, since data assimilation itself becomes more effective with model skill.
In essence, data assimilation into coupled models should be treated as a fully multivariate estimation problem, i.e. one piece of data related to physical or biological properties should impact the whole coupled state vector. However, as most coupled models do not include significant dynamical feedback from biology to physics (i.e. biology can be considered as mostly 'driven' by physics to a large extent), the assimilation can be treated as a decoupled problem. In this case, however, spatial mismatches between assimilated physical and biogeochemical fields could arise (i.e. front position), owing to the different data sets used in the assimilation processes. In addition, owing to the limited amount of information available to data assimilation, controlling many different error sources (in the physical or biological modelling components) at once is often impossible (Solidoro et al. 2003). The vast majority of operational systems provide physical fields with reduced uncertainty regarding horizontal advective fields, which is already very useful to reduce errors from the physics into coupled systems. However, the consistency of assimilated vertical advective or diffusive fluxes (associated with sub-grid scale parameterizations) may be unchanged or even degraded by data assimilation ). This is recognised as a major problem with present assimilative systems that still requires dedicated research effort from the community.
In parallel, the assimilation of ocean colour data has received a growing attention by the research community during the past years. A variety of approaches have been tested [see the review by Gregg et al. (2009)] with some of them expanded today to (pre-)operational status (Ford et al. 2012;Fontana et al. 2013;Teruzzi et al. 2014aTeruzzi et al. , 2014b. A key issue for future operational applications will be to build appropriate observation operators for data-assimilation systems to efficiently link remotely sensed ocean colour radiances with phytoplankton biomass concentrations from models by using optical modules as described in Fujii et al. (2007) and Ciavatta et al. (2014). Another key development will be to ensure consistent merging of ocean colour maps at the surface with subsurface vertical profiles from a sparse Argo network equipped with nitrate, oxygen and fluorescence sensors.
Finally, it is worth noting that exploratory work has been initiated recently by applying the concept of image assimilation to ocean colour data (Gaultier et al. 2013). This could be a promising way to refine the control of eddies and associated mesoscale circulation from the physical space using ocean colour data assimilation into operational systems. It might also help to mitigate the frequent problem of the low persistence of data assimilation innovation impact. As a matter of fact, if the modification of the biogeochemical fields is not supported by a concomitant adjustment of the governing physical fields, often the biological signal quickly fades away (e.g. blooms driven by winter mixing).

Global ocean biogeochemical reanalyses of the ocean colour era
A key challenge for the future will be to develop a capacity to routinely deliver reanalyses of the biogeochemical state of oceans and regional seas, combining ocean colour and other relevant observations from satellites (e.g. SST, altimetry, surface winds, cloud coverage), in situ measurements (e.g. T/S profiles, fluorescence, oxygen, pCO 2 , nutrient concentrations) and coupled physical-biogeochemical models. Such products will be unique to monitor the impact of global change on marine primary production, regional ecosystems, the role of oceans in anthropogenic carbon storage, the occurrence statistics and trends of harmful algal blooms etc.
The ocean colour era started with the CZCS satellite mission in the late seventies, and since 1997 ocean colour has been monitored almost continuously by SeaWiFS, Envisat/MERIS and MODIS missions after an over 10year gap between 198510year gap between and 199710year gap between (Wilson, 2010. In spite of significant progress in sensor technologies, algorithmic tools and spatial resolutions, the production of long-term time series needed to capture interannual to decadal variability still requires a sustained effort to develop and maintain the required ocean colour constellation and intercalibration between instruments (Antoine et al. 2005). Additional opportunities for placing sensors on geostationary platforms in synergy with polar orbiting missions should be promoted (Ruddick et al. 2014) to enable the observation of the fine, mesoscale to submesoscale signature of surface phytoplankton and its associated diel variability. A key scientific challenge is indeed to incorporate the effect of those very fine space and time-scales on the local, basin and large-scale biogeochemical budgets estimated by integrative models (Lévy et al. 2012). In parallel to the space component, a denser network of in situ data from advanced autonomous platforms (e.g. BioArgo) will be required to extrapolate vertically the surface information s180 M. Gehlen et al.
collected from ocean colour satellites into the entire euphotic layer through assimilative models (Claustre et al. 2010). The sustained effort on biogeochemical monitoring has to be paralleled by a targeted development of biogeochemical models in order to allow the best possible utilisation of novel data. For example, current operational biogeochemical and lower trophic ecosystem models use very sophisticated treatments of the hydrodynamics and fairly sophisticated biology, but use grossly oversimplified treatments of the optics. Most models used for global and regional predictions are often driven by daily-averaged or hourly above-surface irradiances and very simple analytical formulations of light penetration through the water column. Some recent models have begun to include feedbacks between biology and hydrodynamics (Danabasoglu et al. 2006;Jochum et al. 2010), but still rely on very simple treatments of the optics. Fujii et al. (2007) incorporated spectrally resolved inherent optical properties (IOPs) into a coupled physical-biogeochemical model (Chai et al. 2002) connected with a radiative transfer model (Eco-Light-S) (Mobley, 2011). This fully coupled physical-biogeochemical-optical model simulates detailed and spectrally resolved underwater light field that not only can improve biological simulations but also can model explicitly biological feedbacks to ocean temperatures. The fully integrated model can be improved with direct comparison of model-predicted and satellite-measured remote-sensing reflectance, without the intermediate step of converting satellite data to chlorophyll. The model output also can be compared and improved with in situ optical measurements on various platforms, such as buoys, Argo floats, gliders or AUVs. Finally, all satellite and in situ measurements of optical variables can be assimilated into a fully coupled physical-biogeochemicalopticals model to constrain ecosystem conditions and improve forecasts.

Ocean carbon cycle
Operational oceanography systems expanded to biogeochemistry have the potential to contribute to ocean carbon cycle research by providing (1) near-real-time monitoring of surface ocean pCO 2 , pH and CO 2 fluxes and (2) multi-year reanalyses of the ocean carbon cycle. These products will ultimately serve a variety of purposes, such as environmental monitoring (e.g. of pH), yearly global carbon cycle assessment (Global Carbon Project, http://www.globalcarbonproject.org/), and academic carbon cycle research. Since the first review of potential contributions of operational systems to carbon cycle research by Brasseur et al. (2009), the overall volume of carbon system observations (ALK, DIC, pCO 2 , CO 2 flux, pH) has increased, and new data compilations have been released. Despite this positive trend, data coverage is still too sparse to consider operational carbon system data assimilation. Since pCO 2 is a function of salinity and temperature, improved distributions of these variables and of their temporal-spatial variability should translate into increased model skill for pCO 2 . Coupled biogeochemical-physical forecast systems could contribute to ocean carbon-cycle research by providing ocean physical-state estimates as forcing for the biogeochemical model component. Modelled pCO 2 is, however, very sensitive to errors in salinity and temperature. Errors introduced by data assimilation into ocean physics will be amplified in carbon-related variables such as pCO 2 ). Shortcomings in current biogeochemical models are an additional source of error of modelled pCO 2 . Despite current limitations, the potential of assimilation of pCO 2 data into coupled physical-biogeochemical models for improving modelled pCO 2 distributions needs to be emphasized. Recent efforts to produce maps of surface ocean pCO 2 by means of statistical techniques (Landschützer et al. 2013) hold promise of an improved data availability in the years to come.

Coupling up the food chain
The current set of modelling systems for the higher trophic levels is described by a wide range of conceptual approaches and purposes. One strong focal point here is the development of systems built around the population of single species that are commercially important, such as the above-mentioned models for tuna and sand eel. On the other side of the spectrum, there are the whole-system approaches that can be largely divided into food-web and size-spectrum approaches (Blanchard et al. 2012) assessing ecosystem services in the wider sense. The food-web-based models (Shin et al. 2004;Christensen et al. 2009) consider extensive and complex food-web structures well beyond the complexity level of current LTE models but are consequently operating at coarser time and space scales than the current operational systems offer and that may be required for some important management applications. Size-spectrum models synthesize the energy transfer across the primary producers to top predators in a much lighter framework but lose information on taxonomical and functional diversity beyond size. A step forward here may be the combination of some of these approaches, similar to Fernandes et al. (2013), who combined a bio-climate envelope model with a size-spectrum approach, or Travers-Trolet et al. (2014), who coupled a regional physical model to an NPZD representation of the LTE and an individual-based multi-species higher trophic ecosystem model to assess the long-term impact of climate change on marine fish distributions.
Some major issues concerning the end-to-end modelling from biogeochemistry to top predators are largely still unaddressed (Rose et al. 2010). A main concern, particularly in the operational context, is that all of the above approaches suffer from a substantial lack of data suitable for independent validation of their outputs in most areas of the global ocean, an essential requirement in order to demonstrate maturity for operational systems.
Another key point is the identification of the zooplankton modules as the weak link in the coupled systems, as these modules were often developed with emphasis on the adequate description of the biogeochemical cycles, but do not give a good representation of their role as prey to higher-level predators.
Future research will consequently have to provide an adequate assessment of these frameworks and approaches as to which is most adequate for each operational context in terms of forecast precision, relevance of the information provided and the feasibility of implementation into operational systems in terms of computational requirements. A significant step in this direction is currently being pursued in the EU-FP7 Operational Ecology project (OpEc, www.marine-opec.eu), where a series of near-operational forecasting systems for the physical, biogeochemical and ecological state of the European Seas are implemented and compared, bringing together the variety of the above-mentioned approaches in one framework that takes a strong focus on the operational delivery of relevant and reliable ecological indicators.

Recommendations
Several demonstrations have been made recently, with Lazzari et al. (2010), Ford et al. (2012) and Fontana et al. (2013) showing that concepts of biogeochemical data assimilation using state-of-the-art BOGCMs are getting mature enough to be transitioned to operational centres. The capacity for producing biogeochemical forecasts, along with reanalyses, is a key requirement for the development of biogeochemical/ecological applications and services. However, several blocking points still need to be considered which require further research: . the identification of the appropriate level of complexity that biogeochemical models need to include in their representation of biological processes to properly 'ingest' the observed data; . the objective characterization of the uncertainty of BOGCMs associated with physical atmospheric and ocean circulation forcings, unresolved space and time-scales, and parameterizations of biogeochemical processes and fluxes; . the coupling of BOGCMs with bio-optical models to assimilate the radiances observed by satellites in the context of 'the situation of the day' instead of derived products that often rely on climatologies (e.g. chlorophyll); . the development of coupled physical-biogeochemical data-assimilation schemes; . the development of ocean colour data-assimilation methodologies that will be able to tackle highly non-linear model responses and non-Gaussian probability distributions associated with ensemble coupled physical-biological models, as well as, in the case of multivariate assimilation schemes, to ensure appropriate projection of observed fields onto unobserved state variables using physically balanced increments ; . the specification of relevant sets of metrics to assess the scientific value of operational products, the development of intercomparison exercises between operational groups and the identification of weaknesses in the observing systems; . the use of Observing System Simulation Experiments to inform on the optimal architecture of an observing system designed for biogeochemical data assimilation and to identify biogeochemical products that can be robustly provided; . improvement in biogeochemical datasets for better initial conditions and model validation.

Supplemental data and research materials
Supplemental data for this article will be hosted alongside the article on tandfonline.com and can be accessed at http://dx.doi.org/10.1080/1755876X.2015.1022350.