Natural ocean carbon cycle sensitivity to parameterizations of the recycling in a climate model

Sensitivities of the oceanic biological pump within the GISS (Goddard Institute for Space Studies ) climate modeling system are explored here. Results are presented from twin control simulations of the air–sea CO2 gas exchange using two different ocean models coupled to the same atmosphere. The two ocean models (Russell ocean model and Hybrid Coordinate Ocean Model, HYCOM) use different vertical coordinate systems, and therefore different representations of column physics. Both variants of the GISS climate model are coupled to the same ocean biogeochemistry module (the NASA Ocean Biogeochemistry Model, NOBM), which computes prognostic distributions for biotic and abiotic fields that influence the air–sea flux of CO2 and the deep ocean carbon transport and storage. In particular, the model differences due to remineralization rate changes are compared to differences attributed to physical processes modeled differently in the two ocean models such as ventilation, mixing, eddy stirring and vertical advection. GISSEH(GISSER) is found to underestimate mixed layer depth compared to observations by about 55% (10%) in the Southern Ocean and overestimate it by about 17% (underestimate by 2%) in the northern high latitudes. Everywhere else in the global ocean, the two models underestimate the surface mixing by about 12–34%, which prevents deep nutrients from reaching the surface and promoting primary production there. Consequently, carbon export is reduced because of reduced production at the surface. Furthermore, carbon export is particularly sensitive to remineralization rate changes in the frontal regions of the subtropical gyres and at the Equator and this sensitivity in the model is much higher than the sensitivity to physical processes such as vertical mixing, vertical advection and mesoscale eddy transport. At depth, GISSER, which has a significant warm bias, remineralizes nutrients and carbon faster thereby producing more nutrients and carbon at depth, which eventually resurfaces with the global thermohaline circulation especially in the Southern Ocean. Because of the reduced primary production and carbon export in GISSEH compared to GISSER, the biological pump efficiency, i.e., the ratio of primary production and carbon export at 75m, is half in the GISSEH of that in GISSER, The Southern Ocean emerges as a key region where the CO2 flux is as sensitive to biological parameterizations as it is to physical parameterizations. The fidelity of ocean mixing in the Southern Ocean compared to observations is shown to be a good indicator of the magnitude of the biological pump efficiency regardless of physical model choice.

Abstract. Sensitivities of the oceanic biological pump within the GISS (Goddard Institute for Space Studies ) climate modeling system are explored here. Results are presented from twin control simulations of the air-sea CO 2 gas exchange using two different ocean models coupled to the same atmosphere. The two ocean models (Russell ocean model and Hybrid Coordinate Ocean Model, HYCOM) use different vertical coordinate systems, and therefore different representations of column physics. Both variants of the GISS climate model are coupled to the same ocean biogeochemistry module (the NASA Ocean Biogeochemistry Model, NOBM), which computes prognostic distributions for biotic and abiotic fields that influence the air-sea flux of CO 2 and the deep ocean carbon transport and storage. In particular, the model differences due to remineralization rate changes are compared to differences attributed to physical processes modeled differently in the two ocean models such as ventilation, mixing, eddy stirring and vertical advection. GISSEH(GISSER) is found to underestimate mixed layer depth compared to observations by about 55 % (10 %) in the Southern Ocean and overestimate it by about 17 % (underestimate by 2 %) in the northern high latitudes. Everywhere else in the global ocean, the two models underestimate the surface mixing by about 12-34 %, which prevents deep nutrients from reaching the surface and promoting primary production there. Consequently, carbon export is reduced because of reduced production at the surface. Furthermore, carbon export is particularly sensitive to remineralization rate changes in the frontal regions of the subtropical gyres and at the Equator and this sensitivity in the model is much higher than the sensitivity to physical processes such as vertical mixing, vertical advection and mesoscale eddy transport. At depth, GISSER, which has a significant warm bias, remineralizes nutrients and carbon faster thereby producing more nutrients and carbon at depth, which eventually resurfaces with the global thermohaline circulation especially in the Southern Ocean. Because of the reduced primary production and carbon export in GISSEH compared to GISSER, the biological pump efficiency, i.e., the ratio of primary production and carbon export at 75 m, is half in the GISSEH of that in GISSER, The Southern Ocean emerges as a key region where the CO 2 flux is as sensitive to biological parameterizations as it is to physical parameterizations. The fidelity of ocean mixing in the Southern Ocean compared to observations is shown to be a good indicator of the magnitude of the biological pump efficiency regardless of physical model choice.

