Explicit planktic calcifiers in the University of Victoria Earth System Climate Model

Introduction Conclusions References Tables Figures


Introduction
Earth system models are incorporating ever larger ecological schemata to represent the growing mechanistic understanding of biological connections to global biogeochemical cycles.In the ocean, "Dynamic Green Ocean Models" (Le Quéré et al., 2005) use multiple plankton functional types (PFTs) to explicitly link marine organisms to global chemical cycles through ecology and physiology.These PFTs are not explicit organisms but are instead conceptual classifications of marine plankton according to their biogeochemical role (Hood et al., 2006).
Pelagic calcifiers not only contribute to deep sea and benthic carbon inventory but also affect the atmosphereocean carbon dioxide (CO 2 ) gradient through a small release of CO 2 during calcification (Eq.(1); Zondervan, Zeebe, Rost, & Riebesell, 2001).This release of CO 2 provides a chemical link between calcification and photosynthesis (Eq.( 2)), during which some of the CO 2 is used for POC production (with a net fixation of carbon).
Including explicit calcifying phytoplankton in a fully interactive ocean-atmosphere-biogeochemical model is warranted given their importance in carbon cycling.Furthermore, inclusion of a ballasting parameterization is desirable given its demonstrated significance for ocean oxygenation (Hofmann & Schellnhuber, 2009).The following describes their application to a climate model of intermediate complexity and assesses the model's performance using available biogeochemical and biomass data.
In this latest version, the general phytoplankton PFT is exactly replicated but given new parameter values to reflect key physiological characteristics of phytoplankton calcifiers, albeit biased towards Emiliania huxleyi.This new model version therefore contains "phytoplankton calcifiers," "diazotrophs," and "general phytoplankton."The general phytoplankton PFT includes diatoms as well as all other autotrophic non-calcifying phytoplankton.Just as the general phytoplankton PFT cannot perfectly describe the physiology or ecology of any of the individual classifications of phytoplankton it represents, the calcifying PFT represents a group of phytoplankton with a common role in the carbonate cycle (calcification) and a few generalized shared physiological traits.In this new model, only calcifying phytoplankton and the zooplankton PFT produce CaCO 3 .The CaCO 3 is calculated prognostically as a model tracer, and dissolution of phytoplankton calcifier and zooplankton export is now dependent on ambient carbonate concentration.The new model schematic is shown in Fig. 1.
In the following model description, notation will generally follow the symbols used in Keller et al. (2012), with additionally P standing for the general phytoplankton PFT, C standing for the phytoplankton calcifier PFT, and Z representing zooplankton when a distinction is necessary.Relevant model parameters are listed in Tables 1-4 with T including all transport terms (advection, diffusion, and convection) and S representing all source and sink terms.
1 PHYTOPLANKTON General phytoplankton and calcifying phytoplankton (X representing either) population source and sink terms are where the growth rate (J), mortality (m), and fast recycling (μ * ) terms are described below, and losses to zooplankton grazing (G) are described in Section 2b2.The diazotroph population sources and sinks follow (5) As in Schmittner et al. (2005), the maximum possible growth rate of general phytoplankton and phytoplankton calcifiers (J max ) is a modified Eppley curve (Eppley, 1972) and is a function of seawater temperature (T), an e-folding temperature parameter T b , and iron availability (u Fe ).Parameter values are listed in Table 3. Phytoplankton calcifiers are assigned a lower maximum growth rate (a) than mixed phytoplankton, an assumption used previously by Le Quéré et al. (2005) but also justified by comparing measured growth rates for a selection of four coccolithophores by Buitenhuis, Pangerc, Franklin, Le Quéré, and Malin (2008) (0.3-1.0 d −1 at 15°C) with the general range for phytoplankton by Eppley (1972) (a maximum rate of about 2.2 d −1 at 15°C).
Iron limitation is calculated from the concentration of iron prescribed in interpolated monthly-mean fields using an iron half-saturation approximation constant (k Fe ) (Galbraith, Gnanadesikan, Dunne, & Hiscock, 2010;Keller et al., 2012).A prognostic iron cycle was recently implemented in the UVic ESCM (Nickelsen, Keller, & Oschlies, 2014), though it adds computational expense.However, accounting for iron limitation on growth rates by means of a limitation mask improves phytoplankton biogeography without additional computational cost (Keller et al., 2012).Calcifying and mixed phytoplankton are assigned different k Fe values that vary the degree of iron limitation and are tuned to produce the best possible PFT distributions, not actual iron affinities.Phytoplankton calcifiers are assigned a lower k Fe value than mixed phytoplankton to simulate the relatively low iron half-saturation constant for phytoplankton calcifiers recommended by Le Quéré et al. (2005): The maximum potential growth rate is then multiplied by a nutrient availability (u) for both nitrate (NO − 3 ) and phosphate (PO 3− 4 ) to calculate growth under nutrient limitation, where k N and k P are half-saturation constants: These equations are applied to obtain maximum possible growth rates as a function of temperature and nutrients.As in Schmittner et al. (2005), the maximum possible growth rate under limited light availability (J I ) is calculated as: where α is the initial slope of the photosynthesis versus irradiance (I) curve.Phytoplankton calcifiers have a lower α than diatoms, though it is similar to non-bloom-forming mixed phytoplankton (summarized in Le Quéré et al., 2005).Therefore, a lower α value for phytoplankton calcifiers is used here.Additionally, light scattering by coccoliths is considered in calculating available irradiance at each depth where PAR stands for the photosynthetically available radiation; k w , k c , k CaCO 3 , and k i are the light attenuation coefficients for water, all phytoplankton (calcifiers, diazotrophs, and general phytoplankton), CaCO 3 , and ice, respectively; z is the effective vertical coordinate; a i is the fractional sea-ice cover; and h i and h s are calculated sea-ice and snow cover thickness.Values for k c and k CaCO 3 come from Balch and Utgoff (2009).
The actual growth rate (J PorC ) of the general phytoplankton and calcifying phytoplankton PFTs is taken to be the minimum of the three growth functions described above: Diazotroph growth is not dependent on NO − 3 concentration and hence follows Two loss terms other than predation (described below) are considered.Mortality from old age or disease is parameterized using a linear mortality rate (m).Temperature-dependent fast remineralization is a loss term used to account for the microbial loop and dissolved organic matter cycling and is parameterized using a temperature dependency multiplied by a constant (μ * 0 ): 2 ZOOPLANKTON Zooplankton population (Z) is calculated as the total available food (POC) scaled with a growth efficiency coefficient (ϖ) minus mortality.In addition to old age and disease, zooplankton mortality also encompasses losses from higher trophic level predation and starvation.
Zooplankton grazing (G) follows Keller et al. (2012).Relevant parameters are listed in Table 4. Grazing of each food source (mixed phytoplankton, calcifying phytoplankton, diazotrophs, zooplankton, and detritus) is calculated using a Holling Type-II function, in which a calculated maximum zooplankton grazing rate (μ max Z ) is reduced by a scaling that is weighted by a food preference (ψ X , where X stands for any of general phytoplankton, calcifying phytoplankton, diazotrophs, zooplankton, or total detritus), the total prey population, and a half-saturation constant for zooplankton ingestion (k z ): Other formulations of grazing exist; see Anderson, Gentleman, and Sinha (2010) for a detailed comparison using two zooplankton size classes.Because the focus of this study is the implementation of a calcifying phytoplankton functional type, the grazing parameterization was not modified.The calculated maximum potential grazing rate is a function of a maximum potential grazing rate at 0°C (μ θ Z ), temperature, and oxygen (sox stands for suboxic), where grazing activity is capped when temperatures exceed 20°C: Grazing is also reduced under hypoxic conditions (r O 2 sox ): where O 2 is dissolved oxygen in micromol.

