Particle cycling rates at Station P as estimated from the inversion of POC concentration data

The ocean influences the concentration of carbon dioxide (CO₂) in the atmosphere, partly due to a network of biologically mediated processes known as the biological carbon pump (BCP; Ducklow et al., 2001). Photosynthetic algae fix CO₂ into biomass (i.e., particulate organic carbon, POC) in the surface ocean. As algal cells collide and coalesce, or are consumed by zooplankton and repackaged into fecal pellets, they produce aggregates and become part of the marine snow that settles through the water column. Particles can also disaggregate as a result of a combination of physical (e.g., shear stress) and biological (e.g., zooplankton feeding) processes (Hill, 1998; Dilling and Alldredge, 2000; Stemmann et al., 2004). Most POC is fragmented and/or solubilized in the upper water column, particularly in the mesopelagic zone between about 100 and 1000 m, which precludes long-term carbon sequestration (Martin et al., 1987). Particles that are larger or denser sink faster and are more likely to reach the deep ocean, where their decomposition products can reside for centuries (Buesseler et al., 2007b; Marsay et al., 2015). Thus, processes that facilitate the aggregation of particles into larger and/or denser units tend to increase the “effectiveness” of the BCP by promoting long-term carbon sequestration, whereas those that enhance disaggregation tend to reduce long-term sequestration. Accordingly, understanding the particle cycling processes is critical for assessing the effectiveness of the BCP in transferring atmospheric CO₂ to the deep ocean and the ocean’s role in the global climate system.

In this paper, we estimate particle cycling rates in the upper ocean by combining POC concentration ([POC]) data with a two-particle size class model that includes a second-order formulation of particle aggregation and an explicit formulation of diel vertical migration (DVM) by zooplankton. Data used in this study originate from the NASA EXport Processes in the Ocean from RemoTe Sensing (EXPORTS) program, which seeks to “develop a predictive understanding of the export and fate of global ocean net primary production and its implications for present and future climates” (Siegel et al., 2016). Specifically, we use data gathered during an extensive campaign that took place in summer 2018 at Station P in the Gulf of Alaska (EXPORTSNP). This campaign presents an unprecedented opportunity to estimate particle cycling rates in the euphotic zone (EZ) and the upper mesopelagic zone (UMZ) of the ocean. Several research groups in EXPORTS are investigating specific processes involved in the cycling of organic carbon, such as export production, bacterial respiration, and zooplankton migration, from concurrent field measurements. These investigations provide independent estimates of particle cycling rates that are inferred here from the inversion of [POC] data. Our estimates of bulk POC budgets for the upper water column should provide a useful context for the interpretation of observations that pertain to specific biogeochemical processes. We expect that the incorporation of results from process-specific studies with the present, integrative study will help to develop a comprehensive understanding of the processes that governed the efficiency of carbon export to the deep ocean during EXPORTSNP.

The organization of this paper is as follows. First, in Section 2 we summarize how particle cycling rates have been estimated in previous studies. In Section 3, we describe the hydrographic context at Station P. We then present the data and the two-particle class model of POC cycling used in this study in Sections 4.1–4.2. In Section 4.3, the inverse method that is applied to infer particle cycling rates from [POC] data from the upper 500 m at Station P is described. Prior knowledge of particle cycling processes is assembled in Section 4.4. Results from inversions are presented in Section 5, and discussed in Section 6. Finally, we provide a summary in Section 7 and outline future research directions in Section 8.

Particle cycling processes in the ocean (e.g., sinking, remineralization, aggregation, disaggregation, zooplankton ingestion and egestion) are difficult to measure directly and likely vary in magnitude on a wide range of temporal and spatial scales, particularly along the water column and across biogeochemical regimes. Amongst the most widely studied of these processes is particle sinking. Laboratory experiments have provided means of quantifying particle sinking speeds in controlled environments. Such experiments have involved the use of settling columns (O’Brien et al., 2006; Ploug et al., 2008b; Bach et al., 2012), roller tanks (Ploug et al., 2010), and vertical flow chambers (Ploug and Jørgensen, 1999; Ploug et al., 2008a; Iversen and Ploug, 2010). While roller tanks and vertical flow chambers allow for less destructive manipulation of particles that are collected in situ (Ploug et al., 2010), all three methods fall short in replicating in situ conditions that may affect sinking speed.

Particle sinking speeds are also measured in situ by using sediment traps. This can be done by measuring the time-lag of particle export events recorded by traps deployed at different depths in the water column (Berelson, 2002; Honjo et al., 2008; Xue and Armstrong, 2009). Sediment traps equipped with an indented rotating sphere can facilitate estimation of sinking speeds by allowing for collection of particles in discrete speed bins, which is made possible by a built-in settling column (Peterson et al., 2005; Peterson et al., 2009; Armstrong et al., 2009; Lee et al., 2009; Alonso-González et al., 2010; Ebersbach et al., 2011). Traps can also be amended with cameras (Asper, 1987; Diercks and Asper, 1997; Iversen et al., 2010; Nowald et al., 2015) and polyacrylamide gel traps (McDonnell and Buesseler, 2010; 2012) to estimate sinking speeds of relatively small particles. However, estimation of particle sinking speeds by traps is hampered by a number of issues such as collection inefficiences, contamination by zooplankton (i.e., “swimmers”), and post-collection solubilization of particles (Buesseler et al., 2007a).

Estimates of particle sinking speed in situ may also be obtained by using other techniques built on the principle of the settling column. The marine snow catcher allows for field collection of a water sample in a long cylindrical tube. Upon retrieval, subsamples are taken along the vertical axis of the tube in order to categorize particles into suspended, slow, and fast-sinking bins (Riley et al., 2012; Cavan et al., 2015; Giering et al., 2016; Cavan et al., 2017). Settling columns can otherwise be deployed in concert with optical instruments to estimate sinking rates in situ (Kineke et al., 1989; Spinrad et al., 1989; McCave and Gross, 1991; Murray et al., 1996; Smith and Friedrichs, 2015), although such methods are usually employed in regions with higher particle loads such as in coastal or near-bottom environments (Giering et al., 2020).

As with particle sinking speed, different techniques are available to estimate rates of particulate remineralization both in the field and in the laboratory. In the field, particulate remineralization rates have been estimated from changes in dissolved oxygen concentration (Boyd et al., 2015), electron transport activity in microplankton (Packard et al., 1988), apparent oxygen utilization (Jenkins, 1998; Pahlow, 2000), and bacterial carbon demand (Smith et al., 1992). These oxygen-based methods require assumptions about the relationship between carbon supply and oxygen consumption, and may be confounded by factors such as changes in ambient temperature (Boyd et al., 2015). Laboratory experiments have used sensors to measure oxygen consumption directly from resuspended marine particles (Ploug and Grossart, 1999; Ploug and Jørgensen, 1999; Ploug et al., 1999; Ploug et al., 2008b; Iversen et al., 2010; Iversen and Ploug, 2010; 2013), but such studies have relied on assumptions regarding how diffusion of oxygen in particles is related to remineralization (Boyd et al., 2015).

One of the most recently developed techniques to measure particle remineralization in situ involves the REspiration of Sinking Particles In the subsuRface ocEan (RESPIRE) particle interceptor, which captures and incubates marine particles (Boyd et al., 2015; McDonnell et al., 2015). This instrument mitigates shortcomings associated with previous techniques by allowing for measurement of oxygen consumption directly from particle environments. However, it is also prone to issues of its own that largely arise from particle transformations in the incubation chamber (Boyd et al., 2015).

Unlike particulate sinking and remineralization, no methods exist (to our knowledge) for directly measuring particulate aggregation or disaggregation rates in the ocean (Giering et al., 2020). Such rates have been measured in laboratory experiments that rely on imaging techniques to observe the coalescence and breakup of particles in settling tanks (e.g., Alldredge et al., 1990; Waite et al., 1997; O’Brien et al., 2004), but these experiments lack the ability to replicate in situ conditions that may influence (dis)aggregation rates (Jackson, 2015). Models have been combined with field observations of particle size distributions to estimate rates of (dis)aggregation in situ (Stemmann et al., 2004; Karakaş et al., 2009). However, such observations utilize imaging techniques with particle diameter detection limits on the order of 100 µm, thus excluding interactions between smaller particles that often make up a considerable fraction of POC in ocean waters (Lam et al., 2018). More recently, rates of in situ particulate fragmentation were estimated from backscattering sensors onboard Biogeochemical-Argo floats (Briggs et al., 2020). As interactions between particles of different sizes represented in backscattering data are better understood and more floats are deployed, these sensors may permit estimation of particulate (dis)aggregation rates across biogeochemical provinces with high vertical resolution.

