Global Biogeochemical Cycles

DINAUER ET AL.

10.1029/2021GB007131

3 of 34 depth z, F(z), using the power-law "Martin curve," i.e., F(z) = F 100 (z/100) -b , where F 100 is the POC export flux at 100 m and b is the attenuation coefficient (Martin et al., 1987). This approach has some important limitations. It has become clear that the value of b cannot be uniformly applied (Henson et al., 2012;Marsay et al., 2015). It is also clear that the POC flux attenuation predicted by a power law does not hold for oxygen minimum zones (Weber & Bianchi, 2020). Importantly, ESMs that employ an empirical flux-vs.-depth relationship are limited in their ability to simulate the past or future response of flux attenuation to changing climate/ocean conditions.

Marine ecosystem models have seen a number of improvements in the representation of the biological pump, with the caveat that they tend to focus on surface-ocean food web processes rather than the details of particle sinking and transformations in the mesopelagic zone. Many ecosystem models now include some formulation of temperature and/or oxygen effects on microbial remineralization along with implicit or explicit parameterizations of zooplankton grazing (e.g., Aumont et al., 2015;DeVries & Weber, 2017;Keller et al., 2012;Quéré et al., 2016;Yool et al., 2013). However, particle settling velocities are either prescribed as a fixed value or calculated as a function of depth, while fragmentation processes are either parameterized in an implicit fashion or neglected altogether. More complex 1D models have been developed to provide insight into the mesopelagic food web structure and associated particle transformations (e.g., Anderson & Tang, 2010;Boyd & Stevens, 2002;Dadou et al., 2001;Giering et al., 2014;Jackson & Burd, 2002;Stemmann et al., 2004), yet these models are often designed and evaluated for a single site. In addition, most particle transformation models have focused on the role of physical aggregation in determining the particle size spectrum (Burd, 2013;Burd andJackson, 2002, 2009;Jackson, 1995;Jackson & Burd, 2015). A mechanistic marine aggregate sinking scheme has recently been developed and deployed in an ESM (Maerz et al., 2020), but a thorough examination of parameter uncertainties remains difficult owing to higher computational costs inherent in higher model complexity. The added complexity of aggregate dynamics can also introduce new sources of model errors and uncertainties by increasing the number of potentially poorly known input parameters. There is still a need for computationally efficient, dynamic marine particle cycling models that (a) are suitable for implementation in the ESMs used in climate research and (b) provide a viable tool for identifying the main drivers of global T eff variability.

The broad objectives of this paper are to (a) provide a modeling framework that enables the simulation of T eff on a global scale, (b) explore the range of model behavior in response to variations in the input parameters believed to be especially important for controlling T eff , and (c) determine the relative importance of various environmental factors in the spatial pattern of T eff . To achieve these goals, a new, dynamic 1D marine particle cycling model has been developed and implemented on a latitude-longitude-depth grid. The model performs well against a compilation of 27 field-derived sinking POC flux and oxygen consumption rate profiles. Our modeling approach has a number of advantages. First, POC fluxes are calculated interactively with reasonably sophisticated treatments of particle sinking and transformation processes, thereby requiring a relatively small number of input parameters. Second, the model's ocean state variables, including temperature, dissolved oxygen, density, viscosity, and calcite and aragonite saturation states, are prescribed using observation-based climatologies, which avoids some of the uncertainties and biases inherent in standard ocean circulation models (e.g., model deficiencies in the volume of oxygen minimum zones). Third, the computational efficiency of the model allows for an extensive assessment of parameter uncertainties.

Methods

Global Ocean Particle Cycling Model

The Model of Sinking Particles and Cycling of Marine Aggregated Matter (MSPACMAM) is a computationally efficient 1D marine particle cycling model that is well suited for implementation in ESMs. The model includes two representative size classes of biogenic particles that are sufficiently large and/or dense to gravitationally settle through the water column. The small particle size class represents nano-and micro-sized sinking particles, whereas particles in the large size class may be marine snow, fecal pellets, and/or phytodetritus. The model implicitly assumes that pico-sized particles are either remineralized or repackaged into larger (≥nano-sized) particles in the surface ocean. Nonsinking (suspended) particles are not resolved. We assume that the sinking particles consist of four explicit components: particulate organic matter (POM), two forms of CaCO 3 (calcite and aragonite), and opal (SiO 2 ). This assumption is important since the various component particles making up aggregated matter each have their own densities and thereby play a role in determining particle settling velocities. The model explicitly represents the respiratory decomposition of organic matter (POC), mineral (calcite, aragonite, DINAUER ET AL.

10.1029/2021GB007131

4 of 34 and opal) dissolution, and fragmentation of large particles into small ones. The impacts of ocean temperature, dissolved oxygen, and calcite and aragonite saturation states on the loss rates are resolved.

The sinking POC flux associated with each size class is equal to the POC concentration, [POC], multiplied by the particle sinking speed, w sink . The bulk sinking POC flux, F POC , on a grid cell face is diagnosed by

where subscripts S and L refer to the small and large particle size classes, respectively.

Model Equations for a Grid Cell

The model carries four different particle tracers in small (S) and large (L) size classes: POC (POC S , POC L ), calcite ( 𝐴𝐴 calc CaCO 𝑆𝑆 3 ; calc CaCO 𝐿𝐿 3 ), aragonite ( 𝐴𝐴 arag CaCO 𝑆𝑆 3 ; arag CaCO 𝐿𝐿 3 ), and opal ( 𝐴𝐴 SiO 𝑆𝑆 2 ; SiO 𝐿𝐿 2 ). The time evolution of the concentration (in mol m -3 ) of particle tracer C in size class X, C X , in a grid cell within the ocean model domain (defined here as 100 m to the ocean bottom) is given by

The first term on the right-hand side of the equation is the advective flux divergence where 𝐴𝐴 𝐴𝐴 𝑋𝑋 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 is the settling velocity of the small or large particle; the second term is the loss term due to remineralization or dissolution where k C is the particle tracer C loss rate in units of time -1 ; and the third term is the source/sink term due to large particle fragmentation where 𝐴𝐴 𝐴𝐴 𝑆𝑆 𝑓𝑓 𝑓𝑓𝑓𝑓𝑓𝑓 = -𝐴𝐴 𝐿𝐿 𝑓𝑓 𝑓𝑓𝑓𝑓𝑓𝑓 is the fragmentation rate constant.

Surface Boundary Fluxes

The prescribed boundary fluxes used to force the model include POC, calcite, aragonite, and opal export fluxes applied to the surface boundary at 100 m. We assume opal and aragonite export each originate from large plankton (diatoms and pteropods, respectively) and thus enter the large particle size class, whereas the POC and calcite export fluxes are partitioned into small and large particles. For simplicity, we assume the fraction of POC or calcite export that enters the large particle size class, σ L , is proportional to the fractional contribution of opal and aragonite to total particle export:

where κ L is the constant of proportionality, and POC e , calc e , arag e , and opal e are the prescribed POC, calcite, aragonite, and opal export fluxes, respectively (see Appendix A2 for more details of these fluxes). The proportionality constant κ L used in this study is one. The remaining fraction of POC or calcite export, 𝐴𝐴 (1 -𝜎𝜎𝐿𝐿) , is routed to the small particle size class. The scientific rationale for this export partitioning scheme is discussed further in Section 3.4.

Particle Settling Velocities

The sinking speeds of small and large particles, 𝐴𝐴 𝐴𝐴 𝑆𝑆 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 and 𝐴𝐴 𝐴𝐴 𝐿𝐿 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 , are diagnosed separately at each time step using Rubey's (1933) equation, which, unlike Stokes' (1851) law, applies in both the laminar (viscous) and turbulent (inertial) regimes. Rubey's equation equates the effective weight of a spherical particle to two different forces: (a) the viscous resistance of the fluid through which it moves and (b) the momentum of the fluid that strikes the particle. For an arbitrary size class, the particle settling velocity 𝐴𝐴 𝐴𝐴 𝑋𝑋 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 (in cm s -1 ) at a tracer grid cell face is given by

where g = 981 cm s -2 is the acceleration due to gravity, ρ sw is the in situ density of seawater (in g cm -3 ), 𝐴𝐴 𝐴𝐴 𝑋𝑋 ptcl is the bulk density of the small or large particle (in g cm -3 ), r X is the prescribed effective radius of the small or large particle (in cm) (see Tables 1 and 2), and μ is the dynamic viscosity of seawater (in g cm -1 s -1 ). DINAUER ET AL.

10.1029/2021GB007131

5 of 34 For an arbitrary size class, the bulk density of the particle, 𝐴𝐴 𝐴𝐴 𝑋𝑋 ptcl , is equal to the density of the solids (POM, CaCO 3 , and opal) plus any porewater (seawater):

where ϕ X is the prescribed porosity of the small or large particle (see Tables 1 and 2) and 𝐴𝐴 𝐴𝐴 𝑋𝑋 solid is the total density of the various debris making up the small or large particle (i.e., solid particle density).

