A skill assessment of the biogeochemical model REcoM 2 coupled to the finite element sea-ice ocean model ( FESOM 1 . 3 )

Introduction Conclusions References


Introduction
Primary production plays a large role in ocean carbon cycling, and understanding the drivers behind primary production is therefore of paramount importance when it comes to understanding the changes that a future warmer climate will bring.Observations, as well as coupled biogeochemical-ocean models, indicate that climate change will decrease the oceanic net primary production (NPP) (Behrenfeld et al., 2006;Steinacher et al., 2010).This would have far-reaching implications, from changes of the carbon cycle to effects on fisheries.
Coupled biogeochemical-ocean models are important tools used to analyse the net primary production in the ocean and the effects of climate change on it (e.g.Le Quéré et al., 2003;Bopp et al., 2013).The biogeochemical results of such models are highly impacted by the mixing and circulation of the ocean model as it controls processes such as horizontal advection and nutrient supply to the surface layer (Doney et al., 2004).Supply of nutrients through upwelling is especially important when it comes to modelling the equatorial Pacific (Aumont et al., 1999) and the Southern Ocean, where production is iron limited and sensitive to new supply.Results from the second Ocean Carbon-Cycle Model Intercomparison Project (OCMIP-2) highlighted the importance of the ocean model; they showed how the representation of the ocean circulation in the Southern Ocean has a large impact on the calculations of present and future uptake of CO 2 (Doney et al., 2004), and reported that the global export production varied between 9 and 28 Gt C yr −1 when the same biogeochemical model was coupled to different OGCMs (Najjar et al., 2007).
Traditionally, global OGCMs employ structured grids with relatively uniform spatial resolution in the entire domain, and local refinement is done by utilizing nested models.
The unstructured-mesh technology is emerging as an alternative to nesting in ocean models, and is gradually becoming more widespread within the ocean modelling community (e.g., Chen et al., 2003;Danilov et al., 2004;Piggott et al., 2008).As solutions for the global ocean state provided by models formulated on unstructured meshes have improved (e.g., Sidorenko et al., 2011), it has become feasible to exploit the advantages offered by such models in biogeochemical modelling by coupling a biogeochemical model to an unstructured-mesh ocean model (Hill et al., 2014).One may then benefit from the possibility of aligning the grid with the bathymetry, or refining it in areas of interest without the loss of accuracy that the nesting introduces at boundaries.This is especially relevant when it comes to modelling features such as mixed layer depth, upwelling and the presence of fronts and eddies that are of vital importance for realistic modelling of ecosystems.
A drawback of the unstructured-mesh technology is that, although computer time is saved by using high resolution in chosen areas only, it still uses a substantial amount of computer time as it is less efficient per degree of freedom as compared to structured models.Furthermore, extra care must be taken for models formulated using the continuous finite elements as their local conservation of volume and tracers is formulated in the cluster-weighted sense.This brings some ambiguity into analyzing fluxes between grid cells, while divergences are well defined (Sidorenko et al., 2009).
Before using a newly coupled biogeochemical-ocean model for scientific studies, the skill of the model must be assessed (e.g.Assmann et al., 2010).Performing a skill assessment is not a trivial exercise, considering both the lack of data, especially for parameters such as dissolved iron and export production, and also the inherent uncertainty of the biogeochemical models, in which complex biochemical processes are described by comparably simple mathematical parameterizations.We have coupled the Regulated Ecosystem Model 2 (REcoM2) to the Finite Element Sea Ice-Ocean Model (FESOM), and in this paper a skill assessment of the coupled model is carried out with emphasis on the Southern Ocean.We show to what extent the results are comparable to observations and discuss how they compare to results from other models.

Ocean model
The ocean component of FESOM solves the hydrostatic primitive equations under the Boussinesq approximations (Danilov et al., 2004;Wang et al., 2008).Elastic viscous plastic (EVP) sea ice dynamics is used together with the thermodynamics adopted from Parkinson and Washington (1979) as described in detail by Timmermann et al. (2009).Currently, FESOM is used for simulation of the three-dimensional global ocean with special focus on the Arctic and the Antarctic (Haid and Timmermann, 2013;Wekerle et al., 2013).The latest FESOM version is comprehensively described in Q. Wang et al. (2014).
FESOM operates on unstructured meshes that permit the main feature of the model: the capability of local grid refinement in an otherwise global set-up without nesting.The model domain is discretized by a horizontally triangulated and unstructured, but vertically stratified, mesh with tetrahedral volumes.Integration is carried out on an Arakawa Agrid, which uses vertical z coordinates for simplicity.The mesh used in this study (Fig. 1) is similar to the one used by Sidorenko et al. (2011), in which the horizontal resolution ranges from 15 km in the polar regions to 180 km in the subtropical gyres.In the vertical it has 32 layers, nine of which are located in the upper 100 m.
The bottom topography of FESOM's grid is constructed using a combination of different data products; the bulk of FESOMS bathymetry, from 60 • S to 64 • N, is based on topography data from the General Bathymetric Chart of the Oceans (GEBCO, 1 min resolution), south of 60 • S, the bottom topography from Timmermann et al. (2010) with a resolution of 1 min (Rtopo-1) is used and north of 69 • N it is based on data from the International Bathymetric Chart of the Arctic Oceans with 2 km resolution (IBCAO, version 2; Jakobsson et al., 2008).Between 64 and 69 • N, a combination of the GEBCO and IBCAO data sets is used.FE-SOM's bottom topography is created using bilinear interpo-lation, whereupon smoothing is performed to remove gridscale noise.The topography data also defines the coastline using bilinear interpolation from the data to the model's grid points.For a further description of the creation of bottom topography for FESOM, please refer to Q. Wang et al. (2014).
The version of FESOM used here utilizes a linear representation on triangles (in 2-D) and tetrahedrals (in 3-D) for all model variables.The same is true for the biological tracers, which are treated similar to temperature and salinity.The temporal discretization is implicit for sea surface elevation and a second order Taylor-Galerkin method together with the flux-corrected transport (FCT) is used for advection-diffusion equations.The forward and backward Euler methods are used for lateral and vertical diffusivities, respectively, and the Coriolis force is treated with a second order Adams-Bashforth method.
The vertical mixing is calculated using the PP-scheme first described by Pacanowski and Philander (1981) with a background vertical diffusivity of 1×10 −4 m 2 s −1 for momentum and 1 × 10 −5 m 2 s −1 for tracers.Redi diffusion (Redi, 1982) and Gent and McWilliams parameterization of the eddy mixing (Gent et al., 1990) are applied with a critical slope of 0.004.
The skill of FESOM has been assessed within the CORE framework (Griffies et al., 2009;Sidorenko et al., 2011;Downes et al., 2014), where several sea ice-ocean models were forced with the normal year (CORE-I) and interannually varying (CORE-II) atmospheric states (Large andYeager, 2004, 2009) and results compared.In these assessments, the full flexibility of FESOM's unstructured mesh was not utilized, but the results from FESOM were still within the spread of the other models, and it was consequently concluded that FESOM is capable of simulating the large-scale ocean circulation to a satisfactory degree.