Diel vertical migration likely contributes a considerable fraction of global carbon export out of the EZ (Archibald et al., 2019), but the amount of POC consumed by zooplankton in the EZ that is subsequently egested as fecal pellets at depth is poorly constrained (Maas et al., 2021). Indeed, the amount of POC exported out of the EZ by DVM is highly variable and depends on many factors including zooplankton community composition as well as food lability and availability (Steinberg and Landry, 2017). While direct measurements and application of allometric relationships have been used to quantify fecal pellet egestion by mesozooplankton in the EZ during EXPORTSNP (Stamieszkin et al., 2021), the egestion flux resulting from a daytime zooplankton biomass maximum below the EZ (Omand et al., 2021) remains unconstrained.

In addition to the above methods, measurements of the concentration of particle-reactive radionuclides in sediment traps and water samples have been used to infer rates of particle cycling in situ. In a number of pioneering studies, R. Murnane and his colleagues have constrained particle cycling rates of a two-particle (size) class model from measurements of thorium (Th) isotope concentration and vertical particle fluxes using inverse methods (Murnane et al., 1990; Murnane, 1994; Murnane et al., 1994; Murnane et al., 1996). Th has a pronounced affinity for marine particles and has long been used as a tracer of particle processes in the ocean (e.g., Bacon and Anderson, 1982; Nozaki et al., 1987; Clegg and Whitfield, 1990; Clegg et al., 1991; Lerner et al., 2017; Black et al., 2018). The use of a two-class model is convenient as the particle classes in the model – including small “suspended” particles and large sinking particles – approximate the size-fractionated particles traditionally sampled in the field. Moreover, a two-class model requires the solution of only a few differential equations, making it more tractable than other models which describe the particle size spectrum and particle processes in more detail but are computationally more expensive (e.g., Burd, 2013). On the other hand, the approximation of the particle size spectrum with only two classes implies that parameters of two-class models should be interpreted with caution.

The studies by R. Murnane and his colleagues have paved the way towards the rigorous interpretation of oceanic data in terms of particle processes using inverse methods, which allow for consideration of both data and model errors. More recently, inverse methods have been applied to analyze particle and Th isotope concentration data from the U.S. GEOTRACES North Atlantic section GA03 (Lerner et al., 2016; Lerner et al., 2017). However, these studies relied on a single-particle class model and did not aim to estimate rates of particle (dis)aggregation. Studies that did report estimates of particle (dis)aggregation rates of a two-class model from the inversion of field data have been limited in their ability to produce precise estimates partly due to low vertical sampling resolution. For example, studies of R. Murnane and colleagues have relied on sediment trap measurements of POC flux from only a few depths, often below 1000 m (e.g., Murnane et al., 1990; Murnane, 1994; Murnane et al., 1994). The extent to which the particle cycling rates estimated in these studies could apply to other oceanic regions and the upper ocean is unclear, as particle cycling processes are expected to show large spatial and temporal variations. Other efforts to estimate cycling rates in the upper ocean (Murnane et al., 1996) have also been limited by poor vertical resolution. The studies mentioned in this paragraph also assumed that particulate aggregation is a first-order process with respect to particle concentration, which is likely not an accurate representation of particle coalescence in nature (Jackson and Burd, 2015). Finally, these studies all lacked a formal representation of the effect of DVM on particle cycling.

Particle samples were collected during EXPORTSNP in a region surrounding Station P from August 14 to September 9, 2018 (Siegel et al., 2021). Station P is located in a high-nutrient, low-chlorophyll region in the Gulf of Alaska at 50°N, 145°W (Peña and Varela, 2007). In order to describe the hydrographic context, cruise-averaged profiles of various properties measured during CTD casts are presented (Figure 1). At the time of survey, the sea surface temperature was about 14.1°C and decreased to about 3.9°C at 500 m, with a thermocline extending from the base of the mixed layer at about 30 m to roughly 120 m. Salinity was about 32.3 in the mixed layer and increased to approximately 34.1 at 500 m. An upper halocline spanned from the base of the mixed layer to roughly 90 m and overlaid a deeper, stronger halocline extending to about 150 m. A strong halocline is a well-known feature of the local hydrography and limits mixing of deep waters with shallower depths (Plant et al., 2016). Similarly to the halocline, the pycnocline was separated into an upper region ranging from the base of the mixed layer to about 90 m, and a lower region ranging from about 90 m to about 150 m. The concentration of chlorophyll (measured from fluorescence) was approximately 0.17 mg m⁻³ in surface waters and increased to roughly 0.29 mg m⁻³ in the deep chlorophyll maximum, which occurred at about 60 m. Deeper in the water column, chlorophyll concentration decreased sharply down to about 120 m. The concentration of dissolved oxygen was about 265 µmol kg⁻¹ in the mixed layer, peaked at roughly 306 µmol kg⁻¹ near 50 m, and decreased to about 39 µmol kg⁻¹ at 500 m. Beam attenuation coefficient for particles (c_p) measured by transmissometry was approximately 0.11 m⁻¹ in the mixed layer, and dropped sharply from the base of the mixed layer to about 100 m, reaching about 0.01 m⁻¹ at 500 m.

For the purpose of this study, the depth of the mixed layer was chosen as 30 m based on the cruise-averaged value of 31 m (Siegel et al., 2021). The base of the EZ was taken as 100 m, where photosynthetically active radiation (PAR) measured by a radiometer was typically between 0.1–1% of PAR at the surface (Estapa et al., 2021; Siegel et al., 2021). The base of the EZ coincided approximately with the depth at which the halocline and pycnocline transitioned from stronger to weaker gradients, indicating that the EZ depth corresponded with an important physical as well as biological boundary.

All particle samples considered for this study were collected via large volume in situ filtration (LVISF; Roca-Martí et al., 2021) during multiple casts at approximately the same six depths (Table 1). For two of the casts, a seventh sample was collected within the mixed layer. The samples were fractionated into three size-classes (1–5 µm, 5–51 µm, and > 51 µm) by using quartz microfiber filters (QMA, 1-µm nominal pore size) and Nitex screens (5- and 51-µm nominal pore size). For this study, POC concentrations from the two smallest size fractions were added to obtain what is referred to here as the small size fraction (SSF), representing particles between 1 and 51 µm. The other size fraction (> 51 µm) is referred to as the large size fraction (LSF). POC concentrations were obtained from the difference between blank-corrected total particulate carbon measured using a high-temperature combustion technique and particulate inorganic carbon measured by coulometry (Roca-Martí et al., 2021).

The uncertainty in the cruise-averaged [POC] at a given depth and for a given size fraction was taken as the standard error of the [POC] data at that depth and for that size fraction from the different casts. Because the mixed layer was only sampled twice by LVISF, the uncertainty in the cruise-averaged [POC] in the mixed layer (for SSF or LSF) is estimated by multiplying the cruise-averaged [POC] in the mixed layer by the coefficient of variation of [POC] data at 50 m divided by the square root of the number of samples (N = 2) in the mixed layer. Note that the uncertainties of cruise-averaged quantities considered in this paper only account for inter-station variability. We expect that such variability, and not the errors associated with sampling, storage, and/or measurement, dominates the uncertainties in the cruise averages.

Note that the 1-µm nominal porosity of the QMA filters used for LVISF likely resulted in poor collection efficiency of smaller particles such as picoplankton and bacteria observed at Station P (Stephens et al., 2020). Indeed, McNair and Menden-Deuer (2020) showed that at the time of sampling, phytoplanktonic cells at 95 m were generally smaller than 1.2 µm and dominated by Synechococcus spp. Thus, [POC] measured on samples collected by LVISF may underestimate [POC] that was present at Station P.