For an arbitrary size class, the solid particle density, 𝐴𝐴 𝐴𝐴 𝑋𝑋 solid , is given by the total mass of particulate matter per unit volume of seawater divided by the total volume of particulate matter per unit volume of seawater:

where

] + are the mass concentrations (in g m -3 ) of POM, CaCO 3 (calcite + aragonite), and opal in the small or large particle size class, and ρ POM , 𝐴𝐴 𝐴𝐴CaCO 3 , and 𝐴𝐴 𝐴𝐴SiO 2 are the prescribed densities of POM, CaCO 3 , and opal in units of g m -3 (see Table 1). The masses of POM, CaCO 3 , and opal per unit volume of seawater are derived from the tracer concentrations POC

X , 𝐴𝐴 calc CaCO 𝑋𝑋 3 , 𝐴𝐴 arag CaCO 𝑋𝑋 3 , and 𝐴𝐴 SiO 𝑋𝑋 2

MW𝐶𝐶 𝛼𝛼𝛼 (7)

Parameter Symbol Value Units References Gravitational acceleration g 9.81 m s -2 Radius of large particle r L 250 μm Alldredge and Silver (1988) Porosity of small particle ϕ S 0 Proportionality constant for export partitioning a κ L 1 Density of POM ρ POM 1.06 g cm -3 Density of CaCO 3 𝐴𝐴 𝐴𝐴CaCO 3 2.71 g cm -3 Density of opal 𝐴𝐴 𝐴𝐴SiO 2 2.10 g cm -3 Molecular weight of carbon MW C 12 g/mol Molecular weight of CaCO 3 𝐴𝐴 MWCaCO 3 100 g/mol Molecular weight of opal 𝐴𝐴 MWSiO 2 60 g/mol Mass ratio of organic matter to carbon α 2.7 Sarmiento and Gruber (2006) Calcite dissolution rate constant b k calc For 0.8 < Ω calc < 1, 4.33 × 10 -6 ; for Ω calc ≤ 0.8, 0.0137 mol m -2 d -1 Naviaux et al. (2019) Calcite dissolution reaction order n calc For 0.8 < Ω calc < 1, 0.11; for Ω calc ≤ 0.8, 4.7 Naviaux et al. (2019) Aragonite dissolution rate constant c k arag For Ω arag < 1, 5.28 × 10 -5 mol m -2 d -1 Dong et al. (2019) Aragonite dissolution reaction order n arag For Ω arag < 1, 1.37 Dong et al. (2019) Specific surface area of small particle calcite 𝐴𝐴 SSAcalc 𝑆𝑆 1040 m 2 mol -1 Subhas et al. (2018) Specific surface area of large particle calcite𝐴𝐴 SSAcalc 𝐿𝐿 430 m 2 mol -1 Subhas et al. (2018) Specific surface area of aragonite SSA arag 217 m 2 mol -1 Honjo and Erez (1978) a Constant of proportionality between the fraction of POC/calcite export entering the large particle size class, σ L , and the relative abundance of opal and aragonite export (see Equation 3). b For Ω calc ≥ 1, k calc = 0. c For Ω arag ≥ 1, k arag = 0. where MW C , 𝐴𝐴 MWCaCO 3 , and 𝐴𝐴 MWSiO 2 are the molecular weights of carbon, CaCO 3 , and opal in units of g/mol and α is the mass ratio of organic matter to carbon, i.e., the POM:POC mass ratio (see Table 1). For simplicity, we neglect the effect of increasing pressure with depth on bulk particle densities.

Microbial Remineralization of POC

The respiratory decomposition of organic matter is parameterized in the model using first-order kinetics. The POC remineralization rate constant, 𝐴𝐴 𝐴𝐴 ′ POC , depends on ocean temperature T and dissolved oxygen O 2 :

where k POC is the base remineralization rate constant (in time -1 ) referenced to temperature

is the modification of k POC by temperature, and

is the oxygen limitation of k POC . The exponent b T (in °C -1 ) in the temperature modification function is trivially related to the remineralization rate increase per 10°C increase, i.e., Q 10 , by

Because they are easier to conceptualize, Q 10 values will be used in this paper instead of b T values. Oxygen limitation is parameterized using Michaelis-Menten kinetics, where

𝐴𝐴 𝐴𝐴𝑂𝑂 2 is the half-saturation constant for oxygen uptake. The parameter values chosen for k POC , Q 10 , and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 are shown in Table 2. Parameter Observed values Mean value Description References r S (μm) 5, 7.5, 10 15.625 Nanoplankton a Durkin et al. (2015), Richardson (2019) 12.5, 15, 20, 25, 30 Small particles b ϕ L (%) 95, 96, 97, 98, 99 97 Aggregates Alldredge and Gotschalk (1988), Ploug, Iversen, and Fischer (2008) k POC (d -1 ) 0.02 (8°C), 0.038 (21°C), 0.066 (21°C), 0.078 (20°C), 0.144 (20°C), 0.2 (18°C) 0.122 (16.47°C) POM c Grill and Richards (1964); Garber (1984), Westrich and Berner (1984), Seiki et al. (1991) 0.03 (4°C), 0.08 (15°C), 0.089 (16°C), 0.11 (17°C), 0.12 (15°C), 0.13 (15°C), 0.13 (20°C), 0.15 (15°C), 0.21 (15°C) Aggregates and fecal pellets Ploug et al. (1999), Ploug and Grossart (2000), Grossart and Ploug (2001), Ploug, Iversen, Koski, and Buitenhuis (2008), Iversen et al. (2010), Iversen and Ploug (2010), Iversen and Ploug (2013) 0.16 (20°C), 0.32 (20°C) Optimized with dataconstrained model DeVries and Weber (2017) Q 10 (unitless) 1.8, 1.88, 1.94, 3.3, 3.63, 3.9 2.29 Measured, or experimentally determined Eppley (1972), White et al. (1991), Seiki et al. (1991); Thamdrup and Fleischer (1998), Iversen and Ploug (2013), Brewer and Peltzer (2016) 1.51, 1.74, 2.01, 2.1, 2.4, 2.5 Determined empirically with data-constrained model DeVries and Weber (2017), Laufkotter et al. (2017), Weber and Bianchi (2020) 𝐴𝐴 𝐴𝐴O 2 (μmol L -1 ) 1, 2, 3, 4 8.25 Experimentally determined Devol (1978) 8, 12, 16, 20 Determined empirically with data-constrained model DeVries and Weber (2017), Laufkötter et al. (2017) Note. Parameters are allowed to take on "realistic" values, which are defined as observed values in the literature. The means of the literature-based parameter values (third column) are used to define the "base-case" point in the parameter space for the purpose of the sensitivity analysis. a Operationally defined by 20 μm diameter filter. b Operationally defined by 51 μm diameter filter. c Reactive marine algal POM, or labile marine planktonic POM.

Dissolution of CaCO 3 Minerals

The rate of CaCO 3 dissolution, 𝐴𝐴 𝐴𝐴CaCO 3 , is given by the empirically derived expression:

where k is the dissolution rate constant in units of time -1 , Ω is the saturation state of seawater with respect to the CaCO 3 mineral of interest, and n is the reaction order. In this expression, (1 -Ω) is the fractional extent of undersaturation where Ω represents the thermodynamic driving force for dissolution of calcite (Ω calc ) or aragonite (Ω arag ). k and n are empirically determined variables. k and n values from Naviaux et al. (2019) and Dong et al. (2019) are used for calcite and aragonite dissolution, respectively (see Table 1). Their rate constants (in mol m -2 time -1 ) are converted to units of time -1 by multiplying by the specific surface area SSA of CaCO 3 particles.

The SSA of small and large calcite component particles is set to that of coccoliths, 𝐴𝐴 𝐴𝐴𝐴𝐴𝐴𝐴calc 𝐴𝐴 = 1,040 m 2 mol -1 , and planktonic foraminifera, 𝐴𝐴 𝐴𝐴𝐴𝐴𝐴𝐴calc 𝐿𝐿 = 430 m 2 mol -1 , respectively (Subhas et al., 2018), whereas the SSA of aragonite is set to that of aragonitic pteropods, SSA arag = 217 m 2 mol -1 (Honjo & Erez, 1978).

Dissolution of Opal

The temperature sensitivity of opal solubility suggests that the rate at which opal dissolves is faster in warm waters than in cold waters. The model incorporates the effect of temperature on opal dissolution using the Arrhenius-type equation from Gnanadesikan (1999):

where k opal is the first-order opal dissolution rate constant in units of d -1 and T K is the absolute temperature (in Kelvin) of seawater. This equation results in k opal values of 0.25 d -1 (equivalent to 1/(4 d)) at 25°C and 0.01 d -1 (equivalent to 1/(100 d)) at 2°C.

Large Particle Fragmentation

The fragmentation of large particles is parameterized in the model as a first-order reaction with rate constant 𝐴𝐴 𝐴𝐴 𝐿𝐿 frag . For simplicity, we assume 𝐴𝐴 𝐴𝐴 𝐿𝐿 f rag (in d -1 ) decreases exponentially with depth following Briggs et al. (2020):

