M EDUSA-2 . 0 : an intermediate complexity biogeochemical model of the marine carbon cycle for climate change and ocean acidification studies

MEDUSA-1.0 (Model of Ecosystem Dynamics, nutrient Utilisation, Sequestration and Acidification) was developed as an "intermediate complexity" plankton ecosystem model to study the biogeochemical response, and especially that of the so-called "biological pump", to anthropogenically driven change in the World Ocean (Yool et al., 2011). The base currency in this model was nitrogen from which fluxes of organic carbon, including export to the deep ocean, were calculated by invoking fixed C:N ratios in phytoplankton, zooplankton and detritus. However, due to anthropogenic activity, the atmospheric concentration of carbon dioxide (CO2) has significantly increased above its natural, inter-glacial background. As such, simulating and predicting the carbon cycle in the ocean in its entirety, including ventilation of CO2 with the atmosphere and the resulting impact of ocean acidification on marine ecosystems, requires that both organic and inorganic carbon be afforded a more complete representation in the model specification. Here, we introduce MEDUSA-2.0, an expanded successor model which includes additional state variables for dissolved inorganic carbon, alkalinity, dissolved oxygen and detritus carbon (permitting variable C:N in exported organic matter), as well as a simple benthic formulation and extended parameterizations of phytoplankton growth, calcification and detritus remineralisation. A full description of MEDUSA-2.0, including its additional functionality, is provided and a multi-decadal spin-up simulation (1860–2005) is performed. The biogeochemical performance of the model is evaluated using a diverse range of observational data, and MEDUSA-2.0 is assessed relative to comparable models using output from the Coupled Model Intercomparison Project (CMIP5).


Introduction
Since the beginning of the industrial era, the atmospheric concentration of carbon dioxide (CO 2 ) has significantly increased above its natural, inter-glacial background concentration.Further increases are predicted by climate models, e.g. to 450-650 ppm by the mid-21st century (Houghton et al., 2001).Rising atmospheric CO 2 is mitigated by uptake on land and in the ocean, with the latter accounting for about 30 % of anthropogenic emissions (Sabine et al., 2004).This uptake by the ocean is driven by what are known as the solubility and biological pumps, the former via dissolution of CO 2 in cold waters that are mixed to depth, and the latter as the sinking and downward mixing of organic matter into the ocean interior (Volk and Hoffert, 1985).Global warming will likely cause significant changes in ocean circulation, ecosystems and carbon export (Doney et al., 2012).Some analyses of phytoplankton suggest that change is already detectable and that abundance is declining in response to rising sea surface temperatures (Boyce et al., 2010).Modelling studies have similarly indicated that increased stratification in response to future CO 2 emission scenarios leads to decreased primary production and associated export of carbon (e.g.Bopp et al., 2001;Steinacher et al., 2010; but see Taucher and Oschlies, 2011).
The potential of the ocean to take up CO 2 from the atmosphere is vast because CO 2 is buffered by the carbonate chemistry of seawater, keeping concentrations low relative to other components (HCO − 3 and CO 2− 3 ).Ocean acidification is a further consequence of the chemical equilibrium in seawater because, as anthropogenic CO 2 invades, it combines with H 2 O to form HCO − 3 and H + .Model hindcasts indicate that surface ocean pH has declined from its preindustrial value of 8.2-8.1 today, an increase in acidity of 30 % (Orr et al., 2005).Forward predictions indicate substantial further decreases, e.g.0.3-0.4pH units, by 2050 depending on future CO 2 emissions (Orr et al., 2005).The chemical impact of ocean acidification has the potential to affect ocean ecosystems and associated biogeochemistry in many ways (Doney et al., 2009).In particular, it leads to decreasing saturation state for the two main forms of calcium carbonate (CaCO 3 ) produced by marine calcifiers, aragonite and calcite.Coccolithophores, foraminiferans and pteropods are thus particularly vulnerable to such changes (Fabry et al., 2008;Gangstø et al., 2011).Acidification and decreasing CaCO 3 production have several consequences for the ocean carbon cycle.Production of CaCO 3 removes twice as much alkalinity as it does CO 2 from seawater (Frankignoulle et al., 1994) such that decreasing CaCO 3 leads to elevated pCO 2 and a negative feedback with the atmosphere.On the other hand, the rain ratio, i.e. the ratio of CaCO 3 : POC in sinking particulate organic carbon (Archer, 1991) will decline and with it carbon export flux to the deep ocean.Furthermore, if the export of organic carbon is closely bound by ballasting minerals including carbonate (Armstrong et al., 2002;Klaas and Archer, 2002), a decrease in CaCO 3 production could lead to a substantial shallowing of the depth scale of remineralisation (Heinze, 2004).Previously, we introduced an "intermediate complexity" plankton model, MEDUSA-1.0:Model of Ecosystem Dynamics, nutrient Utilisation, Sequestration and Acidification (Yool et al., 2011).This model expanded beyond the traditional nutrient-phytoplankton-zooplanktondetritus (NPZD) models by having multiple currencies (N, Si and Fe) and by separating plankton into "small" and "large" size classes, yet incorporated sufficiently few tracers to be readily tractable in global ocean general circulation models.A multi-decadal spin-up simulation was undertaken and results presented for global nutrient fields, primary production, distributions of phytoplankton types and export of detritus.Here, we introduce MEDUSA-2.0,an expanded successor model which represents dissolved inorganic carbon (DIC) and pCO 2 in the ocean, thereby allowing the calculation of air-sea CO 2 fluxes as well as an explicit representation of ocean acidification and its impact on ecosystem processes.The new model includes additional state variables for dissolved inorganic carbon, alkalinity, dissolved oxygen and detritus carbon (permitting variable C : N in exported organic matter), as well as a simple benthic formulation and extended parameterizations of phytoplankton growth and detritus remineralisation.A full description of the additional functionality of MEDUSA-2.0 is provided.A multi-decadal spin-up simulation is described , and this is used to provide a means of evaluating the performance of MEDUSA-2.0.

State variables
MEDUSA-1.0 resolves 11 state variables distributed between the nitrogen (6), silicon (2) and iron (1) cycles.The remaining 2 state variables denote chlorophyll for each of the model's 2 phytoplankton classes.Because of its key role in organising marine productivity, nitrogen is MEDUSA-1.0'sprimary currency.In this framework, the cycling of carbon (and other elements) can only be estimated from the explicitly modelled elemental cycles, and then only if fixed stoichiometric relationships are assumed.
In order to incorporate the carbon and oxygen cycles, MEDUSA-2.0adds a further 4 state variables to the existing framework.These include total dissolved inorganic carbon (DIC), total alkalinity (TA) and dissolved oxygen.The final additional state variable is detrital carbon for the slowsinking component of non-living particulate organic carbon (POC).For simplicity, MEDUSA-2.0retains MEDUSA-1.0'sassumption of fixed C : N ratios for the plankton pools (phytoplankton, zooplankton), but since these pools do not have identical C : N ratios (e.g.zooplankton are assumed to have a lower ratio; Anderson, 2005) the flow of organic material to detrital pools, both slow-and fast-sinking, has a variable C : N ratio depending upon which processes (plankton mortality, zooplankton egestion) contribute to it.In the case of fast-sinking detritus, this is still handled implicitly within MEDUSA-2.0,so can be easily accommodated.Since slowsinking detritus is already represented by an explicit nitrogen state variable, a corresponding carbon variable must be added to accommodate this.Note that, again for simplicity, iron is still coupled rigidly to nitrogen, so there is no corresponding state variable for detrital iron.Figure 1 presents a schematic diagram of MEDUSA-2.0,showing the state variables (pelagic and benthic) and the ecological connections between them.
The full list of 3-D water column state variables for MEDUSA-2.0 is as follows: The smaller boxes at the bottom of the diagram refer to benthic reservoirs of model currencies that are fed by sinking detrital material (slow-and fast-sinking).For reasons of diagrammatic clarity, dissolved oxygen and its connections to other state variables are omitted here.Note that the dissolution of benthic CaCO 3 releases both DIC and alkalinity.
In addition to the state variables for the 3-D water column, 4 further state variables have been added to represent 2-D pools of organic and biogenic material at the seafloor.These pools permit temporary storage of particulate material before it is returned to dissolved pools, and they represent an extremely crude submodel of the benthic ecosystem.This approach contrasts with that in MEDUSA-1.0 in which all particulate material reaching the seafloor is instantaneously remineralised (or dissolved).The primary motivation for this addition is to prevent the unrealistic acceleration of nutrient regeneration at the seafloor caused by such simplified model assumptions.This is particularly an issue in the shelf regions of the World Ocean where shallower water columns and strong vertical mixing can quickly return regenerated nutrients to surface waters and unrealistically fuel extra productivity.The simplicity of the submodel used in MEDUSA-2.0means that it does not resolve the complex interplay between benthic ecosystems and seafloor sediments, but instead serves as a precursor to the inclusion of a more sophisticated treatment (e.g.Rowe and Deming, 2011).As in the case of the detritus (slow-and fast-sinking) that fuels these seafloor pools, iron is rigidly coupled to nitrogen and does not have a separate benthic state variable.In principle, it could alternatively be coupled to carbon, but for parity with MEDUSA-1.0,Similarly to MEDUSA-1.0, the oceanic inventories of nitrogen and silicon are fixed, and biogeochemical processes effectively only move these elements between modelled pools.Processes which act to add or remove these elements to or from the ocean (whether "abiotic" such as rivers or burial, or "biotic" such as nitrogen fixation or denitrification) are ignored here.This approach is adopted partly because these unmodelled processes are assumed to be of limited magnitude relative to modelled processes, in part because they are less well-understood and more difficult to model, and partly to simply limit model complexity.(As an aside, note that the residence time of these elements within the ocean is considerably longer than the duration of all simulations of NEMO-MEDUSA-2.0 to date.)In MEDUSA-2.0,this same fixed inventory also applies to alkalinity, but the remaining elemental cycles -iron, carbon and oxygen -have connections to reservoirs external to the ocean.In the case of A. Yool et al.: A description of MEDUSA-2.0iron, aeolian deposition and dissolution of benthic sediments supply this element to seawater, while scavenging actively removes it.Meanwhile, both carbon and oxygen are actively exchanged with the atmosphere at the ocean's surface.However, oxygen also occupies an unusual station in that -as far as the modelled inventory is concerned -within the water column it is both generated from "nothing" by primary production, and dissipated to "nothing" by respiration.
The inclusion of the cycles of carbon, alkalinity and oxygen introduces a number of features to MEDUSA-2.0 that are relevant for studies of future climate change or ocean acidification.These include: -gas exchange of dissolved CO 2 and O 2 with the atmosphere, -a carbonate chemistry module for calculating properties such as the concentrations of carbonate species (H 2 CO 3 , HCO − 3 , CO 2− 3 ), pCO 2 and pH, -a dynamic lysocline depth calculated from the 3-D saturation state of calcium carbonate (specifically the calcite polymorph).
Alongside these major additions, MEDUSA-2.0has a number of less significant differences from MEDUSA-1.0 that relate to aspects such as parameterization and forcing.These differences include: -forcing field of aeolian iron deposition replaced with that of Mahowald (2005), -parameterization of seafloor supply of dissolved iron added, -phytoplankton growth parameterization extended to include option of Liebig "law of the minimum" functionality, - Martin et al. (1987) and Henson et al. (2011) parameterizations of the remineralisation of fastsinking detritus optionally available, -options to use either fixed or dynamic rain ratios and lysocline depths.
A separate development with bearing on the work described here is the utilisation of surface forcing derived from coupled ocean-atmosphere models.This supplants the observationally derived reanalysis forcing (DFS4.1;DRAKKAR Group, 2007) used previously with MEDUSA-1.0(Yool et al., 2011).As well as permitting forecast simulations, adoption of such model-derived forcing permits spin-up or hindcast simulations of the pre-industrial past prior to the ongoing anthropogenic transient.The specific forcing used here is described in Sect.3.1.

Differential equations
The following partial differential equations describe the biogeochemical tendency terms that operate on MEDUSA-2.0'sstate variables.Abbreviations used in the bracketed descriptions are: "PP" for primary production; "µzoo" for microzooplankton; "mzoo" for mesozooplankton; "non-lin" for non-linear; "remin" for remineralisation of organic material; "diss" for dissolution of inorganic material (e.g.opal or CaCO 3 ).The functional forms and parameters used in these equations are expanded upon in Sects.2.3 and 2.4.
non−lin losses (1) non−lin losses (2) fast Si detritus diss air−sea gas exchange ( 13) air−sea gas exchange (15) The above equations are applied throughout the domain of the physical ocean model, without regard to horizontal or vertical position.This approach is inherited from MEDUSA-1.0 but differs from that of some other models (Popova et al., 2006) where different equations are applied in different volumes of the ocean to account, for instance, for photic and aphotic zones.Note that terms such as air-sea gas exchange, aeolian dust deposition and fluxes from the benthic submodel (see below) obviously only apply in ocean grid cells that are in contact with either the atmosphere or the benthos.
Differential equations ( 16)-( 19) describe the storage and release of biogenic material at the base of each water column in the model.Material enters these reservoirs as slow-and fast-sinking detritus, and is remineralised to DIN, iron, silicic acid, DIC and alkalinity.As with the rest of MEDUSA-2.0,iron is coupled via fixed stoichiometry to the nitrogen cycle and so is handled implicitly.Note that there is no horizontal communication between the benthic reservoirs in MEDUSA-2.0.Since release of material from the benthic reservoirs occurs at fixed specific rates, the above equations are complete.