In this section, we describe the model of POC cycling that is considered in this study. The model domain represents the water column between the sea surface and 500 m and includes therefore both the EZ (0–100 m) and the UMZ (100–500 m). It is divided into seven layers whose boundaries coincide with LVISF sampling depths (except at the sea surface). The upper three layers are in the EZ, while the lower four layers are in the UMZ (Figure 1).

Our model of POC cycling includes two particle size classes (Figure 2) and largely relies on a previous model described by Clegg and Whitfield (1990) and later applied by Murnane and colleagues (Murnane et al., 1990; Murnane, 1994; Murnane et al., 1994; Murnane et al., 1996). The two size classes of the model are intended to represent the size fractions 1–51 µm and > 51 µm sampled during EXPORTSNP via LVISF, and are consistent with similar studies that have used LVISF to collect particles in “small” and “large” size fractions (e.g., Lam and Bishop, 2008; Bishop et al., 2012; Black et al., 2018). The model assumes steady state and neglects the effects of lateral transport processes, i.e., it is a steady state, one-dimensional (1D, vertical) model. Furthermore, it omits the influences of vertical transport by turbulent mixing, upwelling, and downwelling. These various assumptions imply that caution should be exercised when applying the model to interpret data suspected to be influenced by unsteadiness (i.e., departures from steady state), such as during bloom conditions, and physical transport processes, as in the vicinity of strong hydrographic fronts (Mahadevan, 2016). In addition, there may be biogeochemical processes that are not explicitly represented in the model but that may affect the distribution of POC. Following the approach of R. Murnane and colleagues, model errors are represented in our [POC] data analysis via residuals introduced in the model equations.

The balance equations describing the cycling of POC in the SSF and the LSF in our model are, respectively:

Here, P_S (P_L) is [POC] in the SSF (LSF), $P˙S$ is the primary production of P_S, w_S (w_L) is the sinking speed of P_S (P_L), $β−1,S$ (⁠$β−1,L$⁠) is the rate constant for remineralization of P_S (P_L), $β2′$ (⁠$β−2$⁠) is the rate constant for P_S aggregation (P_L disaggregation), β₃ is the rate constant for ingestion of P_S by vertically migrating zooplankton in the grazing zone, H is a Heaviside function (⁠$H(zg−z)=1$ if $z<zg$ and $H(zg−z)=0$ otherwise, where z_g is the depth of the grazing zone and z is the depth), E(z) is a function representing the vertical distribution of zooplankton egestion, and $P¯SEZ$ is the average P_S concentration in the EZ. In this study, “remineralization” includes all processes that transfer organic carbon from the particulate (> 1 µm) phase to the organic or inorganic “dissolved” (< 1 µm) phase. ε_S (ε_L) is the residual or error in the equation for P_S (P_L) due to the omission of unsteadiness, transport, and other processes that are not explicitly represented in the model, as well as the sum of errors in the processes that are explicitly represented in the model.

The interpretation of Equations 1.1 and 1.2 is straightforward. Equation 1.1 states that the sources of P_S from particle production and large particle disaggregation should be balanced (within an error given by ε_S) by the divergence of the SSF sinking flux and by the sinks from small particle remineralization, aggregation, and zooplankton ingestion. Likewise, Equation 1.2 states that the sources of P_L from small particle aggregation and zooplankton egestion should be balanced (within an error given by ε_L) by the divergence of the LSF sinking flux and by the sinks from large particle remineralization and disaggregation. A positive residual (ε_S or ε_L) would indicate that the model is missing a source of POC in the corresponding size fraction, while a negative residual would imply a missing loss.

As evidenced by Equations 1.1 and 1.2, all processes that represent the exchange of material between different size fractions are assumed to be first-order with respect to POC concentrations, except for particle aggregation, which is assumed to be second-order with respect to P_S (Jackson and Burd, 2015). The corresponding rate constants are thus first-order (apparent) rate constants with the exception of $β2′$⁠, which is a second-order (apparent) rate constant.

The production rate of small particles, $P˙S$⁠, is formulated as follows in our model. Net primary production (NPP) as measured from radiocarbon incubations was observed to decrease with depth in a quasi-exponential fashion. Thus, primary production of P_S (⁠$P˙S$⁠) in our model is assumed to decrease with depth below the mixed layer according to:

where $P˙S,ML$ is the production rate of P_S in the mixed layer, $zML$ is the mixed layer depth (30 m), and L_P is a vertical scale describing the attenuation of particle production with depth along the water column.

Ingestion of POC and its subsequent egestion by vertically migrating zooplankton acts as a mechanism of carbon export from shallow to deep waters. We assume that POC is consumed as P_S and egested as P_L, which is consistent with the notion that zooplankton repackage smaller units of algal biomass into larger fecal pellets. We also assume that ingestion occurs only in the EZ, where zooplankton feed on algal biomass at night, and that egestion occurs only in the UMZ, where zooplankton reside during the day. Hence, the depth of the ingestion or grazing zone (z_g) is set to the base of the EZ (100 m).

Omand et al. (2021) observed that daytime zooplankton biomass during EXPORTSNP increased from 100 m to a subsurface maximum at 300–400 m, and then decreased with depth. We model this observation in our formulation of P_L egestion. In Equation 1.2, β₃$P¯SEZ$ is the average P_S ingestion rate from the surface to z_g, where

To model the volumetric egestion rate (Equation 1.2), we multiply $β3$$P¯SEZ$ by the dimensionless function $E(z)$⁠:

Here, z_m is the maximum depth at which zooplankton egest POC and $Emax$ is the maximum value of the egestion function $E(z)$⁠, which occurs at $z=(zg+zm)/2$⁠. The value of $Emax$ is determined from a “metabolic constraint” that sets the vertical integral of the egestion flux in the egestion zone (between z_g and z_m) to be a fraction, α, of the vertical integral of the ingestion flux in the grazing zone (between 0 and z_g), such that

Consistent with our assumption that POC is not egested above z_g, the Heaviside function in Equation 4 implies that $E(z)=H(z−zg)=0$ when $z<zg$⁠. The vertical distribution of the zooplankton egestion rate as parameterized by Equation 4 is shown in Figure 3.

Altogether, our model of POC cycling includes a total of 11 parameters (w_S, w_L, $β−1,S$⁠, $β−1,L$⁠, $β2′$⁠, $β−2$⁠, $P˙S,ML$⁠, L_P, β₃, α, z_m; see Table 2). For simplicity, these parameters are referred to as “model parameters” henceforth. The parameters $β−1,S$⁠, $β−1,L$⁠, $β2′$⁠, and $β−2$ are assumed to be uniform within each of the seven layers, but are allowed to vary between layers. The parameters w_S and w_L are defined at the layer boundaries and are allowed to vary between layer boundaries. Finally, $P˙S,ML$⁠, L_P, β₃, α and z_m take on a single value within the entire model domain.

In this study, Equations 1.1 and 1.2 are integrated over the thickness of each model layer, which yields, respectively,

In Equations 6.1 and 6.2, z_t (z_b) is the depth at the top (bottom) of the layer, $h=zb−zt$ is the thickness of the layer, and an overbar denotes an average in the layer. The average [POC] in a given layer (⁠$PS¯$ or $PL¯$⁠) is taken as the arithmetic mean of the concentrations at the two layer boundaries (e.g., $P¯S=[PSt+PSb]/2$⁠) except in the mixed layer, where it is taken as the concentration at the base of the mixed layer. This approach does not involve the implicit assumption that the POC concentrations are varying linearly with depth within each layer, because the POC concentrations in the source/sink terms are not specifically defined as the concentrations that prevail at depths mid-way between the layer boundaries. Our approach assumes that the POC concentrations that appear in the source/sink terms in the POC balance equations can be approximated from the concentrations at the layer boundaries, but it does not make specific assumptions about the vertical variations of the POC concentrations within each layer. Finally, $εS¯h$ (⁠$εL¯h$⁠) is the layer-integrated residual in the SSF (LSF) POC balance. Note that in the mixed layer, the top boundary is the sea surface and $wStPSt=wLtPLt=0$⁠.