where z is depth (in m). This empirically derived function was obtained by simultaneously tracking changes in the concentrations of large and small particles over depth and time. It is assumed that 𝐴𝐴 𝐴𝐴 𝐿𝐿 f rag encompasses all physical (e.g., turbulence-generated shear stress) and biological (e.g., zooplankton "sloppy feeding" and swimming behavior) processes that fragment large particles into smaller ones.

Model Parameter Values

The various choices of parameter values are summarized in Tables 1 and 2. In this study, r S , ϕ L , k POC , Q 10 , and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 are treated as adjustable parameters (see Section 2.2.1). All other parameters are held constant in the model simulations. r L is assumed to correspond to the minimum radius of macroscopic aggregates, operationally defined as 250 μm (Alldredge & Silver, 1988), whereas ϕ S is assumed to be zero. Laboratory studies show that porosity decreases significantly with decreasing particle size (e.g., Alldredge & Gotschalk, 1988;Logan & Wilkinson, 1990). The treatment of r L and ϕ S as fixed parameters is motivated primarily by the need to reduce the number of model runs required for the full factorial design. These two parameters could be revisited in future work, although new observational constraints are needed to limit the model's dynamical freedom.

Grid Resolution and Model Configuration

The model equations are implemented on a 1° latitude-longitude grid with 26 levels in the vertical. The horizontal grid is the standard 1° grid used by the World Ocean Atlas (WOA) (Levitus, 1982;Locarnini et al., 2013). The vertical grid is based on the WOA's 33 standard depth levels. Here, the depth coordinates on the vertical tracer grid are from 112.5 to 5,250 m with grid cells ranging in thickness from 25 to 500 m. The ocean state variables (temperature, dissolved oxygen, density, viscosity, and calcite and aragonite saturation states) are prescribed using observation-based climatologies (see Appendix A1 for more details). The advective flux DINAUER ET AL. 10.1029/2021GB007131 8 of 34 divergence term, 𝐴𝐴 -𝜕𝜕 ( 𝑤𝑤 𝑋𝑋 sink [𝐶𝐶 𝑋𝑋 ] ) 𝜕𝜕𝜕𝜕

, is determined from a first-order "upwind" tracer advection scheme. The model is run with a prescribed surface forcing until it reaches equilibrium. The surface boundary fluxes used to force the model are POC, calcite, aragonite, and opal export fluxes taken from the Coupled Model Intercomparison Project Phases 5 and 6 (CMIP5 and CMIP6). A multi-model ensemble median is constructed for each of the forcing fields from two CMIP5 and five CMIP6 models (see Appendix A2 for further details). The time step used in this study is 5,400 s (1.5 hr), which was chosen so as not to deplete the smallest vertical grid cell over the time step given the estimated maximum probable particle settling velocity. The model is integrated for at least two sinking small particle residence times (≥1,080 days), during which the model reaches equilibrium. The wall-clock time required to run the model, per horizontal grid point and on a single CPU core, is roughly 12 s.

This study makes use of two different configurations of the MSPACMAM. The default model configuration, also referred to as the std configuration, does not include large particle fragmentation, whereas the impacts of fragmentation are explicitly resolved in a second model configuration (so-called frag configuration). Note that, although model runs are performed using both the std and frag configurations, discussion of the model results will focus mainly on the std configuration, given that the generality of the empirical specific fragmentation rate function of Briggs et al. (2020) is uncertain. Note also that the MSPACMAM can be configured to represent a physically protected fraction of POC that is inaccessible to remineralization due to the association with calcite minerals. However, the present study does not make use of this configuration as there is insufficient evidence to suggest that ballast minerals provide physical protection for labile POC against bacterial activity (Iversen & Ploug, 2010;Ploug, Iversen, & Fischer, 2008, Ploug, Iversen, Koski, & Buitenhuis, 2008).

Sensitivity Analysis and Model Evaluation

We explore the sensitivity of model output to five input parameters (r S , ϕ L , k POC , Q 10 , and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 ) thought to be important for POC flux attenuation with depth. A sensitivity analysis is conducted for three main reasons: (a) to assess how the POC flux profile shape responds to changes in each parameter, (b) to determine which parameters contribute most to the variability (or uncertainty) in T eff , and (c) to find combinations of parameter values yielding model results that agree well with available observations. The identification of the most sensitive parameters in the model serves (a) to provide direction for future research so that the uncertainties of the influential parameters are reduced and (b) to guide future calibration efforts such that computational resources are allocated to constraining important parameters.

Full Factorial Experiment

Model simulations are performed with different choices of parameter values. Each parameter is allowed to take on observed values selected from the literature (see Table 2). A full factorial design is used to explore the 5D parameter space, which required 9 × 5 × 18 × 13 × 9 = 94,770 runs of the model for each of the configurations (std and frag). The model is run only at the 26 horizontal grid points where evaluation data are available (see Section 2.2.3.1 and Figure 1).

Differential Sensitivity Analysis

A widely used approach to sensitivity analysis is differential analysis, which is based on the calculation of first-order partial derivatives of the output variable of interest with respect to the input parameters at a "base-case" point in the parameter space (Hamby, 1994;Iman & Helton, 1988;Pianosi et al., 2016). The partial derivatives are used as sensitivity indices (or importance measures) to discriminate between influential and noninfluential parameters. The purpose of the differential analysis applied here is to identify the input parameters that exert the most influence on the model's mean T eff response, 𝐴𝐴 𝑇𝑇eff . Depending on the analysis, T eff is averaged either over the 26 data locations or over a subset of data locations belonging to a distinct ocean regime. 𝐴𝐴 𝑇𝑇eff is recorded for each model run (parameter set) in the factorial design prior to performing the differential analysis.

The base-case point is defined as the mean values of the input parameters p 1 ,…,p n , 𝐴𝐴 𝑝𝑝 . The partial derivative 𝐴𝐴 𝜕𝜕𝑇𝑇 ef f 𝜕𝜕𝜕𝜕 𝑖𝑖 evaluated at 𝐴𝐴 𝑝𝑝 can be thought of as the local sensitivity of 𝐴𝐴 𝑇𝑇eff to input parameter p i . Here, the partial derivatives are numerically approximated as finite differences. To remove the effect of units, each partial derivative is normalized as 𝐴𝐴 ( 𝜕𝜕𝑇𝑇eff ∕𝜕𝜕𝜕𝜕𝑖𝑖 ) ( 𝜕𝜕𝑖𝑖0∕𝑇𝑇eff 0 ) where p i0 is the base-case (mean) value of p i and 𝐴𝐴 𝑇𝑇eff 0 is the 𝐴𝐴 𝑇𝑇eff predicted DINAUER ET AL. 10.1029/2021GB007131 9 of 34 by base-case values 𝐴𝐴 𝑝𝑝 . These standardized partial derivatives, also referred to as normalized sensitivity coefficients S i , are used for parameter ranking. The variance of 𝐴𝐴 𝑇𝑇eff due to parameter uncertainties, 𝐴𝐴 𝐴𝐴 𝐴𝐴𝐴𝐴 ( 𝑇𝑇eff ) , is estimated using the general error propagation formula. 𝐴𝐴 𝐴𝐴 𝐴𝐴𝐴𝐴 ( 𝑇𝑇eff )

is treated as a function of the input parameters and approximated with a first-order Taylor series expansion at 𝐴𝐴 𝑝𝑝 :

where 𝐴𝐴 Var (𝑝𝑝𝑖𝑖) is the variance of the literature-based values of p i . Here, the input parameters are assumed to be independent of one another. While it is only an approximation, the advantage of using the Taylor series expansion in Equation 15 is that it quantifies the relative importance of each parameter's sensitivity versus variability as well as its percentage contribution to the estimated variance of 𝐴𝐴 𝑇𝑇eff , i.e., 𝐴𝐴

Evaluation of Model Performance

The 94,770 model runs, each with a different parameter set, are evaluated against available observations using a cost function that summarizes the model-data misfit. The combinations of parameter values yielding the lowest cost are referred to as the best performing parameter sets. DINAUER ET AL. 10.1029/2021GB007131 10 of 34 2.2.3.1. Observational Constraints Observations of the sinking POC flux and oxygen consumption rate are used to evaluate model performance at 26 locations across the global ocean (see Figure 1). We have compiled a collection of 22 POC flux profiles (defined here as >2 flux observations in the mesopelagic zone) measured either with free-drifting sediment traps (i.e., surface tethered and neutrally buoyant sediment traps) or with Marine Snow Catchers. The compilation includes data from Martin et al. (1987), Devol and Hartnett (2001), Van Mooy et al. (2002), Buesseler et al. (2007), Owens et al. (2013), Durkin et al. (2015), Marsay et al. (2015), Keil et al. (2016), Engel et al. (2017), and Cavan et al. ( 2017). We have also incorporated results from the recent EXPORTS program, which are available online at https://seabass.gsfc.nasa.gov/. Five observation-based oxygen consumption rate profiles (Sonnerup et al., 2013(Sonnerup et al., , 2015;;Stanley et al., 2012) are selected from the compilation of Brewer and Peltzer (2017). In order to make use of these observations, model-predicted POC remineralization rates are combined with DOC remineralization rates from the data-assimilated biological pump model of DeVries and Weber ( 2017) and multiplied by a O 2 :C stoichiometric ratio of 170:117 (Anderson & Sarmiento, 1994).