Interaction functional forms
The following sections expand on the terms that appear above in MEDUSA-2.0'sdifferential equations.Although MEDUSA-2.0includes a number of new state variables as well as several additional biogeochemical processes, it largely overlaps MEDUSA-1.0with regard to the form and parameterization of shared processes.As such, and since this manuscript aims to provide a complete and standalone description of MEDUSA-2.0,there is repetition with the previously published description of MEDUSA-1.0.Parameter definitions and values are listed in Sect.2.4.

Non-diatom limitation and growth
θ Chl Pn is the scaled chlorophyll to biomass ratio, while αPn scales the initial slope of the photosynthesis-irradiance (P -I ) curve, α Pn , by this ratio so that phytoplankton with a high chlorophyll content have an elevated response to irradiance.
This term calculates maximum phytoplankton growth rate as an exponential function of temperature, T , and base growth rate at 0 • C (Eppley, 1972).
Given the (chlorophyll-related) initial slope of the P -I curve and (temperature-related) maximum phytoplankton growth rate, this function calculates realised growth rate given local irradiance, I (W m −2 ).
Nutrient limitation of phytoplankton growth is specified here via standard, hyperbolic Michaelis-Menten terms that use ambient nutrient concentrations and parameters for the concentration at which phytoplankton growth is half its theoretical maximum.
Light-and nutrient-limitation factors are brought together in a multiplicative term that determines nutrient uptake and, via Redfield coupling, primary production.Yool et al. (2011) investigated the significance of an alternative Liebig law of the minimum scheme for multiple nutrient limitation, and use of this approach is permitted in MEDUSA-2.0via a switch, jliebig.
MEDUSA-2.0 uses a light attenuation submodel derived from the simpler LOBSTER model (Levy et al., 2001).This splits photosynthetically available radiation (PAR) into two wavebands ("red" and "green-blue") that are attenuated separately by seawater and by phytoplankton chlorophyll (i.e.not biomass) from both modelled groups.As such, the model includes self-shading by phytoplankton within the water column.Irradiance above, I , is the sum of these two components of PAR.

Diatom limitation and growth
Diatom phytoplankton growth terms are identical to those of non-diatom phytoplankton.However, because of their obligate requirement for silicon, diatom growth is additionally coupled to the availability of this nutrient, and a submodel of silicon uptake and diatom growth is used to represent these processes (Mongin et al., 2006).This places constraints on growth and nutrient uptake based upon the Si : N ratio of the modelled diatom cells, but allows a degree of plasticity in this ratio depending upon ambient growth conditions.
As noted above, the growth of diatom phytoplankton is additionally limited by the availability of the macronutrient silicic acid.
Silicon is largely used by diatom phytoplankton in the construction of their cell walls, or frustules, which can vary significantly in their ornamentation (e.g.spines, girdle bands; Martin-Jézéquel et al., 2000) depending upon silicon availability.As a result, model diatoms have a degree of plasticity in their requirement for silicon, necessitating a separate state variable, Pd Si , and centred around the stoichiometric ratios, In the above equations, U ∞ is the hypothetical growth ratio at an infinite ambient Si : N ratio, and the uptake of nitrogen (and iron) by diatom cells, PP Pd , is governed by the Si : N ratio.If this falls below a critical value, R 0 Si : N , diatom cells are unable to complete their cell division cycle and growth stops (Martin-Jézéquel et al., 2000).At values above this minimum ratio growth is scaled by a factor of the Si : N ratio, and above 3 times this ratio, growth in diatom biomass is unimpeded.
Silicon uptake, PP Pd Si , occurs at the maximum rate permitted by light and silicon availability whenever the Si : N ratio is below a critical threshold, (3 • R 0 Si : N ) −1 .Above this ratio, silicon uptake is linearly decreased to another threshold value, (R 0 Si : N ) −1 , above which no silicon is taken up by diatom cells -though diatom biomass, Pd, can still increase (and, of course, alter the Si : N ratio).Figure 2 illustrates these equations by showing uptake of nitrogen and silicon by diatoms across a range of biomass Si : N ratios.

Chlorophyll growth scaling factors
As noted already, both phytoplankton groups have separate chlorophyll state variables in addition to those of nitrogen biomass.These allow modelled phytoplankton to alter their chlorophyll content dynamically under different light regimes (e.g. in response to season and depth).The following terms for these processes are taken from Taylor et al. (1997).

Microzooplankton grazing
As part of the size-structuring of MEDUSA, microzooplankton graze on smaller non-diatom phytoplankton and on particles of slow-sinking detritus.The ingestion function that balances the availability of these prey items with the preference microzooplankton have for them is drawn from the classic model of Fasham et al. (1990).
In the above, X is Pn or D.
The above term is repeated for each separate prey item consumed by microzooplankton.The term is based around a sigmoid function in which the "substrate" is composed of the sum of the prey items scaled by the preference that microzooplankton have for them.It is assumed here that microzooplankton prefer non-diatom phytoplankton over detritus since they represent a higher quality food item.
Here, the separate quantities of nitrogen, IN Zµ , and carbon, IC Zµ , ingested by microzooplankton are summed.Parameter φ relates to grazing inefficiency, so-called "messy feeding", that returns a fraction of the grazed material back to dissolved nutrient.For the material actually ingested, the resulting C : N ratio, θ Fµ , can be calculated.
Since grazed material may have a different C : N ratio than that required for microzooplankton growth, the assimilation and metabolism submodel of Anderson and Pondaven (2003) is incorporated here to balance growth, excretion and respiration.The C : N ratio of ingested food calculated above is then compared to the ideal ratio preferred by microzooplankton, θ * Fµ .This makes use of the C : N ratio of microzooplankton biomass, θ Zµ , the assimilation efficiencies of nitrogen, β N , and carbon, β C , as well as the carbon growth efficiency, k C , of microzooplankton.Unlike in MEDUSA-1.0,where an implicit treatment of carbon required all C : N ratios to be identical, here θ Zµ adopts a lower value more consistent with that of zooplankton.
Either C or N limits production depending on whether θ Fµ is greater or lower than θ * Fµ , with any excess carbon respired, and any excess nitrogen excreted.Respiration, R Zµ , growth, F Zµ , and excretion, E Zµ , are calculated as follows.
Fµ then N is limiting and Fµ then C is limiting and Figure 3 of Yool et al. (2011) shows the relative partitioning of carbon and nitrogen grazed by zooplankton depending upon food C : N ratio.In MEDUSA-1.0, the flux of C produced by zooplankton respiration was simply diagnostic, since the biogeochemical cycle of C was not resolved.Here, the loss of C through respiration is explicitly balanced by an increase in DIC in Eq. (13).

Mesozooplankton grazing
Mesozooplankton grazing follows that of microzooplankton with the exception that mesozooplankton have a broader range of prey items: non-diatoms, diatoms, microzooplankton and slow-sinking detritus.Because of this longer list of prey items, Eq. ( 55) below is used to simplify the presentation of mesozooplankton grazing.Note that, though mesozooplankton do not utilise grazed silicon from diatoms, Eq. ( 56) is included below to account for the grazing-induced loss of diatom silicon.For simplicity, parameters φ, β N , β C , and k C are identical to those used for microzooplankton.
where X is Pn, Pd, Zµ or D.
Fm then N is limiting and Fm then C is limiting and Note that grazing by both types of zooplankton in MEDUSA-2.0 is not a function of temperature, in contrast with a number of other studies (e.g.Schartau and Oschlies, 2003;Chen et al., 2012).This decision largely reflects the source of the grazing submodel used here, (Anderson and Pondaven, 2003), as well as the likely computational cost of recalibrating this submodel to include temperaturedependence.Nonetheless, model studies such as that of Taucher and Oschlies (2011) illustrate that the response of plankton ecosystems to future climate warming can be more complex than that typically simulated (e.g.Steinacher et al., 2010).As such, this represents an important aspect for the future development of MEDUSA-2.0.

Plankton loss terms
In addition to losses due to grazing, all four living components of the plankton model incur smaller, secondary losses due to other processes.
The above functions are density-independent loss terms for processes such as metabolism that occur without reference to abundance (i.e. the absolute loss scales linearly with abundance).
The above functions are density-dependent loss terms for processes that occur at rates that depend upon plankton abundance (i.e. the absolute loss disproportionately increases with abundance).These include those such as disease (e.g.viruses), intra-trophic trophic "cannibalism" and predation by implicit, higher trophic level actors.By default, densitydependent losses are represented using a hyperbolic function of plankton concentration (Fasham, 1993), although switches in the model code (Table 6) permit linear, quadratic and sigmoid functions.The best choice of a form for a mortality function is unclear, but can have significant consequences for models (e.g.Steele and Henderson, 1992;Edwards and Yool, 2000;Anderson et al., 2010).As such, Yool et al. (2011) investigated alternative functions for this mortality term.While the simplest form examined -linear mortality -had significant (and unrealistic) impacts on the behaviour of MEDUSA-1.0, the differences between simulations using quadratic, hyperbolic (as here) and sigmoid forms was much more minor, and MEDUSA-2.0retains the same default as MEDUSA-1.0.

Miscellaneous losses
As silicic acid occurs at undersaturated concentrations throughout the modern ocean (Yool and Tyrrell, 2003), the silicon component of diatom phytoplankton is additionally vulnerable to dissolution.This is represented here by a simple linear loss rate (Mongin et al., 2006).
Remineralisation of slow-sinking detrital particles to dissolved inorganic pools occurs at rates dependent on ambient temperature.

Iron supply and removal
Following the submodel of Dutkiewicz et al. (2005), iron is added to the ocean by aeolian deposition of iron-carrying dust at the surface, and removed throughout its volume by scavenging.
F atmos = spatially variable rate (79) The field of iron deposition used in MEDUSA-1.0has been updated for MEDUSA-2.0 to take advantage of a newer climatology, and now makes use of the "present-day" field produced by Mahowald (2005).Figure 3 shows a map of annual average iron deposition.However, as with MEDUSA-1.0,aeolian iron solubility was adjusted such that the total addition of dissolved iron to the open ocean by dust was the same as that of Dutkiewicz et al. (2005).Fig. 3.The top panel shows mean annual aeolian iron input to the ocean (i.e. the quantity of iron that dissolves into seawater from deposited dust).The input is shown on a logarithmic scale in units of µmol m −2 yr −1 , and integrated input is 2.564 Gmol Fe yr −1 .The bottom panel shows the fractionation of total iron between "free" and ligand-bound forms across a logarithmic range of total iron concentrations.
A further difference with MEDUSA-1.0lies in the inclusion of a benthic source of dissolved iron.Such a supply route is already known for iron, most noticeably around islands and other areas of shallow water in regions that are otherwise depleted in iron (e.g. the Crozet Archipelago in the Southern Ocean; Pollard et al., 2009), and some existing models already include it (e.g.Moore et al., 2004).Here, a flux of iron is added to ocean cells immediately above the seafloor wherever the water column is shallower than 500 m.There is considerable uncertainty in the addition rate of iron to the ocean by this route (Moore et al., 2004), and here the rate has been chosen such that aeolian and benthic supply routes are of approximately similar magnitude.
MEDUSA's iron state variable, F , represents total iron, and this is assumed to occur in two fractions: "free", F free ; and that bound to organic ligands, F ligand (Gledhill and van den Berg, 1994).In the ocean, it is estimated that more than 97 % of total iron is complexed with ligands (Boye et al., 2003).
The complexation reactions between iron species and ligands occur rapidly, and it is assumed here that they reach equilibrium in a shorter period than the model time step (Rose and Waite, 2003).In the equations above, L total is the total ligand concentration of seawater, and is assumed to be globally constant; k FeL is the ligand binding strength.Given these equations and parameters, Fig. 3 illustrates the resulting partition between "free" and bound iron over a range of total iron concentrations.
Scavenging of iron occurs at a fixed linear rate, k scav , throughout the full volume of the ocean, but is assumed to only remove "free" iron, F free .