The integrals in Equations 6.1 and 6.2 are evaluated as follows:

The evaluation of the integrals in Equations 7.2–7.3 is facilitated by our choice to coincide the depth of the ingestion zone (z_g) with a layer boundary.

The [POC] data described in Section 4.1 are used to estimate the parameters of the POC model described in Section 4.2 by using an inverse method. Central to this method is the concept of a state vector, x, which is a vector of variables that describe POC cycling in the upper 500 m at EXPORTSNP stations according to the model. The state variables (elements of x) are P_S and P_L at each sampling depth, the particle cycling parameters, and the errors in the model equations (⁠$εS¯h$ and $εL¯h$⁠), in each layer (Table 2).

The inverse method applied in this work proceeds as follows. Equations 6.1 and 6.2 are collected in a single equation $f(x)=0$⁠, where $f(x)$ is a vector of functions and 0 is a vector of zeros. Some of the terms in $f(x)$ contain products of unknowns (e.g., β₋₂P_L), which introduces non-linearity into our estimation problem and renders popular methods such as ordinary least-squares inadequate. Here, we utilize the algorithm of total inversion (ATI) (Tarantola and Valette, 1982), which is an iterative approach to solving a non-linear least-squares problem. The ATI is used to find an estimate of the state vector x that minimizes the objective (or cost) function,

where superscript T denotes the transpose. Here, x_o is a vector including prior estimates of the elements of x and C_o is the error covariance, or uncertainty, matrix for the prior values in x_o: the diagonal elements of C_o are the squares of the uncertainties in the prior estimates in x_o and the off-diagonal elements of C_o are the covariances between these uncertainties. While the uncertainties in some prior estimates may be expected to show some amount of covariance (e.g., the uncertainty for a rate parameter in a given layer may co-vary with the uncertainty for the same parameter in neighboring layers), the covariances between the uncertainties in the elements of x_o are largely unknown and set to zero for simplicity.

In the ATI, the vector x that minimizes the cost (Equation 8) is estimated by iteration. At each iteration, the model equations are linearized about the state estimate from the previous iterative step, so that the original non-linear inverse problem is replaced by a sequence of linear inverse problems. The state estimate, $x^$⁠, at iteration k + 1 is given by

Here F_k is a matrix that contains the derivatives of the model equations with respect to the variables in x evaluated at the $kth$ iteration (i.e., a Jacobian matrix), and f_k is a vector that contains the left-hand side of the model equations $f(x^k)=0$⁠.

The iterative solution (Equation 9) is computed until convergence is reached. At the end of each iteration, the estimated state variables are extracted from $x^k+1$⁠. The iteration continues until every variable in $x^k+1$ changes by less than $10−6$ relative to its estimate at the previous iteration in $x^k$⁠. When this criterion is met, convergence is assumed to have been reached and the iteration is terminated. The error covariance matrix of $x^k+1$ is approximated from (see Tarantola and Valette, 1982)

The uncertainties in the state estimates in $x^k+1$ are the square root of the diagonal elements of $Ck+1$⁠. For convenience, the uncertainties extracted from $Co$ and $Ck+1$ will be generically denoted below as $σo$ and $σk+1$⁠, respectively. Because the second term on the right-hand side of Equation 10 is a positive-definite matrix, the posterior uncertainties in $Ck+1$ are always less than or equal to the prior uncertainties in C_o (i.e., $σk+1≤σo$⁠) (Tarantola and Valette, 1982).

Prior estimates of the unknowns of the inverse problem, which include P_S, P_L, the particle cycling parameters, $εS¯h$⁠, and $εL¯h$⁠, are represented in the vector x_o, and the uncertainties in these prior estimates are represented in the matrix C_o. In this section, we present the prior estimates of the unknowns and their respective uncertainties that are assumed in this study. This prior knowledge is summarized in Table 2.

Consider first the prior estimates of P_S and P_L, as well as their uncertainties. The cruise averages of P_S and P_L at Station P show maximum values in the mixed layer of about 2.5 mmol m⁻³ and 0.12 mmol m⁻³, respectively (Figure 4). Below the mixed layer, the POC concentrations in both size fractions generally decrease with depth. These observations were used as prior estimates of P_S and P_L in x_o, and the squares of their standard errors were taken as the squares of their uncertainties in C_o.

Prior estimates and uncertainties of particle cycling parameters were also collected in x_o and C_o, respectively. Prior estimates of particle production parameters (⁠$P˙S,ML$ and L_P) and remineralization rate constant for P_L (⁠$β−1,L$⁠) were obtained from data collected during EXPORTSNP. The prior estimate of $P˙S,ML$ was calculated from the mean of ¹⁴C incubation-based NPP measurements near the base of the mixed layer (⁠$28≤z≤35 m$⁠, $N=10$⁠), and its uncertainty was set equal to the standard error of the mean. The prior estimate of L_P was derived from the vertical distribution of NPP measurements at and below the base of the mixed layer (⁠$z≥28 m$⁠, $N=24$⁠). The reciprocal of the slope of the ordinary least-squares fit of ln(⁠$P˙S$/$P˙S,ML$⁠) versus z was used as our prior estimate of L_P, and the uncertainty in this estimate was set equal to the standard error of the slope. Note that ¹⁴C incorporation rates were measured on particles collected on filters with a smaller pore size (0.2 µm) than that of filters used to collect POC in the SSF (approximately 1 µm). Accordingly, our prior estimate of $P˙S$ may be an overestimate of the production of P_S as sampled by LVISF, and the associated error is assumed to be represented in the corresponding uncertainty σ_o. Finally, our prior estimate of the remineralization rate constant for P_L (⁠$β−1,L$⁠) was obtained from in situ particle incubation data provided by RESPIRE traps deployed during EXPORTSNP (Collins et al., 2015; Santoro et al., 2020).

Prior information about DVM parameters was derived as follows. The zooplankton ingestion rate constant (β₃) was estimated from cruise-averaged depth-integrated NPP measurements from EXPORTSNP and an empirical relationship between primary production and biomass-specific ingestion rate (Calbet, 2001; their Figure 2). Zooplankton are commonly assumed to use roughly 70% of their food for growth and metabolism (Steinberg and Landry, 2017). The fraction that remains (i.e., 30%) was used as a prior estimate for our egestion fraction (α). For convenience, we assumed that the maximum egestion depth, z_m, is equal to 500 m, the lower boundary of our domain. While there is evidence of zooplankton biomass below 500 m at EXPORTSNP (Omand et al., 2021), we assumed that egestion by migrating zooplankton is confined within the UMZ.

Prior estimates and uncertainties for all other parameters (⁠$β2′$⁠, $β−2$⁠, $β−1,S$⁠, w_S, w_L) were also obtained from the literature. Given the large uncertainties in particle (dis)aggregation rates, we conducted two inversions, and thus report two sets of parameter estimates, based on two different pairs of prior estimates of $β2′$ and $β−2$ from two different studies. One pair of estimates comes from Murnane (1994), who published estimates of (dis)aggregation rate constants with relatively low precision from observations in the deep water column (1000–3800 m) at Station P (i.e., same location as our study, but from a deeper part of the water column). The other pair comes from Murnane et al. (1996), who reported estimates of (dis)aggregation rate constants with higher precision from data collected in shallower waters (150–300 m) during the North Atlantic Bloom Experiment, or NABE (i.e., from a depth range included in our model, but at a different location). Prior estimates and uncertainties of $β−2$ were obtained directly from these studies. A prior estimate of $β2′$ for NABE was obtained by dividing the estimate of the first-order aggregation rate constant (β₂) by the estimate of small POC (Murnane et al., 1996). As Murnane (1994) did not estimate P_S at Station P, we obtained an estimate of $β2′$ for Station P by dividing his estimate of β₂ by the mean of P_S data from the upper 300 m at Station P reported by Bishop et al. (1999). Use of P_S data from Bishop et al. rather than from EXPORTSNP was preferred to avoid error covariances in C_o, which is taken as diagonal in our study. Uncertainties for both prior estimates of $β2′$ were propagated from their respective estimates and uncertainties of β₂ and P_S, neglecting error covariances (Bevington and Robinson, 2003).