Model Performance Metrics

A cost function J is used to evaluate model performance with each of the 94,770 parameter sets. J is expressed as a weighted sum of squared differences between model predictions (P) and observations (O):

where M = 2 is the number of different data types or variables (i.e., sinking POC flux F POC and oxygen consumption rate OUR), N v is the number of model grid cells where observations are available for variable v, and P vi and O vi are the predicted and observed values of variable v at grid cell i. To account for differing units, the weight given to each variable is inversely proportional to the square of the standard deviation of its observations, σ v (Friedrichs et al., 2007;Kuhn et al., 2018;Kuhn & Fennel, 2019;Xiao & Friedrichs, 2014). The weighting factor, C v , is included to increase the weight of F POC since it is the output variable of most interest. Here, 𝐴𝐴 𝐴𝐴𝐹𝐹 𝑃𝑃 𝑃𝑃𝐴𝐴 is set to 3 and C OUR to 1. The POC flux profiles were normalized to the flux at the shallowest depth at which observations are available, the oxygen consumption rate data were log transformed, and all observations were binned to the model grid prior to calculating J.

Best performing parameter sets are defined as all parameter sets yielding cost values within 1% of the minimum cost value. This threshold value of 1% is used to reduce the number of "acceptable solutions" to a manageable size. In addition to cost values, model performance is quantified with four univariate statistical measures (Stow et al., 2009): (a) the standard deviation of the predictions relative to the observed standard deviation, i.e., σ mod /σ obs , (b) the correlation coefficient, r, (c) the root-mean-square error (RMSE), and (d) the normalized RMSE, 𝐴𝐴 NRMSE = RMSE∕𝜎𝜎𝑜𝑜𝑜𝑜𝑜𝑜 .

Analysis of Global Transfer Efficiency Variability

Model Experiments

Global model simulations are performed to provide insight into six environmental factors EF 1 ,…,EF n thought to be important for controlling the spatial pattern of T eff , i.e., ocean temperature, dissolved oxygen, seawater density, seawater viscosity, particle size, and particle density. These factors will be denoted by T, O 2 , ρ sw , visc, σ L , and ρ PM , respectively. Table 3 summarizes the suite of model experiments utilized for the analysis of global T eff variability. The baseline simulation (ctrl) provides a benchmark for the spatial pattern of T eff . Idealized experiments have been designed to remove spatial variations in EF 1 ,…,EF n one at a time. In the case of T, O 2 , ρ sw , and visc, the model is run with 3D fields set to global mean mesopelagic values. Note that the temperature adjustment applies only to the temperature field used to calculate the rate of POC remineralization (see Equation 10). Particle size variability is removed by setting the σ L distribution to the global mean of 0.154 from the ctrl simulation. To remove particle density variations, the bulk particle densities of the small and large size classes, 𝐴𝐴 𝐴𝐴 𝑆𝑆 ptcl and 𝐴𝐴 𝐴𝐴 𝐿𝐿 ptcl , are prescribed rather than calculated interactively (Equation 5).

𝐴𝐴 𝐴𝐴 𝑆𝑆

ptcl is set to the global mesopelagic mean of 1.536 g cm -3 from the ctrl simulation, and 𝐴𝐴 𝐴𝐴 𝐿𝐿 ptcl is likewise set to 1.052 g cm -3 .

DINAUER ET AL.

10.1029/2021GB007131

11 of 34

Deconvolving the Environmental Factors

The contribution of environmental factor EF i to the spatial pattern of T eff is estimated by the difference between the patterns of T eff for ctrl and an experiment in which EF i is forced to be spatially uniform:

where

𝐴𝐴 ∆𝑇𝑇 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f is the component of the T eff pattern attributable to EF i variability, 𝐴𝐴 𝐴𝐴 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 ef f is the T eff pattern from ctrl, and 𝐴𝐴 𝐴𝐴 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f is the T eff pattern from an experiment in which the spatial variations in EF i have been removed (Table 3). It should be noted that there is less confidence in 𝐴𝐴 ∆𝑇𝑇 𝜎𝜎 𝐿𝐿 ef f and 𝐴𝐴 ∆𝑇𝑇 𝜌𝜌 𝑃𝑃 𝑃𝑃 ef f than in the other environmental components, since the particle size and density information used to drive the model equations is diagnosed by the model itself, rather than derived from observation-based climatologies as is the case for seawater temperature, density, viscosity, and dissolved oxygen. Nonetheless, the broad spatial features should be properly manifested in the 𝐴𝐴 ∆𝑇𝑇 𝜎𝜎 𝐿𝐿 ef f and 𝐴𝐴 ∆𝑇𝑇 𝜌𝜌 𝑃𝑃 𝑃𝑃 ef f patterns.

Results and Discussion

Model-Data Comparisons

The performance of the MSPACMAM is evaluated against observed sinking POC fluxes and oxygen consumption rates at 26 locations across the global ocean (see Figure 1). Cost function values are similar across a large portion of the parameter sets (Figure 2), indicating that different combinations of parameter values yield equally plausible model solutions. The shape of the vertical POC flux profile results from the balance between the rates of particle sinking and POC remineralization. Without sufficient observations to constrain one or both of these terms, there are many different ways of reproducing the same balance. As a consequence, there is not a single, unique set of parameter values that outperforms the others.

The ranges of the best performing parameter sets, defined here as parameter sets within 1% of the minimum cost value, are shown in Table 4. The best-fit parameter values span much of the literature-based ranges, with the exception of Q 10 and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 whose low-cost values lie at the lower end of their possible ranges. The parameter combinations within 5% of the minimum cost value span even wider ranges as shown in Table B1. The subset of best-fit parameter combinations with a slower rate of POC remineralization (i.e., k POC below the 50th percentile) exhibit a slower rate of particle sinking, which can be explained by smaller particles in the small size class

Experiment name Abbreviation Description Baseline or control ctrl Model allows all environmental factors to vary in space Constant temperature 𝐴𝐴 𝑇𝑇 Model assumes ocean temperature a is spatially uniform at the global mesopelagic mean of 7.26°C Constant oxygen 𝐴𝐴 𝑂𝑂2 Model assumes dissolved oxygen is spatially uniform at the global mesopelagic mean of 195 mmol m -3 Constant seawater density 𝐴𝐴 𝜌𝜌𝑠𝑠𝑠𝑠 Model assumes seawater density is spatially uniform at the global mesopelagic mean of 1029 kg m -3 Constant viscosity 𝐴𝐴 𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 Model assumes seawater viscosity is spatially uniform at the global mesopelagic mean of 0.0015 kg m -1 s -1 Constant particle size 𝐴𝐴 𝜎𝜎𝐿𝐿 Model assumes POC export allocation fraction σ L is spatially uniform at the global mean of 0.154 Constant particle density 𝐴𝐴 𝜌𝜌PM Model assumes small and large particle densities 𝐴𝐴 𝐴𝐴 𝑆𝑆 ptcl and 𝐴𝐴 𝐴𝐴 𝐿𝐿 ptcl are spatially uniform at their global mean mesopelagic values of 1,536 kg m -3 and 1,052 kg m -3 , respectively Note. The baseline model simulation (or control simulation) provides a benchmark for the spatial pattern of T eff . The idealized experiments target six key environmental factors controlling the large-scale T eff pattern: ocean temperature, dissolved oxygen, seawater density, seawater viscosity, particle size, and particle density (see Section 2.3.1 for more details). a

The temperature field used to calculate the rate of POC remineralization (see Equation 10).

Table 3 Overview of Model Experiments Utilized for the Analysis of Global T eff Variability DINAUER ET AL.

10.1029/2021GB007131 An open question raised by sediment trap studies is whether T eff is higher over high latitudes or higher over low latitudes. About one fourth of the literature-based parameter sets yield low-latitude T eff that is higher than high-latitude T eff , although none are within 1% of the minimum cost value. Higher low-latitude T eff is achieved by running the model with the combination of relatively large r S (22.6 μm), high ϕ L (97.4%), and high 𝐴𝐴 𝐴𝐴𝑂𝑂 2 (12.3 μmol L -1 ). An important caveat is that model evaluation is limited to those locations for which observations are available, with the consequence being many of the low-latitude locations are within oxygen minimum zones (see Figure 1) where the lack of oxygen strongly inhibits microbial remineralization. We conclude that the credible representation of higher low-latitude T eff outside oxygen minimum zones is difficult to achieve with literature-based parameter values.