DETRITUS
Detritus sources and sinks now include contributions from phytoplankton calcifiers and are split into free and ballast pools using a fixed ratio (R bal:tot ).Ballasted detritus is formed of the CaCO 3 -protected portion of phytoplankton calcifier and zooplankton grazing and mortality.For simplicity, the same R bal:tot is used for both phytoplankton calcifiers and zooplankton.This protected portion does not interact with nutrient pools directly but instead transfers from the ballast to free detrital pool at the rate of CaCO 3 dissolution (λ; Eq. ( 33)): where γ is the food assimilation efficiency, R CaCO 3 :POC is a fixed production ratio of CaCO 3 and detritus, R C:N is a Redfield molar ratio, and μ D is the detrital remineralization rate.As in Keller et al. (2012), detritus is exported from the surface with a sinking speed (w DorC ) that increases linearly (in per second units) with depth: The initial surface sinking speeds of POC and CaCO 3 (w DorC0 ) are assigned different values to represent the denser structure of CaCO 3 relative to that of POC.Ballasted detritus sinks at the CaCO 3 speed, but once it enters the free pool it uses the detrital sinking speed and remineralization rate.Any detritus reaching the sediments is dissolved back into the water column.

DISSOLVED BIOGEOCHEMICAL TRACERS
Ocean nutrient sources and sinks follow: where R P:N and R O:N are Redfield molar ratios, and u N is the Michaelis-Menten nitrate uptake rate.In suboxic (sox) water, oxygen consumption is replaced by the oxidation of nitrate and ocean surface dissolved oxygen exchanges with the atmosphere (F sfc ).
The DIC and alkalinity tracer sources and sinks are now also a function of sources and sinks of prognostic CaCO 3 (Section 2b5):