1138
A. Romanou et al.: Ocean biological pump and recycling LeQuéré et al., 2012) through the solubility pump. The net uptake flux of CO 2 by the global oceans is estimated to be about 1.6 ± 0.9 Pg Cyr −1 (Takahashi et al., 2009) with the Southern Hemisphere oceans being the largest CO 2 sink taking up about 1.1 Pg Cyr −1 , while the northern subtropical and midlatitude oceans take up about 0.7 Pg Cyr −1 . The equatorial oceans emit 0.7 Pg Cyr −1 to the atmosphere. Uncertainties in the solubility pump stem mainly from uncertainty in the pCO 2 distributions at the surface and the parameterization of the gas exchange at the air-sea surface. These range from 3 to 30 % with significant regional differences (Signorini and Mcclain, 2009;Feely et al., 2004).
There is greater uncertainty in our knowledge of the ocean biological pump (OBP), which describes the amount of carbon fixed by biology at the surface and then redistributed within the ocean (Denman, 2003;Buesseler et al., 2008;Henson et al., 2011), i.e., those processes, such as detrital sinking and remineralization, that remove surface carbon from the euphotic zone of the ocean and transport it to the deeper layers. It is estimated that this deep ocean carbon export is about 6 Pg C yr −1 Martin et al. (1987) to 10 Pg C yr −1 (Falkowski et al., 1998) usually depending on the type of observations while model results range from about 8 Pg C yr −1 to 11 Pg C yr −1 (Taucher and Oschlies, 2011;Dunne et al., 2007) depending on the model parameterizations.
To be able to reduce the uncertainty range and the biases in the model results with respect to the OBP, we need to understand the type of sensitivities different parameterizations of the OBP produce, and how these sensitivities compare to other model biases, e.g., due to model physics and how much of these sensitivities are model dependent. Several recent studies (Oschlies, 2001;Carr et al., 2006;Kwon et al., 2009;Taucher and Oschlies, 2011;Kriest et al., 2010) showed that the biologial pump has important sensitivities to model parameters such as temperature, remineralization depth, nutrient uptake/growth rate and sinking of particulate matter. These key process parameterizations directly affect the biological pump, and through changes in the surface pCO 2 distributions, they have indirect effects on the solubility pump.
Studies based on simple biogeochemical models with few components (Bacastow and Maier-Reimer, 1991;Kwon and Primeau, 2008) demonstrated that when the remineralization length scale is increased, or equivalently the remineralization rate is decreased, near the surface there are increases in the meso-and bathypelagic nutrient distributions and the nutrient transport, which traces the large-scale deep circulation pathways.
More complex biogeochemistry model studies showed that even small changes in remineralization depth lead to redistribution of recycled carbon from the meso-to the bathypelagic waters and lead to reduction of the atmospheric CO 2 (Kwon et al., 2009).
Preformed nutrients are remineralized in the epipelagic zone within the young North Atlantic Deep Water (NADW) and subsequently stored in the oldest waters near the ocean bottom. Therefore, the NADW is a major pathway of atmospheric CO 2 being sequestered into the deep ocean. Kriest et al. (2012) showed that global biases in observed nutrient distributions depend mainly on parameter values of specific processes related to recycling and sinking of nutrients, rather than on other processes in the biogeochemistry model. In their experiments, simple parameterizations of the remineralization and sinking led to large global biases from observations due to deficits in nutrient concentrations in the intermediate depths of the ocean due to excessively large export of organic matter to greater depths. Moreover, they found that even though global biases are low, this may not represent the overall skill of the parameterization, as regional biases often cancel each other. North Pacific deep waters are very old and therefore all biases, related to biology and physics model deficiencies, slowly accumulate there over time . Changes in these parameterizations account for most of the model biases in the North Pacific. Therefore, the simpler models, models with fewer components, have better skill in the representation of the surface properties in this region, at depth nonlocal errors (i.e., model errors from remote regions) are larger. Kriest et al. (2012) also find that adjusting these parameters may enhance phytoplankton production particularly in the high nutrient low chlorophyll (HNLC) regions (e.g., the Southern Ocean, the equatorial regions or the northern North Pacific regions) but this may compensate for the lack of other nutrient limitations in the model (e.g., iron).
Several studies (Kriest et al., 2010;Kriest and Oschlies, 2011) have drawn attention to the role of physical process parameterizations in the models such as vertical mixing, or even numerical diffusion by the model advection scheme, (e.g., the upstream advection operator) in transporting some of the particulate organic matter to depth, which may either compensate or enhance biases due to biological parameterizations. Similarly, low model vertical resolution leads to increased particle flux via increased numerical mixing, and therefore leads to more recycling at depth, compared to a model with a finer vertical resolution.
Finally, improving the model fidelity to observed nutrient distributions by adjusting the parameterizations of the vertical particle flux may degrade simulations of other fields such as primary production (PP) (Carr et al., 2006;Kriest et al., 2010). In these studies, it is shown that a model that predicts PP skillfully may not perform as well in simulating nutrient concentrations.
In this paper we will focus on one biogeochemistry model and vary the remineralization rate within two different ocean models coupled to the same atmosphere. Since we are able to use two ocean models that represent different classes of ocean modeling (z coordinate vs isopycnal) we are able to assess the relative importance of the biological pump sensitivities to biological and physical parameterizations.
The paper is structured as follows: Sect. 2 includes information on the GISS climate model configurations with the two ocean models and the biogeochemistry module. Results from two sets of twin runs where we vary the remineralization rate are presented in Sect. 3. We distinguish effects of the remineralization rate change on the biological pump (Sect. 3.1), i.e., detritus, nitrates, total chlorophyll, primary production fields and deep carbon export, and effects on the gas exchange pump (Sect. 3.2), i.e., the dissolved inorganic carbon (DIC), pCO 2 distributions and the air-sea flux of CO 2 . Our results are summarized and are placed in the context of similar studies in the Discussion section (Sect. 4) and lastly some conclusions are offered in Sect. 5.

