A new marine ecosystem model for the University of Victoria Earth system climate model

Introduction Conclusions References


Introduction
The oceans have a large effect on the Earth's climate and play an important role in global biogeochemical cycles.Although many of the interactions within the oceans and between the oceans and other systems (i.e., atmosphere, land, etc.) are driven purely by physical or chemical mechanisms, biolog-ical processes play an important, but often less understood role.Marine ecosystems affect the climate primarily through the "carbonate" and "soft tissue" pumps (i.e., the "biological" pump) (Longhurst and Harrison, 1989;Volk and Hoffert, 1985).These pumps work through the biological uptake of carbon in the surface ocean and the subsequent transport (mainly by sinking and zooplankton grazing dynamics) of the small fraction of it that is not recycled along the way (mostly through respiratory processes) to the deep ocean, where it is unable to affect the climate (i.e., as CO 2 ) for hundreds to thousands of years.The biological pump has been estimated to export between 5 and 20 Gt C yr −1 out of the surface layer (Henson et al., 2011;Honjo et al., 2008;Laws et al., 2000).However, as indicated by the large range of estimates, there is great uncertainty in our understanding of the magnitude of carbon export (Henson et al., 2011;Oschlies, 2001) and, thus, it's effect on the Earth's climate.
Marine ecosystems effect global biogeochemical cycles in a number of ways.As mentioned above, they play an important role in the carbon cycle through the uptake, recycling and sequestration of some CO 2 .They also play a large role in the cycling of nitrogen, phosphorus and oxygen.In surface waters nitrogen and phosphorus are major nutrients that are consumed by, and drive, primary production (PP) and, thus, are linked to the carbon cycle.Since these nutrients often limit PP their availability can, therefore, influence the climate system by controlling the magnitude and location of biological pump processes.When heterotrophic members of the ecosystem (i.e., bacteria, zooplankton, etc.) consume the organic material that was fixed during PP, biogeochemical cycles are further effected as oxygen is consumed during respiration and new C-, N-, and P-containing chemical compounds are formed from the processed organic matter.Biological oxygen consumption can have a large impact on the ocean's oxygen concentration when it takes place away from the surface and under some physical conditions may cause areas of the ocean to become hypoxic or anoxic.Most of the organic C, N and P is eventually recycled by heterotrophs back into an inorganic form and once again becomes available for uptake by primary producers when it reaches the surface ocean (if it isn't already there).However, some N may be returned to the atmosphere as N 2 or N 2 O when microorganisms in hypoxic/anoxic waters utilise N containing compounds.This loss of N to the atmosphere is not one-way as nitrogen can also be removed from the atmosphere and added to the ocean in bioavailable form by nitrogen-fixing organisms.Of course these biotically driven cycles do not take place without interactions between chemically and physically driven C, N and P cycles and there are frequent exchanges of material between these cycles.Although the major biogeochemical pathways are known, there is much uncertainty in the magnitude and time-scale of each pathway and the biological communities that drive them.
Ocean biology and the effects that it has on important processes are changing because of naturally and anthropogenically driven environmental change.Increasing anthropogenic greenhouse gas emissions have, and are, altering the Earth's climate and warming the planet (IPCC, 2007;Sarmiento and Gruber, 2002).The oceans have taken up much of the heat that has accumulated in the Earth system and near-surface ocean waters have already warmed between 0.5 and 0.7 degrees over the last 100 yr (Hansen et al., 2010;IPCC, 2007).Since most physical, chemical and biological processes are temperature dependent, these changes have had an effect on many marine ecosystem processes.This has, in turn, affected the ocean's role in the Earth's climate system and global biogeochemical cycles (often through feedback effects).The magnitude of the biological pump is predicted to decline in response to global warming, resulting in reduced carbon sequestration and higher atmospheric CO 2 levels (Sarmiento et al., 1998).Given the magnitude of this carbon sink and its effect on the climate, accurately quantifying and understanding how it may change due to anthropogenic activities is important.Furthermore, since manipulation of the biological pump is being considered as a means of climate engineering (Keith, 2000;Oschlies et al., 2010a), it is essential to have a thorough understanding of potential effects and side effects.Temperature changes are also effecting the ocean's capacity to hold dissolved oxygen and causing the global volume of oxygen held by the ocean to decrease (i.e., ocean deoxygenation) (Gruber, 2011), which may effect biogeochemical cycles if more areas become hypoxic or anoxic.The effects of anthropogenic activities on ocean biology are not limited to temperature since CO 2 emissions are having an effect on ocean chemistry and making marine waters more acidic (i.e., ocean acidification) (Hofmann and Schellnhuber, 2010), to the detriment of some species.Furthermore, anthropogenic activities are causing many coastal waters to become eutrophic and are increasing remote open ocean nutrient concentrations (Doney, 2010), which has a fertilisation effect on marine ecosystems.Since marine ecosystems are a key component of oceanic systems and play a role in other Earth system processes it is essential to understand how they function and how they may function in the future as the world changes.
Many approaches have been used to study the role of marine ecosystems in the Earth's climate system and global biogeochemical cycles.Models are one of the few tools available for understanding these dynamics on a global basis and for predicting how they may change in the future as a result of anthropogenic influences.Here we describe improvements made to the marine ecosystem component of an Earth system model of intermediate complexity (EMICs).EMICs were developed roughly two decades ago for a broad spectrum of purposes (Claussen et al., 2002) and could still be considered as being in the early stages of development.Since these models are much more complex than ocean models without explicit atmospheric or terrestrial components, the marine ecosystems described by them are usually much more simplistic when compared to current state-of-the-art marine ecosystem models that now include multiple "functional groups" or "biogeochemical guilds" of organisms and biogeochemical tracers (Hood et al., 2006).While adding multiple functional groups or additional biogeochemical tracers may not be computationally feasible, or likely to improve model skill (Friedrichs et al., 2007;Kriest et al., 2010), there are many improvements that can be made to EMIC marine ecosystem models to improve their performance and realism.Our improvements to the marine ecosystem component of the University of Victoria Earth System Climate Model (UVic ESCM) are designed to correct previous errors, add more realistic processes, and provide a basic platform for future applications.