The prior estimate for the sinking speed of LSF POC (w_L) was obtained from a previous estimate at Station P reported by Murnane et al. (1990). To our knowledge, previous estimates of the sinking speed and remineralization of SSF POC (w_S and $β−1,S$⁠, respectively) are not available at Station P. Here, we used a prior estimate of w_S that is consistent with the bulk particle sinking speed inferred from radionuclide and particle data collected at 11 open-ocean stations in the subtropical North Atlantic (Lerner et al., 2017; their Table 5), and a prior estimate of $β−1,S$ that is consistent with an estimate of particle remineralization in the upper 100 m inferred from a particle cycling model applied to data from the equatorial Pacific and Station P (Clegg et al., 1991).

Note that model parameters that represent DVM (β₃, α, z_m), particle sinking speeds (w_S and w_L), and remineralization (⁠$β−1,S$ and $β−1,L$⁠) have uncertainties that are not well constrained by available data. We assumed a relative error of 50% for all of these parameters (Table 2).

Finally, consider prior estimates of the model errors, $εS¯h$ and $εL¯h$⁠, and their uncertainties. The error (co)variances of oceanographic models are notoriously difficult to estimate, and there are no clear guidelines for this task. Here the prior estimates of $εS¯h$ and $εL¯h$⁠, which are represented in x_o, were set equal to 0 (i.e., Equations 6.1 and 6.2 are assumed to contain no error a priori). Their uncertainties in C_o were set proportional to the observed cruise-averaged particle production at the base of the mixed layer times the mixed layer thickness, with a proportionality constant γ. This treatment is motivated by the fact that particle production is probably the term that is best constrained (relative to other terms) in Equations 1.1 and 1.2. The case γ = 0.5, for example, corresponds to a situation where the errors in the POC equations have a standard deviation equal to half our prior estimate of $P˙S,ML$ times $zML$⁠. Because particle production is likely a first-order (i.e., leading) term in the P_S budget, at least in the EZ (were this not true, maxima in P_S typically observed in ocean surface waters would be difficult to explain), this specific case appears conservative. As such and unless stated otherwise, γ was set to 0.5 in this study.

In this section we present estimates of POC concentrations, particle cycling parameters, and model residuals from two inversions. One inversion (referred to as NA inversion below) incorporates prior estimates of aggregation (⁠$β2′$⁠) and disaggregation (⁠$β−2$⁠) rate constants for the upper water column during the North Atlantic Bloom Experiment (Murnane et al., 1996). The other inversion (referred to as SP inversion) uses prior estimates of the same parameters for the deep water column at Station P (Murnane, 1994) (Section 4.4).

For both inversions, we verified that the solution (elements of $x^k+1$⁠) is consistent with prior values (elements of x_o) given estimates of their uncertainties (square root of diagonal elements of C_o) (Figures S1 and S2). Note that the choice of γ and relative error do not appear to meaningfully affect the results of the inversions (Text S1, Figures S3 and S4). The accuracy of the inverse method is evaluated from twin experiments (Text S2, Figures S5–S12).

Figure 5 compares posterior estimates of [POC] from both the NA and SP inversions with [POC] data obtained from LVISF. Both inversions capture the original data within the uncertainties of the posterior estimates, and the estimates between the two inversions are similar. The precision in the [POC] estimates is not dramatically improved by data inversion.

Figure 6 compares posterior estimates for the non-uniform rate constants from the NA and SP inversions with their respective prior estimates (see also Table 3). For the NA inversion, rate constant estimates for aggregation (⁠$β2′$⁠), disaggregation (⁠$β−2$⁠), and remineralization of large particles (⁠$β−1,L$⁠) do not differ noticeably from their prior estimates and do not vary much with depth. Particle sinking speeds (w_S and w_L) and the rate constant for small particle remineralization (⁠$β−1,S$⁠) have minima in the mixed layer, but closely resemble the prior estimates with no considerable improvement in precision throughout the rest of the water column.

We see similar results for the rate parameter estimates from the SP inversion, with the exception of the disaggregation rate constant. The posterior estimates of $β−2$ display relatively large vertical variations, with a subsurface maximum between 30 and 50 m. Furthermore, the precision of the posterior estimates of $β−2$ is greatly improved compared to its prior estimate. Both of these results can be explained by the fact that the SP prior estimate has a much higher relative uncertainty (2500%) than the NABE prior estimate (11%; Table 2). The value of $β−2$ is allowed to vary more with depth in the SP inversion where it is not bound by a narrow prior uncertainty. The large reduction in the uncertainty of $β−2$ indicates that the [POC] data from EXPORTSNP lead to a great improvement in our knowledge of particle disaggregation rates at Station P if a large uncertainty in the prior estimate is assumed.

For the uniform particle cycling parameters, the NA and SP inversions produce similar posterior estimates with no obvious improvement in precision compared to their prior estimates (Table 4 and Figure 7). Furthermore, the posterior estimates tend to be very similar to the prior estimates, with the largest relative change appearing in the zooplankton ingestion rate constant (β₃).

The posterior estimates of the residual terms in the layer-integrated POC balance equations (⁠$εS¯h$ and $εL¯h$ in Equations 6.1 and 6.2, respectively) are shown in Figure 8. The profiles of the residuals for both the NA and SP inversions are similar, with maxima between 0–30 m. For both inversions and size fractions, the residuals are positive throughout the entire water column (0 and 500 m). This result suggests that source terms are underestimated, missing, and/or loss terms are overestimated in the balance equations for both small and large POC. Posterior errors in the residuals tend to be considerably smaller than their prior errors, indicating a reduction in the uncertainty in model errors.

The posterior estimates of rate constants for remineralization of large particles (⁠$β−1,L$⁠) and aggregation of small particles (⁠$β2′$⁠) are very similar to their prior estimates in both the NA and SP inversions (Tables 2 and 3). The same is generally true for the uniform parameters (Figure 7). The rate constants that show small adjustments relative to their prior estimates indicate that they are consistent with the data being considered in the inversions ([POC] in SSF and LSF), given the assumptions made about error (co)variances. Below, we focus our discussion on rate constants that exhibit noticeable adjustments in comparison to their prior estimates.

Posterior estimates of particle cycling parameters are sensitive to assumptions about prior estimates. This result is particularly highlighted by the disparate estimates of the disaggregation rate constant (⁠$β−2$⁠) between the NA and SP inversions (Figure 6). The NA inversion, which relies on a relatively precise prior estimate of $β−2$ (0.43 ± 0.05 d⁻¹), results in posterior estimates of $β−2$ that are uniform with depth and identical to the prior value (Tables 2 and 3). In contrast, the SP inversion, which relies on a prior estimate of $β−2$ that is much less precise (1.1 ± 27.4 d⁻¹), allows the posterior estimate of this parameter to vary with depth and thus depart from its prior value. In the latter inversion, $β−2$ is inferred to increase slightly from the mixed layer to the layer below (30–50 m), and subsequently to decrease to 500 m.

Note that the use of different prior $β2′$ and $β−2$ estimates (from NA or SP) leads to noticeable differences in the inferred magnitude of particle aggregation and disaggregation rates at Station P (Figure 6 and Table 3). This result suggests that similar POC profiles consistent with POC measurements could be produced from distinct magnitudes of these processes. It stresses that future work is needed for better constraining the magnitude of particle coagulation and particle break-up, as well as the major controls of these processes, in the upper water column.

The remineralization rate constant is thought to decrease as temperature and dissolved oxygen availability decrease (Cram et al., 2018). For both inversions, however, the posterior estimates of the rate constant for remineralization of small POC (⁠$β−1,S$⁠) increase slightly between the mixed layer and the layer below, and are relatively uniform below the mixed layer (Figure 6) despite (i) the 10°C decrease in temperature between the mixed layer and 300 m, and (ii) oxygen concentrations that approach typical oxygen half-saturation constants (0–30 µM) for aerobic microbial metabolism (e.g., Cram et al., 2018) by 500 m (Figure 1). Our results suggest that either $β−1,S$ is insensitive to both temperature and oxygen, or that [POC] data do not provide much information about remineralization rates below 50 m. Alternatively, and more speculatively, processes other than temperature and oxygen concentration may influence remineralization rates between 50 and 500 m at Station P. Mechanistic interpretations aside, both inversions produce estimates of small and large POC that underestimate the data in the mixed layer (Figure 5). In attempting to produce [POC] estimates that fit the mixed layer data, the inversion leads to estimates of the remineralization rate constant for small POC that are relatively small.