Model description
The carbon cycle runs are described extensively in Romanou et al. (2013) and have been submitted to the Coupled Model Intercomparison Project-Phase 5 (CMIP5) as GISS-E2-R-CC and GISS-E2-H-CC, which are the carbon cycle (CC) versions of model E2 coupled to the Russell (R) or HYCOM (H) ocean models. All carbon cycle runs started from equilibrium conditions of coupled atmosphere-ocean-ice 3000 yr runs, and were subsequently integrated again with fully prognostic carbon cycle to equilibrium for an additional 300 yr.
Here we present only preindustrial control CO 2 simulations where the atmospheric pCO 2 responds to the ocean, but the global average is fixed to 285.2 ppmv (parts per million by volume). A very brief description of the two ocean and the biogeochemical models is given here. The Russell ocean model is a non-Boussinesq, massconserving ocean model with up to 32 vertical levels and 1 • × 1.25 • horizontal resolution. The vertical coordinate is mass per unit area because mass is formally conserved in the model and therefore introduction of freshwater fluxes can be easily handled. Normalizing to grid box area ensures that computations and comparisons of fluxes to/from adjacent boxes are properly handled. The model has a free surface and natural surface boundary fluxes of freshwater and heat whose prognostic variable is potential enthalpy. It uses the Gent and McWilliams (1990) scheme for isopycnal eddy fluxes and isopycnal thickness diffusion. Vertical mixing is done according to the KPP (K-Profile Parameterization) vertical mixing scheme (Large et al., 1994) with nonlocal mixing effects included. Vertical advection is computed using a centered differencing scheme.
The Hybrid Coordinate Ocean Model (HYCOM, Sun and Bleck (2006) variant C) uses pressure coordinates near the surface and isopycnals in the ocean interior. Horizontal resolution is 1 • × 1 • cos(φ), where φ is the latitude, vertical resolution of 26 coordinate layers. The vertical mixing scheme is KPP (Large et al., 1994) near the surface and following McDougall and Dewar (1998) in the isopycnic interior. Deep convection at high latitudes is parameterized through a convective adjustment scheme. Brine rejected during freezing is uniformly distributed in the uppermost 200 m when water depth is greater than 1000 m or the total depth otherwise. The coefficient of lateral T /S diffusion is increased near the Equator to account for the tropical instability waves, which are very diffusive and which are not resolved by the coarse grid. The Gent-McWilliams parameterization is implemented via interface smoothing, and its effect therefore is limited to the isopycnic coordinate domain (which excludes the near-surface waters). In the present implementation, only mass fluxes at even time steps are used in building up the mass flux time integrals. This is to avoid leapfrog-related inconsistencies between layer thickness tendencies and horizontal flux fields, which would generate spurious diapycnal mass fluxes. Their sum is then used at the end of a longer time step (1/2 day) to rigorously reconstruct all terms of the timeintegrated continuity equation. (The other terms are the net change in layer thickness over the long time step and the vertical mass fluxes that are derived diagnostically from the others.) This allows us to solve the flux form of the tracer advection equation over a large time step while rigorously conserving everything in an environment of perpetually shifting layer configurations. The single remaining inconsistency is the one caused by time smoothing of the thickness field. It spawns small diapycnal fluxes (and hence interlayer tracer leakage) even in the absence of physically based vertical mixing processes. The numerical method used for 3-D tracer advection is FCT (flux corrected transport), which maintains positivedefiniteness and monotonicity in the tracer fields. The horizontal fluxes, before being clipped, are 2nd order centered, and the vertical fluxes (again before clipping) are computed using the piecewise parabolic method (PPM). Solving transport equations in flux form does not rigorously conserve total tracer when layer thickness is allowed to go to zero. Conservation errors are currently eliminated by applying a global correction factor to each tracer field following each transport step. With this corrective device in place, tracer is rigorously conserved in the model.
Most relevant to the present study is the fact that the different coordinate system in the two ocean models produces different distributions of seawater properties below the mixed layer depth where numerical and spurious diapycnal mixing can become important (Bleck et al., 1992;Griffies et al., 2000). Other relevant differences between the two ocean models are the extent of implementation of the vertical mixing scheme (KPP); in GISSER KPP extends to the bottom of the water column but nonlocal terms are not included, whereas in GISSEH it only reaches the bottom of the mixed layer. Eddy mixing along isopycnals (Gent and McWilliams, 1990) is only explicitly modeled in GISSER, whereas in GIS-SEH it is approximated by isopycnal thickness diffusion for the regions where these isopycnals do not outcrop. Moreover, there is no tracer horizontal advection in GISSEH. Although the models are different in other ways, the significantly decreased vertical diffusion/mixing below the mixed The interactive ocean carbon cycle model is the NASA Ocean Biogeochemistry Model (NOBM, (Gregg and Casey, 2007;Romanou et al., 2013), which is interactively coupled to the ocean model from which it uses temperature, salinity, mixed layer depth, the ocean circulation fields, horizontal advection and vertical mixing schemes and the atmospheric model from which it uses the surface shortwave radiation (direct and diffuse) and surface wind speed. NOBM includes four phytoplankton groups (diatoms, chlorophytes, cyanobacteria, and coccolithophores), four nutrients (nitrate, silicate, ammonium and iron), three detrital pools (nitrates/carbon, silicate and iron) and a single herbivore component. The carbon submodel includes dissolved organic and inorganic carbon (DOC and DIC, respectively) and surface alkalinity, which is a constant function of salinity at the surface. The underwater irradiance needed to drive phytoplankton growth is calculated explicitly for the entire spectrum of light with the use of a radiative transfer model that takes into account water absorption and scattering based on the optical properties of the phytoplankton groups that were derived from laboratory studies (Gregg, 2002). Detrital settling and phytoplankton sinking are modeled in HYCOM using a onedimensional piecewise parabolic transport method whereas in the Russell ocean model using a piecewise linear method.
In Gregg and Casey (2007) phytoplankton sinking was modeled as an additional advection in the vertical. Sinking rates were specified at 31 • C and derived from Stokes' law using representative phytoplankton sizes for each group from a very thorough survey of the literature (Gregg and Casey (2007) and references therein). The sinking rates w s of each phytoplankton group were adjusted to account for viscous effects, according to Stokes' law as a function of temperature T : (1) Coccolithophores were assumed to be able to sink at speeds that were a function of the growth rate from 0.3 to 1.4 m d −1 based on observations by Fritz and Balch (1996). A linear relationship was assumed: where w s is the sinking rate of coccolithophores (m d −1 ) and μ(high) is the maximum growth rate actually achieved for the previous day. For three of the nutrient groups, i.e., nitrate, silica, and dissolved iron, corresponding detrital pools are defined in order to simulate the remineralization and return to dissolved form. The equations for these detrital pools are given in the Appendix of Gregg and Casey (2007), Eqs. (A7-A9). Only carbon detritus is explicitly modeled in NOBM since the C : N ratio is constant and conversion to the nitrate pool is straight forward. Detritus pools are changed due to advection, diffusion and sinking. Similar to phytoplankton sinking, detrital sinking depends on viscosity and is a function of temperature. Remineralization is also temperaturedependent and varies with the phytoplankton growth rate for each phytoplankton group. The remineralization rate for nitrogen/carbon is 0.1 d −1 , for silica is 0.002 d −1 and for iron is 0.7 d −1 .

Results
Results are presented below for two sets of model runs for each climate model configuration. GISSER-l(h) stands for Russell ocean model with low(high) remineralization rate and GISSEH-l(h) stands for HYCOM and low(high) remineralization rate.
Remineralization rates in the ocean are highly uncertain and vary significantly regionally as the VERTIGO experiment has demonstrated (Buesseler et al., 2008). Tian et al. (2000) showed that regionally (northern Atlantic subtropical gyre) the remineralization rate in the upper 400 m can be as high as 0.8 d −1 . Christodoulaki et al. (2013) suggested a value of 0.45 d −1 . At the same time, little is known about how and how much such large regional differences affect estimates of the biological pump and the air-sea gas exchange. We chose here to present results for two cases: a "low" remineralization rate of 0.1 d −1 , which is the value suggested by Denman (2003), and a "high" rate of 0.5 d −1 . Analysis is extended to cover more cases in order to investigate the robustness of our conclusions, and the parameter range extends from 0.01 to 0.5 d −1 .
Both NASA-GISS ocean-atmosphere climate models are coupled to NOBM and run out to equilibrium under constant preindustrial atmospheric CO 2 concentrations of 285.23 ppmv. Equilibrium is achieved when the gas exchange between the ocean and the atmosphere stabilizes to a value about 0.9 Pg Cyr −1 . The differences in the biological and gas exchange pumps between the two models are then evaluated. Detailed model assessment is given in Romanou et al. (2013). Here we will focus on the differences between the two physical models, particularly in the regions/cases where there are large biological pump differences as well.

Large-scale mixing
To cast the biological pump sensitivity in light of the physical model differences, we review some of the main findings of Romanou et al. (2013) where a more detailed model evaluation was presented. In that paper, the model ocean carbon cycle is shown to be very sensitive to biases at the surface but also the interior of the ocean. Surface thermal, radiative and momentum biases give rise to mixed layer depth biases whereas deeper ocean vertical mixing differences explain the the lower nutrient load that upwells to the surface to sustain the primary production in the GISSEH model. Figure 1 shows a regional comparison and evaluation of the model mixed layer depth (MLD). Both models underestimate MLD in the Southern Hemisphere, particularly in the Southern Hemisphere's midlatitudes, where the temperature bias in GISSER is 2 • C and in GISSEH is 3 • C (Fig. 2). In Romanou et al. (2013) it was shown that wind speeds over the Southern Ocean are underestimated in both models ( Fig. 2) which, compounded by overestimated heat losses in the winter and cloud coverage during summer, lead to shallower mixed layer depths.
Below 100 m and down to 3000 m ( Fig. 2), GISSER exhibits significant warm bias whereas GISSEH is in better agreement with the observations. The GISSER deep warm bias is due to significant eddy transport in the Southern Ocean.

Detritus
In NOBM, detritus, i.e., the nonliving part of the particulate organic carbon (POC) pool, is changed by advection, diffusion, detrital sinking, remineralization, biological production and bacterial degradation according to where, D is the detritus concentration, K is the vertical diffusivity, V is the horizontal current velocity, w is the vertical advection (sinking) speed, R is a temperature dependence term, α is the remineralization rate at 20 • C, is the carbon : chlorophyll ratio, κ is the senescence rate, P i is the vector of phytoplankton species, where i is for diatoms, chlorophytes, cyanobacteria, and coccolithophores, η 1 and η 2 are heterotrophic loss rates, H is the herbivore concentration, is an excretion rate for nitrate/carbon and λ is the detrital breakdown rate at 20 • C. The values of the parameters are given in Gregg and Casey (2007).
As we increase the remineralization rate, the sink term, RαD, on the right-hand side of Eq. (3) increases, and that leads to decreased detritus concentrations at all levels but more surface intensified (Fig. 3). Most of the reduction occurs in the Southern Ocean north of 60 • S, within the winddriven divergence zone and up to the subtropical convergence zone (STCZ) along 40 • S approximately. The distributions of detritus in the Southern Ocean upwelling and divergence zone reflect the effect of changes in remineralization rate and the combined effect of the thermohaline circulation pattern: as we increase the remineralization rate in the North Atlantic convection regions, less detrital material reaches deeper waters, therefore less is remineralized back to inorganic forms along the deep water pathways and eventually the NADW carries less detritus as it upwells in the Southern Ocean divergence zone.
In the Southern Ocean, the GISSEH model shows the largest detritus reduction in the region that lies between the Antarctic Polar Front (APF) at about 50 • S and the STCZ (≈ 40 • S, Fig. 3). Similar reductions are shown in the GISSER model but over a wider region (70 • -40 • S). The reduction region in GISSEH is narrower than in GISSER, because the former model has significantly shallower mixed layer depths in the region 70 • -50 • S since the temperature there is much higher than in GISSER (see Fig. 3 in Romanou et al. (2013)). Higher temperatures in GISSEH lead to higher remineralization rates and therefore lower detritus distributions than in GISSER (see Fig. 17 in Romanou et al. (2013)).
A close inspection of Eq. (3) shows that in NOBM the remineralization sink term scales as the local temperature profile, through the Bissinger et al. (2008) parameterization, which is a refinement to the Eppley curve (Eppley, 1972) for most ocean areas.
If we define "effective remineralization" as α eff = Rα: where α = 0.1 d −1 (the low remineralization case), then at the surface ( Fig. 4 top panels) the remineralization term is lower in the tropical upwelling regions and in the polar regions where the surface waters are colder than in the rest of the ocean. The effective remineralization distributions follow the SST (sea surface temperature) climatologies, shown in Fig. 2 of Romanou et al. (2013). Fronts in temperature correspond to fronts in the remineralization rate. The two models row panels) the differences are larger, up to 50 %, for most of the ocean. GISSEH mixes down heat more vigorously in the deep convection region of the North Atlantic, in the northwestern Pacific and in the northern Indian Ocean. This effect is reversed below 1000 m (Fig. 4e, f) where the GISSER ocean is warmer than the GISSEH ocean except in a narrow latitude band adjacent to Antarctica. Therefore effective remineralization is 20-60 % larger in GISSER than in GISSEH. In the abyssal ocean, at 3000 m (Fig. 4g, h), effective remineralization is about 30-40 % higher in GISSER than in GISSEH, except in the narrow latitude band against the coast in Antarctica where GISSEH effective remineralization is about 10 % larger than GISSER. Depicted in the right-hand column in Fig. 4 are the differences between the effective remineralization in GISSEH and GISSER. Away from polar latitudes, GISSER is warmer and therefore has higher remineralization rates at all depths. Only in the North Atlantic subtropical and midlatitudes is GISSEH warmer at depth (down to 1000 m at least) and therefore has higher remineralization rates. The settling speed term in NOBM is parameterized as follows: for ice-free ocean. The rate wsd is 20 md −1 for particulate nitrogen/carbon.
The vertical flux of detritus (normalized to the flux at 75 m) is very sensitive to the remineralization rate changes. Both model flux profiles follow the Martin curve (Fig. 5) for the lower remineralization rate cases. Agreement with the Martin curve, i.e., the exponential 0.86 form, extends to 400 m depth for the GISSER model and to 200 m depth for GISSEH. Below these depths, however, the models underestimate the Martin curve significantly. This result implies that perhaps our models remineralize detritus much faster than they allow it to sink. However, as described in Buesseler et al. (2008), there may be large variation in the Martin exponent magnitude in the world ocean, at different locations and depths. Other studies have also shown that the exponential may vary anywhere in between 0.3 and 2 (Kwon et al. (2009) and references therein).