Comparisons of modeled and observed normalized POC flux profiles are shown in Figure 3 for 12 of the locations where observations are available (see Figure B1 for all 22 locations). Model results are in good agreement with observations across a range of temperature and oxygen regimes (see Figure B2 and Figure B3 for the temperature and oxygen profiles at the data locations). Model performance metrics, including σ mod /σ obs (a measure of the amplitude of modeled variability relative to observed variability), r (a meas- r S (μm) ϕ L (%) 𝐴𝐴 𝐴𝐴 𝑟𝑟𝑟𝑟𝑟𝑟 𝑃𝑃 𝑃𝑃𝑃𝑃 (d -1 ) a Q 10 (unitless) 𝐴𝐴 𝐴𝐴O 2 (μmol L -1 ) A priori 5-30 95-99 0.033-0.32 1.51-3.9 1-20 Standard configuration (N = 73) Mean 17.0 96.2 0.11 1.51 2.01 Range 10-25 95-99 0.037-0.22 1.51-1.51 1-3 Configuration with fragmentation (N = 57) Mean 16.4 95.4 0.10 1.51 1.79 Range 10-25 95-97 0.033-0.22 1.51-1.51 1-3 Note. The choice of threshold for defining these "acceptable solutions" is a 1% change in the minimum cost value (see Equation 16 for the cost function). For comparison, the ranges of the parameter combinations below a threshold of 5% are shown in Appendix B (Table B1). a The POC remineralization rate constant, k POC , referenced to 20°C. Since the original, literature-based k POC values are referenced to various temperatures T ref (see Table 2), they have been normalized to a common reference temperature of 20°C by multiplying by 𝐴𝐴 𝐴𝐴 𝑏𝑏 𝑇𝑇 (20-𝑇𝑇ref ) where 𝐴𝐴 𝐴𝐴𝑇𝑇 = ln 𝑄𝑄10∕10 . For k POC in the best performing parameter sets, the accompanying Q 10 value is used. For the a priori literature-based values of k POC , a Q 10 value of 1.51 is used.

Table 4 Means and Ranges of the Best Performing Parameter Sets in the 94,770-Run Factorial Experiment DINAUER ET AL.

10.1029/2021GB007131 5. The normalized error for the POC flux profiles (NRMSE = 0.35) is smaller than that for the oxygen consumption rates (NRMSE = 0.93), which is expected since the latter is a reconstruction. There is a strong positive relationship between the predicted and observed POC flux values with a correlation coefficient of 0.94. For oxygen consumption rates, the correlation between predictions and observations is weaker (r = 0.64). The ratio of the modeled to the observed standard deviation is close to 1 for both output variables. The ability of the MSPACMAM to capture important aspects of the observed variability in POC flux profiles, including amplitude and pattern, adds confidence to its generality for global-scale simulations of T eff .

Model Sensitivity to Parameter Variations

Response of the Vertical POC Flux Profile

Figure 4 shows the mean POC flux profiles resulting from varying each input parameter over its literature-based range, while all others are held constant at their base-case values (see Table 2). The shape of the profiles Global Biogeochemical Cycles DINAUER ET AL. 10.1029/2021GB007131 14 of 34 reflects the strength of POC flux attenuation with depth. Larger values of r S and lower values of ϕ L increase bulk particle sinking speeds and therefore decrease POC flux attenuation. Higher k POC values increase POC remineralization rates, thereby increasing flux attenuation. Increases in Q 10 decrease remineralization rates in cold waters deeper in the water column, leading to decreased flux attenuation. Higher values of 𝐴𝐴 𝐴𝐴𝑂𝑂 2 result in weaker remineralization rates in oxygen deficient waters and thus lower POC flux attenuation.

The response of the POC flux profile to perturbed parameter values varies vertically. Changes in k POC have a large impact on the sinking POC flux throughout the water column, although the deep ocean response below ∼2,000 m is noticeably weaker. While the influence of r S on the upper ocean POC flux profile is similar to that of k POC , its impact is more limited in the deep ocean (>1,000 m) where small particles become less prevalent. In contrast to k POC and r S , the influence of ϕ L and Q 10 on the POC flux profile increases with increasing depth. Changes in Q 10 induce changes in the POC remineralization rate via

. Because the base-case value of k POC is referenced to the relatively warm temperature of 16.5°C (i.e., T ref = 16.5°C), the deeper the depth (the colder the ocean temperatures) the stronger the influence of Q 10 on the remineralization rate. In the case of ϕ L , large particles contribute more to the bulk sinking POC flux with depth, and hence changes in ϕ L have a stronger influence on the deep ocean POC fluxes. While influential at depth, ϕ L variations have the weakest influence on the POC flux in shallow waters above ∼500 m, where small particles tend to be generally more abundant than large ones. The sensitivity of the POC flux to changes in 𝐴𝐴 𝐴𝐴𝑂𝑂 2 appears to be very weak except in the shallow water column from ∼500 to 1,000 m, the depth range over which suboxic conditions occur at some of the 26 locations considered in this analysis (see Figure B3).

Sensitivity Analysis of Transfer Efficiency

Relative Influence of Model Parameters on Transfer Efficiency

To identify the most influential input parameters, we first examine the range of 𝐴𝐴 𝑇𝑇eff output, averaged over the 26 data locations, resulting from varying each parameter over its literature-based range, while all others are held constant. Figure 5 2 are referenced to various temperatures and therefore have been normalized to a common reference temperature of 20°C in order to place them on a comparable level for visualization purposes. (d) The remineralization rate increase per 10°C increase, Q 10 , is varied from 1.51 to 3.9. (e) The half-saturation constant for oxygen uptake, 𝐴𝐴 𝐴𝐴𝑂𝑂 2 , is varied from 1 to 20 μmol L -1 . Note that each parameter is varied over its literature-based range, while all others are held constant at their base-case values (i.e., their mean values).

DINAUER ET AL. . Based on the normalized sensitivity coefficients S i , the most sensitive parameters in order are ϕ L , r S , k POC , Q 10 , and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 (see Table 6). The fact that 𝐴𝐴 𝑇𝑇eff is least sensitive to Q 10 and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 might be expected since their impact on the 𝐴𝐴 𝑇𝑇eff output was among the smallest. However, the high sensitivity of 𝐴𝐴 𝑇𝑇eff to ϕ L is a rather surprising result given that changes in ϕ L resulted in the smallest range of 𝐴𝐴 𝑇𝑇eff output. This underlines the importance of examining not only the range of modeled 𝐴𝐴 𝑇𝑇eff behavior in response to p i variations but also the rate at which 𝐴𝐴 𝑇𝑇eff changes with respect to p i . Parameter pi 𝐴𝐴 𝜕𝜕𝑇𝑇 ef f 𝜕𝜕𝜕𝜕 𝑖𝑖 S i 𝐴𝐴 Var (𝑝𝑝𝑖𝑖) 𝐴𝐴 ( 𝜕𝜕𝑇𝑇 ef f 𝜕𝜕𝜕𝜕 𝑖𝑖 ) 2 Var (𝜕𝜕𝑖𝑖) % of 𝐴𝐴 Var ( 𝑇𝑇eff ) r S (μm) 0.95 0.86 66.8 60.4 52.51 ϕ L (%) -1.02 5.70 2.50 2.59 2.25 𝐴𝐴 𝐴𝐴 ref POC (d -1 ) a -78.7 0.74 0.007 40.9 35.55 Q 10 (unitless) 3.18 0.42 0.57 5.80 5.04 𝐴𝐴 𝐴𝐴O 2 (μmol L -1 ) 0.35 0.17 43.7 5.34 4.65 𝐴𝐴 Var ( 𝑇𝑇eff ) 115.0 SD of 𝐴𝐴 𝑇𝑇eff (%) 10.7 Note. The normalized sensitivity coefficients S i shown in the third column are the partial derivatives (second column) normalized by the quotient 𝐴𝐴 𝐴𝐴𝑖𝑖0∕𝑇𝑇eff 0 , where pi 0 is the base-case (mean) value of pi and 𝐴𝐴 𝑇𝑇eff 0 is the 𝐴𝐴 𝑇𝑇eff predicted by the basecase values 𝐴𝐴 𝑝𝑝 . The last column shows each parameter's percentage contribution to the estimated variance of 𝐴𝐴 𝑇𝑇eff , 𝐴𝐴 𝐴𝐴 𝐴𝐴𝐴𝐴 ( 𝑇𝑇eff ) (see Section 2.2.2 for more details). a The POC remineralization rate constant, k POC , referenced to 20°C. Since the original, literature-based k POC values are referenced to various temperatures T ref (see Table 2), they have been normalized to a common reference temperature of 20°C by multiplying by 𝐴𝐴 𝐴𝐴 𝑏𝑏 𝑇𝑇 0 (20-𝑇𝑇ref )

, where b T0 = 0.0827 is the base-case value of b T (which is equivalent to a Q 10 value of 2.29). 15). The estimated standard deviation of 𝐴𝐴 𝑇𝑇eff is ±10.7% for the std configuration and ±13.0% for the frag configuration. r S makes the largest contribution of 52.51% to the estimated 𝐴𝐴 𝑇𝑇eff variance.