A brief description of the UVic ESCM
The UVic Earth System Climate Model (Eby et al., 2009) version 2.9 consists of three components: a three dimensional general circulation model of the ocean, a terrestrial model and a simple one layer atmospheric energy-moisture balance model (Weaver et al., 2001).All model components use a common horizontal resolution of 1.8 • latitude × 3.6 • longitude and the oceanic component has nineteen levels in the vertical with thicknesses ranging from 50 m near the surface to 500 m in the deep ocean.The oceanic circulation model (Modular Ocean Model 2) includes physical parameterisations for diffusive mixing along and across isopycnals, eddy induced tracer advection (Gent and McWilliams, 1990), and, in our configuration, a scheme for the computation of tidally induced diapycnal mixing over rough topography (Simmons et al., 2004).The atmospheric energy-moisture balance model interactively calculates heat and water fluxes to the ocean, land and sea ice.Wind velocities to calculate the momentum transfer to the ocean and to a dynamicthermodynamic sea ice model, surface heat and water fluxes, and the advection of water vapour in the atmosphere are prescribed from NCAR/NCEP monthly climatology data.The terrestrial model of vegetation and carbon cycles (Meissner et al., 2003) is based on the Hadley Center model TRIFFID.Continental ice sheets are assumed to remain constant.The simulation of sea ice has been evaluated and found to be in good agreement with observations (Bitz et al., 2001;Saenko et al., 2002).

Rationale for improving the marine ecosystem model
The most recent UVic ESCM marine ecosystem module of Schmittner et al. (2008) has been used extensively in biogeochemical studies of climate change (Oschlies et al., 2008;Schmittner et al., 2009b), paleooceanography (Schmittner and Galbraith, 2008), climate engineering (Oschlies, 2009;Oschlies et al., 2010a, b), and the nitrogen cycle (Somes et al., 2010a, b).These studies have provided valuable insight on a number of research topics and have helped to answer questions about pressing issues of climate change.However, despite the usefulness of this model and the relatively good skill at which it simulates biogeochemical tracer distributions, there are some aspects of the model that can be improved and some errors that need to be corrected, for future applications.
A first set of errors in the original Schmittner et al. (2008) model had already been identified and corrected by Schmittner et al. (2009a).Additional analysis of various aspects of the model results led to the identification of a few other previously overlooked model errors and aspects that could be improved.Originally, we had focused solely on correcting these errors and improving the equations and code describing light limitation of phytoplankton growth and zooplankton grazing.The improvements made to the equations/code describing light limitation of phytoplankton growth were made because light attenuation due to diazotroph biomass was not included in them.Since some species such as Trichodesmium can form dense blooms at the surface that can significantly attenuate light (Capone and Zehr, 1997;Sellner, 1997) it was viewed as important to simulate this effect on the phytoplankton community.In addition, there was an error in the code (in both versions 2.8 and 2.9) describing light attenuation in the layers below the surface.This error occurred because the biomass of diazotrophs from the layers above was not used in calculations of light attenuation at depth.A second set of changes was made to the equations/code describing zooplankton grazing on phytoplankton and diazotrophs.The original equations were not formulated to provide a multiple-prey functional response, making grazing unrealistic when both types of plankton ("ordinary" phy-toplankton and diazotrophs) were present (i.e., ingestion was incorrectly calculated).Correctly modelling grazing on multiple prey items is important for realistically simulating zooplankton growth and the different components of their diet (Anderson et al., 2010;Gentleman et al., 2003).While correcting zooplankton grazing, we also decided to make their growth and grazing more realistic.However, these changes, along with conclusions drawn from recent studies on the role of iron in phytoplankton growth, required us to also modify phytoplankton growth and mortality in the model.
In Schmittner et al. (2008) zooplankton were limited to grazing on phytoplankton and diazotrophs, excretion was not a function of grazing, and their growth rate was not a function of temperature.Since the zooplankton community grazes on itself and detritus, in addition to phytoplankton (Calbet and Saiz, 2005;Dilling and Brzezinski, 2004;Kleppel, 1993;Sherr and Sherr, 2002), we decided to represent these food web pathways in the model.Furthermore, since the excretion of inorganic nitrogen by zooplankton is related to their ingestion and diet (Davidson et al., 1995;Glibert et al., 1992;Saba et al., 2009), we feel that it is important to allow excretion to be a function of grazing, even if our current formulation is simplistic due to Redfield stoichiometry and the desire to keep computational costs low (i.e., no complex formulations as in Anderson, 1992 or Pahlow andProwe, 2010).Finally, since the growth of a major prey item, phytoplankton, in the model is temperature dependent and zooplankton growth is also known to be dependent on temperature (Hirst and Bunker, 2010), we feel that it is important to allow zooplankton growth to change with temperature.However, making zooplankton growth and grazing a function of temperature had dramatic effects on phytoplankton growth as it was formulated in Schmittner et al. (2008).
In Schmittner et al. (2008) phytoplankton growth at high latitudes in cold water, especially in the Southern Ocean, had been controlled by unrealistically high rates of grazing due to the fixed (temperature independent) growth rate of zooplankton (i.e., this fixed value while reasonable for warmer waters was too high for colder waters).When this top-down control was reduced by our modifications (i.e., making zooplankton growth a function of temperature), phytoplankton biomass and productivity became very high in these nutrient rich waters.Since phytoplankton biomass and productivity are known to be low in these regions due to iron limitation (Martin, 1992;Martin et al., 1991), it was necessary to include some form of iron limitation in the model.Additionally, by including the effects of iron limitation we hoped to improve the spatial distribution of diazotrophs, whose growth is also known to be limited by iron availability (Berman-Frank et al., 2001, 2007).Unfortunately, iron cycling and the interactions between iron and biology are quite complex and not fully understood, making it difficult to model (Galbraith et al., 2010;Moore and Braucher, 2008;Moore et al., 2004;Tagliabue et al., 2009).Somes et al. (2010a, b) successfully used an iron mask, based on measured atmospheric dust deposition, in a nitrogen isotope study with the UVic 2.8 model to constrain diazotroph growth and achieve a more reasonable diazotroph spatial distribution than in Schmittner et al. (2008).Following this simple approach, we also elected to use an iron mask to constrain the growth of both diazotrophic and non-diazotrophic phytoplankton.However, unlike in Somes et al. (2010a, b), the mask that we use is of dissolved iron and, thus, accounts for sources of iron from both aeolian dust sources and sedimentary efflux (Mahowald et al., 2005;Moore and Braucher, 2008).