Nutrients
Similar to detritus, nutrients and in particular nitrates are modeled as follows: with terms analogous as in Eq. (3) and, in addition, C : N is the carbon to nitrogen Redfield ratio, b N is the nitrogen : chlorophyll ratio, μi is the phytoplankton maximum growth rate for each species i (i.e., diatoms, chlorophytes, cyanobacteria, and coccolithophores) and f (NO 3 ) is the nutrient uptake function. The parameter values are the ones used in Gregg and Casey (2007); namely, C : N ratio = 106 / 16, nitrogen : chlorophyll ratio = 50 (C : N ratio), phytoplankton maximum growth rate for diatoms = 2.00 (in units d −1 and at 20 • C), for chlorophytes, diatom growth rate = * 0.84, for cyanobacteria, diatom growth rate = * 0.670 and for coccolithophores (E. huxleyi only), diatom growth rate = * 0.755.
Because detritus at depth is a major source of remineralized nutrients (Eq. 2), we expect that increased remineralization rate leads to decreases in the amount of detritus at depth, and subsequent reductions in the recycled nutrient load of the deep waters. Therefore Fig. 6a and b shows lower concentrations of nitrates at the surface in the deep convection regions of the Northern Hemisphere and in the Southern Ocean, extending from the coast of Antarctica to the SAF, because of mixing with nutrient depleted deep water. Detritus-strapped water masses that ventilate the world ocean produce less recycled nutrients and result in lower nutrient distributions in the main upwelling zone in the Southern Ocean. This is a robust result in both models.
However, unlike detritus, north of the STCZ, nutrient concentrations are higher in the high remineralization rate cases. This is because these waters originate from the tropical ocean upwelled surface and are therefore richer in remineralized nutrients. In fact, the increases in remineralized nutrient supply at the surface along the STCZ and the South American west coast are twice as high as the background values.
Furthermore, the change in nutrients due to recycling changes in the Pacific Ocean is not robust in the two models (Fig. 6a, b). GISSER gives more nutrients in the upwelling zones at the Equator than GISSEH. The reason for that is that GISSER has much stronger coastal upwelling due to stronger winds over the cold tongue region (Romanou et al., 2013), which brings more nitrate to the surface.
As a result, increasing remineralization causes increase of nutrients in the equatorial eastern Pacific in GISSER but not in GISSEH where the distributions are already too low.
Similarly, in the northeastern Pacific, higher remineralization rates in GISSER than in GISSEH, due to higher SSTs in this model, result in higher nitrates in the Russell model but lower in the HYCOM (Fig. 6a, b). The mixing and the mixed layer depth in GISSEH extends deeper than in GISSER and the combination of remineralization and local mixing result in GISSEH showing a reduction in nutrients as we increase the remineralization rate.
In the subpolar North Atlantic, increasing the remineralization rate reduces nutrients at the surface because deep mixing (through seasonal convection) lowers surface values, since less recycled material reaches the deeper ocean ( Fig. 6a, b).
The vertical distributions of nitrates show that increasing the remineralization rate produces a deeper, less steep nutricline because more recycled nutrients are generated in the upper 1000 m due to more rapid detritus remineralization. Vertical mixing of nutrients due to the local nutrient gradient at the base of the nutricline is enhanced. The increase at the surface occurs mostly in the Southern Ocean (Fig. 6c, d) and amounts to about 20 % of the background values.
Therefore, mesopelagic reductions (at depths of 100-1500 m, Fig. 6c, d) are a robust feature in both models, although in GISSEH the vertical extent is shallower and not as pronounced. The largest reductions in nutrients occur in the region between the SAF and the STCZ in each model where the nutricline slopes upward. The GISSEH subsurface nutrient reductions are more pronounced, because in GIS-SEH upwelling extends deeper than in GISSER. The region of the largest sensitivity to changes in the remineralization rate is the mesopelagic "twilight" zone as has been previously pointed out by Buesseler et al. (2008).
These results are robust across the models particularly in the 30 • -50 • S belt. Thus, in this area nutrient distributions are primarily controlled by sensitivities to biological parameterizations rather than differences between the two physical models, which mostly are due to physical parameterizations.