k POC also contributes substantially to the 𝐴𝐴 𝑇𝑇eff variance, accounting for 35.55%. The percentage contributions from Q 10 , 𝐴𝐴 𝐴𝐴𝑂𝑂 2 , and ϕ L to the variance of 𝐴𝐴 𝑇𝑇eff are relatively small. Q 10 and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 make similar contributions of 5.04% and 4.65%, respectively, while ϕ L contributes only 2.25%. It is interesting to note that while ϕ L contributes least to the estimated variance of 𝐴𝐴 𝑇𝑇eff , it is the model parameter to which 𝐴𝐴 𝑇𝑇eff is most sensitive (in terms of the normalized sensitivity coefficients). This demonstrates that a sensitive parameter is not necessarily important for the output variance if the parameter's reasonable range (defined here as its literature-based variance) is relatively small as is the case for ϕ L .

Regional Differences in the Sensitivity of Transfer Efficiency

In this section, we will use 𝐴𝐴 𝑇𝑇eff to denote the model's T eff response averaged over a specific biogeographical region. Figure 6 shows a regional breakdown of the changes in 𝐴𝐴 𝑇𝑇eff output in response to perturbed values of r S , ϕ L , Q 10 , and 𝐴𝐴 𝐴𝐴𝑂𝑂 2 . Variations in r S have a larger impact on 𝐴𝐴 𝑇𝑇eff at low latitudes (i.e., in regions equatorward of 45°) than at high latitudes (i.e., in regions poleward of 45°), whereas the opposite is the case for ϕ L . r S accounts for 55.3% of the estimated 𝐴𝐴 𝑇𝑇eff variance at low latitudes and 38.9% at high latitudes, while ϕ L contributes 8.9% at high latitudes and 1.19% at low latitudes. Figure 7 demonstrates that large particles contribute more to the bulk sinking POC flux in the high-latitude regions than at low latitudes, whereas the contribution from small particles to the total POC flux is larger at low latitudes than in the high-latitude regions (see Figure C1). As a result, the 𝐴𝐴 𝑇𝑇eff output at high latitudes is more sensitive to assumptions about ϕ L and less sensitive to those about r S , and vice versa for the low-latitude regions.

The most striking regional difference in the 𝐴𝐴 𝑇𝑇eff response is its much larger sensitivity to 𝐴𝐴 𝐴𝐴𝑂𝑂 2 in oxygen-deficient waters than in well-oxygenated waters.

𝐴𝐴 𝐴𝐴𝑂𝑂 2 contributes 12.25% to the estimated variance of 𝐴𝐴 𝑇𝑇eff in suboxic regions (defined here as waters where [O 2 ] is < 5 mmol m -3 ), compared to only 0.08% in high-oxygen regions (defined here as waters where [O 2 ] is > 100 mmol m -3 ). This localized sensitivity suggests that the actual value of 𝐴𝐴 𝐴𝐴𝑂𝑂 2 is not critical on a global scale, unless the threshold oxygen concentration at which anaerobic respiration begins is assumed to be very high (e.g., 𝐴𝐴 𝐴𝐴O 2 = 100 mmol m -3 ). The specification of the value of 𝐴𝐴 𝐴𝐴O 2 , however, is crucial for the ocean's most-oxygen-depleted oxygen minimum zone regions, i.e., at oxygen concentrations less than 5 mmol m -3 . As an example, Figure 8 Different temperature regimes are expected to exhibit varying degrees of Q 10 sensitivity. Interestingly, the regional breakdown provided in Figure 6 shows that the 𝐴𝐴 𝑇𝑇eff output is more sensitive to variations in Q 10 in cold, high-latitude regions (defined here as seasonally stratified subpolar regions north of 45°N or south of 45°S) than in the warmer mid-latitudes and tropics (defined here as permanently stratified subtropical regions and equatorial regions). Q 10 contributes 12.47% to the estimated 𝐴𝐴 𝑇𝑇eff variance in subpolar regions, compared to 3.60% in thermally stratified regions. For subtropical and equatorial regions with a sharp vertical temperature gradient, the positive effect of increasing Q 10 on remineralization rates in warm surface waters is counteracted by the opposite effect at deeper depths where cold waters predominate, leading to a dampening of POC flux attenuation over the top ∼1,000 m. In contrast, when Q 10 is increased at high-latitude locations, the rate of remineralization is decreased everywhere in the vertical since high-latitude water columns tend to be uniformly cold. As an example, the response of the POC flux profile to variations in Q 10 is shown in Figure 9 for two locations belonging to different temperature regimes. The range of T eff values is considerably larger in the Northwest Pacific subarctic gyre than in the North Pacific subtropical gyre, reflecting differences in the steepness of their vertical temperature gradients.

Spatial Variability in Mesopelagic Transfer Efficiency

Global Pattern of Transfer Efficiency

Figure 10 shows the global distribution of T eff from the ctrl simulation for model configurations std and frag. Under both configurations, the model simulates a spatial pattern of T eff that is latitudinal in character with a north-south T eff gradient in both hemispheres. Values of T eff are generally larger at high latitudes (18.5 ± 3.76%) than at low latitudes (11.7 ± 2.47%), although there are a few exceptions to this, including high T eff values in the equatorial Eastern Tropical Pacific. The spread in T eff values for the std configuration is similar to that for the frag configuration. The 5th and 95th percentiles of the T eff values for the std configuration are 8.63% and 21.9%, a spread which is only slightly higher than the 9.41% to 20.8% range for the frag configuration. However, the spatial variance of T eff is noticeably smaller for the frag configuration (14.8) than for the std configuration (20.5). This is because the inclusion of large particle fragmentation in the model under the frag configuration acts to reduce the latitudinal gradient in particle size initially established at the surface, leading to a relatively more homogenous pattern of T eff compared to the std configuration. DINAUER ET AL.

Environmental Drivers of Spatial Variations in Transfer Efficiency

The spatial pattern of T eff has contributions from several environmental factors 𝐴𝐴 𝐴𝐴𝐴𝐴 = (𝐴𝐴𝐴𝐴1, . . . , 𝐴𝐴𝐴𝐴𝑛𝑛) , including ocean temperature, dissolved oxygen, seawater density, seawater viscosity, particle size, and particle density. The focus of this section is primarily on the environmental factors whose spatial variability has a significant impact on the large-scale pattern of T eff . A chi-square variance test shows that the predicted spatial pattern of T eff from the 𝐴𝐴 𝜌𝜌sw experiment (see Table 3) is statistically indistinguishable from the T eff pattern in the ctrl simulation. Hence, the impact of seawater density on large-scale T eff will not be discussed further. One important caveat to this analysis is that the relative importance of environmental factors in the spatial pattern of T eff is sensitive to the assumptions underlying the model. The results of the analysis will depend on model structure and formulations, the choice of parameter values and the surface boundary fluxes used to force the model. In this work, the POC, calcite, aragonite and opal export flux fields are each constructed from a CMIP multi-model ensemble median.

Exploring the sensitivity of model output to alternative particle export flux fields will be the focus of future work. DINAUER ET AL. ef f for a particular location indicates that T eff is larger than the global mean T eff due to spatial variability in EF i , while a negative value means T eff is smaller than the mean due to space-varying EF i . Here, we identify the environmental factors associated with the largest positive and negative ef f = -4.4% ) induce the largest negative T eff changes. Interestingly, particle size is less influential in the frag configuration (not shown in Figure 11). When fragmentation is included in the model, ocean temperature ( 𝐴𝐴 ∆𝑇𝑇 𝑇𝑇 ef f = 5.9% ) has the largest positive impact on T eff , followed by particle size ( 𝐴𝐴 ∆𝑇𝑇 𝜎𝜎 𝐿𝐿 ef f = 4.8% ), whereas particle density ( 𝐴𝐴 ∆𝑇𝑇 𝜌𝜌 PM ef f = -5.6% ), seawater viscosity ( 𝐴𝐴 ∆𝑇𝑇 visc ef f = -3.2% ) and ocean temperature ( ∆𝑇𝑇 𝑇𝑇 ef f = -3.1% ) are the environmental factors with the largest negative impacts on T eff . The less important role for particle size can be explained by a tendency for large particle fragmentation to reduce the DINAUER ET AL. 21 of 34 latitudinal gradient in particle size, initially established in the surface ocean, over the global ocean's mesopelagic zone. As large particles fragment into smaller ones as they settle through the water column, the sinking particle pool becomes mostly dominated by small particles everywhere in the horizontal and thus particle size variability becomes less important large-scale T eff .

While the majority of 𝐴𝐴 ∆𝑇𝑇