Configuration of the circulation model
The ocean circulation model described in Sect. 2 and the standard physical settings as set in the version 2.9 download (http://www.climate.uvic.ca/model/)have been modified slightly to have the similar physical dynamics to those in Schmittner et al. (2008).These modifications include turning off the Bryan-Lewis vertical mixing option, turning on the tidal mixing option, increasing the vertical diffusivity parameter in the Southern Ocean, and implementing an anisotropic viscosity scheme in the tropics to improve the simulation of the equatorial currents (see the supplemental mk.in model configuration file).Based on the UVic 2.8 studies by Goes et al. ( 2010) and Schmittner et al. (2009b), the vertical background mixing parameter, K vb , in the Southern Ocean (south of 40 • S) was set to 1.0 cm 2 s −1 in our implementation of UVic 2.9.The sinking of detritus is also different than in Schmittner et al. (2008) as it is not constant below 1000 m, but continues to increase linearly with depth (this is the standard formulation in the downloadable model version).An anisotropic viscosity scheme (Large et al., 2001) is implemented, as in Somes et al. (2010b), to improve equatorial circulation.

New ecosystem model description
As discussed above, the marine ecosystem/biogeochemical model (Fig. 1) is a modified version of the NPZD model of Schmittner et al. (2008).As in the original model, it consists of seven prognostic variables that are embedded within the ocean circulation model described above.The state variables include two phytoplankton classes (nitrogen fixers and other phytoplankton), zooplankton, particulate detritus, nitrate (NO 3 ), phosphate (PO 4 ) and oxygen (O 2 ).Additional biogeochemical tracers include dissolved inorganic carbon (DIC) and alkalinity (ALK).All biological variables and particulate detritus are expressed in units of mmol N m −3 .Constant (∼Redfield) stoichiometry relates the C, N and P content of the biological variables and their exchanges with the inorganic variables (NO 3 , PO 4 , O 2 , ALK, and DIC).Parameters that are new or differ from those of Schmittner et al. (2008)    rameters and variables.The model code is available in the Supplement.
Each variable changes its concentration C according to the following equation where Tr represents all transport terms including advection, isopycnal and diapycnal diffusion, and convection.S denotes the source minus sink terms, which describe the following biogeochemical interactions: Note that in Eq. ( 8) the first term, F sfc , calculates dissolved oxygen exchanges with the atmosphere according to the OCMIP protocol and the second term calculates oxygen production from photosynthesis or consumption due to respiration.The rates at which oxygen production or consumption  occur are equal to the consumption or remineralisation of PO 4 , multiplied by a constant ratio, R O:P , except in suboxic waters (< 5 µM) where oxygen consumption is limited.

Phytoplankton
The maximum potential growth rate of phytoplankton in the top three vertical layers is dependent on both temperature (T ) and dissolved iron (Fe) (illustrated in Fig. 2a).Although the maximum potential growth rate of phytoplankton in models is often calculated using temperature, we include dissolved iron in our calculations because iron is necessary for photosynthesis, the reduction of nitrate to ammonium, and a number of other key cellular processes (Galbraith et al., 2010) (i.e., we assume that iron must be available before photosynthesis or the uptake and utilisation of nitrogen and phosphate is possible).For non-diazotrophic phytoplankton the maximum potential growth rate in these layers is and for diazotrophs it is Temperature is calculated by model and the maximum potential growth rate of non-diazotrophic phytoplankton at 0 • C is set at 0.6 d −1 , based on observations of mixed phytoplankton community growth rates (Le Quéré et al., 2005).The maximum potential growth rate of diazotrophs is set with a handicap, c D , to be 0.4 times that of non-diazotrophic phytoplankton.For the e-folding temperature of biological rates, T b = 15.65 • C (Schmittner et al., 2008), diazotrophs can grow only at temperatures higher than 15 • C. The dissolved iron concentration is determined by a three-vertical layer (∼ the upper 240 m), global mask (Fig. 3) of mean monthly dissolved iron concentration outputs from the BLING model (Galbraith et al., 2010).To facilitate the use of monthly data, daily dissolved iron concentrations are calculated using linear interpolation routines from Press et al. (1992).In addition, since the BLING model was run at a different global resolution than the one here, these differences were accounted for when creating the mask by linearly interpolating BLING output to the UVic ESCM grid.Below the top three vertical layers (i.e., below 240 m depth) where light is likely the most limiting factor of phytoplankton growth, iron limitation is assumed to be negligible and the maximum potential growth rate is dependent only on temperature as in Schmittner et al. (2008).Since this masking approach does not resolve the iron cycle or even allow phytoplankton to take up iron, the half-saturation constants that determine phytoplankton iron limitation were set during the tuning process to achieve a reasonable spatial distribution of phytoplankton and thus, do not reflect any real phytoplankton affinities for this micronutrient.
Once the maximum potential growth rate has been calculated the realized growth rate of "ordinary" non-diazotrophic phytoplankton, J O , is then determined by irradiance (I ) and the concentrations of NO 3 and PO 4 , while the growth rate of diazotrophs, J D , is determined by irradiance (I ) and the concentration of PO 4 .
For both types of phytoplankton light-limited growth, J (O or D)I , is basically calculated as in Schmittner et al. (2005Schmittner et al. ( , 2008) ) using However, these equations now account for the effects of diazotroph biomass on light attenuation.Modifications were also made in the code so that the integrated biomass of phytoplankton and diazotrophs from the layers above is now used in these calculations for depths below the surface.Thus, shortwave radiation at depth z is − 1) ( 14) with I z=0 denoting the downward shortwave radiation at the sea surface, as calculated by the atmospheric model component, and PAR denoting the fraction of that radiation that is photosynthetically active (Schmittner et al., 2005).Lightlimited growth, averaged over depth, with a triangular shaped diurnal cycle becomes with z calculated as in Schmittner et al. (2005Schmittner et al. ( , 2008) ) and for I evaluated at the top of the layer z − z/2, and with d denoting the day length in a fraction of 24 h.Following the correction by Schmittner et al. (2009a) the function is used in Eq. ( 15) to calculate the potential light limited daily averaged growth rate.Non-diazotrophic phytoplankton mortality is no longer set as a quadratic function as in Schmittner et al. (2008).Instead these phytoplankton now die at linear rate of 0.03 d −1 which in combination with the additional loss term representing the temperature-dependant fast remineralisation process, sets their total mortality at ≥ 0.045 d −1 .