Total chlorophyll, primary production and deep carbon export
Chlorophyll distributions are determined by light and nutrient availability and the temperature of seawater. GISSEH predicts higher chlorophyll concentrations than GISSER (Fig. 7) in two main regions: (1) the tropical upwelling regions, where water originates from deeper, generally richer in nitrates, layers and (2)  gyres because of steeper nutriclines (as well as thermoclines/haloclines), which prohibit vertical mixing with the deeper layers, resulting in unfavorable conditions for phytoplankton growth.
In the Southern Ocean, in the region between the SAF and the STCZ, increasing the remineralization rate leads to increased nutrient supply, which in turn results in 30-50 % increase in total chlorophyll in both ocean models (Fig. 7). South of the SAF, however, increased recycling of nutrients reduces nutrient concentrations (Fig. 6c, d) which in turn leads to lower surface chlorophyll concentrations. For similar reasons, in the high latitude Atlantic ocean, phytoplankton is also reduced (see Figs. 6a, b and 7e, f).
In the tropical/subtropical ocean (40 • S-50 • N), chlorophyll changes with remineralization are robust: both models agree quantitatively in the magnitude of chlorophyll as well as the fact that increases in remineralization rate lead to increased chlorophyll (Fig. 8a). In the subpolar regions, however, the two physical models, with different high latitude treatment of deep convection, ice formation and brine rejection, as we have seen in the Model Description section, lead to quite different distributions of chlorophyll. It is found that at high latitudes, the chlorophyll sensitivity to the physi-cal parameterizations is larger than the sensitivity to changes in the remineralization rate. Overall, increasing the remineralization rate tends to increase surface chlorophyll, except in the Southern Ocean south of the SAF, in the region 40-70 • S, because of lower nutrient supply there (Fig. 6a, b) from the upwelling deep water.
Primary production (PP) is very different in the two models. HYCOM has consistently lower PP than the Russell model, because of the much shallower mixed layer depths. However, the response to increased recycling rates is robust in both models, with primary production increasing with the recycling rate in all the regions (Fig. 8b). These increases are important in the Northern Hemisphere's frontal regions 30-60 • N, in the subtropical gyres, because of the more readily available remineralized nutrients there. Little change is found in the Southern Hemisphere's high latitudes, south of 40 • S, where changing the remineralization rate does not impact PP (Fig. 8b). Diatoms are less influenced by the remineralization rate changes in both models. Cyanobacteria and coccolithophores are the most sensitive particularly in the tropical/subtropical regions and chlorophytes mostly in the Northern Hemisphere's subtropical front region. The reason for the disparate response of the chlorophyll species to recycling changes is that growth, which is the process that is affected most by remineralization changes, depends on the geographical distributions of species and nutrient. Diatoms, which are mostly found in the high latitude oceans in our models, benefit less from the increase in nutrients. Coccolithophores and cyanobacteria, which are more abundant than diatoms particularly in low and midlatitudes, benefit more from the increases in nutrients when the remineralization rate increases. Carbon export, i.e., the flux of particulate organic material across 75 m depth, illustrates the strength of the biological pump. The carbon export is very sensitive to both physical and biological parameterizations (Fig. 9a, b). Altering the recycling rate in each model leads to changes in carbon export. These changes are more pronounced in the GISSER model, however, the relative changes between the high and the low remineralization rate runs are of the same order of magnitude in both models (Fig. 9a). We find that an 80 % increase in the remineralization rate leads to about a 50-70 % reduction in the carbon export. Carbon export decreases as the remineralization rate increases, which is due to the reduction of particulate organic matter concentration at depth as more of it remineralizes near the surface. The largest differences between the high and low remineralization runs occur at the frontal regions of the subtropical gyres and at the Equator, which are the same regions where detritus differences are found to be largest (Fig. 3). In these regions, the two physical model responses are closer than the differences in response due to the biogeochemistry. This becomes apparent from the comparison of the vertical fluxes across the 75 m isobath (Fig. 9b). These vertical fluxes are associated with processes such as  vertical mixing (KPP), vertical advection (which is also sensitive to the advection scheme) and the Gent-McWilliams mesoscale eddy mixing. It is found that in the Russell ocean model, all physical circulation related vertical fluxes are 1-2 orders of magnitude smaller than the carbon export at that level.
Additional model runs with different remineralization rates were performed to examine the sensitivity of the ocean biological pump to a wider range of values for this parameter. The results show that both models (GISSER and GIS-SEH) exhibit the same sign of the sensitivity, although the magnitude is different (see Fig. 10a, b) for the primary pro- The bars are estimates from the models and the symbols represent analytic solutions obtained using Eq. (7) for each model, PP from diatoms, PP from chlorophytes, PP from cyanobacteria and PP from coccolithophores. duction and the carbon export: increasing the remineralization rate leads to increased primary production and decreased carbon export. Primary production responses (Fig. 10a) show a tendency towards saturation at remineralization rate values of about 0.3 and higher, whereas carbon export does not (Fig. 10b). The reason for this is that there is a limit of nutrient remineralization increase, beyond which, limitations of other nutrients (iron and silica) start playing a role. Carbon export however continues to decrease as remineralization rate is increased, because particulate organic carbon continues to decrease.
Estimates of the f ratio (or carbon export efficiency, i.e., ratio of export at 75 m to the vertically averaged PP) are shown in Fig. 10c: the carbon export efficiency in both models decreases as we increase the remineralization rate. Model carbon export efficiency estimates (0.1-0.7) fall within observational estimates as reported in Laws et al. (2000) and Henson et al. (2012).
Both models show similar rates of decrease of the f ratio with increasing remineralization although the magnitude of the f ratio itself is very different in the two models. To further investigate the sensitivity of the carbon export efficiency to model parameterizations, both physical and biological, we used parametric curve fitting that describes the functional relationship between the f ratio and the remineralization in the two models. Initially, two parametric equations were obtained, one for each physical model. The two sets of coefficients were then expressed in terms of a single metric of the physical models, here chosen to be the model mixed layer depth normalized by deBoyer Montégut et al. (2004) observed mixed layer depth. The choice of mixed layer depth in expressing the dependence of export efficiency on remineralization is probably not unique, but was dictated by the fact that MLD is an expression of the physical model's internal dynamics and ability to redistribute tracers. Different forms of the mixed layer depth such as the globally averaged mixed layer depth, the North Atlantic mixed layer depth, which is associated with the strength of deep convection in that region, and the Southern Ocean mixed layer depth, which is the region were remineralized nutrients resurface, were tried. The metric that provided us with the best fit in the least square sense to both model estimates of the f ratio at all the remineralization rate values was the normalized Southern Ocean averaged mixed layer depth (see Eq. 7). Equation (7) describes how the biological pump efficiency, K bp , changes with the remineralization rate α: where rMLD is the Southern Ocean normalized mixed layer depth (annual mean, regional average 70 • -40 • S normalized by observations) in each model. Equation (7) implies that although there is a tight relationship between the remineralization rate and the carbon export efficiency, this is modulated by how well the physical model depicts mixing in the Southern Ocean. It would be interesting to see if this or similar functional dependencies of the biological pump efficiency to physical and biological parameterizations are valid in other biogeochemical models.