Biogeochemical model
The Regulated Ecosystem Model 2 (REcoM2) belongs to the class of so-called quota models (Geider et al., 1996(Geider et al., , 1998)), in which the internal stoichiometry of the phytoplankton cells varies depending on light, temperature and nutrient conditions.Uptake of macronutrients is controlled by internal concentrations as well as the external nutrient concentrations, and the growth depends only on the internal nutrient concentrations (Droop, 1983).Iron uptake is controlled by Michaelis-Menten kinetics.An overview of the compartments and fluxes in REcoM2 can be seen in Fig. 2.
The model simulates the carbon cycle, including calcium carbonate as well as the nutrient elements nitrogen, silicon and iron.It has two classes of phytoplankton: nanophytoplankton and diatoms, and additionally describes zooplankton and detritus.The model's carbon chemistry follows the guidelines provided by the Ocean Carbon Model Intercomparison Project (Orr et al., 1999), and the air-sea flux cal- culations for CO 2 are performed using the parameterizations suggested by Wanninkhof (1992).
We do not add external sources to the macronutrient pools since the timescale of the runs is short compared to the residence time of the macronutrients in the ocean (Broecker et al., 1982).
Iron has a much shorter residence time (Moore and Braucher, 2008) and is strongly controlled by external sources as well as scavenging.Dissolved iron is taken up and remineralized by phytoplankton, it reacts with ligands and it is scavenged by detritus in the water column (Parekh et al., 2005).New iron is supplied to the ocean by dust and sedimentary input.For dust input, REcoM2 uses monthly averages (Mahowald et al., 2003;Luo et al., 2003), which have been modified to fit better to the observations from Wagener et al. (2008) (N.M. Mahowald, personal communication, 2011).The model assumes that 3.5 % of the dust field consists of iron and that 1.5 % of this iron dissolves when deposited in the surface ocean.This gives a total aeolian input of 2.65 × 10 9 mol DFe yr −1 (DFe -dissolved iron) to the ocean on average.A flux of iron from the sediment has been added accounting for an input of 2.67 × 10 8 mol DFe yr −1 on average.It is incorporated following Elrod et al. (2004) with the magnitude of the iron concentration released by the sediment being dependent on the rate of carbon remineralization in the sediment.
The model has 1 zooplankton class, which is the model's highest trophic level.Grazing is calculated by a sigmoidal Holling type 3 model with fixed preferences on both phytoplankton classes (Gentleman et al., 2003).
The sinking speed of detritus increases with depth, from 20 m day −1 at the surface, to 192 m day −1 at 6000 m depth (Kriest and Oschlies, 2008).Sinking detritus is subject to remineralization.
REcoM2 has sediment compartments for nitrogen, silicon, carbon and calcium carbonate, which consist of one layer into which the detritus sinks when reaching the lower-most ocean layer.Remineralization of the sunken material subsequently occurs in the benthos, and the nutrients are returned to the water column in dissolved state.
REcoM1 and 2 have previously been used for large-scale simulations with focus on the Southern Ocean in set-ups with the MITgcm (MIT general circulation model) (Hohn, 2009;Taylor et al., 2013;Hauck et al., 2013), and the purpose of the current coupling between REcoM2 and FESOM is likewise studies of the Southern Ocean.
A full description of the model equations can be found in Appendix A along with lists of parameters used in the current run.

Model experiment
We present a numerical hindcast experiment with a newly coupled biogeochemical-ocean general circulation model.The run was forced using the CORE-II data set, which was developed for the use of coupled sea ice-ocean models and gives an interannually varying forcing for the years 1949 to 2008 (Large and Yeager, 2009).As focus here is on evaluating the biological surface processes of a newly coupled model, we follow the example of Vichi and Masina (2009) and Yool et al. (2011) and let the coupled model run for a total of 38 years, from 1971 to 2008.The first 33 years are considered spin-up and we present the results for the years 2004 to 2008.Prior to activating the biogeochemical module, the ocean model had been spun up for 300 years, which is sufficient to reach a quasi-equilibrium state (Fig. 8 in Sidorenko et al., 2011).The length of the time step used throughout the run was 1800 s.
In REcoM2, the tracers for dissolved inorganic nitrogen (DIN) and dissolved silicon (DSi) were initialized with values from the Levitus World Ocean Atlas climatology of 2005 (Garcia et al., 2006), and the dissolved inorganic carbon (DIC) and total alkalinity (TA) tracers were initialized with contemporary values from the Global Ocean Data Analysis Project (GLODAP) data set (Key et al., 2004).Due to scarcity of observations for DFe, the iron field was initialized with an output from the Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES) model (Aumont et al., 2003), which has been modified south of 45 • S with mean observed profiles from Tagliabue et al. (2012).All other tracers were initialized with arbitrary small values.
We used a constant value for the atmospheric CO 2 during the simulation.Because of the duration of the simulation, the carbon cycle is not in equilibrium at the end of the run, and we do not to focus on this part of the model here.