Zooplankton
In contrast to the earlier model of Schmittner et al. (2008), zooplankton are now allowed to graze on detritus and other zooplankton (self-predation), in addition phytoplankton and diazotrophs.Grazing varies with prey density and is characterized by a multiple-prey Holling II functional response that assigns preferences, ψ i , for different types of prey (i).The rate of grazing on each type of prey is where and www.geosci-model-dev.net/5/1195/2012/ In these equations the maximum potential growth rate of zooplankton (Eq.28) is dependent on temperature (T ) up to 20 • C, at which point it is capped, as in other models (Anderson et al., 2010), to avoid inordinately high grazing rates in the tropics.Since zooplankton are represented as a single state variable in the model, a maximum grazing rate of 0.4 d −1 at 0 • C was chosen because it is approximately mid-way between the measured growth rates of micro-and meso-zooplankton at 0 • C (0.6 and 0.24 d −1 , respectively) (Le Quéré et al., 2005).Figure 2b illustrates the resulting annually averaged maximum potential zooplankton growth rate in the surface layer.
In these equations zooplankton grazing is also inhibited in hypoxic waters ( 10 µM O 2 ) with grazing ceasing when O 2 decreases below 5 µM.Inhibition of zooplankton grazing in hypoxic and anoxic waters was based on observations that there is little effect of low O 2 levels on the biomass of zooplankton down to ∼ 10 µM, at which point there are pronounced effects on species distributions and biomass (Seibel, 2011).However, since some species can survive and maintain their filtering activity at levels as low as 6 µM O 2 and most euphausiids and copepods avoid only the core of oxygen minimum zones with less than 4.5 µM O 2 (Ekau et al., 2010), our formulation (Eq.28) allows some grazing to occur between 5 and 10 µM O 2 .
The growth of zooplankton (1st term in Eq. 6) is determined by a growth efficiency term, , that was set at 0.4 based on measured (Kiørboe, 1989;Rivkin and Legendre, 2001) and theoretical (Landry and Calbet, 2004;Mitra, 2006) studies of zooplankton growth efficiencies and production.The amount of C that they respire, and N and P that they excrete, is the difference between their assimilation and growth coefficients (γ -).Their assimilation efficiency (1−γ ), where γ has been set to 0.7 based on measured (Fenchel, 1982;Geider and Leadbeater, 1988;Hasegawa et al., 2001;Landry et al., 1984;Pelegrí et al., 1999) and theoretical (Anderson, 1994) studies of zooplankton assimilation and production, also determines the production of detritus from sloppy feeding, egestion, and fecal pellet production.

Other modifications
Formulations of air-sea gas exchange, oxygen cycling, and carbon chemistry are as described in Schmittner et al. (2008), except for a modification to the equation describing the production of CaCO 3 so that it is

Model evaluation
For evaluation, and to generate a starting point for future modeling research, a model spin up of more than 10 000 yr has been performed using preindustrial boundary conditions for insolation and a fixed atmospheric CO 2 concentration of 283 ppm.Since the other model components have previously been evaluated in detail and published elsewhere (Meissner et al., 2003;Schmittner et al., 2005Schmittner et al., , 2008;;Somes et al., 2010b;Weaver et al., 2001) the evaluation below focuses primarily on marine ecosystem dynamics and biogeochemical cycles.To evaluate how our modifications have changed the model dynamics comparisons are made with a 13 000 yr spin up of the Schmittner et al. (2008) model (this spin up is used as the control run in Oschlies et al., 2008), which is hereafter referred to as the "old model".Making these comparisons to the old model is valuable because in addition to evaluating whether the changes have improved the model results, they also provide us with an idea of how the full Earth system model with the new ecosystem model component may compare to other models since the ESCM with the previous ecosystem formulation has been evaluated in several model inter-comparison projects (Archer et al., 2009;Cao et al., 2009).Comparisons are also made with observations from the World Ocean Atlas 2009 (WOA09) database (Antonov et al., 2010;Garcia et al., 2010a, b;Locarini et al., 2010), the Global Data Analysis Project (GLODAP) (Key et al., 2004), and other observational studies.To make the statistical comparisons and model-observation difference plots with WOA09 and GLODAP data the simulated values were regridded onto the observed data grids in Ferret (software version 6.72) using the nearest function (@NRST).Note that since both model spin ups have been performed using preindustrial conditions, comparisons with WOA09 and GLO-DAP data are expected to show some differences as preindustrial forcing will not account for any recent changes in global biogeochemical cycles or climate (i.e., the addition of anthropogenic CO 2 to the atmosphere).However, since there is little preindustrial data available for comparison and the largest anthropogenically caused changes are likely to have occurred in the surface ocean we feel that it is still valuable to make comparisons to these data sets.Moreover, since most Earth system climate models use preindustrial spin ups as a starting point for experiments we feel that it is important to make these comparisons here before climate change or other experiments are conducted with the model.
In addition to making comparisons to annual observations we also examine how well the model simulates the seasonal cycles of plankton and certain biogeochemical tracers.As far as we are aware this type of analysis has not been conducted for any of the previous UVic ESCM marine biogeochemical models.However, it is valuable to conduct these analyses because analyzing only annual averages, as is commonly done for EMICs, may hide systematic deficiencies in simulating the seasonal cycle.Since the seasonal cycle is driven by changes in environmental forcing (light, stratification, temperature, etc.), the ability to correctly simulate the sensitivity of the marine biogeochemistry to this seasonal forcing might be viewed as an indication of whether a model can correctly simulate the sensitivity to other environmental forcing such as those driven by anthropogenic CO 2 emissions.