DIC and pCO 2
Changes in remineralization impact CO 2 gas exchange between the ocean and the atmosphere through changes to surface DIC and pCO 2 . Similar to nitrates (Fig. 6a, b), in both models increasing the remineralization rate leads to increases of DIC in the region between the SAF and the STCZ (Fig. 11), again due to the fact that carboclines get steeper (as do nutriclines), and therefore mixing with underlying waters is suppressed. In the Southern Ocean divergence zone, DIC increases are small as a percentage of the total DIC (about 10 %, Fig. 11) but they are of the same order of magnitude as the local meridional gradient. South of the SAF, towards Antarctica, there are lower concentrations of DIC in the high remineralization case, due to the upwelling of deep carbondeficient water. As in nitrates, the increase of surface DIC with the remineralization rate in the tropics is found only in GISSER, and is due to the more nutrient rich waters in GISSER than in GISSEH. In the subtropical/equatorial Pacific the two models disagree. GISSER predicts increased DIC at the surface in the runs where recycling is increased while GISSEH shows a weak reduction. The disparity is attributed to two competing processes that set the DIC distribution at the surface when we increase the remineralization rate: remineralization sources vs. vertical mixing and sinking. As the surface remineralized DIC distribution increases as we increase the rate of remineralization, the subsurface profile of DIC changes as well due to sinking. As we increase the remineralization rate, subsurface DIC decreases. Vertical mixing (or upwelling) will then dilute the surface distributions with deeper and lower distributions of DIC. Here, GISSER has higher temperatures than GISSEH and therefore more remineralized DIC at the surface. As we increase the remineralization rate, initially both models will have more remineralized DIC at the surface; but the GISSER depth profile of remineralized DIC will be more surface enhanced and GISSEH will have more remineralized GISSER below the surface. Vertical mixing comes into play then, as GISSEH has deeper mixed layers (see Fig. 7 in Romanou et al. (2013)) and therefore mixing with subsurface DIC is such that GISSER ends up with more DIC whereas GISSEH ends up with less DIC. Similar to nitrates (Fig. 6c, d), increasing the remineralization rate leads to decreased concentrations of DIC at mesopelagic depths (100-2000 m, Fig. 12). As detritus is remineralized back to inorganic carbon, less of it descends through the water column, reducing the source of remineralized carbon at depth. This is a robust behavior, evident in both physical models. In GISSER, increasing the remineralization rates deepens the carbocline and makes it more diffuse. This effect extends down to 1500 m depth. In GISSEH, there is a similar effect but it only extends down to 500 m.
In both models, as we increase remineralization rates, the carbocline gradient at mesopelagic depths decreases and therefore more DIC can reach the surface through vertical mixing, which explains the increased DIC distributions there. This is mostly evident north of 50 • S (i.e., north of the SAF). The changes in partial pressure of CO 2 at the surface of the ocean (pCO 2 ) mirror the changes in the surface DIC with the remineralization rate.