CALCITE PRODUCTION AND EXPORT
The original model fixed CaCO 3 production to POC using a uniform ratio of CaCO 3 production to non-diazotrophic POC (detritus) production (R CaCO 3 :POC ).The CaCO 3 produced then contributed to DIC and alkalinity with a fixed remineralization profile exponentially dependent on depth.In our model, the general phytoplankton PFT no longer contributes to CaCO 3 but is instead replaced with the phytoplankton calcifier PFT.Different R CaCO 3 :POC values for zooplankton and phytoplankton calcifiers can be assigned in case the ballast model is turned off, but a second ballasted detritus tracer would be required for this feature to be used with the ballast model.This second detritus tracer is not yet implemented, so the tuned model presented here includes ballast and a shared R CaCO 3 :POC value for zooplankton and phytoplankton calcifiers.In earlier versions of the UVic ESCM, an R CaCO 3 :POC value of 0.03 was used.In this version, it is increased to 0.04, which places it closer to (but still outside) the low end of the 0.05-0.25 range estimated by others and summarized by Fujii, Ikeda, and Yamanaka (2005).A CaCO 3 :POC production ratio for E. huxleyi is summarized by Paasche (2001) to vary between 0.51 and 2.30, depending upon nutrient status and strain.A lower rain ratio for the model, therefore, indicates that the phytoplankton calcifier PFT cannot be considered to represent calcifiers exclusively, with other non-calcifying phytoplankton sharing the physiological traits also represented by the PFT.Likewise, calcification by the zooplankton PFT must be considered a global zooplankton average, with the zooplankton PFT representing all calcifying and non-calcifying zooplankton.
Production and dissolution of CaCO 3 are now a source and sink of a prognostic CaCO 3 tracer (Eq.( 31)).Calcite held in living tissue is calculated separately as the net source-sink from phytoplankton calcifiers and zooplankton (Eqs (4) and ( 15)), converted to CaCO 3 units: where R C:N is the Redfield ratio (Table 1).New model tracer particulate CaCO 3 (in non-living form) follows the same general model structure as detritus, though the base units are millimoles of carbon per cubic metre rather than millimoles of nitrogen per cubic metre.One critical difference exists, which is that any CaCO 3 consumed by zooplankton is assumed to pass through the digestive system unaltered and is not assimilated into zooplankton growth.The source and sink terms for CaCO 3 include both phytoplankton calcifier and zooplankton sources from grazing and mortality, as well as losses from dissolution (31) The full equation of particulate CaCO 3 is thus: with T including all transport terms (advection, diffusion, and convection).
A CaCO 3 dissolution rate (λ) that allows for supersaturated dissolution (Milliman et al., 1999) is calculated using a fixed dissolution rate parameter (k), following the calculation used in the Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES) model family (Aumont et al., 2003): where δ sat is the deviance of the ambient seawater carbonate concentration from saturation (δ sat = [CO 3 ] − [CO 3 ] sat ) and any negative λ is set to zero.Particulate CaCO 3 that reaches the sediments accumulates in an oxygen-only respiration model following Archer (1996b) that is unchanged from earlier versions of the UVic ESCM.During model spin-up, losses of alkalinity and carbon to the sediment model are exactly compensated for by a terrestrial weathering flux (diagnosed from the net sediment burial rate) that is applied as a flux of alkalinity and DIC to the ocean through river discharge.Once the model is in equilibrium either a constant or a prognostic terrestrial weathering flux can be used (Meissner et al., 2012).
The changes described in the above section improve the mechanistic realism of the UVic ESCM by explicitly including a phytoplankton functional type (phytoplankton calcifiers) that is both uniquely vulnerable to resource competition and uniquely responsible for CaCO 3 production and export.Representation of the CaCO 3 cycle in the UVic ESCM is additionally improved by including thermodynamic dissolution and detrital ballasting.
3 Model tuning a Stabilizing Alkalinity Model tracers, alkalinity and DIC are very sensitive to the prognostic CaCO 3 described in the previous section, which made model tuning a challenge.Stabilizing the model with realistic parameter values required multiple steps.After each step, conservation of global alkalinity and carbon was confirmed before proceeding.The original UVic ESCM ocean chemistry is in fairly good agreement with observations, so the initial goal was to tune the model as closely as possible to the original output.To do this, annual mean CaCO 3 dissolution at a pre-industrial equilibrium was diagnosed from the original model.The output file was then used as input to the new model (but, for structural reasons, without the ballast model turned on) to prescribe CaCO 3 dissolution.Production of CaCO 3 in the new model was not dramatically changed from the original model, but the possibility of greater dissolution than production in any given grid box meant a correction term was required to avoid negative CaCO 3 concentrations.The CaCO 3 tracer was, therefore, calculated from the ocean bottom to the surface in a reverse depth loop, and when the tracer was calculated as negative, a correction term was added to the concentration to set the tracer equal to zero.The correction term was then carried into the tracer calculation in the grid box above, which was likewise adjusted with a correction term if needed.At the surface, the integrated correction was added to the total CaCO 3 production to conserve carbon.In this way, the new model with a prognostic CaCO 3 tracer was able to reproduce the alkalinity and DIC fields of the original instant-export-production model.
b Adjusting Phytoplankton Production The next step was to tune production as closely as possible to average global estimates and PFT distributions in the modern ocean (Table 5).Production parameters have been shown to be highly model dependent (Kriest & Oschlies, 2011).Parameters were adjusted under the constraints that mixed phytoplankton parameters be kept at original model parameter values and that phytoplankton calcifier growth rate, nitrogen and iron uptake, and α values would all be lower than the mixed phytoplankton parameter values (Le Quéré et al., 2005).According to Scott, Kettle, and Merchant (2011), over 10% of model variance in primary production is attributable to four parameters: maximum growth rate (a, in the high and mid-latitudes); the initial slope of the photosynthesis-irradiance curve (α, at all latitudes); mortality (a more model-dependent variable having the largest impact at low latitudes, μ 0 * and μ 0 here); and the carbon to chlorophyll ratio (at low latitudes, but not included in this model).Growth rate was by far the most sensitive of the production parameters in this model, with mortality and α having less of an influence on biomass and biogeography.As has been shown previously (e.g., Cropp & Norbury, 2009), achieving multiple extant PFTs required careful model tuning.Nutrient half-saturation constants for nitrate and iron provided phytoplankton calcifiers a competitive advantage, whereas a lower growth rate and α produced a disadvantage.The tuning of these parameters (a, α, k N , and k Fe ) required an iterative process to sufficiently "balance" the relative advantages with the relative disadvantages such that both general phytoplankton and calcifying phytoplankton populations remained extant and roughly realistically distributed in the surface ocean.As in other multiple-PFT models (e.g., Cropp & Norbury, 2009) similar growth rates for phytoplankton calcifiers and general phytoplankton were required to maintain both populations, but more variable nutrient uptake and grazing parameter values were possible.
c Implementing Prognostic CaCO 3 Dissolution With fixed CaCO 3 dissolution, tuned production, and stable alkalinity and DIC, the next step was to tune the CaCO 3 sinking rate.A sensitivity study across a range of w C0 and R CaCO 3 :POC values was conducted to determine the combination that yielded the best fit to the original model CaCO 3 export (Table 5) and did not substantially alter ocean alkalinity distributions.After these parameter values were determined, the model was integrated for several thousand years to achieve an equilibrated state.The new model CaCO 3 dissolution scheme was run in parallel as a diagnostic only and tuned to approximate the original model dissolution.The new model CaCO 3 dissolution scheme replaced the original dissolution scheme after model equilibrium was achieved.The reverse-loop correction of CaCO 3 was not necessary after this step, so it was turned off and CaCO 3 was treated the same as any other tracer in the model.The CaCO 3 ballasting of detritus was the last component of the model to be turned on.Parameters R bal:tot , w D0 , w C0 , and R CaCO 3 :POC were then re-evaluated to determine the optimal values.The model was tuned to reproduce (as best as possible) Global Data Analysis Project (GLODAP) and World Ocean Atlas (WOA) observations (Figs 2 and 3; Garcia et al., 2009;Key et al., 2004).The Taylor diagrams (Taylor, 2001) shown in Fig. 3 compare the correlation and standard deviation of model-simulated tracers normalized against the standard deviation of GLODAP and WOA observations.This section offers an example of how major structural modifications to existing models may not be   Popova, 1997) and the inclusion of kill-the-winner feeding (Vallina, Ward, Dutkiewicz, & Follows, 2014) have been proposed to address the common multi-PFT model problem of phytoplankton population instability encountered using the Holling Type-II grazing function.These alternatives produce a top-down control on biodiversity that can improve model agreement with phytoplankton diversity and bloom succession observations (Prowe et al., 2012) and reduce the competitive exclusion (Prowe et al., 2012;Vallina et al., 2014).Although these parameterizations might similarly improve UVic ESCM model performance, they have not been included in this model for a number of reasons.The primary objective of this study is to include prognostic CaCO 3 and a phytoplankton type resembling calcifiers and to compare these changes with the previous model.Simultaneously changing the grazing formulation would complicate the comparison and would, furthermore, require considerable additional experimentation to choose and tune the new formulation.Such an ambitious modification would be better suited to a separate study and will be seriously considered for inclusion in the model in the future.It is also important to remember that although prey switching does have a stabilizing effect, this does not necessarily mean that the assumptions behind it, or any grazing formulation for that matter, are correct (Anderson et al., 2010).