Annual results
The model simulates the observed present-day distributions of multiple tracers fairly well (Figs. 4,5,6,7,and 8) with most tracer distributions similar to those of the old model.However, when compared to the old model there are some notable differences with the new formulation slightly overestimating nutrients (phosphate and nitrate) and alkalinity in the deep ocean but better simulating oxygen, apparent oxygen utilisation (AOU), DIC, and C14 in many places.Some of these differences are highlighted (Fig. 4 bottom panels) at different depths with a model-data misfit evaluation (Kriest et al., 2010;Eqs. 32 and 34) that shows the absolute difference between the observations and model results.Simulated surface nutrient concentrations compare reasonably well with observed annual surface nutrient concentrations at lower latitudes and the new model resolves major features such as the Eastern Pacific equatorial upwelling region as well as the old model (Fig. 9, PO 4 data not shown).At higher latitudes the simulated surface nutrient concentrations are often too high or too low when compared to observations.Simulated oxygen concentrations at ∼ 300 m with the new model are slightly better in the eastern equatorial and north Pacific regions than with the old model (Fig. 10), but similar discrepancies exist in other regions when compared to observations.

Seasonal results
For surface nitrate and phosphate concentrations, the new model produces a seasonal pattern that is similar, especially in terms of magnitude, to WOA09 data (Figs.11a, c and 12a, c).However, there is often a slight temporal mismatch between the occurrences of seasonal changes in these nutrients.Nonetheless, the new model does a much better job of simulating these seasonal changes when compared to the old model (Figs.11a, b and 12a, b; Fig. 13).Since surface nutrient concentrations are strongly influenced by biological nutrient uptake and regeneration these results, and the differences between the model formulations, are mostly due to differences in simulated biological seasonal cycles (i.e., the spring bloom) at higher latitudes (see below).