Fast detritus production
Sinking detrital material in MEDUSA-2.0occurs in two forms: -Small particles that are assumed to sink slowly and are modelled explicitly (as D and D C ); these particles remineralise at a temperature-dependent rate and are a food item of both micro-and mesozooplankton.
-Large particles that are assumed to sink quickly and whose attenuation down the water column is modelled implicitly; these particles remineralise exponentially with depth and are not available as a food item.
As in MEDUSA-1.0,fast-sinking detrital particles are remineralised down the water column using a variant of the socalled ballasting hypothesis (Armstrong et al., 2002).This scheme posits a relationship between organic material and associated -and protective -biominerals.As the description in Yool et al. (2011) includes extensive treatment of the scheme used in MEDUSA-1.0,here we give a summary overview and focus on the differences in MEDUSA-2.0.
In the first instance, the components of fast-sinking detrital particles are produced by a series of ecosystem processes.Organic material (N, Fe, C) is derived from losses from diatoms and mesozooplankton, the larger components of the plankton.Note that, as with other processes, iron is again coupled to nitrogen via a fixed Fe : N ratio.
Equations 87-91 below relate to the total quantities of fastsinking detrital components, T X , being transferred downwards between model levels -that is, increasing values of model grid index k.The equations express the amount of material entering, k, and exiting, k + 1, a given model layer and the processes that act as sources ("production") and sinks ("remineralisation") for these quantities.Later equations describe these sink, or loss, terms, LD X .
FD production (87) Inorganic biogenic opal (Si) is derived directly (via cell mortality) or indirectly (as a product of mesozooplankton grazing) from diatom phytoplankton.In MEDUSA-1.0, the fraction of grazed opal that became associated with fastsinking detritus was the same as the fraction of mesozooplankton losses that were similarly channelled, D2 frac .Here, a new parameter, D3 frac , has been introduced to allow the separate specification of this transfer efficiency.
Calcium carbonate, CaCO 3 , is also an important biomineral in the ballast hypothesis, but its production is not modelled explicitly in either version of MEDUSA.This decision to omit calcification in MEDUSA stems from the diversity (phylogenetic and trophic) of organisms that manufacture CaCO 3 and the uncertainty in the ecological factors that regulate it, as is evidenced by the wide range of approaches used to model it (e.g.Tyrrell and Taylor, 1996;Moore et al., 2002;Gehlen et al., 2007;Zahariev et al., 2008;Yool et al., 2010).Instead, MEDUSA adopts an empirical approach in which the only calcification explicitly considered is that associated with sinking material; CaCO 3 that is synthesised and dissolved without significant vertical movement is considered tangential.
Following Dunne et al. (2007), MEDUSA-1.0used a simple triangular function of latitude, fc(lat), to calculate the relative quantity of CaCO 3 associated with fast-sinking detrital particles, the so-called "rain ratio" (highest values at the equator and lowest values at the poles).MEDUSA-2.0retains this functionality as an option, but introduces a further option that instead calculates associated CaCO 3 as a function, fo( calcite ), of the ambient saturation state of the CaCO 3 polymorph calcite.
This is based on the formulation of Ridgwell et al. (2007), and uses the concentrations of calcium (seawater average; scaled by salinity) and carbonate (calculated from DIC) ions to calculate calcite .Options exist in MEDUSA-2.0for the rain ratio to be based on calcite at the ocean surface or at the local position within its interior (via switch jrratio; see Table 6).Parameter r 0 has been scaled in MEDUSA-2.0so that total production of CaCO 3 using Eq. ( 92) approximately matches that in MEDUSA-1.0(see later).Note that in the real ocean a second polymorph of CaCO 3 is also produced, aragonite, but for simplicity calculations are performed as if all CaCO 3 in MEDUSA-2.0 is the more stable polymorph, calcite (though the saturation state of aragonite, aragonite , is calculated as a diagnostic variable).

Fast detritus remineralisation
The ballast hypothesis of Armstrong et al. (2002) posits that a fraction of the sinking organic material is quantitatively associated with sinking inorganic material (here calcium carbonate and biogenic silica), and that this provides "protection" for the organic matter, allowing it to penetrate deeper into the water column than might otherwise be expected.Follow-up work by Klaas and Archer (2002) derived a parameterization of the hypothesis based on a global data set of sediment trap measurements, and this latter study has subsequently been used as the basis for other work.Its implementation by Dunne et al. (2007) was that adopted by MEDUSA-1.0,and this has been retained by MEDUSA-2.0.
By way of summary, the fast-sinking detrital flux of organic carbon is proportioned into so-called "protected", TC protect = (TC bSi +TC bCaCO 3 ), and "excess", TC excess , portions as follows.
Where M Si and M CaCO 3 convert molar silicon and calcium carbonate ballast into mass equivalents that can then be used with mass-based organic carbon protection ratios f Si and f CaCO 3 .The "protected" fraction passes through unscathed to the next level down the water column, while the "excess" fraction is attenuated across a particular level, with a corresponding release of inorganic carbon.Not all "excess" carbon is remineralised in a given level, and the surviving portion, TC survive , is calculated as follows.
Leaving aside that added through production (see Eq. 89), the quantity of fast detritus reaching the next model layer, T C (k + 1), is then as follows.
The flux of remineralised carbon to level k is then simply as shown below.
The remineralisation fluxes of nitrogen and iron follow that of carbon, with the same fraction of sinking material "protected" by ballasting minerals.By contrast, the sinking fluxes of both biogenic silica, T Si (k), and calcium carbonate, T CaCO 3 (k), attenuate with depth independently of organic carbon.
In the case of biogenic silica, this attenuation occurs at all depths because it is globally undersaturated with respect to ambient silicic acid concentrations.The equations governing sinking biogenic silica and its dissolution are as follows.
Unlike biogenic silica, CaCO 3 is generally not soluble in surface waters because of supersaturating concentrations of the carbonate ion.However, at depth, specifically below the lysocline, concentrations become undersaturating and dissolution can occur.if z(k) < lysocline(lat, lon) In MEDUSA-1.0, the depth of the lysocline, lysocline(lat, lon), was precalculated using physical and biogeochemical fields from the World Ocean Atlas and GLODAP climatologies (Locarnini et al., 2010;Antonov et al., 2010;Key et al., 2004).Here, the inclusion of DIC and alkalinity, as well as a carbonate chemistry submodel, allows MEDUSA-2.0 to calculate the saturation state of CO 2− 3 at all depths, and to use this to determine the point in each water column at which biogenic CaCO 3 will begin to dissolve.The dissolution flux calcium carbonate is then simply as follows.
Figure 20 (see later) shows the depth of the simulated lysocline.

Computation
The structure of MEDUSA is such that the production of particles of fast-sinking detritus has a variable vertical distribution that depends upon location-specific details in plankton dynamics.Consequently, production and remineralisation of fast-sinking particles occur in parallel down the water column, unlike the situation in the source model for this part of MEDUSA, Dunne et al. (2007).As described previously in Yool et al. (2011), the computation of the distribution and fate of fast-sinking detritus is performed layer-bylayer down the water column.All fast-sinking detritus produced within one layer is exported to the next layer, together with the fast-sinking detritus from preceding layers that is not remineralised within that one layer.MEDUSA then iterates this process down the water column to the seafloor, at which point all fast-sinking detritus has either been remineralised or is transferred to the benthic reservoirs.The entire procedure is implicit and occurs within a model time step.As such, the fast-sinking detritus component of MEDUSAunlike the slow-sinking component -does not explicitly consider detrital sinking speeds or remineralisation rates, but instead operates in an e-folding length scale approach akin to the canonical empirical export model of Martin et al. (1987).
Additional explanation of the fast-sinking detritus scheme can be found in commentary within the model source code that accompanies this document.

Alternative models
Separate from the ballast model, MEDUSA-2.0includes a code switch, jexport, to permit the use of two alternative remineralisation schemes for the organic components of fastsinking detritus: the classic Martin et al. (1987) curve; and the variant developed by Henson et al. (2012).Both models attenuate organic material using the same power relationship shown below.
Parameterized using the limited data that was available at the time, the Martin et al. (1987) In the work described here, only the ballast scheme is formally used, though the significance of these (and, potentially, other) schemes will be the subject of future work.

Air-sea gas exchange
MEDUSA-2.0 includes gas exchange for two modelled constituents, O 2 and CO 2 .In the case of O 2 , the scheme developed by Najjar and Orr (1999) for the OCMIP-2 project is used.In this, the saturation concentration of O 2 is calculated based on local temperature and salinity, and this is used in conjunction with ocean surface O 2 concentration and wind speed (via standard gas transfer calculations) to calculate airsea exchange.
The case of CO 2 is complicated by the intricacies of carbonate chemistry, which necessitates the iterative calculation of surface ocean pH to determine surface H 2 CO 3 concentration.As with O 2 , this is then combined with atmospheric pCO 2 and wind speed to calculate the air-sea exchange of CO 2 .The numerical scheme used here is that published by Blackford et al. (2007) (and utilised in Artoli et al., 2012).Alongside air-sea exchange, this scheme calculates other carbonate chemistry properties that are utilised by MEDUSA-2.0,such as calcite .It also permits the calculation of all these properties at arbitrary depths down the water column, and is used in MEDUSA-2.0 to determine the location of the CCD.
Surface gas exchange calculations are performed at every model timestep.Carbonate chemistry calculations are only performed for the full water column on a monthly timescale to reduce computational burden.

Oxygen cycle
Since its cycle is tightly coupled to that of nitrogen and carbon, the differential equation for dissolved oxygen, Eq. ( 15), contains a large number of terms.However, these are largely replicated from other differential equations, and scaled by the appropriate stoichiometric ratio, θ nit or θ rem .The oxygen stoichiometry used here follows that of Yool et al. (2010), and is ultimately derived from one of the range of estimates put forward by Anderson (1995) to account for the production by phytoplankton of a suite of organic molecules (including lipids, proteins and nucleic acids) in addition to carbohydrates, (CH 2 O) n .The resulting C : N : O 2 stoichiometry of organic matter production is as follows.
The remineralisation of organic matter in MEDUSA-2.0 is assumed to be be precise reverse of this.This translates to a C : N : O 2 ratio of 106 : 16 : −151 for the organic matter produced, and effectively assumes that primary production is fuelled by nitrate, as well as the complete nitrification of organic nitrogen back to nitrate during remineralisation.
In addition to suggesting a reduced H and O content of organic material relative to conventional carbohydrate synthesis, 106 : 16 : −138, it also results in an increased production of O 2 per mol of fixed carbon.For reference, using similar assumptions to those above, organic matter production using ammonia, and remineralisation back to the same, would have a corresponding C : N : O 2 stoichiometry of 106 : 16 : −119.
To facilitate accounting given variable C : N ratios across MEDUSA-2.0,the terms listed in Eq. ( 15) (and the model code) separate oxygen production and consumption according to whether nitrogen or carbon are being remineralised by a particular process.Given the above stoichiometry, this gives a O 2 : N ratio for θ nit of 2 : 1 (= 2.0), and a O 2 : C ratio for θ rem of 119 : 106 ( 1.1226).
Note that, following Najjar and Orr (1999), dissolved oxygen is consumed down to a minimum concentration, O min , below which remineralisation can still take place (using unspecified and unmodelled oxidants) but without consuming oxygen.

Alkalinity cycle
The partial differential equation for alkalinity, Eq. ( 14), has only three terms: one for CaCO 3 production, and one each for pelagic and benthic dissolution.As described previously, the production of CaCO 3 is a function of the production of fast-sinking detritus and ambient calcite .Dissolution occurs below the calculated CCD (see Fig. 20) and at the seafloor regardless of CCD depth, in order to prevent drift in pelagic alkalinity inventory.This simplicity reflects the aim, in both MEDUSA-1.0and MEDUSA-2.0, of representing the dominant driver of alkalinity distributions, the so-called "hard tissues" component of the biological pump (cf.Najjar and Orr, 1999).A result of this approach is the omission from consideration of secondary processes that would require a more complex treatment (cf.Wolf-Gladrow et al., 2007;Paulmier et al., 2009).For instance, as noted above, bulk DIN is considered in MEDUSA-2.0rather than separate nitrate and ammonia (as well as other species), the differing use (and remineralisation) of which impacts proton consumption and production and, thus, distribution of alkalinity.Similarly, the inclusion of explicit calcifiers with variable abundance, and potentially dynamic calcification, would require a more complex alkalinity cycle.However, the restricted set of actors and processes selected for inclusion in the current version of MEDUSA limit the need for a more sophisticated submodel.

Miscellaneous
In MEDUSA-1.0, the same Redfield C : N ratio of 6.625 was assumed for both phytoplankton and zooplankton so that the pool of detritus was fed C and N at the same ratio regardless of the source.With the inclusion of a separate detrital carbon pool, D C , these ratios no longer need to be identical, and both micro-and mesozooplankton are assumed to have a lower C : N ratio, 5.625 (Anderson and Pondaven, 2003).