Equilibrium flux of CO 2
In the model runs described here, changes in the surface pCO 2 depend on changes in DIC due to varying the remineralization rate. Surface alkalinity is modeled as a constant function of surface salinity and is also affecting the flux. However, since we do not simulate the carbonate pump in the model we do not account for changes in alkalinity at depth. Since both models are at equilibrium, T, S and surface winds do not change within runs of the same model as we vary the remineralization rate. Therefore we are able to examine how much gas exchange depends on changes in DIC and the biological pump. Note that in the following, the flux sign convention is positive down, i.e., positive air-sea flux denotes oceanic uptake of CO 2 .
The Southern Ocean sink (Fig. 13) is substantially reduced in the fast recycling case. This reduction occurs mainly because the uptake region (between 40 • and 60 • S) is significantly reduced while although, the source region south of 40 • S is also reduced, this reduction in outgassing is not sufficient to outweigh the reduction of the actual sink more to the north. Carbon uptake is reduced because increased remineralization leads to higher DIC and pCO 2 at the surface at the divergence zone of the Southern Ocean. South of the SAF, outgassing is reduced because although surface DIC (and pCO 2 ) is reduced, due to carbon deficient NADW, the imbalance with the atmosphere still favors outgassing, but less so than in the low remineralization case. At the same time, tropical outgassing is similar in the two models, although slightly enhanced in GISSER due to more upwelling DIC. The North Atlantic sink strengthens because surface DIC decreases as increased remineralization leads to less carbon at depth. Therefore more CO 2 is sequestered into the ocean. On the other hand, the North Pacific sink of CO 2 gets stronger in the GISSER model when we increase the remineralization rate, but the opposite happens in the GISSEH model.
As a global average, the outgassing flux increases as we increase the remineralization (region Glb in Fig. 14). This increased outgassing occurs mostly in the Southern Ocean (ANT in Fig. 14) and in the equatorial regions (EIN, EPA and EAT in Fig. 14). The North Atlantic sink increases in the GISSER model but slightly decreases in GISSEH because mixing in GISSEH brings much less DIC near the surface. . Red denotes uptake, blue denotes outgassing. In unshaded regions the flux is low. These are results from the two models GISSER and GISSEH using different remineralization scales. Hi refers to high remineralization rate (i.e., fast recycling of nutrients) and lo refers to low remineralization rate.