Annual results
Rates of annual global net primary production (NPP) and export production with the new model are almost the same as with the old model, while the rate of new production (the fraction of integrated NPP that is not supported by nutrients from the fast remineralisation process, zooplankton excretion and the remineralisation of detritus) is lower and nitrogen fixation and denitrification rates are a few percent higher (Table 3).These simulated annual global NPP rates compare well to present day estimates of annual global NPP (44-67 Gt C yr −1 ) derived from satellite measurements (Behrenfeld et al., 2005;Carr et al., 2006;Westberry et al., 2008).However, the spatial patterns of NPP in the models are noticeably different (Fig. 14), although major features such as low NPP in the oligotrophic gyres are resolved in both.With the new model (Fig. 14a) well-defined areas of intense NPP are present in the equatorial eastern Pacific, Atlantic and Indian oceans and the western south Atlantic ocean, while with the old model NPP is highest in broad swaths of the equatorial eastern Pacific and Indian oceans (Fig. 14b).Since PP at higher latitudes varies significantly on a seasonal basis, we have not made further comparisons to other estimates of PP here and instead refer you to Sects.5.2.2 and 6.2 for further  analysis of PP on a seasonal basis.However, we should note that aside from seasonal light limitation at high latitudes and the constraints on growth imposed by temperature and our growth-limiting iron mask (Fig. 3a), PP in both old and new models is almost entirely nitrogen limited (data not shown).
The model rates of N 2 fixation (Table 3) are within the wide range (∼ 100-200 Tg N yr −1 ) of global N 2 fixation estimates (Codispoti, 2007;Deutsch et al., 2007;Gruber and Sarmiento, 1997;Karl et al., 2002).However, as with NPP the spatial patterns of N 2 fixation in the models are quite different (Fig. 15).With the new model (Fig. 15a) N 2 fixation is high in the tropical/subtropical north Pacific, the western tropical/subtropical south Pacific, the western tropical/subtropical south Atlantic Ocean and the Indian Ocean, while with the old model (Fig. 15b) N 2 fixation was high in the central and eastern tropical/subtropical Pacific.Although the distribution of diazotrophs is poorly known in many areas of the ocean due to a lack of data, the patterns of N 2 fixation from the new model are much more consistent with observations (Sohm et al., 2011) than are the patterns from the old model.However, as discussed by Somes et al. (2010a), who produced a similar distribution of N 2 fixation with the old model when they constrained diazotroph growth with an Fe mask, N 2 fixation does not extend as far north as it should in the north Pacific.Moreover, there is too much N 2 fixation in the south Atlantic and not enough in the north Atlantic, where some of the highest rates of N 2 fixation have consistently been measured (Capone et al., 2005;Sohm et al., 2011).In the Indian Ocean, the model also simulates too much N 2 fixation in the Bay of Bengal and not enough in the Arabian Sea.Since iron cycling and availability, which is poorly understood, and phosphorous availability appears to control most N 2 fixation (Monteiro et al., 2011), other coarse-resolution global biogeochemical models have had similar problems in correctly simulating the spatial distribution and intensity of N 2 fixation because they fail to adequately simulate these nutrient cycles (Monteiro et al., 2010;Moore and Doney, 2007).
Recent global estimates of denitrification (defined here as the ensemble of biological processes that can convert fixed-N to N 2 ) vary widely, are poorly constrained and are the subject of considerable debate (Codispoti, 2007;Gruber, 2004),   which makes it difficult to validate this aspect of the model results.Furthermore, the model does not include benthic denitrification and, therefore, in the steady state essentially reached at the end of the spin up, N 2 fixation must be approximately equal to denitrification even in the absence of benthic denitrification.Thus, our model estimate of water-column denitrification is less than optimal.Nonetheless, the model rates of water column denitrification (Table 3) are the same as the 150 Tg N yr −1 suggest by Codispoti et al. (2001) and Codispoti (2007), but are higher than the 50 (±20) and 80 (±20) Tg N yr −1 suggested by Gruber and Sarmiento (2002) and Gruber (2004).As with other model results, the location of water column denitrification with the new model is somewhat different than with the old model (Fig. 16).offshore in a broad area with the most intense denitrification occurring around the Galapagos Islands (Fig. 16b).In both model versions denitrification also occurs off the coast of Namibia, Africa and in the Bay of Bengal.However, much more denitrification now occurs in the Bay of Bengal with the new model when compared to the old one.Additionally, with the new model, denitrification occurs in the Gulf of Aden as well.Validating whether the model simulates water column denitrification in the correct locations is difficult since observations of denitrification are sparse and the spatial extent of denitrification zones is often estimated by using many assumptions, and an often limited set of phosphate and nitrate observations, to calculate if a nitrate deficit occurs (Codispoti, 2007;Kamykowski and Zentara, 1990;Paulmier and Ruiz-Pino, 2009).However, it is generally agreed that the three major quasi  model formulation (see Sect. 6.2 for more discussion on productivity).The new model does especially well in simulating the North Atlantic spring phytoplankton bloom (Fig. 19), which was almost nonexistent with the old model, and which is known to be the most pronounced spring bloom of any open ocean region (Yoder et al., 1993).The seasonal cycle of zooplankton, in terms of biomass, is also more realistic with the new model formulation simulating a clear "classic" lag between the peak in phytoplankton biomass and that of zooplankton during the spring/summer at mid to high latitudes (Fig. 17).Zonal diazotroph biomass also varies seasonally with the new model formulation, which appears to be in contrast to the fairly constant presence of diazotrophs with the old model (Fig. 17c and d).However, Fig. 17d is somewhat misleading due to the zonal averaging.Diazotroph biomass does vary seasonally with the old model, particularly in the North Eastern Equatorial Pacific and Northern Indian Ocean, but with a spatial and temporal pattern that averages these variations out when presented in a Hovemöller figure.Although data on the seasonal abundance of diazotrophs is limited, the available datasets suggest that diazotroph biomass does vary seasonally in many regions (Sohm et al., 2011) as in the model runs.

Annual and seasonal particulate fluxes
Since the export flux of particulate matter out of the euphotic zone is a primary driver of marine biogeochemical cycles and the biological pump, we now discuss how well it is simulated by the models.For both the old and new model configuration, annual estimates of export production and CaCO 3 export out of the euphotic zone (Table 3) are within the observational estimates (5.8-13 and 0.38-1.64Gt C yr −1 , respectively) of a number of studies (listed in Table 4 in Dunne et al., 2007).The spatial patterns of our annual simulated fluxes (data not shown) are, for both models, similar to those of primary production (Fig. 14) with the highest export occurring below or near areas of high productivity.To evaluate  4) to the sediment trap data of Honjo et al. (2008).Although neither model formulation performs particularly well, the new model is slightly better at simulating the flux of POC to the deep ocean (Fig. 20 RMSE value).Globally, the simulated flux of POC at about 2 km depth with the new model is lower than the observationally based estimates, while the flux of CaCO 3 (i.e., PIC) is higher (Table 4).Additionally, our simulated rain ratio at 2 km is lower than observations suggest.The simulated transfer efficiency (the ratio of the POC flux at 2 km to export production) is also lower than estimated by Honjo et al. (2008), while the change in POC (the difference between export production and the POC flux at 2 km) is higher.However, Honjo et al. (2008)   duction using the model of Laws et al. (2000) and obtained a value that is lower than most other estimates (Dunne et al., 2007).If a higher export production value of 9.6 Gt C yr −1 (from Dunne et al., 2007). is used to recalculate the transfer efficiency and the change in POC with the same deep ocean POC flux data of Honjo et al. (2008), then our transfer efficiency is higher and our change in POC is lower than the deep ocean flux data suggests (Table 4).Uncertainties in the observations are relatively large because there exists only a limited amount of flux data, with possible biases in the data, and uncertainty in extrapolating sparse sediment trap data to the global ocean, as well as uncertainty in estimates of export production derived from various models.The fact that our simulated annual global flux values are close to observational estimates, together with the relatively good agreement with mean profiles of biogeochemical tracers, suggests that the model does a reasonable job of simulating the global biological carbon pump with the associated impacts on nutrient and oxygen distributions, despite some regional biases and discrepancies.
Although both the new and old models have very similar annual global POC fluxes at 2 km (Fig. 20 and Table 3) their simulated seasonal fluxes of POC differ considerably (Fig. 22) due to the differences in seasonal ecosystem cycles.However, despite these differences in both models the peaks in POC flux occur approximately one to two months after the seasonal peaks in primary production (Fig. 18a and b), which is consistent with observations that maximum flux timings are generally synchronous with the timing of blooms in temperate and polar waters and have a production-to-flux lag that is typically between 40 and 80 days (Lutz et al., 2007).The formulation of the downward CaCO 3 flux assumes instantaneous export and remineralisation in both model configurations.The seasonal fluxes of CaCO 3 are, thus, at any depth, tightly and with zero phase shift linked to CaCO 3 production in the surface ocean, which is in turn proportional to the production of detritus.As a consequence, deep CaCO 3 fluxes with the new model (Fig. 23a) occur approximately a month earlier than the highest POC fluxes, which results in a seasonally variable rain ratio at depth (Fig. 23b).While this is consistent with the fact that calcite and aragonite are some of the heaviest biominerals found in ocean water and, thus, sink more rapidly than POC, it is doubtful that in re-   ality CaCO 3 fluxes and lag times are as uniform as in the simulations since the size of the source material (i.e., coccoliths versus foraminiferals tests or pteropod shells) and its production, which is determined by the composition of the community and, thus, not accounted for in the model, will play a large role in determining its sinking speed and export to depth (Honjo et al., 2008).Note also that the significance of the very high rain ratios in polar regions during the winter must be interpreted with caution because PP is very low or   nonexistent during these times of year and, therefore, there is not a large flux of material to deeper waters.