Parameter values
Tables 1-6 list model parameters, a brief description of each, and their respective units and default values.For ease of use, Fe nutrient uptake half-saturation constants 0.33, 0.67 µmol Fe m −3 the ordering of parameters closely reflects their appearance in the namelist.trc.smsfile in which they are specified (see Appendix A and accompanying model code).
Almost all parameter values in MEDUSA-2.0are identical to those from MEDUSA-1.0,though there are a small number of minor changes, and several additional parameters that relate to new state variables.Regarding parameters with reassigned values, the diatom half-saturation concentration for silicic acid uptake, k Si , has been increased (0.75 → 3.0) to a value more congruent with studies such as Fasham et al. (2006).Small detritus sinking velocity, w g , has been slightly decreased (3.0 → 2.5) to favour shallower remineralisation and near-surface nutrient retention.Reflecting the addition of the carbon cycle, the assimilation efficiencies of both zooplankton types are now specified separately for the nitrogen (0.77) and carbon (0.64) ingested during grazing (from MEDUSA-1.0'scommon value of 0.69; Anderson and Pondaven, 2003).New parameters include a separate remineralisation rate for detrital carbon, a series of oxygen stoichiometry parameters, a minimum concentration for dissolved oxygen consumption, and a series of remineralisation/dissolution rate parameters for the benthic reservoirs.
In addition to the parameters above, MEDUSA-2.0includes a number of control parameters that allow the model to switch between different functional forms for a small number of processes.These appear in namelist.trc.sms and are listed in Table 6.As noted above, the control parameters available in MEDUSA-1.0have been augmented by several new options including export submodel, jexport, rain ratio calculation, jrratio and CCD calculation, jocalccd.

Default simulation
The following section describes a simulation and evaluation of MEDUSA-2.0using the default equations, functional forms and parameter values described previously.Evaluation is performed against observational data, but also with MEDUSA-1.0itself.
Both NEMO and MEDUSA-2.0were initialised at the timepoint of midnight on 1 January 1860.This is a standard point in HadGEM2-ES simulations for CMIP5.The model was then run out to 30 December 2005.Note that this is the final day of the year in the 360 day calendar of the atmospheric forcing used here.

Physical model
The underlying physical model used in this simulation is version 3.2 of NEMO (Madec, 2008).This is comprised of an ocean general circulation model, OPA9 (Madec et al., 1998;Madec, 2008), coupled with a sea-ice model, Louvainla-Neuve Ice Model version 2 (LIM2; Timmermann et al., 2005).This physical framework is configured at approximately 1 • × 1 • horizontal resolution (292 × 362 grid points), with a focusing of resolution around the equator to improve the representation of equatorial upwelling.Vertical space is divided into 64 levels, which increase in thickness with depth, from approximately 6 m at the surface to 250 m at 6000 m.To improve the representation of deep water circulation, partial level thicknesses are used in the specification of bottom topography.Vertical mixing is parameterized using the turbulent kinetic energy (TKE) scheme of Gaspar et al. (1990), with modifications by Madec (2008).
The sea-ice submodel used here, LIM2, is based upon viscous-plastic ice rheology (Hibler, 1979) and three layer (two layers of sea ice, one layer of snow) thermodynamics (Semtner, 1976), with a number of updated physical processes (see Timmermann et al., 2005; and references therein).Model sea ice is coupled to the ocean every 5 ocean timesteps through the non-linear quadratic drag law of the shear between sea ice and ocean surface velocity (Timmermann et al., 2005).Freshwater exchange between the ocean and sea ice is calculated from precipitation and ice formation/melting (Fichefet and Morales Maqueda, 1997), where sea-ice salinity is assumed to be 4 psu and rain/snow are assumed fresh.The heat flux between the sea ice and ocean is proportional to the departure in temperature from salinity-dependent freezing point and the friction velocity at the ice-ocean interface.Solar radiation can penetrate sea ice not covered by snow, and is dissipated by brine pockets within the ice where it increases latent heat storage (Fichefet and Morales Maqueda, 1997).
In Yool et al. (2011), NEMO was forced at the ocean surface for the period 1966-2005 using DFS4.1 fields developed by the European DRAKKAR collaboration (DRAKKAR Group, 2007).As MEDUSA-2.0includes the ocean's carbon cycle, and since this is currently undergoing secular change driven by increasing atmospheric concentrations of CO 2 , simulations running over a longer period of time are necessary.There are a number of approaches to achieve this including, for instance, the use of a climatological average or "normal year" (e.g.Najjar et al., 2007), or the repeated cycling of historical forcing (e.g.Yool et al., 2010).These have the advantage of using actual observationally derived forcing, but also assume that the recent past from which they are derived is representative of earlier periods of time (in spite of ongoing climate change).An alternative approach is to utilise forcing derived from either atmospheric models or coupled ocean-atmosphere models.These are routinely run in long duration simulations that span pre-industrial or pre-20th century periods when there was comparatively little change in climate or the carbon cycle.They also offer the opportunity to forecast biogeochemical cycles into the future with a significantly different climate from that of the present-day.
Here, NEMO is forced following this latter approach, using output from a simulation of the HadGEM2-ES Earth system model run by the UK Meteorological Office (UKMO).HadGEM2-ES is a development of the physical climate model, HadGEM1 (Johns et al., 2006), that includes representations of the terrestrial and oceanic carbon cycles, atmospheric chemistry and aerosols (Collins et al., 2011).The HadGEM2-ES simulation used here, identifier AJKKH, was performed as part of the UKMO's input (Jones et al., 2011) to the Coupled Model Intercomparison Project 5 (CMIP5) and Assessment Report 5 (AR5) of the Intergovernmental Panel on Climate Change (IPCC).Operationally, HadGEM2-ES output was processed into the same forcing fields as that provided by the DFS4.1 forcing previously used with MEDUSA-1.0.The frequency of the output fields also matched that of DFS4.1, namely monthly for precipitation (rain, snow, runoff), daily for radiation (downwelling shortand long-wave) and 6-hourly for the turbulent variables (air temperature, humidity and wind velocities).Note that the reference height of forcing in HadGEM2-ES differs from that of DFS4.1, but that NEMO's bulk formulae allow this height to readily be changed to accommodate HadGEM2-ES.
For maximum congruence with the surface forcing, temperature and salinity fields are initialised here using output from HadGEM2-ES valid for the same time as the forcing.To prevent excessive drift, sea surface salinity (SSS) is relaxed towards that derived from HadGEM2-ES.Unlike simulations under DFS4.1,where an invariant monthly mean climatology of SSS values is used, here the SSS target consists of a monthly time series running across the forcing period.The relaxation timescale is approximately 30 days for the open ocean, and 12 days under sea ice.The freshwater budget is also monitored for imbalances between integrated downward and upward fluxes, and a correction term applied between years (i.e. an imbalance in year X is corrected for in year X + 1).
Further details concerning physical model configuration can be found in Barnier et al. (2006), Penduff et al. (2007) and Penduff et al. (2010), but note that these describe higher resolution instances of NEMO.

Biogeochemistry
MEDUSA-2.0'sfields of DIN, silicic acid and oxygen were initialised using January values from the World Ocean Atlas 2009 (Garcia et al., 2010a, b).Similarly to MEDUSA-1.0,total iron was initialised using an iron field derived from a longduration simulation of a lower resolution GCM (Parekh et al., 2005;Dutkiewicz et al., 2005).DIC and alkalinity were initialised using a modified form of the GLODAP climatology (Key et al., 2004).It was assumed that GLODAP's preindustrial DIC field is approximately valid for the 1860 start of this simulation, though this approach has known issues concerning the ocean's anthropogenic CO 2 inventory in 1860 (e.g.Yool et al., 2010).
The GLODAP fields used here required modification to account for large regional lacunae including the Arctic Ocean, the Caribbean Sea, the Mediterranean Sea and the Malay Archipelago.These were filled through an approach utilising multiple linear regression (MLR) together with the more complete WOA 2009 fields of temperature (Locarnini et al., 2010), salinity (Antonov et al., 2010), DIN, phosphate, silicic acid and oxygen.For each missing region, values of these tracers in immediately adjacent areas were used to construct a unique MLR.The calculated MLR was then used to fill the lacuna using field values from the WOA 2009.As biogeochemical tracers frequently show strong vertical gradients, separate MLRs were constructed for a series of intervals down the water column (0-50, 50-100, 100-200, 200-500, 500-1000, 1000-2000, below 2000 m).This procedure was used first with alkalinity, and then the resulting alkalinity field was added to the list of input fields for the construction of MLRs to fill DIC lacunae.While extrapolating in this fashion is likely to introduce some spurious values, particularly where WOA 2009 fields are already uncertain (e.g. the Arctic Ocean), it resulted in fields of DIC and alkalinity that appeared more credible than extrapolation by simple floodfilling was able to achieve.
All other model tracers (plankton and detritus) were initialised to arbitrary small values.Benthic reservoirs of nutrients, carbon and CaCO 3 were set to zero.Note that, unlike in MEDUSA-1.0,no coastal relaxation fluxes were applied to nutrients (N, Si) in MEDUSA-2.0.This change reflects both the switch to forcing periods outside the "present-day", and the finding in Yool et al. (2011) that this relaxation scheme did not universally emulate the riverine addition of nutrients as originally intended.

Results
In this section, a selection of model results are presented with the aim of providing an overview of MEDUSA-2.0'sperformance.In the first instance, model outputs that can be compared to observational fields are presented.These are followed by Taylor diagrams that aim to provide a quantitative  This format of presentation and analysis is generally repetitive of that for MEDUSA-1.0as described in Yool et al. (2011).However, since the simulation of MEDUSA-2.0here is of considerably longer duration than than analysed for MEDUSA-1.0(146 yr versus 41 yr), the results are of particular interest because they permit evaluation of the model's longer-term behaviour and stability.To extend the utility of this analysis, it concludes with an intercomparison of MEDUSA-2.0with a selection of CMIP5 models.Observational fields used in comparison with MEDUSA-2.0are comprised of WOA 2009 nutrients (Garcia et al., 2010b), SeaWiFS chlorophyll (O'Reilly et al., 1998), estimated primary production (Behrenfeld and Falkowski, 1997;Carr et al., 2006;Westberry et al., 2008), GLODAP carbon and alkalinity (Key et al., 2004) and air-sea CO 2 exchange (Takahashi et al., 2009).Because of its biogeochemical importance, and the diversity in estimates of it, observational primary production is drawn here from three empirical models: VGPM (Behrenfeld and Falkowski, 1997); Eppley-VGPM (Carr et al., 2006); and CbPM (Westberry et al., 2008).The observational fields of chlorophyll and productivity used here represent averages over the same 5 yr period from 2000 to 2004 inclusive, and this same period is used throughout the following analysis as a standard interval except where noted.
The philosophy behind selecting these fields for the purpose of model validation has several facets.Firstly, they are ocean properties that have been observed at the global scale that MEDUSA-2.0 is simulated at.In the case of surface chlorophyll, this is now estimated by remote sensing at a fine spatial scale on a continuous basis.Dissolved tracers are much less well-sampled, but coverage has still been sufficient for high quality climatologies of each to be assembled.Secondly, they generally represent quantities that are believed to be the foundation of biological oceanography.Nutrient distributions, for instance, play a critical role in structuring ocean communities in both space and time, while primary production is the overwhelming route by which energy-rich organic carbon enters the marine food web.Thirdly, their measurement is well-defined and open to relatively little ambiguity.Properties that are more directly related to biological entities or processes can be more difficult to measure in the field, and more difficult to marry with model "equivalents".That said, synoptic estimates of primary production -a property examined here -still carry relatively high uncertainty, as evidenced by the range in estimates produced from the same inputs.Finally, and this is in part a corollary of the above, they are properties which, if modelled poorly, can cast legitimate doubt over the utility of a biogeochemical model as a whole.Models will always have discrepancies with observations, but any systematic failure to capture at least qualitative aspects of these quantities in particular will strongly suggest model weakness.
Note that, as well as from these geographically synoptic fields, MEDUSA-2.0 is also compared with more sparse observations of modelled quantities such as zooplankton and with globally integrated estimates of quantities such as biogenic opal production.

Validation
Figures 4 and 5 compare MEDUSA-2.0'sperformance in representing, respectively, surface concentrations of the macronutrients DIN and silicic acid (note that here, and subsequently, "surface concentration" refers to concentration within the uppermost model level, 0.00-6.06m).In the case of DIN, MEDUSA-2.0shows generally good agreement in the Northern Hemisphere, but with noticeably higher concentrations in both equatorial upwelling regions and in the Southern Ocean.A similarly strong Southern Ocean bias was found with MEDUSA-1.0,though equatorial waters there showed a slight bias in the opposite direction.Silicic acid concentrations are very similar between both MEDUSA versions, and show the very same patterns of bias.Most noticeably, markedly elevated Southern Ocean concentrations, uniformly too-low equatorial concentrations, and concentrations in the Northern Pacific lower than those observed in this HNLC region.Figures 6 and 7 show corresponding, basinaveraged Hovmöller diagrams of DIN and silicic acid for the Atlantic and Pacific Oceans.
Focusing on the deep ocean, Figs. 8 and 9 show zonally averaged sections of DIN and silicic acid down the Atlantic and Pacific basins (the Atlantic includes the Arctic Ocean; both basin sections include the Southern Ocean).In both cases, most large-scale structure has persisted in MEDUSA-2.0across the run duration.However, there are some important differences, of which the Southern Ocean is the most extreme.In this region, excessive ventilation acts to homogenise horizontal and vertical gradients, most notice-  ably those of silicic acid.A similar problem in the Southern Ocean was noted by Yool et al. (2011) and ascribed to a deficiency in NEMO, but the problem here is somewhat worse and that this may stem from the change in surface forcing.
An examination of the large-scale circulation of the run finds that the Antarctic Circumpolar Current (ACC) is significantly stronger (220 Sv) in this simulation compared to that used with MEDUSA-1.0(160 Sv), and toward the high end of other models (CMIP5 range of 90-264 Sv; Meijers et al., 2012).This is associated with stronger Antarctic Bottom  (and similarly for DIC and alkalinity; Figs. 17 and 18).However, as noted above, much of the zonal structure in the rest of the World Ocean is maintained, even in the case of dissolved oxygen.So while an improved circulation state would certainly be preferred, the impacts for MEDUSA-2.0 of NEMO's "robust" Southern Ocean ventilation are somewhat restricted.
Leaving aside these significant circulation-driven changes in tracer distributions, the section-averages of Figs. 8 and 26 indicate that MEDUSA-2.0itself may also be a source of model-observation discrepancy.In both the real and modelled oceans, the remineralisation of sinking detrital material results in concentrations of DIN being generally elevated with depth, while those of dissolved oxygen broadly decline.However, in the case of MEDUSA-2.0, the highest DIN concentrations, and lowest oxygen concentrations, occur at noticeably shallower depths than in the WOA (2009).This is most pronounced in the Pacific panels of Figs. 8 and 26, but this mismatch is also apparent in the Atlantic panels.While deficiencies in circulation will also play a role, the remineralisation of sinking detrital material at depths that are too shallow plays a part.There are a number of parameters in MEDUSA-2.0that may be implicated in this including the remineralisation rate of slow-sinking detrital particles, µ D , the corresponding sinking velocity, w g , and the remineralisation length scale of fast-sinking detritus, d excess .Ironically, one of these parameters, w g , was changed from its MEDUSA-1.0value to improve model behaviour, and was altered in exactly the direction that would cause the current mismatch.Note, however, that the convolution of biogeochemical processes and circulation (as well as the long spinup periods required to detect model-observation mismatch) prevents a clean separation or quantitative evaluation of the cause or causes of discrepancies such as this.Returning to the surface ocean, Figs. 10 and 11 compare MEDUSA-2.0'ssimulated total chlorophyll (non-diatom plus diatom) to corresponding SeaWiFS fields (note that a logarithmic colour scale is used to best represent the large range in ocean colour).Not uncommonly for ocean models, and similarly to MEDUSA-1.0, the representation of chlorophyll exhibits significant discrepancies with observations.MEDUSA-2.0shows much less pronounced seasonality, particularly at higher latitudes in the Northern Hemisphere, and spatial boundaries that are significantly more sharply defined and consistently lower "background" chlorophyll concentrations in the ocean gyres.While the latter regions are not productive areas of the ocean, they represent a significant fraction of its total area.This was also noted with MEDUSA-1.0,and speculatively attributed to the assumption of geographically invariant nutrient kinetics.This prevents model phytoplankton from adapting to oligotrophic conditions when, in the real world, nutrient uptake kinetics are more plastic (e.g.Smith et al., 2009).However, given the globally uniform parameterization of ecosystem actors in MEDUSA, it may be difficult to resolve this deficiency without more fundamental changes to the model framework.For instance, the addition of further phytoplankton types with parameter values more "at home" in oligotroph conditions.
Figures 12 and 13 compare MEDUSA's simulated total primary production (non-diatom plus diatom) to a simple average of the estimates of the VGPM, Eppley-VGPM and CbPM models.The average estimated production has been used here both to simplify intercomparison and because, while sharing inputs, the separate estimates disagree significantly with one another (to the extent that model-observation difference is comparable with the range of observational estimates; Fig. 48).In broad terms, MEDUSA-2.0captures some   of the spatial and seasonal patterns in productivity, though it does show significant systematic differences as well.These include: consistently low subtropical gyre productivity; elevated productivity in iron-limited regions including the Southern Ocean, equatorial Pacific and (seasonally) North Pacific; and a weaker bloom across the North Atlantic.In terms of total oceanic primary production, MEDUSA-2.0predicts 41.6 Pg C yr −1 , a value slightly below the bottom of the Figures 14-15 show the corresponding model-observation comparisons using Taylor diagrams.These illustrate both the correlation between (circumference axis) and relative variability (radial axis) of model and observations.For each comparison two plots are shown.The first uses annually average fields, but separates the analysis between major ocean regions; the second uses globally average fields, but separates the analysis between months.In all cases, modelobservation agreement is greater the closer plotted data are to the red/black bullseye on the horizontal axis.
Similarly to MEDUSA-1.0, the best agreement occurs with nutrient fields, particularly DIN.While there remains significant scatter, MEDUSA-2.0generally shows good correlation with World Ocean Atlas 2009 fields, and comparable magnitudes of variability.In the case of surface silicic acid, there is considerable variability between basins with the Pacific performing very poorly, and the Indian exhibiting significantly elevated variability.Much as with MEDUSA-1.0,agreement is still very weak in the case of chlorophyll, where the model both correlates poorly and shows much less variability than the observed SeaWiFS fields.Although estimated productivity is based on the same SeaWiFS chlorophyll fields, MEDUSA-2.0'sagreement with the three productivity models is actually much greater, particularly the VGPM and CbPM models (results not shown), although correlations are still relatively weak.
Extending beyond MEDUSA-1.0,Fig. 16 compares MEDUSA-2.0'ssurface fields of 1990s average DIC and alkalinity concentration with those from the GLODAP climatology.As with preceding fields, there is broad agreement, but a number of notable differences.In the case of DIC, regions of high concentration are typically slightly higher (< 50 mmol C m −3 ) in MEDUSA-2.0,with the greatest difference occurring adjacent to Antarctica.In the case of alkalinity, the Southern Ocean again shows elevated concentrations caused by the increased ocean ventilation in this region mentioned earlier, but in the Pacific discrepancies are generally downwards.While, again, circulation deficiencies are in part responsible, excessive production in the Pacific, particularly in the equatorial region, is responsible for this discrepancy.In this region, MEDUSA-2.0'srain ratio formulation means that elevated primary production and export also drives a stronger export of CaCO 3 -with lower surface alkalinity a principle result.As Fig. 19 shows, this general pattern of elevated DIC but decreased alkalinity also manifests in other surface properties.Here, modelled pH is globally lower than that estimated from GLODAP, with the result that calcite is also lower than that estimated, though insufficiently to significantly impact Eq. ( 92) and the rain ratio of export.
Following on from this, Fig. 20 shows a comparison of the observationally based lysocline (calculated from WOA and GLODAP fields) with that simulated by MEDUSA-2.0for the 1990s (which corresponds to GLODAP's "presentday").In both cases, the same carbonate chemistry routine used in MEDUSA-2.0 is used to calculate calcite , and the CCD shown is the shallowest depth at which calcite has a value below 1.0 (i.e. the first depth at which the concentration of CO 2− 3 is undersaturated).While the two maps broadly agree in terms of overall pattern, there are several regions where there are significant differences and -globally averaged -modelled CCD is noticeably shallower, 2646 m compared to 2864 m from observations.The Arctic (model vs. observations: 1347 m vs. 1347 m) and Indian (3322 m vs. 3308 m) basins have very close average CCDs, but the Southern (3318 m vs. 3166 m), Atlantic (3165 m vs. 3341 m) and, especially, the Pacific (2008 m vs. 2482 m) oceans show much larger discrepancies.Throughout much of the World Ocean, these differences are caused by relatively minor rearrangements of the vertical gradients of DIC and alkalinity by modelled circulation and biogeochemistry.In the case of the eastern Pacific, the large errors in CCD depth are a direct result of the excessively shallow remineralisation of sinking organic material previously identified in the zonal sections of dissolved inorganic nitrogen, carbon and, especially, oxygen (Figs. 8,17 and 26).Between 150 and 1500 m depth, this region shows significantly elevated concentrations of DIC that depart markedly from observations, with the result that values of calcite drop below 1.0 much closer to the surface than in reality.Near-surface values of calcite are naturally less than 2.0 in this region, so model errors of this kind can easily impact CCD.They also, in part, reflect positive     feedback driven by the ballast submodel: too-shallow remineralisation increases near-surface DIC; this decreases local calcite ; this shallows the depth at which CaCO 3 dissolves; this decreases the "protection" offered to sinking organic material by ballast; and, in turn, this shoals remineralisation and increases near-surface DIC.It is difficult, however, to determine whether such feedback is the original cause of this error or merely a downstream consequence of some other problem with NEMO or MEDUSA-2.0(e.g.physical circulation; biogeochemical parameter changes).Changing to the ocean carbon cycle's interaction with the atmosphere, Figs.21-24 compare MEDUSA-2.0 to observationally derived fields of pCO 2 and air-sea CO 2 flux for year 2000 (Takahashi et al., 2009).The former is simply the localised difference between surface ocean pCO 2 and that of the atmosphere (assumed a globally uniform but timevarying quantity in the model).The latter is an estimate of the actual net exchange of CO 2 between the ocean and the atmosphere (where positive values indicate net air-to-sea flux), based on pCO 2 , air pressure, piston velocity and sea-ice concentration.(June-July-August; top) and northern (December-January-February; bottom).∆pCO2 in ppm.probably an underestimate because of undersampling (which they suggest would increase it to 1.6 Pg C yr −1 ).MEDUSA-2.0 also exchanges oxygen with the atmosphere but, as Fig. 25 illustrates, its solubility is both greater than that of CO 2 and is strongly temperature-dependent.As a result, though there are discrepancies in MEDUSA-2.0'sprediction of surface dissolved oxygen concentrations, the modelled seasonal distributions are very similar.What differences there are occur in the over-ventilated Southern Ocean, and in areas where large-scale circulation produces strong sea  surface temperature features that are displaced in NEMO.For example, in the Atlantic, where the position of the Gulf Stream determines SST gradients and, thus, those of surface dissolved oxygen.As previously noted, Fig. 26 shows larger discrepancies at depth, caused in part by circulation deficiencies (Southern Ocean; all depths) and partly by MEDUSA-2.0itself (Pacific; midwater).In terms of ocean anoxia, suboxic regions (< 20 mmol O 2 m −3 ) occupy 13.1 × 10 6 km 3 in the WOA but are more than double this in MEDUSA-2.0at 30.4 × 10 6 km 3 .
Switching to ecosystem properties for which observations are less synoptic, Figs.27-35 show seasonal and geographical plots for a range of model fields.
Figures 27 and 28 respectively show the split between surface biomass and integrated production for MEDUSA-2.0'stwo phytoplankton groups (shown on the same colour scales to facilitate intercomparison).Much as with MEDUSA-1.0,non-diatoms are dominant across most of the World Ocean, and particularly in the oligotrophic gyres, where diatom abundance and productivity is extremely low.However, diatom biomass can seasonally exceed that of the non-diatoms in regions such as the North Atlantic, and they still contribute modestly to total primary production (15.9 %; 17.0 % in terms of total nitrogen biomass).Observational estimates of this fraction at the global scale are rare.While a survey by Mann (1999) suggested 40-45 %, this is much greater than that estimated by either MEDUSA-2.0or localised observations (13-34 %; Nelson and Brzezinski, 1997;Blain et al., 1997;Brzezinski et al., 1998).The left panels of Fig. 29 show the fraction (0-1) of total primary production that MEDUSA-2.0predicts for the upper mixed layer, with the remainder occurring deeper in the water column, often in simulated deep chlorophyll (or productivity) maxima.In general, this fraction is lower in the summer, when nutrients are more limiting than light, and higher in the winter, when light limits production more.Patterns are less clear in the tropics and upwelling regions where the interplay of nutrient and light availability is more complex.In the case of northern latitudes, the ratio generally shifts between 0.5 and 1.0, but in the Southern Ocean the seasonal range is 0.7-1.0,reflecting this basin's all-year macronutrient  availability.Integrating to the global scale, slightly more than two-thirds (67.3 %) of production occurs in the mixed layer in MEDUSA-2.0,and while there is geographical variation in this fraction between basins, the Southern Ocean is the most different at 85.8 %.
The right panels of Fig. 29 instead show ocean productivity from the perspective of the benthic communities that ultimately rely on them.The panels are shown on a log scale because the geographical variability of export production is compounded exponentially by variability in the seafloor depth that sinking material needs to reach.It shows strong seasonality at high latitudes and low seasonality in the tropics, with the actual magnitude of supply to the benthos strongly tied to seafloor depth (e.g.compare North Sea and Patagonian Shelf regions with adjacent deep water regions).Aside from this depth effect, the seafloor supply of organic material mirrors that of its source, primary production, and so discrepancies in the former (e.g. the excessive seafloor flux around Antarctica) are simply symptomatic of problems in the latter.Predictably, the seafloor distribution of organic material (results not shown) largely follows this input of sinking detrital matter.Globally integrated, MEDUSA-2.0has an inventory of 62.27 Mt C of benthic organic carbon, with a C : N ratio, 6.72 : 1, higher than that of surface plankton.This higher ratio is driven by the (marginally) preferential remineralisation of nitrogen in slow-sinking detritus.Dynamically, after an immediate spike following initialisation at zero, this total inventory is quickly approached (within 5 yr), and the long-term, interannual behaviour of this reservoir tracks that of overlying primary production (see later).The reservoirs of inorganic benthic material, silica and CaCO 3 , behave similarly, though with slightly greater stability since the modelled production of both of these biominerals has lower variability with time.
Following up on ocean productivity, Fig. 30 shows the patterns of limitation by nutrients for both modelled phytoplankton groups.The leftmost panels show overall phytoplankton growth limitation by nutrients (separate from light), while the rightmost panels indicate which nutrient provides the strongest limitation.In the case of diatoms, nitrogen and iron limitation are joined by silicon limitation.In broad outline, nitrogen is most limiting for both phytoplankton groups in oligotrophic gyres, while iron plays a more significant role in high latitude regions, particularly the Southern Ocean, and the equatorial Pacific.Note that, although iron is indicated as most-limiting in both the north Atlantic and Pacific, its impact is greater in the Pacific, particularly in the eastern region.For diatoms, the boundaries between regions of N-and Fe-limitation are typically where Si-limitation occurs, though in the North Atlantic in particular, the scarcity of silicon almost completely displaces iron stress.The geographical patterns in MEDUSA-2.0generally parallel those of MEDUSA-1.0(and other models; Moore et al., 2004), though the change in dust deposition forcing means that the equatorial Pacific experiences a greater degree of iron limitation.As this is one of the more productive regions in MEDUSA, this iron-mediated change also decreases total oceanic primary production.
Switching from the production of organic material, Fig. 31 shows the seasonal production of the biominerals opal and CaCO 3 .As with the preceding plots, production of both is highly seasonal at high latitudes, and more constant at low latitudes.Because of the differential availability of silicic acid, opal production is highest in the North Pacific and Southern Ocean, higher even than that in the tropics, though the latter's annual constancy leads to greater overall production.Patterns of CaCO 3 production -technically its export in MEDUSA-2.0-are similarly seasonal, though the northern Atlantic and Pacific basins swap from the patterns shown with opal.However, within the northern reaches of both basins, opal and CaCO 3 production show different, nonoverlapping geographical patterns.Global total opal production in MEDUSA-2.0 is 194 Tmol Si yr −1 , around 20 % lower than that estimated by Tréguer et al. (1995) (and lower than that in MEDUSA-1.0).Total CaCO 3 export in MEDUSA-2.0 is 0.41 Pg C yr −1 , at the bottom end of the broad 0.4-1.8Pg C yr −1 range estimated by Doney et al. (2009).
Switching again, this time to the consumption of organic material, Fig. 32 shows the seasonal distributions of surface concentrations of both zooplankton groups.Unsurprisingly, both show the same strong seasonality at high latitudes already seen.Though they have slower maximum growth and are less efficient at lower prey concentration, the wider range of available prey types provides a wider base for mesozooplankton and, coupled with their role as predator, makes them dominant in terms of biomass over microzooplankton.Figure 33 compares the spatio-temporal distribution of surface mesozooplankton in MEDUSA-2.0with observational data from the COPEPOD database (O'Brien, 2005).Because of the relative scarcity of zooplankton data, annual means and seasonal zonal means are used.In general, though MEDUSA-2.0has similar patterns of geographical distribution (which ultimately relate to those of phytoplankton), Fig. 33 shows that the model broadly overestimates abundance.This overestimation is greater at higher latitudes during summer months, and MEDUSA-2.0underestimates winter abundance at the same latitudes.Noticeably,  MEDUSA-2.0'ssouthern hemisphere summer peak of mesozooplankton is displaced northwards by more than 15 • .
Figures 34 and 35 show the production of the two size classes of detritus in MEDUSA-2.0,and the export of this material to the deep ocean.As with MEDUSA-1.0, the production of small particles dominates in the surface ocean (70.6 %), but this dominance declines down the water column as these particles are quickly remineralised, such that, by 100 m, small particles are the minority component of the  export flux (38.8 %).By 1000 m, small particles are of almost no importance to abyssopelagic or benthic communities (2.2 %).As an aside, since large, fast-sinking particles have a narrower spatial distribution of production than do slowsinking particles (per Fig. 34), deep water benthic communities in MEDUSA-2.0experience greater variability in supply than do shallow water communities.
In terms of model performance, Fig. 36 compares the total flux of detritus in MEDUSA-2.0at 2000 m with the synthesis data set of Honjo et al. (2008) (both expressed here in car- bon units).As with zooplankton earlier, such observational data is relatively rare, and is plotted here as individual sites on a global map.Model output is plotted as the field of detrital flux at 2000 m overlaid with 3 • × 3 • averages of this field at the same locations as the observations.In general, MEDUSA-2.0agrees relatively well with data, though there are a number of stations at which modelled fluxes are noticeably lower than those observed (e.g. in subtropical stations fringing the oligotrophic gyres).And in the equatorial Pacific region, MEDUSA-2.0'sexcessive productivity is responsible for a mismatch towards higher detrital fluxes than observed.In passing, note that observations are biased towards locations where fluxes are expected to be appreciable, and that large regions of less productive ocean are much less wellsampled.
To illustrate longer-term trends in the performance of MEDUSA-2.0,Figs.37-39 show basin-average vertical profiles of several of the model's major nutrient elements.Since the model is being simulated for a period during which climate change is comparatively limited (though anthropogenic CO 2 is increasing), and since these directly influence the behaviour of MEDUSA-2.0'secological actors (unlike DIC, alkalinity and oxygen; results not shown), they illustrate the degree to which the model has equilibrated.
In the case of nitrogen, per Fig. 37, while globally there is a steady rise in near-surface concentrations -and a steady decrease in deep (2000-5000 m) concentrations -this is largely driven by changes in the Pacific Ocean, with the other basins showing much weaker trends.The pattern of Pacific dominance in the global signal continues with silicon, per Fig. 38, but surface changes are broadly much less significant.Both elements show rapid and significant changes in the Southern Ocean that are consistent with the circulation and water mass changes described earlier.Figure 39 shows the corresponding situation for iron, where the situation is complicated by large removal (scavenging) and addition (aeolian/benthic) fluxes.Here, changes are greatest at depth, where continual scavenging removes deep iron, but there is also a slight general decrease in the surface ocean.An exception lies with the Indian Ocean, which shows almost static surface concentrations.
Complementing these profiles, Fig. 40 shows annual time series of surface nitrogen, silicon and chlorophyll, and integrated primary production for the duration of the simula-  tion.Note that the vertical scales have been focused to emphasise change across the simulation.Consistent with the profile plots above, both surface DIN and silicic acid show an increase during the simulation, but while DIN continues to gradually rise throughout its duration, silicic acid saturates relatively quickly (by 1880).In the case of DIN, this global trend generally reflects that at the surface of the Pacific Ocean, but for silicic acid the global trend is driven by the large increase in surface concentrations in the Southern Ocean.The lower two panels show trends in biological variables that are similar to that of DIN, and this similarity extends to its source, with, again, changes in the Pacific Ocean driving the wider global trend.In the case of primary production, almost all of the increase during the simulation (5 Pg C yr −1 ) is driven by the corresponding increase in the Pacific.
Finally, Tables 7 and 8 intercompare MEDUSA-1.0and MEDUSA-2.0for a range of common properties on a basin average basis.Table 7 focuses on the surface concentrations of major model components, plankton and nutrients.Generally, the two models show very similar patterns, though this is unsurprising given the relatively minor differences between their core nutrient dynamics.However, there are several notable differences in nutrient concentrations.For  instance, surface nitrogen (+9.1 %), silicon (+5.2 %) and iron (+5.8 %) are elevated globally in MEDUSA-2.0,but there are strong regional biases.In the case of the Pacific Ocean, nitrogen increases by +74.4 %, while iron falls by −25.8 %.And some of the largest differences between the models occur in the Arctic Ocean (increased N, decreased Si).In part, the longer duration of MEDUSA-2.0'ssimulation (146 yr versus 40 yr) appears responsible for these differences, but the change in iron deposition forcing, especially in the Pacific, also appears a key factor in the change between the otherwise very similar models.
In terms of major biogeochemical fluxes, Table 8 shows, again, much congruence between the two versions of the model.As mentioned previously, the simulation of MEDUSA-2.0exhibits lower productivity (−8.1 %), with knock-on consequences in phytoplankton biomass (−5.3 %) and opal production (−12.3%).CaCO 3 production is more substantially impacted (−31.3 %), reflecting the compounded declines in both organic production and rain ratio in MEDUSA-2.0.Via the ballast submodel, decline in the production of both opal and CaCO 3 has a further impact on export production: in MEDUSA-1.0,5.3 % of the 100 m flux of organic matter reached 1000 m, while in MEDUSA-2.0,only 4.4 % did.

CMIP5 intercomparison
The preceding agreement between the behaviour of MEDUSA-1.0and MEDUSA-2.0 is to be expected given the parental role of the former in the latter.Their similar disagreement with observations is also to be expected since most of the evolution leading to MEDUSA-2.0relates to extending the model to include additional elemental cycles rather than alterations to the core plankton ecosystem model.Of greater importance, however, is the performance of MEDUSA-2.0relative to comparable biogeochemical models.To examine this, a series of standard outputs already examined above (DIN, silicic acid, chlorophyll, primary production, DIC and alkalinity) were collated from simulations of nine CMIP5 models.Table 9 lists the models examined, together with the number of biogeochemical tracers in each, which serves as a crude measure of their complexity and computational cost.In each case, model output was averaged for the same time period in the preceding analysis, and was compared -together with that from MEDUSA-2.0-against the same observational data sets (see also Appendix B).Note that these models do not share a common physical model framework, and it is not possible to deconvolute the performance of the biogeochemical models from that of their underlying physical models.As such, any weaker performance here may actually be driven by deficiencies in modelled physics rather than in modelled biogeochemistry.
Table 8.Mean annual (2000Mean annual ( -2004) ) biogeochemical properties in MEDUSA-1.0(upper row) and MEDUSA-2.0(lower row).Units indicated for each property.Figure 41 shows the resulting intercomparison as a series of Taylor diagrams based on the annual means of each ocean property.In broad outline, the models tend to show similar patterns of variability and correlation with observations between the considered fields.Fields that MEDUSA-2.0previously correlated well with are generally also well correlated with for CMIP5 models, and vice versa.DIN, DIC and alkalinity, for instance, show strongly clustered, good correlation across all models, while all models show very poor correlation for chlorophyll and intermediate correlation for primary production.Though there are clearly differences between the models, none stand out as being noticeably better or worse than their rivals, in spite of a range of model complexity from 13 to 30 prognostic tracers.Nonetheless, for the fields examined, MEDUSA-2.0 is typically among the best models in terms of model-observation correlation, and exhibits spatial variability that is comfortably within the range of models assessed.
Figures 42-47 show corresponding global fields of the same properties for all of the CMIP5 models examined as well as MEDUSA-2.0and observations.Almost all of the models capture the broad zonal patterns in DIN distribution of moderate concentrations in equatorial upwelling regions, very low concentrations in the subtropics and high concentrations at high latitudes, particularly in the Southern Ocean (Fig. 42).The deficiencies already noted in MEDUSA-2.0adjacent to South America are not unique, though they are the most extreme.Significant deficiencies in other models include elevated gyre concentrations in one model, and excessive Arctic nutrient concentrations in two others (see also Popova et al., 2012).Regarding surface silicic acid, while MEDUSA-2.0'sSouthern Ocean concentrations are too high,  Lengaigne et al. (2007) and those of its North Pacific too low, the range in behaviour of the other CMIP5 models examined is noticeably more divergent (Fig. 43).For instance, the two HadGEM2 model runs systematically overpredict silicic acid throughout much of the World Ocean.Of the other models, they split between over-/under-predicting Southern Ocean silicic acid, and almost all underpredict North Pacific conditions to a stronger degree than MEDUSA-2.0.As Fig. 41 has already suggested, none of the models examined does well at predicting the geographical pattern of surface chlorophyll (Fig. 44).Several models do better than MEDUSA-2.0 in predicting higher (but still low) gyre conditions, but they do so with significantly larger gyre regions than observed.And several models predict large regions with chlorophyll concentrations markedly higher than either observations or MEDUSA-2.0-particularly within the Southern Ocean.Interestingly, in spite of the foregoing issues with surface chlorophyll, the models examined are noticeably better in terms of primary production (Fig. 45).No one model stands out as globally better than the others, but several show much better Southern Ocean productivity than MEDUSA-2.0,though most exhibit even lower productivity in the North Atlantic.Finally, Figs.46 and 47 compare surface DIC and alkalinity across the models.As previously noted, excessive ventilation in MEDUSA-2.0'sSouthern Ocean is responsible for the largest discrepancy with observations for both fields.However, most of the CMIP5 models examined exhibit larger discrepancies, most notably in the North Atlantic, where concentrations of both tracers are elevated by more than 100 mmol m −3 (or meq −3 ) in several models.Nonetheless, for both fields, it is noticeable that the GFDL-ESM2M model (Dunne et al., 2013) performs best.In passing, that observed fields are fairly consistently represented better or worse across different models is suggestive of a gradation of difficulty in modelling aspects of oceanic biogeochemistry.Bulk properties, such as nutrients, which place a strong constraint on all biological activity, and for which there are relatively good and widespread measurements, are tightly constrained, whereas primary production, for which there are global estimates but far fewer direct observations, is less well-constrained.Chlorophyll, by contrast, is something of a confounding factor here, since its surface values have become extremely well-observed since the beginning of the satellite era.This unexpected disparity likely stems from the strong plasticity of chlorophyll : C ratios in phytoplankton, which acts to decouple its concentrations both from well-observed ambient factors such as nutrients and more poorly observed aspects such as phytoplankton biomass.

Discussion
Despite the ever-increasing processing power of supercomputers, incorporating marine ecosystem models into global GCMs is a computationally expensive business when it comes to undertaking climate simulations, especially if high resolution is desired and/or the ocean is coupled to an atmospheric model.MEDUSA-1.0,the precursor to the version of the model described herein, was explicitly developed with this consideration in mind as an "intermediate complexity" plankton ecosystem model for global biogeochemical modelling.It simulates primary production, grazing and export of detritus to the deep ocean, the sinking particles containing both organic and inorganic C with the latter via a latitudinally dependent "rain ratio".The base currency of MEDUSA-1.0 is nitrogen, the selection of a nutrient element (N or P) for this purpose being a necessity given the role of nutrients in limiting primary production in the ocean.Although MEDUSA-1.0predicts the fluxes and cycling of organic carbon, simulating the full carbon cycle in the ocean, including ventilation of CO 2 with the atmosphere and the resulting impact of ocean acidification on marine ecosystems, requires additional tracers.Here, we describe MEDUSA-2.0,an expanded successor model which includes dissolved inorganic carbon, alkalinity, dissolved oxygen and detrital carbon as additional state variables, as well as a simple representation of the benthos.
In principle, the two versions of the model ought to give similar predictions given that phytoplankton, at the base of the food chain, are not limited by the availability of dissolved inorganic carbon (at least in the model -work such as Riebesell et al., 2007, suggests that the availability of DIC species may actually affect carbon fixation).Differences do occur, however, for several reasons.Principally, the simulation described here used forcing from output of the HadGEM2-ES coupled model and was run for 145 yr , rather than observationally derived DFS4.1 forcing and a simulation length of only 40 yr  as used with MEDUSA-1.0(Yool et al., 2011).There were also minor parameter tweaks to adjust near-surface nutrients (see later), as well as parameter changes and additions to accommodate the carbon and oxygen cycles.
However, in general, the performance of the two models is very similar, in both cases successfully reproducing major features such as the oligotrophic gyres and the seasonal progression of plankton blooms at high latitudes.At the global scale, predicted primary production of 45.3 and 41.6 Pg C yr −1 for MEDUSA versions 1.0 and 2.0 respectively, are in line with, although slightly below, observationally derived estimates of 46.3-60.4Pg C yr −1 (Behrenfeld and Falkowski, 1997; Carr et al., 2006;Westberry et al., 2008).The lower primary production in MEDUSA-2.0 is in    part a consequence of changes to productivity in the Pacific that result from the switch in MEDUSA-2.0 to a more modern aeolian deposition field (Mahowald, 2005).The resulting deficiency of iron leads, in part, to excess DIN in the surface waters of the equatorial Pacific, though this is convoluted with excessively shallow remineralisation in this region which more generally increases near-surface concentrations of DIN (and DIC).Both models do a reasonable job at capturing the spatial and seasonal patterns of productivity, although various discrepancies with observations are seen including lower primary production in the subtropical gyres and elevated productivity in iron limited high-nutrient-lowchlorophyll (HNLC) regions including the Southern Ocean, equatorial Pacific and subarctic North Pacific.Predicted concentrations of DIN and, especially, silicic acid are too high in the Southern Ocean, a result of excessive ventilation in this basin which acts to homogenise horizontal and vertical gradients.The problem is somewhat worse in MEDUSA-2.0because ocean circulation in this region is strong under the HadGEM2-ES forcing used here.Note, however, that the longer duration of the MEDUSA-2.0simulation allows for any deficiencies in either physics or biogeochemistry to more obviously manifest themselves.
A new feature of MEDUSA-2.0 is the inclusion of a simple benthic model.This serves as a series of four reservoirs for detrital material (slow-and fast-sinking) that reaches the seafloor -nitrogen, silicon, organic carbon, CaCO 3 ; but not iron, which is coupled to nitrogen.In MEDUSA-1.0such material was instantaneously remineralised (or dissolved) upon reaching the seafloor.While this latter, simplistic approach has limited consequences in the deep ocean where the recycled dissolved inorganic nutrients cannot be consumed by phytoplankton growth, in shallower regions such as the shelves it has the potential to unrealistically enhance production.Patterns in the supply of organic matter to the seafloor closely mirror those of primary production in the surface ocean in the model.This supply thus shows strong seasonality at high latitudes and low seasonality in the tropics.However, given the turnover of sinking particles as they descend through the water column, the magnitude of benthic supply is closely tied to seafloor depth.
As was the case with MEDUSA-1.0, the modelling of iron is still problematic.Aeolian deposition balances uneasily with scavenging, with the result that iron distributions diverge from those of the initial condition (admittedly modelderived; Dutkiewicz et al., 2005) at both the surface and, especially, at depth.The latter discrepancy has limited impact on the simulations here and in Yool et al. (2011) but it does illuminate gaps in understanding of this elemental cycle.While understanding of iron in the ocean has progressed in recent years (e.g.Boyd and Ellwood, 2010;Breitbarth et al., 2010), accurately representing iron in ecosystem models remains difficult for a number of reasons.For instance, accurate estimation of the iron supply to the ocean is hampered by our ignorance of both dust supply and dust solubility once in the ocean (Schulz et al., 2012).Furthermore, even once in the ocean, iron's bioavailability is influenced by its various speciation and redox states, biological cycling and the various uptake strategies of phytoplankton and bacteria.And though increasingly complex representations of iron are being developed and incorporated into ecosystem models (e.g.Weber et al., 2007;Ye et al., 2009), the current generation of ocean biogeochemical GCMs typically only include a single iron pool and so cannot account for the roles of ligand complexation and nonbiological processes (light and temperature) in controlling bioavailable Fe and therefore the extent of phytoplankton limitation (Tagliabue et al., 2009).
Moving on to the carbon cycle, predicted patterns of pCO 2 and air-sea CO 2 exchange throughout the world ocean generally compare favourably to maps based on observations (Takahashi et al., 2009).Some areas, such as the northeast leading to a diverse range of approaches in models (e.g.Tyrrell and Taylor, 1996;Moore et al., 2002;Gehlen et al., 2007;Ridgwell et al., 2007;Gregg and Casey, 2007;Zahariev et al., 2008;Yool et al., 2010).Anderson (2005) used calcifiers as an example of how difficult it is to reliably parameterize complex models for use in forecast projections.For example, the ecology of coccolithophores is poorly understood including the relative roles of bottom-up (via different nutrients) and top-down (grazing, viral lysis) controls on their dynamics.Further, calcifiers are a diverse group of organisms, including hundreds of species of coccolithophores, as well as foraminiferans and pteropods.Grouping them into a single model state variable and then parameterizing based on, for example, the well-known species Emiliania huxleyi is a potentially hazardous strategy.Ecosystem models used in global biogeochemical modelling studies have therefore adopted relatively simple approaches to the representation of calcification.
In MEDUSA-1.0,where nutrient cycles were of greater concern, the rain ratio of CaCO 3 : C org was made a simple empirical function of latitude (following Dunne et al., 2007).This approach captured some of the first order features of the rain ratio (e.g.equator-pole gradients) but prevented any sensitivity to physico-chemical changes (though changes in productivity would still impact the absolute quantity of calcification).For MEDUSA-2.0, the parameterization of calcification was therefore improved to permit dynamic change under the influence of ambient marine chemistry (Riebesell et al., 2000;Zondervan et al., 2001).As noted above, there are good reasons why a representation of a CaCO 3 production via a dedicated state variable ("coccolithophorid phytoplankton", "pteropod zooplankton") may be problematic.To this end, MEDUSA-2.0adopts a calcification parameterization which straightforwardly replaces that in MEDUSA-1.0,and which was developed for, and optimised to, the global scale (Ridgwell et al., 2007).Though developed within the framework of a low resolution Earth System Model, GE-NIE, and coupled to a simple "nutrient-restoring" biogeochemical framework, this parameterization serves the same purpose there as in MEDUSA-2.0-the production of exported CaCO 3 .Of course, the relationship that it assumes between calcite and CaCO 3 export is known to be diverse (e.g.Buitenhuis et al., 1999;Iglesias-Rodriguez et al., 2008;Langer et al., 2006), but it serves here as an obvious stepping stone in complexity for MEDUSA-2.0,and the potential impacts of adopting it are explored in a separate study (Yool et al., 2013).
In terms of future developments for MEDUSA-2.0,a number of avenues suggest themselves.The performance of the chlorophyll submodel remains somewhat problematic (though see the comparable CMIP5 results), with the model failing to simulate spring bloom concentrations as high as those observed, while having much lower concentrations in the unproductive oligotrophic gyres.The literature contains more sophisticated treatments of phytoplankton physiology than that used here (up to long-standing submodels such as Flynn, 2001) and the adoption of such a submodel may improve this aspect of MEDUSA-2.0.With the inclusion of the oxygen cycle, and the simulation of suboxic regions, MEDUSA-2.0'somission of denitrification could also be addressed (Deutsch et al., 2007).At the other end of the nitrogen cycle, the factors regulating the distribution of nitrogen fixation are increasingly well-understood (e.g.Moore and Doney, 2007;Monteiro and Follows, 2012), and this process both interacts with denitrification (Deutsch et al., 2007;Fernandez et al., 2011) and is expected to change into the future (Levitan et al., 2007;Barcelos e Ramos et al., 2007).It should be noted that such changes would have additional implications beyond the nitrogen and oxygen cycles including, for instance, impacts in the distribution of ocean alkalinity (Wolf-Gladrow et al., 2007).Though much has been made above of MEDUSA-2.0'smore sophisticated treatment of CaCO 3 , it is clear that this remains just one "solution" for this aspect of the ocean's carbon (and alkalinity) cycle.The broad range of ongoing research into the impacts of ocean acidification on calcifiers will continue to inform the modelling of CaCO 3 , and will hopefully provide a more "universal" understanding and formulation -for instance, a consensus on the ecophysiological factors that govern calcifier abundance.On a related point, the role played by CaCO 3 in the export of organic material to the deep ocean has been questioned (Passow and De La Rocha, 2006;Wilson et al., 2012), and MEDUSA-2.0'sutilisation of the ballast hypothesis may require revisiting.And there are further omissions of MEDUSA-2.0,less immediately pressing, that could be considered.For example, CO 2 -enhanced carbon fixation (Riebesell et al., 2007) or DOM production (Engel, 2002), a more thorough treatment of elemental ratios (Burkhardt et al., 1999), the importance of food quality in grazing interactions (Mitra and Flynn, 2005) or phytoplankton mixotrophy (Hartmann et al., 2012).
Another route to model improvement lies with parameterization.As remarked upon earlier, MEDUSA-2.0features several minor parameter changes to amend problems that appeared following its evolution from MEDUSA-1.0.One set of such changes aimed to increase nutrient retention in near-surface waters by decreasing the export of sinking detritus (Sect.2.4).However, as noted subsequently (Sect.4), MEDUSA-2.0'smidwater distributions of nutrients and oxygen in the Pacific, as well as its CCD depth, feature errors that arise from remineralisation of organic material occurring at depths that are too shallow; errors that, in part, are caused by this very "fix".These, and other, changes were made to MEDUSA-2.0'sparameters in a relatively ad hoc fashion, being both "best guess" solutions to (usually) nutrient distribution problems, and were "resolved" by short duration test runs.This inefficient approach is driven by the expense (in compute and wall-clock time) of 3-D simulations, which favours limited numbers of limited duration "nudges" to resolve errors.Limiting the spatial domain of simulations to 0-D (Fasham and Evans, 1995) or 1-D (Schartau and Oschlies, 2003) does permit a more time-efficient "solution" but at the expense of critically neglecting the horizontal transports that are convoluted with biogeochemical processes to create important gradients throughout the ocean (but see Hemmings and Challenor, 2012).However, in recent years, techniques such as that of Khatiwala (2007) have been developed to permit efficient, offline simulation of ocean biogeochemistry in 3-D, and recent studies such as Kriest et al. (2010Kriest et al. ( , 2012) ) have used this approach to assess the performance and evaluate the parameter sensitivity of simple NPZD models.Development of such approaches offers the potential in the future of optimising the parameterization of models like MEDUSA-2.0,such that the full consequences of changes to parameters (or even model structure) can be properly investigated and accounted for.
However, notwithstanding the considerable room for improvement -or expansion -outlined above, the further development of MEDUSA-2.0runs counter to the stated intention that the model occupies the "intermediate complexity" niche of biogeochemical modelling.Furthermore, while potentially extending the reach -and utility -of MEDUSA-2.0 on several fronts, they present no method for expanding the model in a systematic or quasi-objective fashion.Piecemeal additions to model complexity, however warranted and justifiable, run the risk of creating a succession of "specialist models", that while individually useful may not sit within a consistent hierarchy of complexity.As such, it may be difficult to fine-tune the biogeochemical complexity to suit a particular task (with particular resources) to hand.
Nonetheless, despite the limitations outlined above, MEDUSA-2.0still represents an efficiently sized tool for realistically simulating the ocean's major biogeochemical cycles.
6 Conclusions -MEDUSA-2.0builds traceably on MEDUSA-1.0by adding carbon, alkalinity and oxygen cycles, a simple benthos submodel and options for CaCO 3 production and export remineralisation.
-Calcification submodel permits dynamic response to ambient seawater chemistry allowing investigation of ocean acidification feedbacks at an appropriate level of additional complexity.
-MEDUSA-2.0performace evaluated at the global scale using observational nutrient, chlorophyll and carbon cycle fields following a century-scale simulation .
-Similarly to its predecessor model, MEDUSA-2.0has excessive nutrient concentrations in Southern Ocean and low chlorophyll and productivity in oligotrophic gyres.
-Global productivity is slightly lower than in MEDUSA-1.0, in part due to changes in aeolian iron deposition that decrease Pacific productivity and increase excess surface DIN.
-Excessive circulation-driven ventilation and tooshallow remineralisation of sinking organic material introduce discrepancies to interior concentrations of biogeochemical tracers.
-MEDUSA-2.0 has generally good agreement on surface carbon cycle properties ( pCO 2 and air-sea flux) and CaCO 3 production within observational range, but discrepancies created in CCD field by too-shallow remineralisation.
-Though not without discrepancies, intercomparison with CMIP5 models of similar (or greater) complexity places MEDUSA-2.0'sperformance among the best of those examined. A.
-trcoxy_medusa.F90 this routine is called by trcbio_medusa.F90 to perform calculations associated with saturation concentration and air-sea O 2 flux; while modified to interface with MEDUSA-2.0, it is derived from Najjar and Orr (1999) As with MEDUSA-1.0, the above routines are included in the Supplement that accompanies this article.

Appendix B CMIP5 intercomparison
As part of the evaluation of MEDUSA-2.0,an intercomparison was performed for a number of surface ocean biogeochemical properties against a selection of nine models drawn from the 5th round of the Coupled Model Intercomparison Project (CMIP5; Taylor et al., 2012).Output from these models was accessed and downloaded via the CMIP5/OCMIP5 FileFinderAR5 Web Application (Brockmann, 2012; http://ocmip5.ipsl.fr/FileFinderAR5/).
The fields downloaded were global-scale annual averages of surface dissolved inorganic nitrogen (= nitrate + ammonium for some CMIP5 models), surface silicic acid, surface chlorophyll, vertically integrated primary production, surface dissolved inorganic carbon and surface alkalinity.In each case, the output used in the intercomparison was generated from model output covering the same temporal periods described in 4. Note that not all of the models selected were able to supply all of the fields examined (NorESM1, chlorophyll; GISS-E2, alkalinity).
The CMIP5 models used in the intercomparison, their underlying ocean biogeochemistry submodel, together with the groups that performed the simulations and provided output, are listed in Table 9.The bracketed numbers indicate the total number of tracers in each model and serves as a proxy of model complexity and computational burden (MEDUSA-2.0has 15 tracers).Edited by: R. Marsh

Fig. 1 .
Fig. 1.Schematic diagram of the components and interactions in the MEDUSA-2.0model.Boxes with solid borders indicate explicitly modelled state variables, while boxes with dashed borders indicate implicitly modelled components.Overlapping boxes indicate components for which multiple currencies are modelled (e.g.different elements, chlorophyll).The smaller boxes at the bottom of the diagram refer to benthic reservoirs of model currencies that are fed by sinking detrital material (slow-and fast-sinking).For reasons of diagrammatic clarity, dissolved oxygen and its connections to other state variables are omitted here.Note that the dissolution of benthic CaCO 3 releases both DIC and alkalinity.
its fate remains bound to that of nitrogen.The full list of 2-D state variables represented are:B N Benthic organic nitrogen mmol N m −2 B C Benthic organic carbon mmol C m −2 B Si Benthic inorganic silicon mmol Si m −2 B Ca Benthic inorganic CaCO 3 mmol C m −2 Fig. 2. Diatom uptake of nitrogen (top) and silicon (bottom) against the Si : N ratio of diatom biomass.

Fig. 2 .
Fig. 2.The top panel shows mean annual aeolian iron input to the ocean (i.e. the quantity of iron that dissolves into seawater from deposited dust).The input is shown on a logarithmic scale in units of µmol m −2 y −1 , and integrated input is 2.564 Gmol Fe y −1 .The bottom panel shows the fractionation of total iron between "free" and ligand-bound forms across a logarithmic range of total iron concentrations.
curve uses a fixed value of −0.858 for parameter b in Eq. (105).Using a more modern data set of thorium-derived POC export, Henson et al. (2012) developed a variant scheme in which parameter b is instead a function of local surface temperature.b = −1.06+ (0.024 • T )

Fig. 4 .
Fig.4.Observational (World Ocean Atlas, 2009; left)  and simulated (right) surface dissolved inorganic nitrogen for northern summer (June-July-August; top) and northern winter (December-January-February; bottom).Concentrations in mmol m −3 .In passing, note that in this and subsequent geographical plots, the Mollweide equal area projection has been preferred in order that ocean regions are presented without undue emphasis.

Fig. 12 .
Fig. 12. Observational (left) and simulated (right) integrated primary production for northern summer (June-July-August; top) and northern winter (December-January-February; bottom).The observational field shown here is an average of the VGPM, Eppley-VGPM and CbPM estimates.Production in g C m −2 d −1 .

Fig. 13 .Fig. 13 .
Fig. 13.Observational (left) and simulated (right) integrated primary production for northern summer (June-July-August; top) and northern winter (December-January-February; bottom).The observational field shown here is an average of the VGPM, Eppley-VGPM and CbPM estimates.Production in g C m −2 d −1 .

Fig. 14 .Fig. 14 .
Fig. 14.Taylor diagrams of spatial (left) and temporal (right) model-observation comparisons for surface dissolved inorganic nitrogen (top) and silicic acid (bottom).In the leftmost panels, simulated annual means for different regions are compared to corresponding observational fields.In the rightmost panels, simulated global average means for different months are compared to corresponding observational fields.

Fig. 15 .
Fig. 15.Taylor diagrams of spatial (left) and temporal (right) model-observation comparisons for surface chlorophyll (top) and integrated primary production (bottom).In the leftmost panels, simulated annual means for different regions are compared to corresponding observational fields.In the rightmost panels, simulated global average means for different months are compared to corresponding observational fields.

Fig. 15 .
Fig. 15.Taylor diagrams of spatial (left) and temporal (right) model-observation comparisons for surface chlorophyll (top) and integrated primary production (bottom).In the leftmost panels, simulated annual means for different regions are compared to corresponding observational fields.In the rightmost panels, simulated global average means for different months are compared to corresponding observational fields.

Fig. 17 .
Fig. 17.Intercomparison of observational (left) and model (right) fields of zonally averaged dissolved inorganic carbon for the Atlantic (top) and Pacific (bottom) basins.Note that the GLODAP climatology has no values within the Arctic.Concentrations in mmol C m −3 .

Fig. 18 .
Fig. 18.Intercomparison of observational (left) and model (right) fields of zonally averaged alkalinity for the Atlantic (top) and Pacific (bottom) basins.Note that the GLODAP climatology has no values within the Arctic.Concentrations in meq m −3 .

Fig. 20 .Fig. 19 .
Fig. 20.Calcite compensation depth (CCD) as calculated from observations (left) and MEDUSA-2.0(right).The CCD is defined here as the depth at which carbonate ion concentration falls below the local saturation concentration, that is, where calcite falls below a value of 1. l et al.: A description of MEDUSA-2.033

Fig. 29 .
Fig. 29.Simulated mixed layer primary production fraction (left) and seafloor detrital flux (right) for northern summer (June-July-August; top) and northern winter (December-January-February; bottom).Production fraction is dimensionless; seafloor detrital flux in mg C m −2 d −1 , and shown on a logarithmic scale.

Fig. 31 .
Fig. 31.Simulated vertical profiles of dissolved inorganic nitrogen concentration averaged for the World Ocean (top left) and 5 major regions.Concentrations in mmol N m −3 .Note that depth is shown on a logarithmic scale.

Fig. 37 .
Fig. 37. Simulated vertical profiles of dissolved inorganic nitrogen concentration averaged for the World Ocean (top left) and 5 major regions.Concentrations in mmol N m −3 .Note that depth is shown on a logarithmic scale.

Fig. 32 .
Fig. 32.Simulated vertical profiles of silicic acid concentration averaged for the World Ocean (top left) and 5 major regions.Concentrations in mmol Si m −3 .Note that depth is shown on a logarithmic scale.

Fig. 38 .
Fig. 38.Simulated vertical profiles of silicic acid concentration averaged for the World Ocean (top left) and 5 major regions.Concentrations in mmol Si m −3 .Note that depth is shown on a logarithmic scale.

Fig. 33 .
Fig. 33.Simulated vertical profiles of iron concentration averaged for the World Ocean (top left) and 5 major regions.Concentrations in µmol Fe m −3 .Note that depth is shown on a logarithmic scale.

Fig. 39 .
Fig. 39.Simulated vertical profiles of iron concentration averaged for the World Ocean (top left) and 5 major regions.Concentrations in µmol Fe m −3 .Note that depth is shown on a logarithmic scale.

[Fig. 34 .
Fig. 34.Globally averaged surface dissolved inorganic nitrogen (top left), surface silicic acid (top right), surface chlorophyll (bottom left) and integrated primary production (bottom right).Solid black lines are annual averages/integral; individual points are individual months.Note that individual monthly primary production values have been normalised so that they appear on the same scale as annual integrals.

Fig. 40 .
Fig. 40.Globally averaged surface dissolved inorganic nitrogen (top left), surface silicic acid (top right), surface chlorophyll (bottom left) and integrated primary production (bottom right).Solid black lines are annual averages/integral; individual points are individual months.Note that individual monthly primary production values have been normalised so that they appear on the same scale as annual integrals.

Fig. 41 .Fig. 41 .
Fig. 41.Taylor diagrams of the performance of MEDUSA-2.0and a selection of CMIP5 models for a range of biogeochemical properties: dissolved inorganic nitrogen (top left), silicic acid (top right), chlorophyll (middle left), primary production (middle right), dissolved inorganic carbon (bottom left) and alkalinity (bottom right).In each case, diagrams show spatial model-observation comparisons based on annual average fields.All diagrams share a common model key.

Fig. 42 .
Fig. 42.Intercomparison of annual mean surface dissolved inorganic nitrogen for MEDUSA-2.0and a range of comparable CMIP5 models.DIN concentration in mmol N m −3 .

Fig. 44 .
Fig. 44.Intercomparison of annual mean surface chlorophyll for MEDUSA-2.0and a range of comparable CMIP5 models.Chlorophyll concentration in mg chl m −3 and is shown on a logarithmic scale.

Table 3 .
Plankton and detritus loss parameters.

Table 9 .
CMIP5 models used for intercomparison with MEDUSA-2.0.Each entry lists the component biogeochemical submodel, the number of state variables (in brackets) and the host institution for each model.

Yool et al.: A description of MEDUSA-2.0 -trclsm_medusa
.F90 this routine initialises the parameters to the values specified in namelist.trc.sms-trcsms_medusa.F90 this routine is called by the NEMO model during a simulation and in turn calls the MEDUSA-2.0routines that calculate biogeochemical sources and sinks -trcopt_medusa.F90 this routine calculates the submarine light field -trcbio_medusa.F90 this is the main model routine and includes (almost) all of the ecosystem equations used for the biogeochemical sources and sinks for tracers -trcsed_medusa.F90 this routine both initialises the aeolian iron deposition and calcite CCD fields (if required) and (for historical reasons) calculates the sinking of the slow detritus tracer -trcco2_medusa.F90 this routine is called by trcbio_medusa.F90 to perform calculations associated with carbonate chemistry and air-sea CO 2 flux; while modified to interface with MEDUSA-2.0, it is derived from Blackford et al.

to this article is available online at http://www.geosci-model-dev.net/6/ 1767/2013/gmd-6-1767-2013-supplement.zip.
's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.