EF 𝑖𝑖 ef f patterns exhibit broad latitudinal structures, there exists disagreement on the sign of 𝐴𝐴 ∆𝑇𝑇 EF 𝑖𝑖 ef f over low versus high latitudes. Particle size has a positive impact on T eff at high latitudes and a negative impact at low latitudes. This can be explained by the fact that the opal-dominated high latitudes exhibit the greatest fractions of POC export going into large, rapidly sinking particles, resulting in enhanced vertical POC transfer. Ocean temperature also induces positive changes in T eff at high latitudes and negative T eff changes at low latitudes. The positive impact of temperature on T eff in cold, high-latitude regions is consistent with the rate of chemical reactions decreasing with decreasing temperature, following the classical Arrhenius equation. In contrast to particle size and ocean temperature, both particle density and seawater viscosity have a positive impact on T eff in low-latitude regions and a negative impact in high-latitude regions. The positive impact of particle density on T eff at low latitudes can be explained by the presence of denser ballast minerals (i.e., CaCO 3 ) that act to increase bulk particle sinking speeds. Seawater viscosity positively affects T eff at low latitudes since sinking particles encounter DINAUER ET AL. A few simple statistical measures are used to determine the relative importance of each environmental factor to global T eff variability (see Table 7). One measure of factor importance involves calculating the root-mean-square error (RMSE) between 𝐴𝐴 𝐴𝐴 ctrl ef f and 𝐴𝐴 𝐴𝐴 EF 𝑖𝑖 ef f , which may be thought of as the error occurring when spatial variations in factor EF i are not considered. The errors obtained in this fashion are utilized to calculate a "factor importance index" (FI), which is equal to the RMSE for EF i divided by the greatest RMSE among all environmental factors 𝐴𝐴 𝐴𝐴𝐴𝐴 = (𝐴𝐴𝐴𝐴1, . . . , 𝐴𝐴𝐴𝐴𝑛𝑛) . Particle size, particle density, and ocean temperature emerge as the most important determinants of global T eff variability. For the std configuration, the spatial variability in T eff can be attributed largely to the contribution from particle size variability (FI = 1.0) and in smaller part to particle density (FI = 0.50) and ocean temperature (FI = 0.37). For the frag configuration, the spatial pattern of T eff is attributable in more equal measure to contributions from ocean temperature (FI = 1.0), particle size (FI = 0.83), and particle density (FI = 0.73). The error originating from the removal of spatial variability in seawater viscosity is nontrivial, especially for the frag configuration (FI = 0.54), whereas the errors due to neglecting variations in dissolved oxygen and seawater density are very small or negligible on a global scale. Table 7 introduces two other importance measures that both come to the same ranking of factor importance.

Basic Assumptions Underlying the Model

Although it provides a useful tool for understanding what controls large-scale T eff , the MSPACMAM makes several simplifying assumptions. The treatment of the partitioning of POC export into small and large particles has obvious simplifications. The straightforward approach used in the model assumes that the fraction of POC export entering the large size class, σ L , is proportional to the fraction of total particle export consisting of opal and aragonite, both of which are biogenic minerals associated with large plankton (i.e., diatoms and pteropods, respectively). For simplicity, the constant of proportionality, κ L , between σ L and the relative abundance of opal and aragonite export is assumed to be 1. There are four main reasons for the POC export partitioning scheme used in the MSPACMAM. First, most ESMs prognostically calculate POC, CaCO 3 and opal export fluxes, so the MSPACMAM modeling framework can be readily deployed in such models. Second, the use of the scheme avoids the need to employ some version of a marine ecosystem and/or particle aggregation model in the surface production zone, both of which can introduce new sources of model errors and uncertainties via additional uncertain parameters. Third, the partitioning of POC export into differently sized particles based on the relative abundance of opal and aragonite export has the advantage of being able to respond to changing forcing, specifically any forcing that affects the prognostic calculation of export production, which is not the case for, e.g., a power law-based particle size distribution. Fourth, the resultant spatial pattern of σ L (Figure A1) is consistent with global distributions of size-partitioned carbon biomass derived from the satellite-based particulate backscattering spectrum (Kostadinov et al., 2016), showing that large particles make a larger contribution in the opal-dominated high latitudes than in the CaCO 3 -dominated low latitudes. Future studies may benefit from improvements in the choice of κ L .

The model also makes an assumption about the partitioning of calcite export between small and large particles. It is straightforward to treat the fraction of calcite export consisting of large particles in the same manner as the

Environmental factor SD of 𝐴𝐴 ∆𝑇𝑇 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f (%) Ratio of SD of 𝐴𝐴 ∆𝑇𝑇 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f to SD of 𝐴𝐴 𝐴𝐴 ctrl ef f RMSE between 𝐴𝐴 𝐴𝐴 ctrl ef f and 𝐴𝐴 𝐴𝐴 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f (%) Factor importance Standard configuration Temperature 1.70 0.37 1.74 0.37 Oxygen 0.41 0.09 0.43 0.09 Seawater density 0.01 0.002 0.01 0.002 Viscosity 1.14 0.25 1.14 0.24 Particle size 4.70 1.04 4.72 1.00 Particle density 2.31 0.51 2.35 0.50 Configuration with fragmentation Temperature 3.11 0.81 3.23 1.00 Oxygen 0.28 0.07 0.29 0.09 Seawater density 0.02 0.005 0.02 0.006 Viscosity 1.71 0.44 1.74 0.54 Particle size 2.63 0.68 2.66 0.83 Particle density 2.33 0.61 2.37 0.73 Note. Rows and columns represent individual factors and measures, respectively. 𝐴𝐴 𝐴𝐴 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 ef f is the spatial pattern of T eff from the baseline model simulation in which all factors vary spatially. 𝐴𝐴 𝐴𝐴 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f is the spatial pattern of T eff from a model experiment in which EF i is forced to be spatially uniform (i.e., its spatial variability has been removed). 𝐴𝐴 ∆𝑇𝑇 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f is the component of the spatial pattern of T eff attributable to space-varying EF i , calculated as the difference between the 𝐴𝐴 𝐴𝐴 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 ef f and 𝐴𝐴 𝐴𝐴 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f patterns. The standard deviation SD of the global 𝐴𝐴 ∆𝑇𝑇 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f pattern is shown in the second column. The ratio of the standard deviations of 𝐴𝐴 ∆𝑇𝑇 𝐸𝐸𝐸𝐸 𝑖𝑖 ef f to 𝐴𝐴 𝐴𝐴 ctrl ef f (third column) reveals whether the two patterns have similar amplitudes of variability. Shown in the fourth column is the RMSE between 𝐴𝐴 𝐴𝐴 ctrl ef f and 𝐴𝐴 𝐴𝐴 𝐸𝐸𝐸𝐸 𝑖𝑖

ef f , which is the root-mean-square error occurring when spatial variations in EF i have been neglected. This measure of importance is used to determine the factor importance index (fifth column), defined as the factor's RMSE divided by the maximum RMSE of all factors.

Table 7 Statistical Measures of the Relative Importance of Each Environmental Factor, EF i , in the Spatial Pattern of T eff DINAUER ET AL. 23 of 34 fraction of POC export going into large particles. This approach was motivated in part by the results of laboratory studies, which suggest that diatom aggregates are effective at scavenging coccoliths and other small calcite particles (De La Rocha et al., 2008;Iversen & Ploug, 2010;Passow & De La Rocha, 2006). While it is sensible to assume that some fraction of calcite export participates in the large particle size class, the constant of proportionality associated with this fraction could be revisited in future work. For example, information from food web models could be combined with satellite-based phytoplankton functional type distributions to refine this constant. Another potential target for model improvement is the treatment of refractory and/or protected POC, which has been hypothesized to be associated with the particles sinking out of the surface ocean in the CaCO 3 -dominated low latitudes (Henson et al., 2012). Rubey's (1933) equation is used to calculate the sinking speeds of small and large particles, 𝐴𝐴 𝐴𝐴 𝑆𝑆 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 and 𝐴𝐴 𝐴𝐴 𝐿𝐿 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 . Since it was derived by equating the effective weight of a spherical particle to both viscous and inertial forces, Rubey's equation represents an improvement over Stokes' law of viscous resistance. Its application is suitable for the turbulent flow regime in which inertial forces dominate over viscous forces, i.e., when the particle Reynolds number Re p is outside the Stokes' range (Re p > 0.5). In contrast, Stokes' law applies only in the laminar (or viscous) flow regime where Re p < 0.5. The particle settling velocities given by the two equations are similar at low Re p but diverge for particles sinking outside the Stokes' range, with the velocities given by Rubey's equation lower than those given by Stokes' law. While a number of possibly important factors (e.g., the effects of shape and surface roughness) have been neglected in its derivation, the advantage of using Rubey's equation in the model is that it allows for the calculation of environmentally dependent, space-varying 𝐴𝐴 𝐴𝐴 𝑆𝑆 sink and 𝐴𝐴 𝐴𝐴 𝐿𝐿 sink . The representation of the effects of zooplankton migration, such as large particle fragmentation, should be a priority target for model improvement. The MSPACMAM does not include diel and/or seasonal vertical migration of zooplankton and only implicitly represents zooplankton-mediated fragmentation by assuming that the rate constant 𝐴𝐴 𝐴𝐴 𝐿𝐿 frag encompasses all physical and biological processes that fragment large particles into smaller ones. While it is likely that neglecting active transport of carbon by vertically migrating zooplankton has a relatively small impact on the climatological annual mean spatial pattern of T eff since POC export at a large spatial scale is thought to be dominated by the gravitational settling of particles (Boyd et al., 2019;Resplandy et al., 2019), recent field work has demonstrated that particle fragmentation is a key process governing POC flux attenuation in a regional context. For example, the breakup of particles by biotic and/or abiotic mechanisms may account for almost half of the observed large particle flux loss in the North Atlantic and Southern Ocean (Briggs et al., 2020), whereas the general absence of zooplankton in oxygen minimum zones contributes to the preservation of sinking POC in these regions (Cavan et al., 2017). A major innovation in the MSPACMAM would be the incorporation of mechanistic representations of fragmentation mechanisms, such as physical disaggregation by turbulence-generated shear stress (Alldredge et al., 1990) and zooplankton "sloppy feeding" (Anderson & Tang, 2010;Lampitt et al., 1990) and/or swimming behavior (Dilling & Alldredge, 2000;Goldthwait et al., 2004). However, a more complex treatment of particle fragmentation will be constrained by available observations and limited by current knowledge. Future research efforts could focus on experimental approaches and field studies in order to generate a more detailed understanding of the individual fragmentation processes and their relative importance in the net vertical transport of carbon.