Primary production
A closer examination of PP with the new model shows that while seasonal PP at higher latitudes compares well with some estimates (i.e., Fig. 19), simulated PP is too high in the tropical upwelling regions along the equator and the northern Indian ocean with both the new and old models when compared to satellite-based estimates (Fig. 18) and observations (data not shown) of PP such as those made by Balch et al. (2011) in the equatorial Pacific and Marañón et al. (2000) in the equatorial Atlantic.Furthermore, with either model there is too little productivity in oligotrophic regions when compared to the satellite-based models (Fig. 18) and observations (i.e., such as those made by Marañón et al. (2000) and  et al. (2008) in the south Pacific).However, it is worth noting that the satellite-based estimates of PP do not agree very well with each other (Figs.18c and d, 19 right panels) and, thus, caution must be taken when making comparisons to these models.Furthermore, satellite-based estimates of PP are inherently based on models, so that this type of assessment essentially makes model-to-model comparisons, which cannot be considered true validation.Nonetheless, there are a number of systematic differences between our model and reported estimates of PP.First, the model parameters were set so that the total annual global PP was within the range of estimates (i.e., 44-67 Gt C yr −1 ; from  et al., 2004;Vichi and Masina, 2009;Yool et al., 2011) also simulate higher than observed PP in many regions (especially the equatorial upwelling regions) to achieve estimated annual PP within the range of 40 to 67 Pg C yr −1 .Furthermore, this problem is not limited to estimates of global PP as Dunne et al. (2007) noted that a similar problem exists in calculations of particle export and burial which use satellite-based PP data that is regridded onto coarser grids excluding coastal areas.The second reason for model and observational differences is the use of invariant half-saturation constants for nutrient uptake, which prevent phytoplankton from adapting to low nutrient concentrations in oligotrophic regions and, thus, severely limits their growth in these areas (e.g., Smith et al., 2009).Finally, the dissolved iron mask that we use to constrain phytoplankton growth may be inaccurate as it is generated by another model (i.e., from the BLING model, Galbraith et al., 2010) and not well validated due to a lack of observations in many ocean regions.Moreover, we do not explicitly model the cycling of iron or it's uptake by phytoplankton, which may be necessary to correctly simulate the effects of iron limitation on phytoplankton growth.

Grazing
Grazing related processes are quite different in the two model versions (Fig. 24).The annual percentage of primary production that is grazed daily is much higher in most regions with the new model when compared to the old one (Fig. 24a  and b).The notable exception to this trend is in the Southern Ocean south of 60 • S where it is not surprising that grazing is lower in the new model since the maximum potential grazing rate is now a function of temperature and much lower in this region (Fig. 3b) than in with the old model where it was set globally at 1.5 d −1 .The differences between the model versions are also particularly striking in the oligotrophic gyres where almost no grazing occurred with the old model because there were essentially no zooplankton present in these regions (data not shown).Validating how much primary production is lost to grazing is difficult because there are few observations to constrain this trophic transfer on a global scale and it can vary significantly regionally and seasonally.Nonetheless, global estimates suggest that mesozooplankton consume ∼ 10-15 % (Behrenfeld and Falkowski, 1997;Calbet, 2001) and microzooplankton consume ∼ 59-75 % (Calbet and Landry, 2004) of oceanic PP.Grazing with the new model, where an average of 71 % of PP is consumed by zooplankton, falls within the range of these estimates, while with the old model an average of only 52 % of PP was consumed which is below the lowest estimate for microzooplankton grazing alone.Furthermore, with the new model the ratio of zooplankton to phytoplankton biomass (integrated over the upper 500 m), which is strongly influenced by grazing, is much higher (0.94) than with the old model (0.67) and more consistent with observations of high heterotrophic to autotrophic biomass ratios in the open ocean (Gasol et al., 1997).

Nutrient regeneration
The annually averaged remineralisation rates of inorganic nutrients by zooplankton (i.e., excretion) are quite different in the two model versions (Fig. 24c and d).Most of the remineralisation with the old model occurs in the tropics, which is not surprising since with this model formulation zooplankton excretion is dependent only on biomass and temperature, whereas ingestion is temperature independent.Remineralization with the new model, where it is now a function of grazing, occurs primarily in areas where grazing is high, with the spatial pattern of remineralisation in Fig. 24c looking much like the pattern of primary production in Fig. 14a.The differences between the model versions are also due to the effects of the new grazing formulation on the pathways of nutrient remineralisation.With the new grazing formulation more nutrients are remineralized through excretion than with the old model formulation where the fast recycling loop (i.e., implicit remineralisation of organic matter by bacteria) was a more important pathway (Fig. 25).The relative importance of one pathway versus another in the euphotic zone is difficult to quantify on a global scale since it depends on many factors that vary both temporally and spatially.Furthermore, neither model explicitly includes bacteria or dissolved organic matter, which may be necessary to effectively simulate bacterial remineralisation.Nonetheless, there is considerable evidence that zooplankton play a more important role in nutrient regeneration than the old model formulation allows for (Banse, 1995;Ferrier-Pages and Rassoulzadegan, 1994;Glibert, 1998).