Discussion
In this paper we assess the sensitivity of the oceanic carbon cycle and the air-sea CO 2 exchange to changes in the recycling rates within the ocean. We have also compared these changes with model uncertainties due to varying parameterizations of physical mechanisms such as vertical advection due to large-scale flow convergence/divergence, mesoscale eddy stirring and turbulent mixing. For this we ran the GISS climate model coupled to two different ocean models and we varied the remineralization rate gradually from 0.01 to 0.5 d −1 . The two ocean models used in this study differ in the vertical coordinate system and therefore simulate diapycnal and isopycnal processes differently. The main result obtained is that carbon cycle sensitivities to the biological pump parameterizations are of the same order of magnitude or larger than the sensitivities due to physical model parameterizations. Consequently, equal attention should be placed in assessing the fidelity of parameterizations such as remineralization scales, sinking rates and others that may affect the biological pump as to those that aim to improve the physical ocean skill. Specifically, organic particulate sinking reductions due to increased recycling in the epi-and mesopelagic ocean are larger than uncertainties in parameterizations of physical processes, particularly in subtropical frontal regions in the Northern and Southern Hemi-spheres. This leads to decreased nutrient remineralization at these depths and nutrient depletion of deep water masses such as the NADW. Primary production increases as we increase the remineralization rate while carbon export at 75 m decreases. CO 2 outgassing is enhanced as we increase the remineralization rate although CO 2 uptake increases in the SAF region.
In this study, the Southern Ocean CO 2 sink emerges as particularly sensitive to the parameterizations of the biological pump, and this is a robust feature in both models. In particular, waters that upwell in the Southern Ocean and which originally mixed and subducted in the Northern Hemisphere (NADW) carry less nutrients and DIC in the runs where regeneration rates are higher. The Southern Ocean between APF and SAF supports less chlorophyll and outgasses less CO 2 to the atmosphere in these runs as well. North of the SAF in the Southern Ocean, between SAF and the subtropical front (STF), surface DIC is increased and the sink region there is reduced. Changes in the biological pump in the Southern Ocean's and the Northern Hemisphere's frontal regions lead to substantial changes in the surface DIC distributions and therefore in the air-sea exchanges. This result is robust in both ocean models and therefore independent of a model's physical parameterizations.
The sub-Antarctic zone (SAZ, 35-50 • S; region between SAF and STF), is an area of mode water formation and the  most vigorous atmospheric sink of CO 2 (Metzl et al., 1999;Sabine et al., 2004;Fletcher et al., 2006;Takahashi et al., 2009). The SAZ is also the location of enhanced phytoplankton blooms as opposed to the HNLC regions in the rest of the Southern Ocean (Lourantou and Metzl, 2011). The Southern Ocean, southward of 30 • S, is a major uptake region for atmospheric CO 2 , in atmospheric and ocean inversion models (Friedlingstein and co authors, 2006;Gruber et al., 2009) as well as in observational studies (Metzl et al., 1999;Takahashi and co authors, 2002;Takahashi et al., 2009) although the magnitude of this uptake is argued (Gruber et al., 2009;Denman and Brasseur, 2007) .
Mixed layer depths in this region where mode waters form, e.g., the Sub-Antarctic Mode Water (SAMW) and the Antarctic Intermediate Water (AAIW), increase particularly along the SAF, but the uncertainty in the observations is large since this region is severely undersampled, especially during austral winter. Our models show that increasing remineralization rate leads to modest decreases in pCO 2 of about 6-10 %, and the Southern Ocean CO 2 flux sensitivity to the biological pump parameterizations is greater than model uncertainty due to physical processes.
As in previous studies (Kriest and Oschlies, 2011;Marinov et al., 2008), we conclude that in the Southern Ocean, the efficiency of the biological pump determines the atmospheric pCO 2 distributions. Surface nutrients are not being efficiently utilized in the Southern Ocean and this leads to deep water masses (which have been ventilated in the Southern Ocean) that have more preformed nutrients (i.e., nutrients that have not been biologically utilized, originate from the surface and are redistributed through mixing, "unremineralized nutrients") than NADW. As the deep water masses, rich in preformed nutrients and ventilated in the Southern Ocean, upwell at low latitudes, they promote increased bio-utilization of surface nutrients and more carbon fixation at the surface. However, increasing remineralization depth (or decreasing remineralization rate) depletes surface nutrients in the NADW more than in the Southern Ocean, because in the latter the nutricline is less steep. Consequently, deep water masses ventilated from the Southern Ocean are not as efficient as NADW at storing atmospheric CO 2 when the remineralization rate decreases, a result in agreement with Kwon et al. (2009).
Carbon export efficiency is determined by the fraction of carbon fixed by primary production that is transported through recycling and sinking to the deep ocean. This carbon flux is important because it removes carbon that is in contact with the atmosphere and sequesters it to greater depths and therefore increases the ocean's ability to uptake more CO 2 . The global mean export ratio is estimated to be between 10 and 40 % (Henson et al., 2012) and both GISS models lie within this range.
The globally integrated carbon export is found to be about 5 ± 2.2 Pg C yr −1 (Henson et al., 2012), 6 Pg C yr −1 (Martin et al., 1987) or even 10 Pg C yr −1 (Dunne et al., 2007). We find carbon export ranging from 1 to 7 Pg C yr −1 according to remineralization rates and model physics. GISSEH has typically less carbon export across the 75 m depth and as we increase remineralization rate, the carbon export decreases. The carbon export efficiency (i.e., the ratio of carbon export to primary productivity, Sarmiento et al. (2004)) is about 0.5 for the low remineralization rate in GISSER and drops to 0.17 for the higher rates we simulated. GISSEH is generally half as efficient in exporting carbon from the euphotic zone. Our remineralization profiles agree with Martin et al. (1987), who found that about half of the carbon exported from the euphotic zone towards greater depths is regenerated at depths shallower than 300 m; whereas two-thirds is regenerated in the upper 500 m and 90 % in the upper 1500 m. This very fast recycling in the euphotic zone implies that the majority of the carbon exported from the surface will be again in contact with the atmosphere within tens of years or less (Martin et al., 1987).
In addition to carbon, nutrient distributions are also sensitive to changes in remineralization rate, particularly in the SAZ (Pollard et al., 2002;Lovenduski and Gruber, 2005), the formation region of the SAMW, which supplies most of the world ocean surface water with nutrients (Sarmiento et al., 2003). This area has been shown to be very sensitive to changes in the Atlantic Meridional Overturning Circulation (AMOC; Schmittner (2005)) and in our results is shown to be very sensitive to biogeochemical changes as well. The fact that HYCOM has shallower penetration of ventilated NADW than the Russell model affects remineralization levels, acting to reduce the remineralization depth.
Finally, although the scope of this paper is not to evaluate the quality of the model predictions against observations (since the model reflects preindustrial equilibrium) a brief comparison to observed fields is given in Table 1. In particular, regional correlations of the annual mean surface distributions in total chlorophyll, diatoms, chlorophytes, cyanobacteria, coccolithophores as well as nutrient fields (nitrate, silicate and iron) and DIC are shown in Table 1. GISSER yields equal or better regional correlations with observations in the low remineralization case for all these fields whereas GISSEH distributions seem to be much less affected by the change in the remineralization rate. While this is not a rigorous model evaluation, it is evident that the biological sensitivity of the two models is different with the GISSER model being more transparent to these changes. The observational data for the estimation of correlation coefficients come from the data set described in Gregg and Casey (2007), Sect. 2.3. This data set includes 469 surface distributions of phytoplankton abundances and is available at http://polar.gsfc.nasa.gov/research/oceanbiology/index.php.

Conclusions
Simulations of the natural carbon cycle using ocean models that parameterize differently physical processes such as convection, diffusion and advection of thermodynamic properties and tracers, can be very instructive. They help set the range of uncertainty in the oceanic CO 2 uptake, nutrient and carbon distributions, both in global but also in regional terms. Varying a biogeochemical parameterization of the ocean biological pump, such as the recycling rate, can help us assess how important biogeochemical processes are in controlling the air-sea flux of CO 2 , the regional distributions of CO 2 sources and sinks as well as the carbon export efficiency, compared to physical model uncertainties.
We found that uncertainties due to biology become as important and occasionally more important than uncertainties in physics and this may suggest that improvements to physical models and biogeochemical models should be carried out in parallel. In our study, the Southern Ocean emerges as particularly sensitive to changes in the recycling rates and less so to physical model differences. The magnitude of the global biological pump depends on the effectiveness of mixing in the Southern Ocean but its sensitivity to recycling is remarkably robust, i.e., independent of physical parameterization choice. Similar analyses to more ocean/biogeochemical models could further generalize these results.