Evaluation of an operational ocean model configuration at 1/12° spatial resolution for the Indonesian seas (NEMO2.3/INDO12) – Part 2: Biogeochemistry

. In the framework of the INDESO (Infrastructure Development of Space Oceanography) project, an operational ocean forecasting system was developed to monitor the state of the Indonesian seas in terms of circulation, biogeochemistry and ﬁsheries. This forecasting system combines a suite of numerical models connecting physical and biogeochemical variables to population dynamics of large marine predators (tunas). The physical–biogeochemical coupled component (the INDO12BIO conﬁguration) covers a large region extending from the western Paciﬁc Ocean to the eastern Indian Ocean at 1 / 12 ◦ horizontal resolution. The NEMO-OPA (Nucleus for European Model of the Ocean) physical ocean model and the PISCES (Pelagic Interactions Scheme for Carbon and Ecosystem Studies) biogeochemical model are running simultaneously (“online” coupling), at the same resolution. The operational global ocean forecasting system (1 / 4 ◦ ) operated by Mercator Océan provides the physical forcing, while climatological open boundary conditions are prescribed for the biogeochemistry. This paper describes the skill assessment of the INDO12BIO conﬁguration. Model skill is assessed by evalu-ating a reference hindcast simulation covering the last 8 years (2007–2014). Model results are compared to satellite, climatological and in situ observations. Diagnostics are performed on nutrients, oxygen, chlorophyll a , net primary production and mesozooplankton.Themodelreproduces large-scale distributions of nutrients, oxygen, chlorophyll a , net primary production and mesozooplankton biomasses. Modelled vertical distributions of nutrients and oxygen are comparable to in situ data sets although gradients are slightly smoothed. The model simulates realistic biogeochemical characteristics of North Paciﬁc tropical waters entering in the archipelago. Hydrodynamic transformation of water masses across the Indonesian archipelago allows for conserving nitrate and oxygen vertical distribution close to observations, in the Banda Sea and at the exit of the archipelago. While the model overestimates the mean surface chlorophyll a , the seasonal cycle is in phase with satellite estimations, with higher chlorophyll a concentrations in the southern part of the archipelago during the SE monsoon and in the northern part during the NW monsoon. The time series of chlorophyll a anomalies suggests that meteorological and ocean physical processes that drive the interannual variability of biogeochemical properties in the Indonesian region are reproduced by the model.


Introduction
The "coral triangle" delineated by Malaysia, the Philippines, New Guinea, Solomon Islands, East Timor and Indonesia is recognised as a global hotspot of marine biodiversity (Allen and Werner, 2002;Mora et al., 2003;Green and Mous, 2004;Allen, 2008).It gathers 20 % of the world's species of plants and animals, and the greatest concentration and diversity of reefs (76 % of the world's coral species; Veron et al., 2009).