In both inversions, the posterior estimates of sinking speeds of small and large particles (w_S and w_L, respectively) also exhibit minima in the mixed layer with little variation from the prior estimate below. Some phytoplankton groups may be able to reduce their settling speeds via buoyancy control (e.g., Arrieta et al., 2015; Gemmell et al., 2016; Borgnino et al., 2019) which could contribute to these mixed-layer minima. As is the case for $β−1,S$⁠, the inversion leads to relatively small estimates of sinking speeds in the mixed layer to better fit the [POC] data.

We compared our posterior estimates of rate constants for remineralization of small particles (⁠$β−1,S$⁠), disaggregation (⁠$β−2$⁠), and aggregation (⁠$β2′$⁠) to previous estimates for the upper ocean (Figure 9). For this comparison, estimates of a first-order rate constant for aggregation (β₂) were derived by multiplying our estimates of the second-order rate constant (⁠$β2′$⁠) by the average P_S concentration in each layer, with error propagated appropriately. Figure 9 illustrates the large range of estimates reported for the rate constants shown and demonstrates that our estimates are generally within the range of previously reported estimates; various factors such as differences in estimation methods and in oceanic environments complicate the direct comparison of our estimates to those from previous studies (see Section 8).

In this section, we discuss POC fluxes calculated for the small and large size fractions by combining posterior estimates of P_S, P_L, and the particle cycling parameters. The uncertainties in the POC fluxes were calculated by propagating the uncertainties in the posterior estimates of POC concentrations and cycling parameters, considering error covariances (e.g., Bevington and Robinson, 2003).

The POC sinking fluxes are shown in Figure 10 (left column). The sinking flux of P_S (w_SP_S) tends to be greater than that of P_L (w_LP_L) in the EZ, but the fluxes in the two size fractions converge to similar values in the UMZ. This result holds for both inversions, and is a reflection of the fact that while w_L is about 10 times greater than w_S throughout the water column (Table 3), P_S/P_L ratios are greater than 10, which leads to larger sinking fluxes from small particles. While previous models with two particle size classes assumed that small particles do not sink (e.g., Nozaki et al., 1987; Clegg and Whitfield, 1990; Murnane et al., 1990; Clegg et al., 1991; Murnane, 1994; Murnane et al., 1994; Murnane et al., 1996), our results suggest that if allowed to sink, small particles may contribute significantly to carbon export. An important contribution of small particles to total POC flux was indeed observed from EXPORTSNP sediment trap data (Durkin et al., 2021; Estapa et al., 2021).

We compared our POC flux estimates to those obtained from other concurrent studies conducted as part of EXPORTSNP. Figure 10 (right column) shows our estimates of the sinking flux of total POC (w_SP_S + w_LP_L), along with estimates of POC export derived from ²³⁴Th disequilibrium and from sediment trap data that were also collected during the cruise (Buesseler et al., 2020; Estapa et al., 2021). For both inversions, our estimates of total POC sinking flux are consistent (within uncertainties) with estimates of POC export based on measurements of ²³⁴Th/²³⁸U disequilibria in the upper water column using a steady state model and POC/²³⁴Th ratios of particles sampled by LVISF, as well as with POC fluxes measured by sediment traps. Discussions of the ²³⁴Th- and sediment-trap based estimates of POC flux at Station P can be found in Buesseler et al. (2020) and Estapa et al. (2021), respectively.

Figure 11A shows the sinking flux divergences for small and large POC (⁠$[wSbPSb−wStPSt]/h$ and $[wLbPLb−wLtPLt]/h$⁠, respectively, where h is the thickness of the considered layer). Particle sinking can supply POC to, or remove POC from, a given layer. A negative flux divergence (i.e., convergence), in which more POC would be sinking through the top of the layer than sinking through the bottom of the layer, would indicate a supply, while a positive flux divergence would indicate removal, in which more POC would be sinking out the bottom of the layer than sinking in from the top. A positive divergence of POC flux for a given size fraction and a given layer must be balanced by, or imply, a net flux of POC into this size fraction and in this layer. This flux can be explicitly represented in our model, or it could be implicitly captured in the residuals. In the mixed layer, the sinking flux divergence is positive for both size fractions, reflecting POC export at the base of the layer, so that particle sinking tends to reduce POC concentrations there. As P_S and P_L attenuate sharply with depth in the euphotic zone, sinking flux divergences become negative, indicating a supply of POC as more particles sink in from above than are removed from below. Reduced attenuation of [POC] in the UMZ implies that the sinking flux divergences for both size fractions are very small and almost vanish at 200 m.

Comparing the magnitude of remineralization rates with the magnitude of (dis)aggregation rates in different parts of the water column is instructive. Figure 11B compares the fluxes of remineralization and aggregation of P_S, both processes that remove POC from the small size fraction. For both inversions, the remineralization flux is greater than the aggregation flux throughout the water column, indicating that more P_S is lost to the “dissolved” (< 1 µm) pool than to the large particulate size fraction (> 51 µm). This difference is more pronounced for the NA inversion than the SP inversion, and aggregation flux estimates from the SP inversion have larger uncertainties presumably due to the larger uncertainty in the prior estimate of $β2′$⁠. In comparison, Figure 11C shows the fluxes of remineralization and disaggregation of P_L, both of which remove POC from the large size fraction. There is a clear difference in the precision of the posterior estimates of disaggregation rates between the two inversions because of the large differences in the errors of their prior estimates (Table 2). Nonetheless, both inversions suggest that disaggregation fluxes consistently remove more P_L from the water column than does remineralization. This result suggests that particle disaggregation is an important process to consider in ocean biogeochemical models (e.g., Aumont et al., 2015; Niemeyer et al., 2019).

The rate of P_S production (⁠$P˙S$⁠) at different depths is estimated from the posterior estimates of the particle production in the mixed layer (⁠$P˙S,ML$⁠) and of the vertical scale (L_P) using Equation 2 (Figure 11D). Our posterior estimates of particle production are consistent with NPP estimates based on ¹⁴C incubation experiments, a result that could be expected from the fact that our prior estimates of production parameters are based on NPP estimates. $P˙S$ dominates all P_S fluxes at shallow depths in the EZ, becomes comparable to or lower than P_S removal fluxes near the bottom of the EZ, and vanishes in the UMZ as expected from low light levels. Albeit small, it is still non-zero there due to the exponential formulation (Equation 2), indicating that particle production in the UMZ may be overestimated.

Finally, we estimated volumetric fluxes related to DVM, which include zooplankton ingestion of P_S in the EZ and zooplankton egestion of P_L in the UMZ (Figure 12). For both inversions, the ingestion flux decreases with depth in the EZ and represents a smaller loss of P_S than remineralization (the largest loss term for P_S in both zones). In the UMZ, the egestion flux exhibits a maximum between 200 and 330 m, which reflects our formulation of this flux (Equation 4). This subsurface maximum approximately coincides with the maximum in daytime zooplankton biomass that was observed between 300 and 400 m during EXPORTSNP (Omand et al., 2021). Interestingly, the egestion flux supplies more P_L below 150 m than does particle settling. Some models suggest that DVM may provide a smaller flux of organic carbon into the UMZ compared to the passive sinking of particles on a global scale (Aumont et al., 2018). Our results, however, suggest that DVM may have been the dominant source of large particles in the UMZ during EXPORTSNP.

In this section, we discuss the budgets of POC in the EZ and UMZ estimated by summing, over the thickness of each zone, the fluxes in the POC balance equations (6.1 and 6.2) from the NA and SP inversions. The uncertainty of each term in the budgets was calculated by propagating the uncertainties in the fluxes from different layers, accounting for error covariances. The estimated budgets are depicted in Figures 13 and 14, with numerical values listed in Table 5. We highlight the most notable features of these budgets here.