Model assessment
For the purpose of evaluation, the NOCAL and CAL versions of the model (from Keller et al., 2012, and as described here) were first brought to pre-industrial equilibrium using a fixed atmospheric CO 2 concentration of 283 ppm and integrated over ten thousand years.In each case the same physical parameters were used, and the sediment model was applied to both integrations (it was not used by Keller et al., 2012).
Model CAL biogeochemical tracers averaged globally and by ocean basin reveal generally improved performance with respect to NOCAL in reproducing GLODAP and WOA observations (Figs 2 and 3;Garcia et al., 2009;Key et al., 2004).This may be partly due to the application of a parameter set in NOCAL that was tuned by Keller et al. (2012) to a model that did not include sediments, whereas the parameter set in CAL is tuned to achieve the best fit including sediments.Both CAL and NOCAL perform well globally and in the Pacific and Southern oceans, with larger differences in the Indian and Atlantic basins.The NO − 3 fields show a clear improvement between CAL and NOCAL model versions whereas the PO 3− 4 fields are more mixed.Note that the objective was to improve the mechanistic realism of the model without sacrificing model performance with respect to the biogeochemistry, and overall this is achieved.Globally integrated biogeochemical properties for CAL and NOCAL (Table 5) reveal that although both model versions calculate global net primary production (NPP) within observational range, much of the production occurs in the eastern Pacific and Indian oceans (Figs 4h and 4i).High production in these regions is primarily from the general phytoplankton PFT in both CAL and NOCAL (Figs 4a and 4e), though phytoplankton calcifiers offer an important contribution in the CAL model (Fig. 4b).High production in the Indian basin can explain generally low surface nutrient, DIC, and alkalinity concentrations in this region (Figs 5 to 7).In the Atlantic, the model performs well with respect to observations of PO 3− 4 and DIC.As with earlier model versions (e.g., Eby et al., 2009)   Satellite data must be used with caution because they have a seasonal bias, do not distinguish between living and dead CaCO 3 (Tyrrell & Merico, 2004), and can overestimate CaCO 3 by two to three times (Balch et al., 2011).Furthermore, Brown and Yoder (1994) estimate subpolar blooms captured by satellite might only represent 0.3% of the total global calcification, with the majority of coccoliths appearing in sediments having a source that is not detectable by satellites.In situ CaCO 3 and POC concentration data are more reliable but sparser.Model living and detached CaCO 3 and POC (detritus and PFT biomass) are used to compare simulated organic and inorganic carbon with in situ samples (Fig. 9).Regression of CAL concentrations with the data compilation of Lam, Doney, and Bishop (2011) show good agreement in simulated POC concentrations in the uppermost 1000 m.Simulated concentrations of living and detached CaCO 3 are underestimated with respect to Lam et al. (2011) for values indicative of blooms (greater than 0.5 mmol C m −3 ).Simulated CaCO 3 concentrations less than 0.2 mmol C m −3 are overestimated, which is consistent with higher simulated biomass in the low latitudes.
A comparison of detached CaCO 3 concentration in CAL with the pre-industrial control concentration of the Coupled Model Intercomparison Project, Phase 5 (CMIP5) multi-model ensemble, normalized to the Lam et al. ( 2011) dataset shows CAL is within the range of CMIP5 models for which these data are available (Fig. 10).
Annually averaged global CaCO 3 export fluxes (Table 5 and Fig. 11) are low compared with sediment trap data from Honjo, Manganini, Krishfield, and Francois (2008), though both CaCO 3 and POC fluxes in CAL agree better with observations than those in the NOCAL version (CaCO 3 root mean square error (RMSE) of 147.14 in CAL compared with 188.02 in NOCAL, POC RMSE of 97.98 in CAL compared with 100.65 in NOCAL).Improved fluxes are likely a result of the addition of the variable dissolution scheme, which calculates a global average dissolution rate of 0.40 Pg C per year (also low but within the range of error when compared with independent estimates, Table 5).Although application of a ballasting scheme was found to improve POC fluxes, the tuned ballasting parameter R bal:tot yields only a small ballasted POC pool that contributes only 2.6% of the POC reaching the sediments, compared with 80-83% estimated by Klaas and Archer (2002).Spatial biases in export fluxes follow those found in CaCO 3 concentration, with too much export in the low latitudes and too little poleward of 60°compared with Honjo et al. (2008).Sarmiento et al. (2002) and Dunne, Hales, and Toggweiler (2012) both concluded, from simple box models, that the major contribution of CaCO 3 to global export must come from low-latitude, non-bloom-forming phytoplankton calcifiers or zooplankton, so perhaps the CAL model is performing better than direct comparison with trap and satellite data suggest.Six percent of global total carbon export flux at 50 m depth is CaCO 3 , compared with the Jin et al. (2006) estimate of 4% of the total carbon flux leaving the euphotic zone  (75 m depth).Model CAL rain ratio follows the pattern calculated in Sarmiento et al. (2002) of a small POC:CaCO 3 export ratio in the low latitudes that increases poleward.Simulated sediment composition in CAL varies only slightly from NOCAL, with overall lower contributions from CaCO 3 (Fig. 12).Although lower CaCO 3 concentrations represent an improvement compared with observational estimates (Archer, 1996a), concentrations are still too high because of the overproduction of phytoplankton calcifiers relative to total production.As in Dunne et al. (2012), the highest CaCO 3 fluxes to the sediments in the UVic ESCM correspond to regions with the highest CaCO 3 surface export production (not shown).
Unlike earlier versions of the UVic ESCM that used instant export and dissolution, CaCO 3 export now peaks about two months after phytoplankton calcifier biomass reaches seasonal maxima (Figs 13a and 13d).The CaCO 3 export is also now lower than in the NOCAL version (Table 5 and Fig. 13e).Model phytoplankton calcifiers bloom too early (March-May rather than June-July; O' Brien et al., 2013) in the northern latitudes.Zooplankton populations in the northern hemisphere high latitudes peaks about three months after phytoplankton calcifier biomass, with the seasonal progression being phytoplankton calcifiers first, then general phytoplankton, then zooplankton (Figs 13a to 13c).The model biomass succession is in contrast to the observed diatom to non-diatom progression (e.g., Joint, Pomroy, Savidge, & Boyd, 1993;Riebesell et al., 2007) though without explicit diatoms in the model it is expected that the model ecology could not replicate the behaviour of this keystone PFT.A previously noted correlation between Bering Sea Shelf E. huxleyi blooms and seasonal peaks in carbonate ion concentration (Merico, Tyrrell, & Cokacar, 2006) is also not seen in the CAL model because the proposed mechanism (precursor drawdown of DIC by a diatom bloom) is missing.Implementing an explicit dependence for phytoplankton calcifier growth on high CO 2− 3 would likely shift the phytoplankton calcifier biomass peak several months later in the season and move the general phytoplankton biomass peak forward, possibly improving model performance.Such a dependence might also improve CaCO 3 distributions by reducing the production and export in the low latitude upwelling zones.While increasing calcification correlates with increasing CO 2− 3 concentrations, no significant correlation between coccolith mass and chlorophyll or cell  Planktic Calcifiers in the University of Victoria Earth System Model / 343 abundance is apparent in global sampling of surface water and sediment core samples (Beaufort et al., 2011).In the southern hemisphere, zooplankton seasonality is the primary driver of CaCO 3 fluxes because of the absence of a model phytoplankton calcifier population south of 40°S.
Model phytoplankton calcifiers in CAL are reported as a molar concentration, whereas actual coccolithophores have cell biovolumes (in typical units of cubic micrometres) that are taxonomically variable (summarized in O'Brien et al., 2013).Hence, predicted phytoplankton calcifier concentrations in the CAL model are more indicative of the presence or absence of the PFT and cannot be expected to reasonably quantify abundance.Phytopankton calcifiers in CAL can be compared with the recent Marine Ecosystem Data (MAREDAT; Buitenhuis et al., 2013a;O'Brien et al., 2013) sample data synthesis.Because the CAL model does not resolve coastal processes, globally integrated total phytoplankton PFT concentrations (Table 5) are lower than the MAREDAT estimate.Phytoplankton calcifiers, however, are overrepresented by a factor of ten.This overestimate is primarily due to the low number of PFTs in the model, which requires that the phytoplankton calcifier PFT use parameter values (i.e., the growth rate factor) more similar to the general phytoplankton PFT than data support if it is to avoid extinction.The sparseness of the MAREDAT dataset limits conclusions to noting that the CAL model phytoplankton calcifiers have a far greater distribution than supported by in situ sampling and have the highest concentrations in low latitudes, in contrast to MAREDAT.The discrepancy mostly results from the overestimate of total production in this region coupled with the necessary overrepresentation of phytoplankton calcifiers to maintain an extant population.It may also be partly due to the likely sampling bias towards high-latitude blooms in the MAREDAT synthesis, with lower latitude open ocean regions having relatively few sample points (O'Brien et al., 2013).The CAL phytoplankton calcifier biomass maxima in the mid-latitudes (40°-60°N, 40°S) are generally consistent with observed high concentrations at 60°N  (2006).Also consistent with MAREDAT is the lack of much seasonality in the Southern Ocean phytoplankton calcifier population.
Regional models with multiple PFTs (e.g., Litchman et al., 2006;Tyrrell & Taylor, 1996) or models using nutrient-restoring methods (e.g., Jin et al., 2006) are better able to represent coccolithophore abundances and community composition with data-based (rather than model-based) parameter values; Jin et al. (2006) estimate coccolithophores contribute only 2% of NPP, which is in better agreement with the MAREDAT relative abundance estimate for coccolithophores (Buitenhuis et al., 2013a).Models with fully prognostic PFTs applied in ocean general circulation models (OGCMs) have a more difficult time reproducing phytoplankton calcifier biogeography and proportionality.Coccolithophores in NASA's ocean biogeochemical model (NOBM; five PFTs) (Gregg et al., 2003;Gregg & Casey, 2007) show an overall positive correlation with in situ data, though fail to appear in the North Pacific and Antarctic regions.As with the CAL   Buitenhuis et al. (2013a), with diazotrophs having the lowest concentration, followed by general phytoplankton.Diazotroph concentration is substantially lower than the general phytoplankton PFT and is able to remain extant because of its critical advantage of not being nitrogen limited.In CAL, as found by Buitenhuis et al. (2013a), zooplankton concentrations are higher than total phytoplankton concentrations.
Overall the CAL model advances UVic ESCM biogeochemistry by improving the mechanistic realism without sacrificing model performance with respect to nutrient and carbon distributions.As with any model, however, this one is not without caveats regarding its application.Collapsing complex and poorly understood natural biogeochemical cycles into a rigid artificial model structure introduces uncertainty into the parameter space of the constructed equations.The degree of underdetermination of the model equations is sufficiently large that a priori assumptions and optimization methods have been shown to influence results, with "optimal" parameter values comprising a broad range, each performing equally well with respect to independent data (Ward, Friedrichs, Anderson, & Oschlies, 2010).It is important to note that although this model has been tuned manually to reduce the model-data error in global state variables, it cannot be considered optimized (Kriest, Khatiwala, & Oschlies, 2010).Furthermore, nutrients, to a degree, and PFT distributions are especially sensitive to model structure and parameter choice (Anderson et al., 2010;Manizza, Buitenhuis, & Le Quéré, 2010;Sailley et al., 2013), as well as physical biases in any given ocean model (Doney et al., 2004;Najjar et al., 2007;Sinha et al., 2010).Similarly, models can perform comparatively well for very different structural reasons (Hashioka et al., 2013).It is, therefore, often difficult to determine whether the model is producing the right answer for the wrong reason (e.g., Friedrichs et al., 2007;Sinha et al., 2010).More specifically, biological parameter choice for the calcifying PFT is biased towards E. huxleyi, which necessarily biases model results.One must, therefore, be careful to interpret model results appropriately, given these limitations.

Conclusions
Calcifying phytoplankton and zooplankton are key components of the ocean carbon cycle and thus their representation in coupled climate models is important for understanding systemic response to change.This model is a unique attempt to include phytoplankton calcifiers as an explicit PFT alongside a general phytoplankton and diazotroph PFT in an intermediate complexity model and to make the phytoplankton calcifiers and zooplankton responsible for CaCO 3 production and prognostic export, as well as detrital ballasting.The UVic ESCM now fills a niche in Earth system modelling that was previously unoccupied in that it is relatively inexpensive to run, yet resolves the complete Earth system carbon cycle including prognostic CaCO 3 and a separate phytoplankton calcifier PFT.Because the UVic ESCM includes ocean sediments and calcite compensation it is now a model that is particularly well suited to reducing the uncertainty of the fate of emissions over the long term.It is now also well suited to testing the parameter space of feedbacks between the carbonate and carbon cycles and the climate system as transient simulations.The modifications maintain the UVic ESCM's performance with respect to nutrient distributions and carbon fluxes and make the model mechanistically more realistic.Primary production, export production, POC and CaCO 3 fluxes at various depths all fall within independent estimates.Though the model is able to reasonably reproduce observed patterns of mid-latitude maximum phytoplankton calcifier concentrations, it also shares biases common to other phytoplankton calcifier PFT models coupled to OGCMs: calcifiers are overrepresented in total biomass and in low latitudes, and underrepresented in high latitudes compared with satellite and sample data (Gregg & Casey, 2007;Sinha et al., 2010;Vogt et al., 2013).In the CAL model, failure to resolve coastal processes results in NPP, CaCO 3 , and POC export-production fluxes that are necessarily too high in the low latitudes in order to match global estimates.With possibly 48% of total global POC flux occurring in water depths less than 50 m (Dunne, Sarmiento, & Gnanadesikan, 2007), lacking any sort of parameterization for these regions imposes a significant bias to the model.In other phytoplankton calcifier multi-PFT models, exact regions of bias are model dependent and attributable to physical and ecosystem differences, but the systematic overrepresentation of phytoplankton calcifiers in low latitudes may have some physical justification.Previous studies have shown global export budgets require high CaCO 3 export in this region (Sarmiento et al., 2002), and high latitude bloom CaCO 3 is underrepresented in sediments (Brown & Yoder, 1994).Vogt et al. (2013) noted the similarity of phytoplankton calcifier model biogeography to observed picophytoplankton biogeography, so inclusion of additional picophytoplankton PFTs might improve phytoplankton calcifiers in models.
The UVic ESCM is now approaching a level of carbonatecarbon cycle complexity in which recent hypotheses regarding internally driven feedbacks in glacial-interglacial atmospheric CO 2 concentration changes can be tested in an Earth system model.Reduction of the Si:N export ratio in the Southern Ocean during glaciation leading to the expansion of diatoms at the expense of coccolithophores (Matsumoto, Sarmiento, & Brzezinski, 2002) can be tested using the prognostic iron cycle of Nickelsen et al. (2014) once diatoms and silicate are implemented (currently underway).Increased calcifier concentration with increased ocean alkalinity driving sawtooth-shaped global CO 2 time series (Omta, van Voorn, Rickaby, & Follows, 2013) can be tested with the implementation of increased calcifier growth or advantage with increasing carbonate saturation state.The role of temperatureenhanced phytoplankton growth (Fowler, Rickaby, & Wolff, 2013) in glacial-interglacial transitions can also be tested.This exercise reiterates the difficulty of simulating realistic CaCO 3 distributions because production and export depend on many physical, physiological, and ecological factors.There are five potential improvements to the CAL model that have not yet been addressed.Simulated phytoplankton calcifiers are wholly dependent on relative competitive advantage and can easily become extinct or cause the general phytoplankton PFT to become extinct with only small adjustments to production parameter values, especially the growth rate.Because their niche is so poorly defined with respect to the general phytoplankton PFT, additional PFTs (particularly diatoms) are expected to improve their population biogeography, stability, and seasonal behaviour and may allow phytoplankton calcifier parameter values to become less model-and more data-dependent (although this assumption has not been tested).Secondly, the ballast model does not include a parameterization for particle aggregation, which would increase the fraction of ballasted POC ending up in the sediments.Thirdly, static stoichiometric ratios in the model ignore their dependence on remineralization processes (Schneider, Schlitzer, Fischer, & Nothig, 2003), carbonate chemistry (Riebesell et al., 2007), biogeography (Weber & Deutsch, 2010), and taxonomy (Arrigo et al., 1999).Including a parameterization of flexible stoichiometric ratios would have a significant influence on the ecology (Flynn, 2010), nutrient distributions, and carbon uptake (Kortzinger, Koeve, Kahler, & Mintrop, 2001;Schneider et al., 2004).In a similar vein, experimentation with the parameterization of zooplankton (numbers of PFTs with variable grazing preferences, prey-switching, assimilation of consumed CaCO 3 , a variable rain ratio, unique CaCO 3 dissolution parameters, etc.) would also likely produce insight into model sensitivity to zooplankton assumptions.Lastly, the model does not account for decreasing calcification with increasing CO 2 concentration (Riebesell et al., 2000), which would doubtless affect simulated tracer distributions and biogeography.Sustained declines in calcification are questionable (Lohbeck, Riebesell, & Reusch, 2012) and were therefore omitted.Furthermore, using a single dissolution parameterization for zooplankton and phytoplankton calcifier CaCO 3 ignores the likely significant contribution of aragonite dissolution to global alkalinity (Gangstø et al., 2008).These last two were considered for inclusion in this model, but the current code structure is not amenable to flexible or multiple rain ratios and will require a significant restructuring should these changes be implemented in the future.
, with the Keller et al. (2012) model being referred to as NOCAL and this version being referred to as CAL.The model description here covers only the most relevant equations and equations that have changed in this newest version; please see Keller et al. (2012), Schmittner et al. (2005), and Schmittner et al. (2008) for a complete description of the other equations.b Model Description Tracer concentrations (C) vary according to:

Fig. 2
Fig. 2 Averaged biogeochemical simulated tracers (CAL, solid red line; NOCAL, dashed blue line) compared with observations (black line).DIC and alkalinityobservations are the standard GLODAP product(Key et al., 2004).Phosphate and nitrate observations are annual averages from the World Ocean Atlas (WOA;Garcia et al., 2009).Bottom row shows globally averaged model-data misfits.

Fig. 3
Fig. 3 Normalized Taylor diagrams of simulated tracers (CAL, red symbols; NOCAL, blue symbols) compared with observations (black circle).Observations used are as in Fig. 2. Ocean basins are denoted as global average (star), Atlantic (square), Indian (triangle), Pacific (plus), and Southern Ocean (diamond).The distance to the origin represents the normalized standard deviation.Normalized correlation with the observations is read from the azimuthal position.Perfect agreement with observations is a normalized standard deviation of 1 and a normalized correlation of 1.
, the most notable discrepancy between Atlantic observations and model results is in surface alkalinity concentrations (Fig 7), in which model alkalinity is too low in the northern hemisphere mid-latitudes and tropics.Surface DIC in the western Pacific Ocean is improved in the CAL version compared with earlier versions (not shown), though DIC in this region remains too low with respect to observations (Fig. 6) because of high model NPP.The concentration of CAL CaCO 3 peaks in latitudinal bands centred on 50°N, the equator, and 40°S (Fig 8).Limited data exist for this key model variable.Comparison with the Aqua

Fig. 8
Fig. 8 Zonally averaged CaCO 3 concentration by ocean basin (left and middle panels).Model CAL CaCO 3 concentration, including living CaCO 3 attached to phytoplankton calcifiers and zooplankton, in the surface grid box (to 50 m depth, top right panel).Bottom right is the standard CaCO 3 product from Aqua MODIS, accumulated over the entire mission (2002-2013; NASA, 2013).

Fig. 9
Fig. 9 Regression between modelled and observed CAL CaCO 3 concentration (left panel) and CAL POC concentration (right panel) from the ocean surface to 1000 m depth.Data are in situ measurements from Lam et al. (2011).

Fig. 12
Fig. 12 Percentage CaCO 3 sediment composition.CAL is shown in upper left panel; NOCAL is shown in bottom left panel, and gridded sample data from Archer (1996a) is shown in top right panel.