E. Gutknecht et al.: Evaluation of an operational ocean model for the Indonesian seas -Part 2
The Indonesian archipelago is located at the centre of this ecologically rich region.It is characterised by a large diversity of coastal habitats such as mangrove forests, coral reefs and sea grass beds, all of which shelter ecosystems of exceptional diversity (Allen and Werner, 2002).The archipelago's natural heritage represents an important source of income and employment, with its future critically depending on the sustainable management of ecosystems and resources (e.g.Foale et al., 2013;Cros et al., 2014).
The wider coral triangle and its sub-region, the Indonesian archipelago, are facing multiple threats resulting from demographic growth, economic development, change in land use practices and deforestation, as well as global climate change (http://www.metoffice.gov.uk/media/pdf/8/f/Indonesia.pdf;FAO, 2007).Human activities cause changes in the delivery of sediments, nutrients and pollutants to coastal waters, leading to eutrophication, ecosystem degradation, as well as species extinctions (Ginsburg, 1994;Pimentel et al., 1995;Bryant et al., 1998;Roberts et al., 2002;UNEP, 2005;Alongi et al., 2013).Surveys report an over 30 % reduction of mangroves in northern Java over the last 150 years and an increase of coral reef degradation from 10 to 50 % in the last 50 years (Bryant et al., 1998;Hopley and Suharsono, 2000;UNEP, 2009), leading to 80 % of the reefs being at risk in this region (Bryant et al., 1998).These changes not only damage coastal habitats, but also propagate across the whole marine ecosystem from nutrients and the first levels of the food web up to higher trophic levels, along with concomitant changes in biogeochemical cycles.
There is thus a vital need for monitoring and forecasting marine ecosystem dynamics.The INDESO (Infrastructure Development of Space Oceanography; www.indeso.web.id/indeso_wp/index.php) project, funded by the Indonesian Ministry of Marine Affairs and Fisheries, aims at the development of sustainable fishery practices in Indonesia, the monitoring of its exclusive economic zone (EEZ) and the sustainable management of its ecosystems.The project addresses the Indonesian need for building a national capability for operational oceanography.The model system consists of three models deployed at the scale of the Indonesian archipelago: an ocean circulation model (NEMO-OPA; Madec, 2008), a biogeochemical model (PISCES; Aumont and Bopp, 2006) with a spatial resolution of 1/12 • , as well as an intermediate trophic level/fish population dynamics model (SEAPODYM -Spatial Ecosystem and Populations Dynamics Model; Lehodey et al., 2008).Since mid-September 2014, the chain of models is fully operational in Perancak (Bali, Indonesia) and delivers 10-day forecast/2-weeks hindcast on a weekly basis (see http://www.indeso.web.id).
The regional ocean dynamics is fully described in Tranchant et al. (2016; hereafter Part 1).The physical model reproduces main processes occurring in this complex oceanic region.Ocean circulation and water mass transformation through the Indonesian Archipelago are close to observations.eddy kinetic energy displays patterns similar to satel-lite estimates, tides being a dominant forcing in the area.The volume transport of the Indonesian ThroughFlow is comparable to INSTANT data.Temperature-salinity diagrams diagrams highlight the erosion of South and North Pacific subtropical waters while crossing the archipelago.
The present paper (Part 2) focuses on ocean biogeochemistry.It is organised as follows.The next section presents an overview of the area of study with emphasis on main drivers of biological production over the Indonesian archipelago.The biogeochemical component of the physical-biogeochemical coupled configuration is described in Sect.3. Satellite, climatological and in situ observations used to evaluate simulation results are detailed in Sect. 4. Section 5 presents the evaluation of the skill of the coupled model to reproduce main biogeochemical features of Indonesian seas along with their seasonal and interannual dynamics.Finally, discussion and conclusion are presented in Sect.6.

Area of study
The Indonesian archipelago is crossed by North and South Pacific waters that converge in the Banda Sea, and leave the archipelago through three main straits: Lombok, Ombaï and Timor.This ocean current (Indonesian ThroughFlow; ITF) provides the only low-latitude pathway for warm, fresh waters to move from the Pacific to the Indian Ocean (Gordon, 2005;Hirst and Godfrey, 1993).On their way through the Indonesian archipelago, water masses are progressively transformed by surface heat and freshwater fluxes and intense vertical mixing linked to strong internal tides trapped in the semi-enclosed seas as well as upwelling processes (Ffield and Gordon, 1992).The main flow, as well as the transformation of Pacific waters is correctly reproduced by the physical model, with a realistic distribution of the volume transport through the three major outflow passages (Part 1).In the Indian Ocean, this thermocline water mass forms a cold and fresh tongue between 10 and 20 • S, and supplies the Indian Ocean with nutrients.These nutrients impact biogeochemical cycles and support new primary production in the Indian Ocean (Ayers et al., 2014).
Over the archipelago, complex meteorological and oceanographic conditions drive the distribution and growth of phytoplankton and provide favourable conditions for the development of a diverse and productive food web, extending from zooplankton to pelagic fish (Hendiarti et al., 2004(Hendiarti et al., , 2005;;Romero et al., 2009).The tropical climate is characterised by a monsoon regime and displays a well-marked seasonality.The south-east (SE) monsoon (April to October) is associated with easterlies from Australia that carry warm and dry air over the region.Wind-induced upwelling along the southern coasts of Sumatra, Java and Nusa Tenggara islands (hereafter named Sunda islands) and in the Banda Sea is associated with high chlorophyll a levels (Susanto et al., 2006;Rixen et al., 2006).Chlorophyll a maxima along Sunda is-lands move to the west over the period of the SE monsoon, in response to the alongshore wind shift and associated movement of the upwelling centre (Susanto et al., 2006).From October to April, the north-west (NW) monsoon is associated with warm and moist winds from the Asian continent.Winds blow in a south-west direction north of the Equator and towards Australia south of the Equator.They generate a downwelling and a reduced chlorophyll a content south of the Sunda islands and in the Banda Sea.The NW monsoon also causes some of the highest precipitation rates in the world.Increased river runoff carries important sediment loads (20 to 25 % of the global riverine sediment discharge; Milliman et al., 1999), along with carbon and nutrients to the ocean.These inputs are a strong driver of chlorophyll a variability and play a key role in modulating the biological carbon pump across Indonesian seas (Hendiarti et al., 2004;Rixen et al., 2006).High levels of suspended matter decrease the water transparency in coastal areas and modify the optical properties of waters, which in turn interferes with ocean colour remote sensing (Susanto et al., 2006).Although several Indonesian rivers are classified among the 100 most important rivers of the world, most of them are not regularly monitored.It is thus currently impossible to estimate the impact of river runoff on the variability of chlorophyll a in the region (Susanto et al., 2006).
Indonesian seas are also greatly influenced by modes of natural climate variability owing to its position on the Equator between Asia and Australia and between the Pacific and Indian oceans.Strength and timing of the seasonal monsoon are modulated by interannual phenomena that disturb atmospheric conditions and ocean currents.A significant correlation between the variability of the ITF and the El Niño-Southern Oscillation (ENSO) was reported (e.g.Meyers, 1996;Murtugudde et al., 1998;Potemra et al., 1997), with ENSO modulating rainfall and chlorophyll a on interannual timescales (Susanto et al., 2001(Susanto et al., , 2006;;Susanto and Marra, 2005).ENSO can be monitored using a multivariate ENSO index (MEI; Wolter andTimlin, 1993, 1998; http://www.esrl.noaa.gov/psd/enso/mei/).In the eastern Indian Ocean, large anomalies off the Sumatra and Java coasts are associated with the Indian Ocean dipole (IOD) mode monitored via the dipole mode index (DMI; Saji et al., 1999).A strong positive index points to abnormally strong coastal upwelling and a large phytoplankton bloom near Java Island (Meyers, 1996;Murtugudde et al., 1999).Inside the archipelago, effects of ENSO and IOD climate modes are more difficult to discriminate as they both influence ITF transport.There is, however, evidence for Indian Ocean dynamics to dominate over Pacific Ocean dynamics as drivers of ITF transport variability (Masumoto, 2002;Sprintall and Révelard, 2014).

The coupled model
In the framework of the INDESO project, a physicalbiogeochemical coupled model is deployed over the domain from 90-144 • E to 20 • S-25 • N, widely encompassing the whole Indonesian archipelago, with a spatial resolution of 1/12 • .The physical model is based on the NEMO-OPA 2.3 circulation model (Madec et al., 1998;Madec, 2008).Specific improvements include time splitting and non-linear free surface to correctly simulate high-frequency processes such as tides.A parameterisation of the vertical mixing induced by internal tides has been developed especially for NEMO-OPA (Koch-Larrouy et al., 2007, 2010) and is used here.The physical configuration called INDO12 is described in detail in Part 1.
Dynamics of biogeochemical properties across the area are simulated by the PISCES model version 3.2 (Aumont and Bopp, 2006).PISCES simulates the first levels of the marine food web from nutrients up to mesozooplankton.It has 24 state variables.PISCES considers five limiting nutrients for phytoplankton growth (nitrate and ammonium, phosphate, dissolved silica and iron).Four living size-classified compartments are represented: two phytoplankton groups (nanophytoplankton and diatoms) prognostically predicted in carbon (C), iron (Fe), silica (Si) (the latter only for diatoms) and chlorophyll content, and two zooplankton groups (microzooplankton and mesozooplankton).Constant Carbon / Nitrogen / Phosphorus (C / N / P) Redfield ratios are supposed for all species.While internal Fe / C and Si / C ratios of phytoplankton are modelled as a function of the external availability of nutrients and thus variable, only C is prognostically modelled for zooplankton.The model includes five non-living compartments: small and big particulate organic carbon and semi-labile Dissolved Organic Carbon (DOC), particulate inorganic carbon (CaCO 3 as calcite) and biogenic silica.PISCES also simulates Dissolved Inorganic Carbon (DIC), total alkalinity (carbonate alkalinity + borate + water), and dissolved oxygen.The CO 2 chemistry is computed following the OCMIP protocols (http://ocmip5.ipsl.jussieu.fr/OCMIP/).Biogeochemical parameters are based on the standard PISCES namelist version 3.2.Please refer to Aumont and Bopp (2006) for a comprehensive description of the model (version 3.2).
PISCES is coupled to NEMO-OPA via the TOP component that manages the advection-diffusion equations of passive tracers and biogeochemical source and sink terms.In our regional configuration, called INDO12BIO, physics and biogeochemistry are running simultaneously ("online" coupling), at the same resolution.Particular attention must be paid to respect a number of fundamental numerical constraints.(1) The numerical scheme of PISCES for biogeochemical processes is forward in time (Euler), which does not correspond to the classical leap-frog scheme used for www.geosci-model-dev.net/9/1523/2016/Geosci.Model Dev., 9, 1523-1543, 2016 the physical component.Moreover, the free surface explicitly solved by the time-splitting method is non-linear.In order to respect the conservation of the tracers, the coupling between biogeochemical and physical components is done every second time step.As a result, the biogeochemical model is controlled by only one leap-frog trajectory of the dynamical model.The use of an Asselin filter allows for keeping the two numerical trajectories close enough to overcome this shortcoming.The advantage is a reduction of numerical cost and a time step for the biogeochemical model twice that of the physical component, i.e. 900 s. (2) As this time step is small, no time splitting was used in the sedimentation scheme.
(3) The advection scheme is the standard scheme of TOP-PISCES, i.e. the monotonic upstream-centred scheme for conservation laws (MUSCL) (Van Leer, 1977).No explicit diffusion has been added as the numerical diffusion introduced by this advection scheme is already important.

Initial and open boundary conditions
The simulation starts on 3 January 2007 from the global ocean forecasting system at 1/4  1. Nitrate, phosphate, dissolved silica, oxygen, DIC and alkalinity are derived from climatological data sets.For tracers for which this information is missing, initial and open boundary conditions come either from a global-scale simulation, are estimated from satellite data or are build using analytical values.The global-scale model NEMO-OPA/PISCES has been integrated for 3000 years at 2 • horizontal resolution, until PISCES reached a quasisteady state (see Aumont and Bopp, 2006).A monthly climatology was built for dissolved iron and DOC based on this simulation.A Dirichlet boundary condition is used to improve the information exchange between the OBC and the interior of the domain.

External inputs
Three different sources are supplying the ocean with nutrients: atmospheric dust deposition, sediment mobilisation and rivers.Atmospheric deposition of iron comes from the climatological monthly dust deposition simulated by the model of Tegen and Fung (1995), and that of silica follows Moore et al. (2002).Yearly means of river discharges are taken from the global erosion model (GEM) of Ludwig et al. (1996) for DIC, and from the global news 2 climatology (Mayorga et al., 2010) for nutrients.An iron source corresponding to sediment reductive mobilisation on continental margins is also considered.For more details on the external supply of nutrients, please refer to the Supplement of Aumont and Bopp (2006).The improved representation of the contribution of local processes to external nutrient supply, as well as of the seasonal variability of river nutrient delivery is hampered by the lack of in situ observations.
In PISCES, external input fluxes are compensated by a loss to the sediments as particulate organic matter, biogenic silica and CaCO 3 .These fluxes correspond to matter definitely lost from the ocean system.The compensation of external input fluxes through output at the lower boundary closes the mass balance of the model.While such equilibrium is a valid assumption at the scale of the global ocean, it is not reached at regional scale.For the INDO12BIO configuration, a decrease of the nutrient and carbon loss to the sediment was introduced corresponding to an increase in the water column re-mineralisation by ∼ 4 %.This slight enhancement of water column re-mineralisation leads to higher coastal chlorophyll a concentrations (about +1 mg Chl m −3 ) and enables the model to reproduce the chlorophyll a maxima observed along the coasts of Australia and East Sumatra (not shown).

Simulation length
The simulation starts on 3 January 2007 and operates up to present day, as the model currently delivers ocean forecasts.For the present paper, we will analyse the simulation up to 31 December 2014.The spin-up length depends on the biogeochemical tracer (Fig. 1).The total carbon inventory computed over the domain (defined as the sum of all solid and dissolved organic and inorganic carbon fractions, yet dominated by the contribution of DIC) equilibrates within several months.To the contrary, DOC, phosphate (PO 4 ) and iron (Fe) need several years to stabilise (Fig. 1).The annual mean for year 2011 is used for comparison to satellite products (chlorophyll a, net primary production).For comparison to climatologies (zooplankton, nutrients, oxygen) and analysis of the seasonal cycle, we use years 2010 to 2014.Interannual variability is assessed over the whole length of simulation except the first year (2008 to 2014).

Satellite, climatological and in situ data
Model outputs are compared to satellite, climatological and in situ observations.These observational data are detailed and described in this section.

INDOMIX cruise
The (цmolC L -1 ) ) ) at the entrance of the archipelago (Halmahera Sea), in the Banda Sea and in the Ombaï Strait (three of them are used for validation; cf.stations on Fig. 4).This data set provides an independent assessment of model skill.To co-localise the model and observations, we took the closest simulated point to the coordinates of the station; 2-day model averages were considered as measurements were performed during 2 consecutive days at the stations selected for validation.

Chlorophyll a
The ocean colour signal reflects a combination of chlorophyll a content, suspended matter, coloured dissolved organic matter (CDOM) and bottom reflectance.Singling out the contribution of phytoplankton's chlorophyll a is not straightforward in waters for which the relative optical contribution of the three last components is significant.This is the case over vast areas of the Indonesian archipelago where river discharges and shallow water depths contribute to optical properties (Susanto et al., 2006).The interference with optically absorbing constituents other than chlorophyll a results in large uncertainties in coastal waters (up to 100 %, as compared to 30 % for open ocean waters) (Moore et al., 2009).Standard algorithms distinguish between open ocean waters/clear waters (case-1) and coastal waters/turbid waters (case-2).The area of deployment of the model comprises waters of both categories and the comparison between modelled chlorophyll a and estimates derived from remote sensing can be only qualitative.Two singlemission monthly satellite products are used for model skill evaluation.MODIS-Aqua ((Moderate Resolution Imaging Spectroradiometer, EOS mission, NASA)") level-3 standard  Due to the large uncertainty in production models, here we compare the simulated NPP to NPP derived from the three aforementioned models using MODIS ocean colour data.

Mesozooplankton
MAREDAT, MARine Ecosystem DATa, (Buitenhuis et al., 2013) is a collection of global biomass data sets for major plankton functional types (e.g.diatoms, microzooplankton, mesozooplankton).Mesozooplankton is the only MAREDAT field covering the Indonesian archipelago.The     3).The flow across the Indonesian archipelago and the transformation of water masses simulated by the model result in realistic vertical distributions of nutrients and oxygen concentrations in the Banda Sea.The ITF leaves the archipelago and spreads into the Indian Ocean with a biogeochemical content in good agreement with the data available in the area.
However, simulated vertical structures are slightly smoothed compared to data (Fig. 3).The vertical gradient of nitrate is too weak over the first 2000 m depth of the water column (North Pacific and Timor), and the area of minima oxygen concentrations is eroded (especially in North Pacific box).This bias is even more pronounced on the vertical gradient of dissolved silica (Fig. 3).The smoothing of vertical structures results from the numerical advection scheme MUSCL currently used in PISCES, which is known to be too diffusive (Lévy et al., 2001).

Chlorophyll a and NPP
The simulation reproduces the main characteristics of the large-scale distribution of chlorophyll a, a proxy of phytoplankton biomass (Fig. 4 ).The bias (as model -observation) is almost positive everywhere, except around the coasts (discussed later) and in the Sulawesi Sea.As mentioned in the preceding section, optical characteristics of waters over the Indonesian archipelago are closer to case-2 waters (Moore et al., 2009).Simulated chlorophyll a concentrations are indeed closer to those derived with an algorithm for case-2 waters (MERIS) and its mean value of 0.48 ± 1.4 mg Chl m −3 .The model reproduces the spatial distribution, as well the rates of NPP over the model domain (Fig. 5).However, as mentioned before, NPP estimates depend on the primary production model (in this case, VGPM, CbPM and Eppley) and on the ocean colour data used in the production models.For a single ocean colour product

Mesozooplankton
Mesozooplankton links the first level of the marine food web (primary producers) to the mid-and, ultimately, high trophic levels.Modelled mesozooplankton biomass is compared to observations in Fig. 6.While the model repro- duces the spatial distribution of mesozooplankton, it overestimates biomass by a factor 2 or 3.This overestimation is likely linked to the above-described overestimation of chlorophyll a and NPP.

Mean seasonal cycle
The monsoon system drives the seasonal variability of chlorophyll a over the area of study.Northern and southern parts of the archipelago exhibit a distinct seasonal cycle (Figs.7, 8 and 9).In the southern part, the highest chlorophyll a concentrations occur from June to September (Banda Sea and Sunda area in Figs. 8 and 9) due to upwelling of nutrient-rich waters off Sunda islands and in the Banda Sea triggered by alongshore south-easterly winds during the SE monsoon.The decrease in chlorophyll levels during the NW monsoon is the consequence of north-westerly winds and associated downwelling in these same areas.In the northern part, high chlorophyll concentrations occur during the NW monsoon (South China Sea in Fig. 7) when moist winds from Asia cause intense precipitation.A secondary peak is observed during the NW monsoon in the southern part and during the SE monsoon in the northern part due to meteorological and oceanographic conditions described above.
The annual signal of chlorophyll a in each grid point gives a synoptic view of the effect of the Asia-Australia monsoon system on the Indonesian archipelago.A harmonic analysis is applied on the time series of each grid point to extract the annual signal in model output and remote sensing data (MODIS).The results of the annual harmonic analysis are summarised in Fig. 10 and highlight the month of maximum chlorophyll a and the amplitude of the annual signal.The timing of maximum chlorophyll a presents a north-south distribution in agreement with the satellite observations.The simulation reproduces the chlorophyll a maxima in July in the Banda Sea and off the south coasts of Java-Nusa Tenggara.Consistent with observations, simulated chlorophyll a maxima move to the west over the period of the SE monsoon, in response to the alongshore wind shift.North of the Nusa Tenggara islands, maxima in January-February are due to upwelling associated with alongshore north-westerly winds.In the South China Sea, maxima spread from July-August in the western part (off Mekong River) and gradually shift up to January-February in the eastern part.
The temporal correlation between modelled chlorophyll a and estimates derived from remote sensing is 0.55 over the entire INDO12BIO domain, but reaches 0.78 in the South China Sea, 0.81 in the Banda Sea and 0.93 in the Indian Ocean (Figs. 7,8,9 and 11).These high correlation coefficients are associated with low normalised standard deviations (close to 1) in the Banda Sea and in the Indian Ocean (Fig. 11) as well as large amplitudes in simulated and ob- served chlorophyll a (Fig. 10).Normalised standard deviations are higher in the south-east China Sea, Java and Flores seas as well as in the open ocean due to larger amplitudes in simulated chlorophyll a.The offshore spread of the high amplitude reflects the too weak cross-shore gradient of simulated chlorophyll a (Sect.5.1.2),and leads to an increase of the normalised standard deviation with the distance to the coast.For semi-enclosed seas, however, this result has to be taken with caution as clouds cover these regions almost 50-60 % of the time period.
The model does not succeed in simulating chlorophyll a variability in the Pacific sector (Figs. 10 and 11).This area is close to the border of the modelled domain and is influenced by the OBCs derived from the global operational ocean general circulation model.Analysis of the modelled circulation (Part 1) highlights the role of OBCs in maintaining realistic circulation patterns in this area, which is influenced by the equatorial current system.Part 1 points, in particular, to the incorrect positioning of Halmahera and Mindanao eddies in the current model, which contributes to biases in simulated biogeochemical fields.
Finally, correlation is low close to the coasts and the temporal variability of the model is lower than that of the satellite product, with normalised standard deviation < 1 (Fig. 11).The model does not take into account seasonal variations in river discharges.Driven by the monsoon system, seasonal input of river runoff is an important driver of chlorophyll a variability at local scale.

Interannual variability
Figures 7, 8 and 9 present interannual anomalies of surface chlorophyll a concentrations between 2008 and 2014 for model outputs and MODIS ocean colour averaged over three regions: South China Sea, Banda Sea and the Sunda area.Simulated fields and satellite-derived chlorophyll a are in good agreement in terms of amplitude and phasing, with temporal correlation coefficients of 0.56 for the South China Sea and Banda Sea and 0.88 for the Sunda area.The model simulates a realistic temporal variability suggesting that processes IOD drives the chlorophyll a interannual variability in the eastern tropical Indian Ocean, with a correlation coefficient of 0.74 (Fig. 9).IOD index and anomalies of chlorophyll a from satellites give a similar correlation coefficient of 0.7.A positive phase of IOD indicates negative SST anomaly in the south-eastern tropical Indian Ocean associated with zonal wind anomaly along the Equator (Meyers, 1996).The abnormally strong coastal upwelling near the Java Island stimulates a large phytoplankton bloom (Murtugudde et al., 1999).In the Banda Sea and in the South China Sea, no clear impact of ENSO or IOD is detected on the first level of the food chain (Figs. 7,8).Inside the archipelago, both climate modes affect the variability of the ITF transport and it is not straightforward to separate their individual contribution (Masumoto, 2002;Sprintall and Révelard, 2014).
While it is established (see references cited in Sect.2) that ENSO and IOD climate modes play a key role in the Indonesian region, their impact on the marine ecosystem remains poorly understood.The length of simulation is too short for a rigorous assessment of the role of these drivers and a direct relationship is only evident in the Indian sector.However, interannual anomalies of simulated chlorophyll a compare well to satellite observations, which suggests that interannual meteorological and ocean physical processes are satisfyingly reproduced by the model.

INDOMIX cruise
Model results are compared to INDOMIX in situ data at three key locations: (1) the eastern entrance of Pacific waters to the archipelago (station 3, Halmahera Sea), (2) the convergence of the western and eastern pathways (station 4, Banda Sea) where intense tidal mixing and upwelling transforms Pacific waters to form the ITF and (3) one of the main exit portals of the ITF to the Indian Ocean (station 5, Ombaï Strait).
The vertical profile of temperature compares well to the data in the Halmahera Sea (Fig. 12).Simulated surface waters are too salty and the subsurface salinity maximum is reproduced at the observed depth, albeit underestimated compared to the data.Waters are more oxygenated in the model over the first 400 m.The model-data bias on temperature, salinity and oxygen suggests that Halmahera Sea thermocline waters are not correctly reproduced by the model in July 2010.The model tends to yield too smooth vertical profiles.Vertical profiles of nitrate and phosphate are well reproduced, while dissolved silica concentrations are overestimated below 200 m depth.It should be noted, however, that 2010 was a strong La Niña year with important modifications in zonal winds, rainfall, river discharges and ocean currents.While interannual variability is taken into account in atmospheric forcing and physical open boundary conditions, external nutrient inputs from rivers are constant, and biogeochemical OBCs come from climatologies.However, dissolved silica profiles computed from the monthly WOA2009 climatology are close to simulated distributions (not shown), suggesting non-standard conditions during the time of the INDOMIX cruise.
Despite the bias highlighted for Halmahera Sea station, an overall satisfying correspondence between modelled and observed profiles is found at the Banda Sea (Fig. 13) and Ombaï Strait stations (Fig. 14).The comparison of modelled profiles and cruise data along the flow path of waters from the Pacific to the Indian Ocean (from Halmahera to Ombaï Strait) suggests that either the Halmahera Sea had no major influence for the ITF formation during the time of the cruise, or that vertical mixing and upwelling processes across the archipelago are strong enough to allow for the formation of Indonesian water masses despite biases in source water composition.Alternatively, it could reflect the weak impact of ENSO on biogeochemical tracer distributions inside the archipelago compared to its Pacific border and the dominant role of Indian ocean dynamics on the ITF (Sprintall and Révelard, 2014).

Discussions and conclusions
The INDESO project aims to monitor and forecast marine ecosystem dynamics in Indonesian waters.A suite of numerical models were coupled for setting up a regional configuration (INDO12) adapted to Indonesian seas.A forecasting oceanographic centre has been fully operational in Perancak In coastal waters, chlorophyll a concentrations are influenced by sedimentary processes (i.e.re-mineralisation of organic carbon and subsequent release of nutrients) and riverine nutrient input.The slight disequilibrium explicitly intro-duced between the external input of nutrients and carbon and the loss to the sediment is sufficient to enhance chlorophyll a concentrations along the coasts and to make it comparable with observations.The sensitivity of the model to the balancing of carbon and nutrients at the lower boundary of the domain ("sediment burial") highlights the need for an explicit representation of sedimentary reactions.
In order to further improve modelled chlorophyll a variability along the coast, time-variant river nutrient and carbon fluxes are needed.According to Jennerjahn et al. (2004), river discharges from Java can be increased by a factor of ∼ 12 during the NW monsoon as compared to the SE monsoon.Moreover, the maximum fresh-water transport and the peak of material reaching the sea can be out of phase depending on the origin of discharged material (Hendiarti et al., 2004).The improved representation of river discharge dynamics and associated delivery of fresh water, nutrients and suspended matter in the model is, however, hampered by the availability of data as most of the Indonesian rivers are currently not monitored (Susanto et al., 2006).Systematic misfits between modelled and observed biogeochemical distributions may in part also reflect inherent properties of implemented numerical schemes.Misfits highlighted throughout this work include too much chlorophyll a and NPP on the shelves, with too weak cross-shore gradients between shelf and open waters, together with noticeable smoothing of vertical profiles of nutrients and oxygen.Currently, the MUSCL advection scheme is used for biogeochemical tracers.This scheme is too diffusive and smooths vertical profiles of biogeochemical tracers.As a result, too many nutrients are injected in the surface layer and trigger high levels of chlorophyll a and NPP.Another advection scheme, OUICKEST (Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms; Leonard, 1979) with the limiter of Zalezak (1979), already used in NEMO for the advection scheme of the physical model, has been tested for biogeochemical tracers.Switching from MUSCL to QUICKEST-Zalezak accentuates the vertical gradient of nutrients in the water column and attenuates modelled chlorophyll a and NPP.This advection scheme is not diffusive and its use would be coherent with choices adopted for physical tracers.However, it would result in an overestimation of the vertical gradient of nutrients, and the nutricline would be considerably strengthened.Neither tuning of biogeochemical parameters, nor switching the advection scheme for passive tracers fully resolved the model-data misfits.Hence, improving the vertical distribution of nutrients and oxygen, as well as chlorophyll a and NPP in the open ocean and their cross-shore gradient first requires improving the model physics.
Finally, monthly or yearly climatologies are currently used for initial and open boundary conditions.Biogeochemical tracers are thus decorrelated from model physics.In order to improve the link between modelled physics and biogeochemistry, weekly or monthly averaged output of the global ocean operational system operated by Mercator Océan (BIOMER) will be used in the future for the 24 tracers of the biogeochemical model PISCES.BIOMER will couple the physical forecasting system PSY3 to PISCES in offline mode.The biogeochemical and the physical components of INDO-BIO12 will thus be initialised and forced coherently, on the base of the PSY3 forecasting system.

Figure 1 .
Figure 1.Temporal evolution of total carbon (a), plankton (b), DIC and DOC (c) and nutrient (d, e) content averaged over the whole three-dimensional INDO12BIO domain.

Figure 2 .Figure 3 .
Figure 2. Annual mean of nitrate (mmol N m −3 ; left) and oxygen concentrations (mL O 2 L −1 ; right) at 100 m depth from CARS (a, d) and WOA (b, e; statistical mean) annual climatologies, and from INDO12BIO as 2010-2014 averages (c, f).Three key boxes for water mass transformation (North Pacific, Banda, and Timor; Koch-Larrouy et al., 2007) were added to the bottom-right figure.

Figure 4 .
Figure 4. Left: annual mean of surface chlorophyll a concentrations (mg Chl m −3 ) for year 2011: MODIS case-1 product (a), MERIS case-2 product (b) and INDO12BIO simulation (c).Right: bias of log-transformed surface chlorophyll (model-observation) for the same year.The model was masked as a function of the observation, MODIS case-1 (d) or MERIS case-2 (e).Location of three stations sampled during the INDOMIX cruise and used for evaluation of the model in Sect.4.4 (f).

Figure 5 .
Figure 5. Annual mean of vertically integrated NPP (mmol C m −2 d −1 ) for year 2011: VGPM (a), Eppley (d) and CbPM (b) production models, all based on MODIS ocean colour, as well as for INDO12BIO (e).Standard deviation of the three averaged production models (PM) (c), and bias between INDO12BIO and the average of PM (f).

Figure 6 .
Figure 6.Annual mean of mesozooplankton biomass (µg C L −1 ) from MAREDAT monthly climatology (left) and from INDO12BIO simulation averaged between 2010 and 2014 (right), for distinct depth interval: from the surface up to 40 m (a, e), 100 m (b, f), 150 m (c, g), and 200 m depth(d, h).Simulated fields were interpolated onto the MAREDAT grid, and masked as a function of the data (in space and time).
Nitrate and oxygen distributions at 100 m depth are presented on Fig.2for CARS, WOA and the model.Dissolved silica has the same distribution as nitrate (not shown).The marked meridional gradient, seen in observations of the Pacific and Indian oceans, is correctly reproduced by the model.Low nitrate and high oxygen concentrations in the subtropical gyres of the North Pacific and South Indian oceans are due to Ekman-induced downwelling.Higher nitrate and lower oxygen concentrations in the equatorial area are associated with upwelling.Maxima nitrate concentrations associated with minima oxygen concentrations are noticeable in the Bay of Bengal and Adaman Sea (north of Sumatra and west of Myanmar).They reflect discharges by major rivers (Brahmaputra, Ganges and other river systems) and the associated

Figure 7 .
Figure 7. (a) Mean surface chlorophyll a concentrations and (b) its interannual anomalies (mg Chl m −3 ) over the South China Sea.INDO12BIO is in black and MODIS case-1 in red.Temporal correlation (r) between both time series is in black.(c) ENSO (blue) and IOD (green) phenomena are respectively represented by MEI and DMI indexes.Indexes were normalised by their maximum value in order to be plotted on the same axis.Interannual anomalies of simulated chlorophyll a are reminded in black.Temporal correlation (r) between the simulated chlorophyll a and ENSO (IOD) is indicated in blue (green).

(
here MODIS), NPP estimates display a large variability (Fig. 5).Mean NPP over the INDO12BIO domain is 34.5 mmol C m −2 d −1 for VGPM with a standard deviation over the domain of 33.8, 40.4 ± 22 mmol C m −2 d −1 for CbPM and 55 ± 52.7 mmol C m −2 d −1 for Eppley.NPP estimates from VGPM are characterised by low rates in the Pacific (< 10 mmol C m −2 d −1 ) and a well-marked cross-shore gradient.The use of CbPM results in low coastal NPP and almost uniform rates over a major part of the domain including the open ocean (Fig. 5).The Eppley production model is the most productive one with rates about 15 mmol C m −2 d −1 in the Pacific and higher than 300 mmol C m −2 d −1 in the coastal zone.The large uncertainty associated with these products precludes a quantitative evaluation of modelled NPP.Like for chlorophyll a, modelled NPP falls within the range of remote sensing derived estimates, with maybe a too weak cross-shore gradient inherited from the chlorophyll a field.The mean NPP over the INDO12BIO domain is, however, overestimated (61 ± 41.8 mmol C m −2 d −1 ).

Figure 10 .
Figure 10.Timing of maximum chlorophyll a (a, c) and amplitude (b, d) for a monthly climatology of surface chlorophyll a concentrations between 2010 and 2014: MODIS case-1 (left) and INDO12BIO (right).The model was masked as a function of the data.

Figure 11 .
Figure 11.Temporal correlation (a) and normalised standard deviation (b; SD(model)/SD(data)) estimated between the INDO12BIO simulation and the MODIS case-1 ocean colour product.Statistics are computed on monthly fields between 2010 and 2014.The model was masked as a function of the data.

Table 1 .
Key et al. (2004)an Mixing program) cruise on-board Marion DufresneRV (Koch-Larrouy et al.,  2016)crossed the Indonesian archipelago between 9 and 19 July 2010, and focused on one of the most energetic sections for internal tides from the Halmahera Sea to Ombaï Strait.Repeated CTD profiles over 24 h, as well as measurements of oxygen and nutrients, were obtained for six stations Initial and open boundary conditions used for the INDO12BIO configuration.From World Ocean Atlas (WOA, 2009) monthly climatology, with increased nutrient concentrations along the coasts (necessary adaptation due to crucial lack of data in the studied area).bKeyetal. (2004).c From SeaWiFS monthly climatology.Phytoplankton is deduced using constant ratios of 1.59 and 122/16 mol C mol N −1 , and exponential decrease with depth.d Low values offshore and increasing concentrations onshore.