As suggested by the flux profiles shown in the previous section, remineralization is the largest loss term for P_S throughout both the EZ and the UMZ (Figures 13 and 14). In the EZ, DVM results in a greater loss of P_S than does aggregation. Note that in our model, “aggregation” represents processes (biological and physical) by which small particles are packaged into large particles within the same layer, whereas “DVM” represents a biological process whereby small particles in the EZ are packaged and egested as large particles in the UMZ. In the UMZ, disaggregation is the largest loss of P_L and the largest source of P_S.

Note that in both inversions and in both zones, the sum of the residuals in the POC budgets are of the same order of magnitude as that of (and in some instances, greater than) other fluxes that are represented explicitly in our model. This result indicates that temporal variability and/or processes not explicitly included in Equations 1.1 and 1.2 may have influenced the POC concentrations observed during EXPORTSNP, and/or that some processes that are explicitly included may be misrepresented. The fact that the summed residuals are all positive suggest that Equations 1.1 and 1.2 miss sources of P_S and P_L in both the EZ and the UMZ. One plausible explanation for the apparent missing source is under-sampling of submicron POC by LVISF near the surface (J. Graff, personal communication). Under-sampling of submicron POC may, for example, force a higher remineralization rate to account for the missing stock of P_S in the EZ, because the prior estimates for particle production were based on measurements using a smaller pore-sized filter and thus included production of small particles not captured in P_S (Section 4.4). Furthermore, particle breakage during pumping or inefficient rinsing of filters from which P_L was measured may have led to under-collection of P_L and contributed to the residuals in P_L (Roca-Martí et al., 2021).

Factors other than particle under-collection by LVISF could contribute to the estimated residuals in the POC budgets. One process that is not explicitly represented in our model and that could contribute to the estimated residuals is primary production of large POC. Our model assumes that primary production is a source of POC only in the SSF (Equations 1.1 and 1.2). This assumption is based on the expectation that photosynthesis leads to the production of primarily small particles (i.e., < 51 µm). However, some phytoplanktonic species may produce cells larger than 51 µm, thus representing a source of large POC that is not accounted for in our model. Takahashi (1986) found that although most diatoms sampled at Station P were small (< 20 µm), large species such as Rhizosolenia styliformis (which secrete frustules 500–1500 µm in length) contributed considerably to opal fluxes in certain months, including July and August. During EXPORTSNP, diatom cells with diameter up to 107 µm were observed in gel traps (Bodel et al., 2020), but measurements of size-fractionated NPP from ¹³C-based incubations suggest that primary production in the > 5 µm size fraction was only about 10% of total NPP (M. Meyer, personal communication). NPP in the > 51 µm size fraction (i.e., P_L) was less than that, which indicates that the inclusion of primary production of P_L in our model may be insufficient to close the P_L budgets.

There is also a possibility that the ¹⁴C-based estimates of NPP are underestimates: ¹³C-based estimates of NPP were roughly two times greater than the ¹⁴C-based estimates used as prior constraints in our inversions (M. Meyer, personal communication). Inversions with higher prior estimates of particle production in the mixed layer (not shown here) result in smaller P_S (P_L) residuals in the EZ (UMZ). Physical transport and temporal variability may also contribute to budget residuals because our model is 1D, neglects vertical transport by advection and mixing, and assumes steady state conditions. For example, a seven-year time series at Station P suggests that net community production tends to peak in the spring and decline slightly over the summer months (Fassbender et al., 2016). The consideration of production rates which are higher than those measured by ¹⁴C incubation during EXPORTSNP and which may have characterized the system in the weeks preceding the cruise may help to close the POC budgets.

Residence and turnover times of P_S and P_L can provide insight into the potential influence of temporal variability on P_S and P_L concentrations at Station P. We used our results from the NA and SP inversions to estimate these timescales in both the EZ and the UMZ. The residence time is the average time a carbon atom in the SSF (LSF) spends in a particular zone, and was calculated as the inventory of P_S (P_L) in that zone (mmol m⁻²; Table 6) divided by the sum of the fluxes in or out of the P_S (P_L) reservoir in that zone (mmol m⁻² d⁻¹; Table 5). In contrast, the turnover time of P_S (P_L) is the average time a carbon atom spends in the SSF (LSF) with respect to a particular process, and was calculated from the ratio of the inventory of P_S (P_L) in the EZ or UMZ to the sum of the fluxes for this process in the same zone (mmol m⁻² d⁻¹). P_S and P_L residence times and turnover times in each of the zones are listed in Tables 7 and 8, respectively.

Broadly, the turnover times (Table 8) reflect the relative importances of the fluxes in each of the budgets shown in Figures 13 and 14. For example, remineralization is the most important loss term for P_S in both the EZ and the UMZ, and could remove the entire inventory in 14 and 11 days, respectively. In comparison, grazing by vertically migrating zooplankton would take 25 days to remove all P_S in the EZ. In the UMZ, DVM would resupply the P_L reservoir in about 7 days, but the entire P_L reservoir can be removed by disaggregation in only 2–3 days.

Note that the turnover time of POC (in the SSF or LSF) in the EZ with respect to POC sinking flux divergence (SFD; Table 8) is equal to the ratio of POC inventory in the EZ to the sinking flux of POC at the base of the EZ. In other words, these turnover times can be interpreted as the residence time of POC (in the SSF or LSF) with respect to POC settling at the base of the EZ. Eppley et al. (1983) reported an empirical relationship between residence time of suspended POC in the EZ and total ¹⁴C production based on data from the Southern California Bight and the central North Pacific (their Figure 3). Applying this relationship to our cruise-averaged ¹⁴C NPP integrated over the EZ (165.8 mg m⁻² d⁻¹), we obtained a residence time of about 120 d, which is close to our estimates of turnover of P_S due to sinking at the base of the EZ (roughly 115 d). These timescales are one to two orders of magnitude larger than the corresponding residence times of P_S shown in Table 7, indicating that fluxes from processes other than sinking are important in governing the supply and removal of P_S in the EZ. The residence times of P_L in both the EZ and UMZ are very short (< 2 days, Table 7), indicating that this reservoir of POC renews itself very quickly.

In this study, we applied an inverse method for estimating particle cycling rates in the euphotic zone (EZ, 0–100 m) and upper mesopelagic zone (UMZ, 100–500 m) at Station P from size-fractionated [POC] data, which were collected in summer 2018 as part of EXPORTSNP. Two inversions were conducted with distinct prior estimates of the disaggregation and aggregation rate constants – the NA inversion considers relatively precise estimates from depths similar to those from which our samples were collected but from a different study site, while the SP inversion considers less precise estimates from the same study site as ours, but from the deep ocean. We find that the posterior estimates of particle cycling parameters are sensitive to the magnitude and precision of the prior estimates of the (dis)aggregation rate constants. The prior estimates for disaggregation and aggregation rate constants in the NA inversion, which are relatively precise (approximately 10% relative error), lead to posterior estimates of these rate constants that are identical to their prior values. The SP inversion uses prior estimates for (dis)aggregation rate constants that are not only greater in magnitude, but also much less precise. This inversion is not only able to greatly reduce the relative error in the posterior estimate of $β−2$ (from 2500% to 100–300%), but the lower precision in the prior estimate allows for greater adjustment in the posterior estimate, revealing a decrease in the disaggregation rate constant with depth at Station P. In this sense, the SP inversion is thus more successful in extracting information about particle cycling rate parameters from [POC] data, at the cost of slightly greater uncertainty in the model residuals.

Despite the sensitivity of posterior estimates to assumptions about prior values, we nonetheless find robust trends in particle cycling processes during EXPORTSNP: for both the NA and SP inversions, remineralization removes more POC from the small size fraction throughout the water column than does aggregation; in contrast, disaggregation is a more important sink of large POC than remineralization. Furthermore, our estimates of the sinking fluxes of total POC (> 1 µm) are consistent with estimates derived from sediment traps and measurements of ²³⁴Th/²³⁸U disequilibria in the upper water column paired with POC/²³⁴Th ratios in particles collected during EXPORTSNP. We estimate that POC export at the base of the EZ averaged 2.2 ± 0.8 mmol m⁻² d⁻¹ during the cruise, and that both small and large particles contributed considerably to the total export flux throughout the water column.