Summary and conclusions
Overall, the new ecosystem model formulation has achieved our goals of improving the realism of simulated biological processes and biogeochemical cycles.When evaluated against observational data the majority of the results show an improvement over the old model.Many of these improvements can be attributed to a better simulation of higher latitude ecosystem seasonality and the effect that it has on biogeochemical cycles.The important modifications include: -The inclusion of diazotroph biomass in calculations of light attenuation and phytoplankton light limitation.
-Improved zooplankton growth and grazing formula that add more realistic dynamics.Notably, the inclusion of a temperature dependent zooplankton growth rate,  which is responsible for most of the improved ecosystem seasonality because it eliminates the restrictively high higher latitude top-down control on phytoplankton growth that is present in the old model formula and, thus, allows phytoplankton to bloom in response to seasonal forcing.
-The use of a dissolved iron mask to constrain phytoplankton growth in iron-limited regions of the ocean.
-Changes to the parameters governing phytoplankton and zooplankton growth (i.e., maximum growth rates) to achieve more realistic rates and better responses to seasonal forcing.
Despite these modifications there are several aspects of the model that might be improved in the future.One of the most important issues that needs further consideration is the discrepancy between modeled PP and observational estimates.A number of factors would need to be addressed to improve simulated PP.First, there is the decision of what levels of annual PP should be accounted for in a model that does not resolve coastal waters and the significant amount of PP that occurs in them.Perhaps models should simulate lower levels of global primary productivity instead of trying to simulate global rates that are as high as satellite estimates.Including coastal areas or accounting for the flux of material from them may also be necessary, although this would involve a significant modification of the model grid, physics and boundary conditions.Perhaps including variable half-saturation constants for nutrient uptake or introducing more phytoplankton types may also be necessary to better simulate PP (e.g., Smith et al., 2009).However, any of these options introduces other problems, which have been discussed extensively by others (Anderson, 2005;Flynn, 2005;Hood et al., 2006), and may not necessarily improve the simulation biogeochemical cycles.Furthermore, the mask of dissolved iron concentrations that we use to constrain phytoplankton growth is derived from a model and could be improved to help address issues of where and how much primary productivity occurs.This limiting factor could be improved either by including a dynamic iron cycling model within the ecosystem model or by obtaining a better mask as more measurements of dissolved iron become available.Regardless, these issues are beyond the goals we set out to achieve with the improvements made here and can be addressed in the future.Our modifications have significantly improved the marine ecosystem model and provided a useful platform for future research with the UVic Earth System Climate Model.

Fig. 1 .
Fig. 1.Ecosystem model schematic which illustrates the flux (arrows) of material between model variables (squares).See text for a detailed description of these fluxes.

Fig. 5 .Fig. 6 .
Fig. 5. Zonally averaged ocean basin comparisons of PO 4 between the new and old models and the World Ocean Atlas 2009 data set.

Figure 10 Fig. 10 .
Figure 10 -permanent sites of water column denitrification are within the suboxic portions of the thermocline in the Arabian Sea, the Eastern tropical North Pacific and the Eastern tropical South Pacific off the coast of Chile (Codispoti, 2007).Of these major denitrification sites, the model is only able to reasonably simulate the eastern tropical North Pacific one.In the Eastern tropical South Pacific, modelled denitrification is too close to the equator and does not occur off the coast of Chile where it should.In the Indian Ocean, Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec 0

Fig. 13 .
Fig. 13.Root mean squared error (RMSE) comparisons between the new and old models and WOA09 monthly (a) nitrate and (b) phosphate data.Zonal statistics were calculated only for grid cells where observations were present and then averaged longitudinally to give a global value.

Fig. 15 .
Fig. 15.Annual vertically integrated rates of N 2 fixation with the (a) new and (b) old model formulations.

Fig. 17 .
Fig. 17.Hovmöller diagrams of the zonally averaged seasonal integrated phytoplankton (a, b), diazotroph (c, d) and zooplankton (e, f) biomasses as simulated by the new (left panels) and old (right panels) models.

Fig. 24 .
Fig. 24.Annual vertically integrated percentage of primary production grazed per day (top panels) and zooplankton excretion (bottom panels) with the (a and c) new and (b and d) old model formulations. Figure25

Fig. 25 .
Fig. 25.Comparison between model versions of the pathways of nitrate regeneration.Values are from globally integrated data.
are listed in Table 1.Table 2 defines additional pa-

Table 2 .
Definitions of parameters and variables that are not specifically mentioned in the text.
i Calculated sea ice thickness Bitz et al. (2001) a i Fractional sea ice cover Bitz et al. (2001) h s Calculated snow cover thickness Bitz et al. (2001)

Table 3 .
Comparison of globally important properties between model versions.CaCO 3 was set at 4500 m which differs from the value of 3500 m used inSchmittner et al. (2008).

Table 4 .
Comparison of critical properties that drive the global biological pump.

Table 3
Westberry et al., 2008)8).Setting the model parameters to achieve this level of annual global PP may be problematic because these satellite-based estimates include coastal PP, which can be quite high and accounts for a significant amount of the estimated total.Since the UVic ESCM model has a course grid resolution (1.8 × 3.6 degrees) the coastal regions are poorly represented, if at all in many places, in the marine ecosystem www.geosci-model-dev.net/5/1195/2012/Geosci.Model Dev., 5, 1195-1220,

2012 1214 D. P. Keller et al.: A new marine ecosystem model for the University of Victoria's ESCM model
. Thus, simulated open ocean productivity becomes much too high when the model is parameterised to simulate levels of annual global PP that include coastal PP.This problem does not appear to be limited to the UVic ESCM as a number of other global marine ecosystem models (Moore