Synthesis of Transfer Efficiency Sensitivity Results

In Section 3.2, we showed that reasonable changes in model parameters can significantly alter the shape of the vertical POC flux profile. Similar to Omand et al. (2020), we find that the flux profile is very sensitive to the choice of the base remineralization rate constant k POC . T eff ranged from 11.5% to 39.0% when k POC referenced to 20°C was varied from 0.035 to 0.32 d -1 . Omand et al. (2020) showed that the export efficiency (defined as export at 100 m divided by net primary production) ranged from 40% to 82% when their model's remineralization rate was varied from 0.01 to 0.16 d -1 . Cram et al. (2018) and Omand et al. (2020) found that the slope of the particle size distribution, which determines the proportion of small, slow-sinking particles versus large, fast-sinking particles, has a significant influence on the bulk POC flux profile. While our model structure does not rely on a power law-based particle size distribution, varying the small particle radius r S , while holding all else constant, has the effect of varying bulk particle sinking speeds. As expected, our sensitivity analysis demonstrated that the choice of the small particle radius r S strongly influences the bulk POC flux profile. As in Cram et al. (2018), the influence of Q 10 on vertical POC flux attenuation varies DINAUER ET AL.

10.1029/2021GB007131

24 of 34 regionally between subtropical and high-latitude regions, which can be attributed to regional differences in the steepness of the vertical temperature gradient. Also consistent with the findings of Cram et al. (2018), the modeled POC flux profile is rather insensitive to the oxygen dependence of remineralization, 𝐴𝐴 𝐴𝐴𝑂𝑂 2 , except in suboxic regions where oxygen is effectively absent from seawater.

As highlighted by DeVries and Weber (2017), the factors controlling the spatial pattern of T eff are the subject of active debate. In experiments with the MSPACMAM, three environmental factors emerged as the most important determinants of the spatial pattern of T eff (see Section 3.3). The largest contributors to the simulated T eff pattern are particle size, particle density, and ocean temperature, with the caveat that these results depend on our model structure/formulations, parameter values, and boundary conditions. The MSPACMAM predicts that vertical POC transfer is more efficient in the cold, opal-dominated high latitudes, which can be explained by the presence of large, rapidly sinking particles and the reduction in the temperature-dependent rate of remineralization, both of which compensate for the absence of dense ballast provided by CaCO 3 and the viscous resistance offered by cold waters at high latitudes. Maerz et al. (2020) reported similar results. Cram et al. (2018) have highlighted the importance of ocean temperature and particle size in the spatial pattern of T eff , whereas DeVries and Weber (2017) concluded that the general pattern of T eff is shaped by temperature and modulated by low oxygen levels. Likewise, our model results indicate that, although it contributes little to the global T eff variability, dissolved oxygen has a strong local impact in the global ocean's three major oxygen minimum zones, i.e., the Eastern Tropical North Pacific, the Eastern Tropical South Pacific, and the Arabian Sea. Cram et al. (2018) also find that oxygen explains very little of the global T eff variability outside oxygen minimum zones. Maerz et al. (2020) appear to overestimate the extent and intensity of T eff in the equatorial Eastern Tropical Pacific, which may be attributed to a tendency for ESMs to overestimate the volume of oxygen minimum zones (Bopp et al., 2013;Cabré et al., 2015;Cocco et al., 2013;Kwiatkowski et al., 2020).

Summary and Conclusions

The spatial pattern of the mean state of T eff simulated by the MSPACMAM is consistent with several recent model studies (Cram et al., 2018;DeVries & Weber, 2017;Maerz et al., 2020;Weber et al., 2016), showing that vertical POC transfer is more efficient at high latitudes than at low latitudes. This correspondence between five structurally different models suggests that the models correctly simulate the climatological annual mean spatial pattern of T eff . However, it is unclear whether high-latitude T eff is higher than low-latitude T eff on all timescales. The simulation of the seasonal cycle of T eff requires high-quality time-series measurements of sinking POC fluxes at strategic geographical locations. The ability of models to simulate the seasonality in POC fluxes is a necessary prerequisite for improving our predictive understanding of seasonal variability in the efficiency of the ocean's biological carbon pump.

The spatial pattern of T eff in a warming climate is not known. Experiments with the MSPACMAM suggest that particle properties (i.e., size and density), ocean temperature, and to a lesser extent, seawater viscosity drive most of the T eff variability over the contemporary global ocean, whereas seawater density has a negligible influence and dissolved oxygen only a local impact. The direct effect of elevated ocean temperatures will be to enhance remineralization rates and therefore reduce T eff , whereas the smaller viscous resistance offered by warmer waters will lead to higher T eff . The indirect effects of ocean warming on particle properties are more complex. Changes in the species composition (floristic shifts) and geographical distribution (faunistic shifts) of the phytoplankton community are expected to accompany climate change (Basu & Mackey, 2018;Boyd, 2015;Boyd et al., 2010;Passow & Carlson, 2012). There is a tendency for models to simulate a shift toward small-celled phytoplankton in response to climate change, along with a poleward shift in coccolithophorids in a warming ocean. However, in addition to uncertainties in future projections of phytoplankton distributions, the associated changes in the biophysical properties of sinking particles are not well established. As a consequence, it is difficult to speculate about whether changes in particle size and density under a warming climate will enhance or reduce T eff . A model such as the MSPACMAM could be deployed in ESMs in order to quantify the cumulative first-order effects of climate change on T eff . Recommendations for future modeling efforts include (a) assessing the reliability of projected future changes in the phytoplankton community structure and (b) investigating how spatial patterns of particle properties will change in response to climate change. DINAUER ET AL.

10.1029/2021GB007131 27 of 34

The performance of the MSPACMAM under various parameter settings was evaluated against observed normalized sinking POC flux and log-transformed oxygen consumption rate data. Figure B1 compares the modeled and observed normalized POC flux profiles at 22 locations for which trap-derived sinking POC flux measurements were available. The data locations represented a wide range of temperature and oxygen regimes (Figure B2 and Figure B3). The model parameter combinations yielding cost function values within 5% of the minimum cost value are summarized in Table B1. r S (μm) ϕ L (%) 𝐴𝐴 𝐴𝐴 𝑟𝑟𝑟𝑟𝑟𝑟 𝑃𝑃 𝑃𝑃𝑃𝑃 (d -1 ) a Q 10 (unitless) 𝐴𝐴 𝐴𝐴𝑂𝑂 2 (μmol L -1 ) A priori 5-30 95-99 0.033-0.32 1.51-3.9 1-20 Standard configuration (N = 1,521)

Range 10-30 95-99 0.033-0.32 1.51-2.09 1-4 Configuration with fragmentation (N = 1,606) Range 10-30 95-98 0.033-0.32 1.51-2.29 1-4 a The POC remineralization rate constant, k POC , referenced to 20°C. Since the original, literature-based k POC values are referenced to various temperatures T ref (see Table 2), they have been normalized to a common reference temperature of 20°C by multiplying by 𝐴𝐴 𝐴𝐴 𝑏𝑏 𝑇𝑇 (20-𝑇𝑇ref ) where 𝐴𝐴 𝐴𝐴𝑇𝑇 = 𝑙𝑙𝑙𝑙 𝑙𝑙10∕10 . For k POC in the model parameter sets, the accompanying Q 10 value is used. For the a priori literature-based values of k POC , a Q 10 value of 1.51 is used.

Table B1

Ranges of the Parameter Sets Within 5% of the Minimum Cost Value