We constructed a budget of POC for each size fraction (1–51 µm and > 51 µm) and in each zone (EZ and UMZ). In both size fractions and in both the EZ and UMZ, the model residual is of the same order as that of the largest fluxes that are explicitly represented in the model. The residuals, which are all positive, indicate that POC sources are underestimated, and/or missing, and/or that loss terms are overestimated in our POC balance equations. Such residuals may be a consequence of (i) biogeochemical processes that are not considered in the model, (ii) unsteadiness, (iii) physical transport, and/or (iv) methodological issues such as the under-collection of picoplankton by filters and the underestimation of NPP. With these caveats in mind, remineralization represents the largest loss of P_S throughout the water column, and DVM is the largest source of P_L in the UMZ, after the residual term. For the NA inversion, the residence of time of POC in the EZ is estimated to be 8.0 ± 1.8 d and 1.5 ± 0.2 d for, respectively, the SSF and the LSF. For the SP inversion, these numbers amount to 6.7 ± 1.7 d and 0.8 ± 0.6 d, respectively.

The work presented in this paper is a first step in the estimation of particle cycling rates from measurements of the concentration of multiple chemical tracers in different size fractions. Besides the fact that this study includes only a single tracer (POC), a few other important limitations are worth noting. The first is that we used a one-dimensional steady state model, which thus neglects temporal variability and lateral transport processes. Vertical transport by advection and diffusion is also omitted. Second, the [POC] data that are inverted include only two observations from the mixed layer for each size fraction. This small sample size lowers our confidence in our estimates of POC cycling rates in the mixed layer. Future work should incorporate more measurements in this portion of the water column. Finally, constraining the error in the POC balance equations is difficult; here it is assumed to be proportional to the observed cruise-averaged NPP at the base of the mixed layer. We applied a conservative value for this proportionality constant. Confidence in this method was obtained through sensitivity tests, which show that the estimates of most particle cycling parameters do not change appreciably with γ (Figures S3 and S4). Alternative approaches to constrain model errors, however, would be desirable.

While an inversion that uses [POC] data seems to provide reasonable estimates of particle cycling rates, the addition of other chemical tracer data generated from EXPORTS or other field programs may help to further constrain the accuracy and precision of these rates. Examples of such tracers include total and dissolved ²³⁴Th (Nozaki et al., 1987; Clegg and Whitfield, 1990; Murnane et al., 1990; Clegg et al., 1991; Murnane, 1994; Murnane et al., 1994; Murnane et al., 1996; Lerner et al., 2017; Black et al., 2018), as well as particulate lithogenic metals, which may provide important constraints on particle aggregation and disaggregation (Ohnemus and Lam, 2015).

Tracer data with higher vertical resolution and from different biogeochemical regimes may also improve inference about the vertical and geographical distribution of particle cycling parameters and particle fluxes in the upper ocean. For example, the same method could be applied to analyze data from the GEOTRACES GP15 transect, which spanned from the Aleutian shelf to Tahiti. Transect GP15 includes profiles of POC and trace metal concentrations with similar vertical resolution measured in the eastern subarctic Pacific, the transition zone chlorophyll front (Polovina et al., 2001), the North and South Pacific subtropical gyres, and the equatorial Pacific upwelling region. Analysis of these data would help elucidate how differences in environmental variables such as nutrient concentration and plankton community structure may affect particle cycling rates, and disentangle the effect of method versus environment in the large range of particle cycling rate estimates that exist in the literature (Figure 9).

In spite of the limitations mentioned above, the estimates of particle cycling rates reported in this work have the benefit that they are based on [POC] measurements in different size fractions obtained from multiple casts, they incorporate prior knowledge about particle cycling gathered from previous studies, and they consider both data and model errors. Our estimates of particle cycling rates, which include POC production, remineralization, (dis)aggregation, settling, and ingestion and egestion due to DVM, are internally consistent (within relatively large residuals). Particle processes are estimated collectively, not separately. Accordingly, we hope that the results of our study, which considers a relatively large number of POC cycling processes, will provide context for the interpretation of data generated from process-specific studies undertaken during EXPORTSNP and help to place the results from these studies into a larger biogeochemical perspective.

The following datasets were used as model input:

POC concentrations: Buesseler, KO, Benitez-Nelson, CR, Resplandy, L, Roca-Martí, M, Pike, S, Umhau, B, Horner, T, Masqué, P. 2019. EXPORTS-EXPORTSNP_IN_SITU_PUMPS_SURVEY, SeaBASS repository. DOI: 10.5067/SeaBASS/EXPORTS/DATA001; https://seabass.gsfc.nasa.gov/archive/WHOI/buesseler/EXPORTS/EXPORTSNP/archive/.

Net primary production: Fox, J, Behrenfeld, MJ, Halsey, KH. 2019. EXPORTS-EXPORTSNP_NPP_14C_process_R0, SeaBASS repository. DOI: 10.5067/SeaBASS/EXPORTS/DATA001; https://seabass.gsfc.nasa.gov/archive/OSU/behrenfeld/EXPORTS/EXPORTSNP/archive/.

The code and input data used to generate the results within this paper can be found at https://github.com/amaralvin7/pyrite.

The supplemental files for this article can be found as follows:

Figure S1. Probability density functions of residuals for NA inversion.

Figure S2. Probability density functions of residuals for SP inversion.

Text S1. Sensitivity to prior uncertainties.

Figure S3. Sensitivity of POC flux estimates in the euphotic zone to prior uncertainties.

Figure S4. Sensitivity of POC flux estimates in the upper mesopelagic zone to prior uncertainties.

Text S2. Accuracy of model estimates.

Figure S5. Estimates of [POC] in the small (P_S) and large (P_L) size fractions from a twin experiment (NA inversion).

Figure S6. Estimates of [POC] in the small (P_S) and large (P_L) size fractions from a twin experiment (SP inversion).

Figure S7. Estimates of non-uniform particle cycling parameters from a twin experiment (NA inversion).

Figure S8. Estimates of non-uniform particle cycling parameters from a twin experiment (SP inversion).

Figure S9. Estimates of uniform particle cycling parameters from a twin experiment (NA inversion).

Figure S10. Estimates of uniform particle cycling parameters from a twin experiment (SP inversion).

Figure S11. Estimates of residuals in the POC balances from a twin experiment (NA inversion).

Figure S12. Estimates of residuals in the POC balances from a twin experiment (SP inversion).

We would like to thank the crews of the R.V. Roger Revelle and R.V. Sally Ride during the 2018 North Pacific EXPORTS cruise for conducting ship operations. In particular, we are grateful for the work done by Steve Pike, Claudia Benitez-Nelson, and the rest of the “pump team” in collecting samples. We are grateful to Alyson Santoro for contributing particle remineralization data from the RESPIRE traps and Hilary Close for contributing [POC] data from the mixed layer. We would also like to thank Margaret Estapa and Amy Maas for sharing their insights regarding sediment traps and zooplankton migration, and Ken Buesseler for providing comments on an early version of this manuscript. Finally, we thank the three anonymous reviewers whose comments help guide the revision of this manuscript.

This study was supported by the National Aeronautics and Space Administration (NASA) program award 80NSSC17K0555, NSF-OCE 1829614 to PJL, and NSF-OCE-1829790 to OM. VJA was supported by the NSF Graduate Research Fellowship Program and the UC Eugene Cota-Robles Fellowship. MRM was supported by the Ocean Frontier Institute International Postdoctoral Fellowship Program and the Woods Hole Oceanographic Institution’s Ocean Twilight Zone study.

The authors have no competing interests, as defined by Elementa, that might be perceived to influence the research presented in this manuscript.

Contributed to conception and design: VJA, OM, PJL.

Contributed to acquisition of data: VJA, MRM, JF, NBN.

Contributed to analysis and interpretation of data: VJA, OM, PJL, MRM.

Drafted and/or revised the article: VJA, OM, PJL, MRM, JF, NBN.

Approved the submitted version for publication: VJA, OM, PJL, MRM, JF, NBN.