Data and skill metrics
The focus of this skill assessment is on the key parameters of the physical, chemical and biological surface fields, for which we examine the model behaviour on the global scale and in the ocean regions shown in Fig. 3.We have a special interest in the Southern Ocean and therefore also look further into the production and its drivers there.On the temporal scale we primarily focus on annual climatologies of the modelled fields, but also show the seasonal development for certain parameters.
The performance of the model regarding sea surface temperature (SST), sea surface salinity (SSS), mixed layer depth (MLD), DIN, DSi, chlorophyll a (Chl), NPP and export production (EP) is summarized in Taylor diagrams, one displaying the spatial agreement between the modelled surface multi-year means and the observations, and one also taking the seasonal variations into account.Due to a lack of monthly data, EP is only plotted in the former Taylor plot.The spatial distribution of the modelled and observed surface climatologies of these same fields, except SST and SSS, have additionally been plotted to show the bias which is not captured by the Taylor diagrams.SST and SSS were omitted as they have been evaluated elsewhere (Sidorenko et al., 2011).Global climatologies of dissolved iron concentrations do not exist, making it impossible to evaluate the modelled iron fields in the manner described above.Instead, the modelled surface climatology of the iron concentration is plotted on its own, and mean values of the ocean basins are compared to observation compilations in table form.NPP is the model's end result, and we have therefore additionally plotted the mean seasonal cycle of the NPP for each of the ocean basins defined in Fig. 3, for the model result as well as satellite-based observations.Due to the large effect the modelled mixed layer has on the prediction of the NPP, we have also chosen to illustrate the seasonal cycle of the modelled and observed mean MLD in these same ocean basins.
The Taylor diagrams (Taylor, 2001) show the correlation (r), the normalized root mean square error (RMSE) and the normalized standard deviation (SD) between the model results and the observations.The correlation between the model and the observations show whether the two data sets increase and decrease simultaneously, the SDs tells us about the magnitude of the changes in the data, but not when these changes occur and the centered RMSE reflects differences in the overall patterns of the two fields after the bias has been removed.The perfect fit between model and observations will have a correlation and a SD of 1 and a RMSE of 0.
A full list of the observations used can be seen in Table 1.For the NPP we use the Vertically Generalized Productivity Model (VGPM) product from the ocean productivity web page (http://www.science.oregonstate.edu/ocean.productivity/index.php),which is based on the Sea-viewing Wide Field-of-view Sensor (SeaWIFS, 2012) chlorophyll measurements and the VGPM NPP model (Behrenfeld and Falkowski, 1997).We have downloaded monthly values from the web page, and from these calculated the spatial and seasonal means.The EP fields from Schlitzer (2004), Siegel et al. (2014) and Laws et al. (2000) were provided as climatologies.The field from Laws et al. (2000) and Siegel et al. (2014) are satellite based, whereas the field from Schlitzer (2004) comes from an inverse model.
Satellite-based estimates of chlorophyll a, NPP and EP provide detailed spatial and temporal data, but obtaining  them is not trivial.Remotely sensed global ocean colour values are first converted to chlorophyll a, and under a number of assumptions about, for instance, mixed layer depth, temperature and light, NPP (Behrenfeld and Falkowski, 1997) and finally EP (e.g.Laws et al., 2000;Siegel et al., 2014) can be estimated.Increasing uncertainty is introduced during the process, and the satellite-based estimates are not as such observations, but rather another way of modelling the chlorophyll a, NPP and EP.The spread between the different satellite-based estimates of NPP is large.Carr et al. (2006) showed that estimates of the global NPP differed by a factor 2 between 24 models, with the largest discrepancies occurring in the high-nutrient low-chlorophyll and extreme temperature areas.The SeaWIFS (2012) algorithms have further been shown to significantly underestimate chlorophyll a concentrations in the Southern Ocean (Gregg and Casey, 2004), and one must consequently be aware of this when using satellitebased estimates from the Southern Ocean.The Arctic Ocean is likewise an area in which observations are scarce, and for the seasonal Taylor diagrams the modelled results have consequently been removed in both of the polar regions when comparable observational values did not exist.As is the case in the Southern Ocean, it is especially the satellite-based observation of NPP that is affected by this during the winter months.The missing seasonal data has likewise led us to cut off the northern basins at 70 • N in Fig. 3.
In addition to comparing our results to available observations, we discuss them in relation to those of other biogeochemical models.For Dfe, Chl and NPP we have plotted our bias plots to the scale used by Schneider et al. (2008).They present the results of the IPSL (Institut Pierre Simon Laplace) model, the MPI (Max Planck Institute) model and the NCAR (National Center for Atmospheric Research) model, providing a range of results to hold our model against.
We additionally compare our model to the results from Assmann et al. (2010), who present the results of another nontraditional ocean model coupled to a biogeochemical model, and to Yool et al. (2011).For the export production, Moore et al. (2004) provided a thorough discussion of the export of particulate organic matter as well as opal, and our spatial plots of export production is consequently plotted to the scale used by them.

Physics: mixed layer depth, salinity and temperature
The fit between the spatial distributions of modelled and observed surface temperature and salinity is very good for both spatial (Fig. 4a) and monthly spatial fields (Fig. 4b), with the correlations being higher than 0.99 and the normalized SDs close to 1 for both fields.As is general practice in oceanonly models (Griffies et al., 2009), FESOM's surface salinity is weakly restored towards the Polar Science Center Hydrographic Climatology (PHC) (Steele et al., 2001) with a piston velocity of 20 m yr −1 .
In both FESOM and the observations (de Boyer Montegut et al., 2004), the mixed layer depth is defined as the depth at which the difference between the potential density at 10 m depth and the MLD is greater than 0.03 kg m −3 .The spatial distribution of the mean MLD has a correlation of 0.68 and a normalized SD of 0.85 when compared to the data-based estimates (Fig. 4).The seasonal variability of the MLD leads to entrainment of water with high nutrient concentrations to the surface water during winter, and the maximum depth of the mixed layer during the year (MLD max ) is therefore especially important from a biological point of view.Overall, the modelled MLD max fits well with the observations, but it is generally too shallow in the Southern Ocean (Fig. 5), with the consequence that limiting nutrients are not adequately replenished during winter.This may lead to a too small NPP in the area as well as a dominance of nanophytoplankton over diatoms, as the former needs a lower iron concentration for production.This will be further discussed in Sect.3.5.
For the monthly fields, the correlation between the modelled MLD and the observations is above 0.6 and the SD equals 1 (Fig. 4b).We investigate this further by plotting the mean depth of the mixed layer in different ocean regions defined in Fig. 3.All basins have correlations above 0.9, except the northern Indian and equatorial basins (Fig. 6), leading us to conclude that the seasonal change in the MLD is well predicted by FESOM.

Nutrients and nutrient limitation
The annual mean surface distribution of DIN and DSi have correlations between model results and observations of 0.91 and 0.86, respectively (Fig. 4a).In the Southern Ocean, the surface DIN concentrations have a negative bias for DIN (Fig. 7) and a positive for DSi (Fig. 8) when the spatial distribution of modelled and observed values are compared.The DIN concentration additionally becomes too high in the sub- tropical Pacific Ocean (Fig. 7), something we will later argue happens due to a strong iron limitation in the area.
The correlation between model results and observations for the spatial-seasonal distribution of DIN and DSi is close to 0.75 for both fields (Fig. 4b).For both nutrients, the seasonal cycle has the best agreement with the observations in the polar regions (not shown).
Iron has been shown to play a large role as limiting nutrients for phytoplankton in the Southern Ocean, as well as the equatorial and subarctic Pacific (Martin et al., 1991), and is therefore a key parameter in the model.We compare the model's surface iron concentration to compilations of observations (Moore and Braucher, 2008;Tagliabue et al., 2012) and to other biogeochemical models (i.e.Schneider et al., 2008).It must be mentioned here that the model is not independent of the observations from Tagliabue et al. (2012) as they are also used for initialization of dissolved iron.But as we only compare surface values, and the residence time of iron in the Southern Ocean is much shorter than the model run, the surface iron concentrations at the end of the model run should not be affected by the initialized values.
The pattern of the surface iron concentration in the Atlantic Ocean (Fig. 9) fits well with the observations in Table 2 as well as with the results from the MPI and NCAR models in Schneider et al. (2008), with relative high concentrations in the equatorial region fed by the dust plume from the Sahara, and concentrations decreasing towards the poles.The   (2010) are somewhat lower in the equatorial and southern part of the Atlantic Ocean than our result.
In the Indian Ocean, our surface iron concentrations agree well with the IPSL and NCAR models as well as with the results from Yool et al. (2011), with values higher than 1 nM in the Arabian Sea and falling towards 0.3 nM in the main In- dian Ocean.Our values also fit well with the observations in the northern Indian Ocean (Table 2), and this, along with the agreement between the models using varying magnitudes of sedimentary iron input, indicates that the coastal upwelling in the Arabian Sea is well captured in these models, and that this upwelling is responsible for the high surface iron concentrations in the area.The lower surface iron concentration in South East Asia is on the other hand evident in all of these models with the exception of the IPSL model, indicating that the sediment source plays a larger role in this area.Here, unfortunately we do not have observations to validate the models.In the Pacific Ocean, our result is closest to the one from Assmann et al. ( 2010), though they have a higher iron concentration along the North and South American west coast, indicating a stronger coastal upwelling in their ocean model.We have a lower surface iron concentration than all models presented by Schneider et al. (2008), even though they all have low concentrations locally.The observations in Table 2 indicate that all models underestimate the surface iron concentrations in the Pacific, especially in the equatorial region where the upwelling plays the largest  role.In the Southern Ocean, Table 2 shows that our surface iron concentration is too low, but the spatial distribution of the surface iron fits well with observed values, with the highest values found in the vicinity of the Antarctic and east of the Patagonian shelf (Tagliabue et al., 2012).This distribution can mainly be attributed to the sediment and dust sources of iron and the seasonal ice coverage impeding iron uptake by phytoplankton near the Antarctic continent.These factors are also responsible for maintaining the relatively high surface iron concentration in the Arctic, which becomes iron limited in the absence of the sediment source of iron.Assmann et al. (2010) and Moore and Braucher (2008) also experienced this, with the latter mentioning that the missing sediment source has a modest impact on productivity and iron concentrations away from the Arctic.In the Southern Ocean, our result agrees the best with the IPSL model, with a relatively high coastal concentration, which then falls towards the north.Both the MPI and NCAR models have relatively constant higher values in the area south of 40 • S of around 0.5 and 0.3 nM, respectively.
The comparison to observations is, however, hindered by the different definitions of the ocean basins.The value of the North Atlantic from Moore and Braucher (2008), for instance, roughly covers the North Atlantic as well as the northern central Atlantic of our definition (Fig. 3).For the equatorial Pacific, Moore and Braucher (2008)  To plot the distribution of the mean surface limitation we follow the example of Schneider et al. (2008), where the nutrient with the lowest MM in a given place is seen as limiting and it is assumed that other factors, such as temperature and light, are limiting when all Michaelis-Menten coefficients are above 0.7 (Fig. 10).
When looking at the yearly mean, iron limits nutrient uptake for both nanophytoplankton and diatoms up to 45 • S and in most of the Pacific.Nanophytoplankton is mainly nitrogen limited in the Atlantic and Indian oceans, concurring with the result by Assmann et al. (2010), Yool et al. (2011) and the IPSL model in Schneider et al. (2008).For diatoms, silicon is limited in the Atlantic and Indian oceans as well as in the Arctic, a feature that we only share with Yool et al. (2011).
In the high latitudes, the modelled nanophytoplankton become light limited during the respective winter months.For the Arctic this is most pronounced in February, where the light limitation reaches down to 45 • N in the Atlantic Ocean.For the Southern Ocean, the highest degree of light limitation occurs in August when the area south of 55 • S is affected.The higher nutrient demand by diatoms means that they are co-limited by iron and light during winter (not shown).

Chlorophyll and net primary production
Global NPP sums up to 32.5 Pg C yr −1 in the model (Table 3), which is lower than the satellite-based estimate of 47.3 Pg C yr −1 (SeaWIFS, 2012; Behrenfeld and Falkowski, 1997) and also slightly below the estimate range of 35-70 Pg C yr −1 given by Carr et al. (2006), but higher than the modelled values ranging from 23.7 to 30.7 Pg C yr −1 reported by Schneider et al. (2008).
On a global scale, diatoms account for 25.9 % of all production in the model.In the subtropical gyres, we see close to zero percent of NPP from diatoms, whereas it constitutes close to 100 % in the Arctic Ocean (not shown).
The correlations between the spatial distribution of modelled results and satellite data are 0.75 for both chlorophyll a and NPP (Fig. 4a), with the mean surface chlorophyll a concentration being somewhat overestimated as compared to the satellite-based estimates in the high latitudes, while the equatorial regions have concentrations that are too low (Fig. 11) and the extent of the subtropical gyres is too large.Yool et al. (2011) have a higher equatorial chlorophyll a concentration in the equatorial regions of the Atlantic and Pacific oceans as compared to our model, but their concentration in the Southern Ocean is even higher than ours.Furthermore, when we compare our model to the IPSL model (Schneider et al., 2008), we again see that our equatorial chlorophyll a concentrations are lower, whereas the concentrations in the North Atlantic and Southern Ocean are fairly similar to our result.(Behrenfeld and Falkowski, 1997) EP glo [Pg C yr −1 ] 6.1 5.8-13.0(Dunne et al., 2007) 6 (Siegel et al., 2014) Opal glo [Tmol Si yr −1 ] 74.5 69-185 (Dunne et al., 2007) NPP SO [Pg C yr −1 ] 3.1 1.1-4.9(Carr et al., 2006) EP SO [Pg C yr −1 ] 1.1 1 (Schlitzer, 2002;Nevison et al., 2012) Opal The spatial distribution of NPP (Fig. 12), follows the same pattern as chlorophyll a, with low production in the oligotrophic gyres along with a higher production in the temperate regions.Production in the gyres is on the low side compared to the satellite-based estimate (Fig. 13), and as they are known to underestimate production here (Friedrichs et al., 2009), our result is most likely much too low here.An explanation may be that the nanophytoplankton in the model does not represent the smallest phytoplankton classes, such as Prochlorococcus and Synechococcus, which are important in the gyres.Even though adaption of the modelled intracellular N : C ratio is possible, this is not enough to increase production here to the level seen in satellite-based estimates.
The missing coastal primary production along the west coast of Africa and South America (Fig. 12) along with a positive temperature bias in these areas (not shown) indicate that the upwelling is too weak here.FESOM has a coastal resolution of 40 km, which is relatively high, but this resolution only covers a narrow path along the coast (Fig. 1), which may not be sufficient for the upwelling zones to be resolved properly; additionally, the low resolution further out in the subtropical gyres could play a role.Furthermore, the driving force for the upwelling is the coastal winds, and the missing upwelling may partially result from a too low resolution of the atmospheric forcing; moreover, higher resolution allows strong surface winds closer to the coast, thereby increasing the strength of the upwelling (Gent et al., 2010).
Another explanation for the low coastal NPP is the missing riverine input of macronutrients, which at least in the case of silicon plays a role locally in places like Amazonas and the Arctic (Bernard et al., 2011).Yool et al. (2011) deal with the missing riverine nutrient input by restoring the nutrient fields along the coasts.Although Yool et al. (2011) have a larger coastal production in their model (especially along the coast of West Africa), they show that the nutrient restoring only has a small influence on this.
When comparing the mean spatial distribution of NPP with other models, our result is the closest to the NCAR model presented by Schneider et al. (2008), with a relatively high production rate in the North and equatorial Atlantic as well as the Indian Ocean.Moderate production in the area of the polar front in the Southern Ocean is a feature that our model shares with the satellite-based estimate (Fig. 12) and with the NCAR and IPSL models from Schneider et al. (2008).Despite our model's strong iron limitation in the Pacific Ocean, our results fit well with the IPSL model, though we have a smaller production rate in the southern Pacific.
Taking seasonal variations into account, NPP and surface chlorophyll a have correlations of 0.66 and 0.57, respectively, when comparing the model and the satellite-based estimate (Fig. 4b), and the normalized SDs are equal to 1.47 and 1.94 for chlorophyll a and NPP, respectively.Both are on the same order as the values presented by Doney et al. (2009).
The timing of the seasonal cycle of NPP is well captured in the majority of the ocean regions defined in Fig. 3 (Fig. 13).The Pearson's correlation coefficients (R) range from 0.31 in the equatorial Atlantic to 0.93 in the southern central Atlantic (Fig. 13), with significant correlations in eight of the fourteen basins.In general, the modelled seasonal cycle is closest to the satellite-based estimate between 10 and 45 • N and S, where the modelled NPP is low, but the magnitude of the seasonal variations fits well with the satellite-based estimate.This is the case for all basins in the mentioned area, except the northern Indian basin.
The Southern Ocean stands out as it has a modelled NPP of the same magnitude as the satellite-based estimate, but the spring bloom occurs too early here, compared to the satellitebased estimate.This will be further discussed in Sect.3.5.

Export production
The export of organic carbon out of the euphotic zone (export production -EP), is calculated at a reference depth, which in our case is set to the standard 100 m (e.g.Schneider et al., 2008;Doney et al., 2009).Here, we regard EP as the organic matter that sinks due to the effect of gravity, whereas the total EP also entails the vertical movement of POC by advection and diffusion plus a contribution from semi-labile DOC.
The global export production sums up to 6.1 Pg C yr −1 in the model (Table 3), close to the satellite-based estimate of 6 Pg C yr −1 from Siegel et al. (2014).It is also within the range of estimates presented by Dunne et al. (2007), but on the low side and closer to modelled estimates than to estimates based on observations or inverse models.
The modelled EP constitutes 20 % of NPP on a global scale, which is similar to the ratio predicted by Laws et al. (2000).The EP field presented by Laws et al. (2000) is calculated at 100 m depth and is based on satellite observations of ocean colour, whereas the EP field calculated by Schlitzer (2004) is based on an inverse model and is calculated at 133 m.Comparing our model to these fields can consequently be argued to be more of a model-model comparison than a model-observation comparison.
The correlation is 0.28 and 0.48, when comparing the spatial distribution of EP in our model to the fields by Schlitzer (2004) and Laws et al. (2000), respectively (Fig. 4a), indicating that our spatial distribution is closer to the field from Laws et al. (2000), and the normalized SDs are 0.90 and 0.60, respectively.
The EP fields from Schlitzer (2004) and Laws et al. (2000) both have high export along the Equator, in the upwelling regions and along 45 • N and S (Fig. 14a and b).In the Southern Ocean, Schlitzer (2004) has a comparably higher export in the Indian and Pacific sector and in the North Atlantic Schlitzer (2004) has less than Laws et al. (2000).
REcoM2 captures the overall pattern with high EP around 45 • N and S and along the Equator (Fig. 14c and d), and the elevated EP in the North Atlantic is a feature that REcoM2 shares with the field from Laws et al. (2000).Turning to the differences, our EP is lower in the North Atlantic, slightly lower in the gyres and higher south of 45 • S when compared to the field by Laws et al. (2000) (Fig. 14e).Compared to the climatology presented by Schlitzer (2004) (Fig. 14b), our EP is generally lower in the Pacific and in the upwelling regions along West Africa and western North and South America, whereas it is higher in the North and South Atlantic (Fig. 14f).In the Southern Ocean, the differences between  the fields are especially visible in the Indian and Pacific sectors, where Laws et al. (2000) have very low export, Schlitzer has a rather high export and REcoM2's export lies in between the two.Schlitzer argues that the satellites do not capture the deep blooms that occur in the area, thereby explaining the lack of EP in the satellite-based estimate.
Both the spatial distribution and magnitude of EP in our model are very similar to what was found by Moore et al. (2004).
Vertical export of opal is similarly calculated across a reference depth of 100 m.On a global scale we have a total opal export of 74.5 Tmol yr −1 .Previous estimates of global export of opal vary widely (Table 3), and our value is in the lower end of the estimates, as are our global values for NPP and EP.

The Southern Ocean
The coupled model FESOM-REcoM2, as a first step, will be used for studies regarding biogeochemical processes in the Southern Ocean south of 50 • S, and we are therefore especially interested in its performance here.The reasons for this focus are further discussed in Sect. 4.
The model's surface salinity and temperature as well as the nutrient fields are well represented in the spatial domain of the Southern Ocean, with all of them having correlation coefficients above 0.9 when compared to observations (Fig. 15).The chlorophyll a and NPP fields both have somewhat lower correlations, with the correlation for NPP and chlorophyll a being equal to 0.75 and 0.48, respectively.
For the MLD, the correlation between the model results and the observational-based estimate (de Boyer Montegut et al., 2004) is 0.63 in the Southern Ocean (Fig. 15).The MLD max is too shallow in the Indian and Pacific sections of the Southern Ocean, especially in the area of the polar front (Fig. 5), causing this low value.Furthermore, FESOM simulates a too deep MLD max in the convection area of the Weddell Sea associated with deep-water formation.This is a common feature in sea ice-ocean models (e.g.Griffies et al., 2009) and should in itself not have a large impact on the production in the Southern Ocean.
The modelled NPP south of 50 • S sums up to 3.1 Pg C yr −1 (Table 3).Carr et al. (2006) summarize previous studies of NPP based on ocean colour and report an average NPP of 2.6 Pg C yr −1 for the Southern Ocean.They also show that the largest uncertainties in satellite-based estimates regarding NPP are found in the Southern Ocean and that biogeochemical models generally predict higher NPP in the area than satellites.
The model's export production equals 1.1 Pg C yr −1 in the Southern Ocean, close to the 1 Pg C yr −1 found by both Schlitzer (2002) and Nevison et al. (2012), and the EP : NPP ratio equals 36 % in the Southern Ocean, similar to what was found by Nevison et al. (2012).
Considering the spatial distribution of EP in the Southern Ocean (Fig. 14), the model is closer to the estimate from The fraction of the the total biomass comprised by diatoms in the Southern Ocean defers between studies (Alvain et al., 2005;Hirata et al., 2011).In the present study, the diatoms are responsible for 25 % of the NPP south of 50 • S, varying from 0 % in the very iron limited waters of the South Pacific to 100 % in the iron replete regions of the Weddell Sea and on the Patagonian shelf (Fig. 16).Vogt et al. (2013) compared the results of four ecosystem models and showed that the percentage of diatom biomass in the Southern Ocean differed significantly between them, ranging from 20 to 100 %.Our diatom percentage is accordingly within the spread of other models.
Production of the silicon-containing diatoms creates a sinking flux of biogenic silica that sums up to 21.5 Tmol Si yr −1 south of 50 • S in the model (Table 3), which is close to the satellite-based estimate of 25 ± 4 Tmol Si yr −1 calculated south of 45 • S (Dunne et al., 2007), but estimates vary significantly between studies (e.g.Moore et al., 2004;Jin et al., 2006;Holzer et al., 2014).
In REcoM2, the opal export in the Southern Ocean accounts for 29 % of the global opal export (Table 3).This number similarly varies widely between studies, with ours being lower than the 70 % suggested by Jin et al. (2006) and  High export fluxes of biogenic silica (Fig. 17) naturally occur in places with a corresponding high percentage of diatom production (Fig. 16).The largest values are found in the temporary ice zone between 60 and 70 • S, as well as in the area east of Patagonia (Fig. 17), where dust and sediments supply iron to the surface water.A band of relatively high opal export is also present in the polar front in the Atlantic and Indian sectors of the Southern Ocean (Fig. 17).In most of the Southern Ocean, the modelled opal flux falls within a range of 0.4-2.5 mol Si m −2 yr −1 .This is slightly higher than the values given by Moore et al. (2004) and lower than the values of up to 9 mol Si m −2 yr −1 in Jin et al. (2006).
The absence of diatom production in the Pacific sector of the Southern Ocean (Fig. 16), leading to a low opal export in the area (Fig. 17), is notable and can be explained by the pronounced iron limitation of the Pacific, which also reaches into the Southern Ocean and limits production here.*

Control of bloom in the Southern Ocean
We will now examine the roles of MLD and iron concentration in explaining the seasonal variability of NPP.For this purpose we define R 2 * as the temporal coefficient of determination multiplied by the sign of the regression slope.R 2 * is calculated for each spatial point in the domain south of 30 • S between NPP and the MLD (Fig. 18a) and between NPP and DFe (Fig. 18b), using the average iron concentration over the top 100 m of the ocean.The R 2 * values show that the Southern Ocean is roughly divided into two zonal bands; one north of 60 • S, in the area of the polar front (Moore et al., 1999), and one south of 60 • S (Fig. 18).The general picture north of 60 • S is that the concentration of dissolved iron and the mixed layer depth both correlate positively with NPP (Fig. 18), indicating that production in the area mainly is iron controlled, and that production starts when the mixed layer deepens and brings iron and other nutrients to the surface.For the mean seasonal cycle of MLD and NPP north of 60 • S (Fig. 19a), the magnitude of the modelled bloom fits nicely with the one from the satellite-based estimate, but the maximum occurs two months earlier.The mean MLD is well predicted by FESOM in the area (Fig. 19a), but it is consistently shallower than what is observed.This has the effect that the modelled phytoplankton receives a larger light intensity than is the case in the ocean, something that may affect the timing of the bloom.The mean iron concentration in the surface water is highly correlated with the depth of the mixed layer north of 60 • S (Fig. 20a).The phytoplankton concentration starts increasing in July, when the iron concentration is high, and reaches a maximum in October, after which a combination of high grazer concentration and decreasing iron concentrations most likely causes the bloom to decline.Under nutrient and light replete conditions, the increase in biomass is a result of the balance between phytoplankton's maximum growth rate and the grazing (Behrenfeld, 2010;Hashioka et al., 2013).On the one hand, this indicates that the model's timing of the bloom could be changed by a smaller maximum growth rate, something that would change the phyto-plankton dynamics on a global scale.On the other hand, the modelled zooplankton concentration is tightly coupled to the increase in phytoplankton concentration (Fig. 20a), and increasing the maximum grazing rate is another way of keeping the growth in biomass down.As modelled grazers are set to prefer nanophytoplankton over diatoms, this may further increase the diatom percentage in the Southern Ocean (Hashioka et al., 2013).
The NPP and MLD fields are negatively correlated south of 60 • S, whereas correlation between NPP and DFe is close to zero here (Fig. 18).This indicates that light is the main limiting factor in this area and that iron is less important as a controlling factor.The intensity of the incoming light decreases with latitude, and is further decreased or blocked by the presence of sea ice during parts of the year, south of 60 • S. The role of the sea ice for the timing of the spring bloom was highlighted by Taylor et al. (2013), who argued that the sea ice melting induces a shallower and more stable mixed layer, increasing the average light intensity received by the phytoplankton, thereby instigating growth.In our study, the modelled bloom is larger than what is estimated by the satellites south of 60 • S, but the timing fits well (Fig. 19b).The difference can be explained by the aforementioned underestimation of NPP by the satellites.NPP starts increasing when the iron concentration is high and decreasing again when the iron concentration is low and the grazer concentration high (Fig. 20b).
It is worth noticing that the increase in production begins at the correct time in both areas, but that the rate of biomass increase is too high.
The sparse observations make it difficult to assess the validity of the modelled seasonal cycles of iron and zooplankton.Tagliabue et al. (2012) presented a seasonal cycle of DFe from the SR3 transect south of Tasmania.Their results indicate that the highest iron concentrations occur in January and February suggesting that our seasonal change in iron concentration, which peaks in September, is off.Our results, however, fit well with the model result from Hoppema et al. (2003), who also see a peak in the iron concentration in September.

Discussion
We have presented a skill assessment of an initial coupling between the biogeochemical REcoM2 and the Finite Element Sea Ice-Ocean Model (FESOM).FESOM's capability of locally increased resolution has not been fully utilized in the current study, with the smallest distance between neighbouring grid points being only a factor of about 10 smaller than the largest, and as with most biogeochemical models (e.g.Yool et al., 2011) we do not resolve eddies.The model run presented here can thus be regarded as a baseline run, from which future work that further explores the capabilities of the new coupling can proceed.Using the current resolution, the advantage of the new FE-SOM-REcoM2 coupling over the older MITgcm-REcoM2 coupling is not obvious, especially taking FESOM's larger demand for computer time into account.The strength of the new coupling will become clearer when new studies on specialized meshes have been carried out.We are currently looking into the effect of the ocean model on the biogeochemistry of the Southern Ocean by comparing two model runs between which only the ocean model differs.We assess the differences of the supply of iron to the surface mixed layer from below, and how it affects the NPP in the area (Schourup-Kristensen et al., 2014).
In regional models of the Southern Ocean, fixed boundary conditions must be added along the northern boundary where the Southern Ocean is connected to the Atlantic, Indian and Pacific oceans, introducing an extra potential for error in such models.Running the model in a global configuration has the advantage that the boundary conditions can be omitted.
As the dynamics of the iron supply is something we plan to study further in the future, and it is a feature we would like to improve, especially in the Pacific Ocean, we will now discuss how changes to the iron cycle may change the results of the model.
In the Southern Ocean, the spatial distribution of iron in the model is reasonable, but it tends towards low values (Table 2), something that may be explained by the crudely constrained external iron sources in the area.The strength of the sediment source varies widely between models (e.g.Moore and Braucher, 2008;Aumont and Bopp, 2006), and in REcoM2 we have an input of iron from the dust and the sediment source, which is in the smaller end of the range.Increasing the strength of the sediment source would especially impact the iron concentration locally in the Southern Ocean, though it has been shown that iron from the sediments can be carried far from the source region (Lam and Bishop, 2008).The Atlantic sector of the Southern Ocean would have an especially large input due to the presence of the Patagonian shelf and the Antarctic Peninsula (e.g.Lancelot et al., 2009), but it would most likely not change the supply to the pelagic areas in the Indian and Pacific Southern Ocean substantially.For the Southern and Pacific Oceans it would be especially important to further explore the influence of the aeolian and sedimentary iron sources as well as the input from ice in the polar areas.
In the remote parts of the Southern Ocean, the input of iron to the mixed layer from below plays a large role (De Baar et al., 1995;Löscher et al., 1997), andTagliabue et al. (2014) showed that the entrainment of iron during deepening of the mixed layer was especially important.FESOM's MLD max is too shallow in the Southern Ocean, especially in the region of the polar front in the Indian and Pacific sectors (Fig. 5), something likely affecting the degree of iron limitation in these two areas.The tight coupling between iron concentrations, NPP and MLD in the polar frontal area (Fig. 20a), further confirms the importance of entrainment as a supply mechanism of iron.
The lower iron input favors the smaller nanophytoplankton, which have a lower iron half-saturation constant, and thereby a lower requirement for iron.A larger input of iron would probably change the species composition towards more diatoms, but would not necessarily increase primary production (e.g. S. Wang et al., 2014).A higher percentage of diatoms would also possibly decrease the models surface silicon concentration, which tends towards high values in the Southern Ocean (Fig. 8).The effect on the silicon concentration is however complicated by the fact that the model's carbon and silicate cycles are decoupled under iron limitation, leading to a higher silicate uptake when the phytoplankton is iron stressed (Hohn, 2009).
In the Pacific Ocean, the surface iron concentration is very low (Fig. 9), inducing iron limitation (Fig. 10), which leads to a build-up of DIN in the surface water (Fig. 7) and low NPP (Fig. 12).The external input of iron from dust and sediment in the model is marginal in the equatorial and southern part of the Pacific and input from upwelling is consequently important here.FESOM produces a reasonable upwelling of 40 Sv along the Equator (Johnson et al., 2001), whereas upwelling is small along the west coast of South America.We have a low iron flux in the upwelled water along the equatorial Pacific in our model (∼ 10 µmol m −2 yr −1 ) compared to the values suggested by Gordon et al. (1997) and Aumont et al. (2003), who reported 44 and 68 µmol m −2 yr −1 , respectively.Our result is however higher than the 5.1 µmol m −2 yr −1 suggested by Fung et al. (2000).The Fe : C ratio in the upwelled water in the equatorial area is 0.0015 µmol Fe mmol C −1 (6.66 × 10 6 mmol C mmol Fe −1 ), which is significantly lower than the prescribed constant intracellular ratio of 0.005 µmol Fe mmol C −1 (2 × 10 6 mmol C mmol Fe −1 ) in the modelled phytoplankton.The upwelled water consequently contains too little Fe to sustain growth, explaining why biological production is not able to utilize the upwelled DIN.
Allowing the model's phytoplankton to adapt to the conditions in the water with a varying intracellular Fe : C ratio would be a possible way to increase production, as the ratio would then decrease in areas with low iron concentrations.Variable intracellular Fe : C stoichiometry, as found by Sunda and Huntsman (1995) and Wilhelm et al. (2013), is used in other models (e.g.Moore et al., 2002;Aumont and Bopp, 2006).The intracellular Fe : C ratio in diatoms ranged from 0.002 µmol Fe mmol C −2 in the equatorial Pacific to 0.007 µmol Fe mmol C −2 in the subtropical gyres in Moore et al. (2002) and the former value fits remarkably well with our Fe : C ratio in the upwelled water in the Pacific.A lower intracellular Fe : C ratio in our model would lower the intracellular Fe : N ratio and bring it closer to the observed nutrient ratio in the upwelled water, indicating that implementing varying ratios would indeed improve the model's performance in the Pacific.
Other features that could potentially improve the iron cycle are spatially varying solubility of iron in the water, spatially varying ligand concentration and scavenging of iron onto dust particles in the water.The latter is present in the iron cycle used by Moore and Braucher (2008) and would likely counter the relatively high iron concentrations in the equatorial Atlantic and Indian oceans that are present in our model (Fig. 9).

Conclusions
In the current study we show that the newly coupled model FESOM-REcoM2 reproduces the large-scale productivity and surface nutrient patterns, with the main deficiency being the strongly iron limited Pacific Ocean.The totals NPP and EP are within the range of previous estimates, but in the lower end, mainly due to the low productivity in the Pacific.The ratio between EP and NPP is 20 %, similar to the result from Laws et al. (2000).
In the Southern Ocean, the modelled spatial mean fields are likewise reasonable, though the comparison here is hindered by the scarcity of observed data.South of 50 • S, the totals NPP and EP agree well with previous estimates, as does the EP : NPP ratio of 36 %.Production is iron and light limited in the Southern Ocean making the external input of iron important as a controller of production.
On a global scale, the model provides reasonable seasonal variations of the NPP, but the main deficiency in the Southern Ocean is the early onset of the spring bloom in the area between 40 and 60 • S.
Overall, the model results at the present resolution are comparable to those of other non-eddy-resolving biogeochemical models and it is well suited for studies of surface processes in the Southern Ocean on a timescale similar to the one used here.

Appendix A: Equations
In biogeochemical models, the biological state variables are subject to change by the ocean circulation through advection and turbulent mixing as well as by biological processes.Detritus further sinks vertically through the water column due to gravity, and exchange occurs across the surface and bottom boundaries for certain variables.
For a given volume of water, the change in concentration of a given biological state variable C can be expressed as follows: Here, the term −U •∇C represents the change in C due to advection, and U = (u, v, w) denotes the velocity of the water in the x, y and z directions, respectively.For sinking state variables, the speed of vertical sinking (V det = (0, 0, w det )) is added to water's velocity in the advection term.
The turbulent motion is taken into account through the term ∇ • (κ • ∇C) where κ is the diffusivity tensor.
The term SMS(C), where SMS stands for sources minus sinks, represents the changes due to biological processes.This is the term that comprises the main body of biogeochemical models.
Certain state variables are subject to fluxes across the boundaries of the ocean model.For these, the flux between the ocean and the benthos (BenF) is calculated at the bottom of the ocean and the flux between the ocean and the atmosphere (AtmF) is calculated for the surface layer of the ocean.
In the following, the equations that make up the source minus sink code in the biogeochemical model REcoM2 are described.
The quota approach makes it necessary to have more tracers than in a model based on fixed ratios, as we need to know the intracellular concentration of each of the modelled elements.REcoM2 has a total of 21 oceanic state variables (Table A1) and four benthos compartments (Table A2).
The state variables DON, PhyC nano , PhyC dia and DetSi are listed in Table A1.The value of the remineralization rate (ρ N ) is listed in Table A3.The temperature dependency of remineralization (f T ) is calculated in Eq. (A54) and the nitrogen and silicon assimilation rates (V N nano , V N dia and V Si , Table A4)  are calculated in Eqs.(A48) and (A49), respectively.ρ T Si will now be explained.
Silicon remineralization: the temperature dependent remineralization rate of silicon (ρ T Si , Table A4) is calculated following Kamatani (1982) up until a set maximum value: T is the local temperature (Table A4).The remineralization rate (ρ Si ) is listed in Table A3 and the temperature dependency (f T ) is calculated in Eq. (A54).
Input from benthos: the bottom grid point of the water further receives remineralized inorganic matter from the benthos: BenF DIN and BenF DSi (Table A5) denote the fluxes of DIN and DSi into the bottom layer of the ocean.ρ ben N and ρ ben Si (Table A3) are constant remineralization rates.BenthosN and BenthosSi denote the vertically integrated benthos concentration of dissolved nitrogen and silicate, respectively (Table A2).

A1.2 DFe
The intracellular iron concentration is connected to the intracellular carbon concentration through an assumed constant ratio (q Fe : C ; Table A6).Biological uptake and release of iron is likewise connected to uptake and release of carbon.

SMS(DFe
The state variables PhyC nano , PhyC dia , DetC and ZooC are listed in Table A1.The value for the constant Fe : C ratio (q Fe : C ) is listed in Table A6 and the DOC excretion rates from phyto-and zooplankton ( phy C and zoo C ) and the degradation rate for detritus C (ρ DetC ) are listed in Table A3.The phytoplankton respiration (r nano and r dia ) is calculated in Eq. (A45), the photosynthesis (P nano and P dia ) in Eq. (A44), the limitation by intracellular nitrogen (f N : Cmax lim ) is described in Sect.A6.1, and the temperature dependency (f T ) is calculated in Eq. (A54).The respiration by zooplankton (r zoo ) is calculated in Eq. (A46) and the scavenging will now be explained.
Scavenging: the calculation of the scavenging in REcoM2 is based on Parekh et al. (2004), case III.Here, the total concentration of dissolved iron (Fe T ) is divided into iron bound to ligands (Fe L ) and free iron (Fe , Table A4): Iron complexed with organic ligands is protected from scavenging.The total ligand concentration (L T ) can be written: Here L denotes the free ligands.
We assume that the reaction between free iron and free ligand (L + Fe Fe L ) is fast enough to be in equilibrium: By prescribing the value of the conditional stability constant (K Fe L ) as well as the assumed constant total ligand concentration (L T ) and combining Eqs.(A8), (A9) and (A10), we can calculate the concentration of free iron (Fe ).This is then used to calculate the scavenging of Fe , which is assumed to be correlated with the concentration of detritus carbon (Eq.A7).The values for K Fe L and L T are listed in Table A6.
The value for the scavenging rate (κ Fe , Table A6) is an important controller of deep-water iron concentrations.
Iron input from dust: the surface layer of the ocean receives an input of iron from aeolian dust deposition.Dust is assumed to contain 3.5 % iron of which 1.5 % is instantaneously dissolved in the ocean.Sea ice blocks dust, and the dust falling here is lost from the system.
Iron input from the benthos: the release of iron to the bottom layer of water is assumed to be proportional to the release of inorganic carbon.This parameterization is based on the work by Elrod et al. (2004).It is calculated as follows: Here BenF Fe (Table A5) is the flux of iron into the bottom layer of the ocean.ρ ben C (Table A3) is the remineralization rate for the benthos carbon and q Fe : C sed (Table A6) is the iron : carbon ratio for the flux.BenthosC (Table A2) denotes the vertically integrated carbon concentration in the benthos compartment.

SMS
The state variables PhyC nano , PhyC dia , DOC, ZooC and Det-Calc are described in Table A1.Respiration by nanophytoplankton (r nano ), diatoms (r dia ) and zooplankton (r zoo ) is calculated in Eqs.(A45) and (A46) and the photosynthesis terms (P nano and P dia ) in Eq. (A44).
The value of the remineralization rate ρ C is listed in Table A3 and the temperature dependency (f T ) is calculated in Eq. (A54).
The dissolution of calcite from detritus (Diss calc ) is calculated in Eq. (A34), and the value of the calcite production ratio (ψ) is listed in Table A7.ψ denotes the percentage of the nanophytoplankton that are calcifiers, and their PIC : POC ratio.
Atmospheric input: the DIC concentration of the surface grid point is affected by the air-sea flux of CO 2 .It is calculated according to the guidelines provided by the Ocean Carbon Model Intercomparison Project (Orr et al., 1999).In the calculations the surface water CO 2 concentration, alkalinity, temperature and salinity are taken into account.
Input from benthos: the bottom grid point of the water further receives remineralized inorganic carbon from the benthos: BenF DIC (Table A5) denotes the flux of DIC into the bottom layer of the ocean and ρ ben C (Table A3) is a constant remineralization rate.The calcite dissolution rate (Diss calc ) is calculated in Eq. (A34).BenthosC and BenthosCalc (Table A2) denote the vertically integrated carbon and calcium carbonate concentration in the benthos compartment.

A1.4 Total alkalinity
The model's total alkalinity is changed by phytoplankton uptake of nutrients (nitrate and phosphate), precipitation and dissolution of calcium carbonate and remineralization of organic matter (Wolf-Gladrow et al., 2007)  constant P : N ratio of 1 : 16.
The state variables PhyC nano , PhyC dia , DON and DetCalc are described in Table A1.The N assimilation (V N nano and V N dia ) is calculated in Eq. (A48).The remineralization rate (ρ N ) can be found in Table A3 and the temperature dependency (f T ) is calculated in Eq. (A54).
Dissolution of calcium carbonate from detritus adds CO 2− 3 to the water and thereby increases the alkalinity with two moles for each dissolved mole calcium carbonate.Diss calc is calculated in Eq. (A34).The parameter ψ, specifying the calcifying fraction of the nanophytoplankton, is listed in Table A7 and the photosynthesis (P nano ) is calculated in Eq. (A44).
Input from benthos: the alkalinity of the bottom grid point of the water is affected by the remineralization of DIN, and thereby also DIP as well as dissolution of calcite from the benthos: BenF Alk (Table A5) denotes the flux of alkalinity into the bottom layer of the ocean.The dissolution rate (Diss calc ) is calculated in Eq. (A34), and the remineralization rate ρ ben N is listed in Table A3.BenthosN and BenthosC (Table A2) denote the vertically integrated nitrogen and carbon concentration in the benthos compartment.The state variables PhyC nano , PhyN nano , PhyC dia and PhyN dia are described in Table A1.
The nitrogen assimilation (V N nano and V N dia ) is calculated in Eq. ( A48) and the constant excretion rate ( phy N ) is listed in Table A3.When the N : C ratio becomes too high, excretion of DOC is downregulated by the limiter function (f N : Cmax lim ) described in Sect.A6.1.A further loss term is phytoplankton aggregation (Agg), which transfers N to the detritus pools (Eq.A27).
The grazing loss (G nano and G dia ) is calculated in Eqs.(A52) and (A53), respectively.The state variables PhyC nano and PhyC dia are described in ) is listed in Table A3, the ratio q CaCO 3 : N nano = PhyCalc / PhyN nano and the calcite dissolution rate will now be explained.
Calcite dissolution: as the detritus calcite sinks through the water column it is subject to dissolution (Diss calc , Table A4) occurring on a length scale of 3500 m (Yamanaka and Tajika, 1996).

Diss calc =
w z 3500 m (A34) w z denotes the sinking speed at depth z and is calculated in Eq. (A28).
Loss to benthos: when the sinking detritus reaches the bottom grid point it is assumed that it continues sinking into the benthic compartment with the speed w det (Eq.A28).This leads to a detrital flux (BenF DetCalc , Table A5) from the water column to the benthos: The state variable DetCalc is described in Table A1.The state variables are described in Table A1.The values for excretion of nitrogen and carbon from phyto-and zooplankton ( phy N , phy C , zoo N and zoo C ) are listed in Table A3 along with the degradation rates for detritus (ρ detN and ρ detC ) and remineralization rates of DON and DOC (ρ N and ρ C ).The limitation factors (f N : Cmax lim, nano and f N : Cmax lim, dia ) are described in Sect.A6.1 and the temperature dependency (f T ) is calculated in Eq. ( A54).

A2 Sources minus sinks, benthos
The model has a benthos compartment which consists of one layer.Matter is supplied to this layer through sinking detritus, and it hence has pools of nitrogen, silicon, carbon and calcite.When sinking detritus reaches the bottom it continues sinking into the benthos with the speed, w det , calculated by Eq. (A28) and is thus lost from the water column.In the benthos, the detritus is assumed to be remineralized to dissolved inorganic matter.This is then re-released to the water's pools of dissolved inorganic matter (DIN, DIC, Alk and (A41) The state variables are described in Table A1 (DetN, DetSi, DetC and DetCalc) and Table A2 (BenthosN, Ben-thosSi, BenthosC and BenthosCalc).The remineralization rates (ρ ben N , ρ ben Si and ρ ben C ) are listed in Table A3 and the calcite dissolution rate (Diss calc ) is calculated in Eq. (A34).

A3.1 Photosynthesis
The rate of the C-specific photosynthesis is calculated for both nanophytoplankton and diatoms (P nano and P dia , Table A4).
The calculation is based on the work by Geider et al. (1998) and differs between nanophytoplankton and diatoms in the nutrient limitation; nanophytoplankton is limited by iron and nitrogen while diatoms are additionally limited by silicon.
P nano max = µ max C • min f Fe lim, nano , f N : Cmin lim, nano • f T (A42) Nutrient limitation is calculated using the Liebig law of the minimum, in which the most limiting nutrient limits production (O'Neill et al., 1989).The value of µ max C can be found in Table A8.The limitation terms (f N : Cmin lim , f Si : Cmin lim and f Fe lim ) differ somewhat from the formulation in Geider et al. (1998) and are described in Sect.A6.2 and the temperature dependency (f T ) is calculated in Eq. (A54).
The actual C-specific photosynthesis rate depends on how much photosynthetically available radiation (PAR, Table A4) the cell can harvest.This is controlled by the light harvesting efficiency (α) and the intracellular Chl : C ratio as well as the

Figure 2 .
Figure 2. The pathways in the biogeochemical model REcoM2.

Figure 3 .
Figure 3. Map of the ocean regions used to examine the model results on a basin scale.

Figure 4 .
Figure 4. Taylor diagrams(Taylor, 2001) showing correlation, normalized SD and the normalized root mean square error between values of the model results and observations (Table1), weighted by area.(a) Spatial distribution; (b) spatial-seasonal distribution.All values are surface values, except the mixed layer depth and the vertically integrated NPP.(a) Uses the yearly mean calculated over 2004-2008 and (b) uses the monthly means of the same years.All fields have been interpolated to a 1 • × 1 • grid, using linear interpolation.

Figure 6 .
Figure 6.Mean MLD over the year in the ocean basins depicted in Fig. 3.The correlation coefficient is written in each plot, and the statistically significant correlations (p values < 0.05) are marked with * .

Table 2 .
Modelled mean surface iron concentrations (0-100 m) in the different ocean basins shown in Fig. 3. Observed values are from Moore and Braucher (2008), except those marked with * , which are from Tagliabue et al. (2012),

Figure 9 .
Figure 9. Spatial distribution of the mean surface concentration of dissolved iron.Plotted to the scale used by Schneider et al. (2008).
report the value 0.84 nM for the whole ocean, and 0.11 nM for the open ocean, where our value is closer to the latter due to the missing coastal processes in the model.Nutrient uptake limitation is described by Michaelis-Menten kinetics in the model.The Michaelis-Menten coefficient (MM) is computed as MM = [Nut]/([Nut] + K Nut ), with [Nut] being the nutrient concentration, and K Nut a nutrient and phytoplankton dependent half-saturation constant.

Figure 13 .
Figure 13.Mean NPP over the year in the ocean basins depicted in Fig. 3.The correlation coefficient is written in each plot, and the statistically significant correlations (p values < 0.05) are marked with * .

Figure 15 .
Figure 15.Taylor diagram for the Southern Ocean south of 50 • S showing correlation, normalized SD and the normalized root mean square error between the spatial distribution of the model results and observed data sets, weighted by area.All values are surface values, except the mixed layer depth and the vertically integrated NPP.The fields have been linearly interpolated to a 1 • ×1 • grid, similar to the World Ocean Atlas. he model SD and root mean square error have been normalized by the observational SD.

Figure 16 .
Figure 16.Contribution of diatoms to NPP in the model.

Figure 17 .
Figure 17.Modelled export of opal across 100 m depth plotted to the scale used by Moore et al. (2004).

Figure 18 .
Figure 18.Maps showing coefficients of determination for crosscorrelation between model results of (a) NPP and MLD and (b) NPP and DFe.DFe has been averaged over the upper 100 m of the water for the calculation.R 2 * is defined as the temporal coefficient of correlation multiplied by the sign of the regression coefficient.

Figure 20 .
Figure 20.Seasonal change in the modelled NPP, MLD, DFe and zooplankton concentration for (a) the Southern Ocean from 40 to 60 • S and (b) south of 60 • S. DFe and zooplankton concentrations are averaged over the top 100 m of the ocean.NPP and MLD are normalized by the maximum of the monthly values.

Table 1 .
List of the observational data sets used for the skill assessment.

Table 2 .
The latter is the mean of the values given for the Antarctic and Subantarctic regions.
*iron concentrations in the IPSL model and Assmann et al.

Table 3 .
Net primary and export production for the global domain and the Southern Ocean south of 50 • S, for REcoM2 and from literature.

Table A3 .
Degradation parameters for sources minus sinks equations.

Table A5 .
Benthos variables.Si from benthos to bottom water [mmol Si m −2 day −1 ] BenF Fe Flux of Fe from benthos to bottom water [µmol Fe m −2 day −1 ] BenF DetCalc Flux of detritus calcite from the water to the benthos [mmol CaCO 3 m −2 day −1 ] BenF DetC Flux of detritus C from the water to the benthos [mmol C m −2 day −1 ] BenF DetN Flux of detritus N from the water to the benthos [mmol N m −2 day −1 ] BenF DetSi Flux of detritus Si from the water to the benthos [mmol Si m −2 day −1 ]

Table A6 .
Parameters for iron calculations.

Table A7 .
Parameters for sources minus sinks equations.