Disentangling Carbon Concentration Changes Along Pathways of North Atlantic Subtropical Mode Water

North Atlantic subtropical mode water (NASTMW) serves as a major conduit for dissolved carbon to penetrate into the ocean interior by its wintertime outcropping events. Prior research on NASTMW has concentrated on its physical formation and destruction, as well as Lagrangian pathways and timescales of water into and out of NASTMW. In this study, we examine how dissolved inorganic carbon (DIC) concentrations are modified along Lagrangian pathways of NASTMW on subannual timescales. We introduce Lagrangian parcels into a physical‐biogeochemical model and release these parcels annually over two decades. For different pathways into, out of, and within NASTMW, we calculate changes in DIC concentrations along the path (ΔDIC), distinguishing contributions from vertical mixing and biogeochemical processes. The strongest ΔDIC is during subduction of water parcels ( + 101 μmol L − 1 in 1 year), followed by transport out of NASTMW due to increases in density in water parcels ( + 10 μmol L − 1 ). While the mean ΔDIC for parcels that persist within NASTMW in 1 year is relatively small at + 6 μmol L − 1 , this masks underlying dynamics: individual parcels undergo interspersed DIC depletion and enrichment, spanning several timescales and magnitudes. Most DIC enrichment and depletion regimes span timescales of weeks, related to phytoplankton blooms. However, mixing and biogeochemical processes often oppose one another at short timescales, so the largest net DIC changes occur at timescales of more than 30 days. Our new Lagrangian approach complements bulk Eulerian approaches, which average out this underlying complexity, and is relevant to other biogeochemical studies, for example, on marine carbon dioxide removal.


Introduction
The ocean is an integral component of the natural carbon cycle, as well as a large sink for anthropogenic carbon emissions.Since 1850, it has taken up 26% of anthropogenic CO 2 from the atmosphere (Friedlingstein et al., 2022).To understand the ocean carbon sink, now and in the future, it is important to understand how the ocean moves carbon from the upper ocean mixed layer through the permanent thermocline, from where it can be further sequestered on timescales of years, decades, or centuries.A major conduit through which this ocean surface-interior exchange occurs is the North Atlantic subtropical mode water (NASTMW).It links the interior to the surface on an annual basis during winter convective events and is responsible for 20% of the carbon uptake by the solubility pump in the 14-50°N latitude band of the North Atlantic (Bates, 2012).
NASTMW, also referred to as Eighteen Degree Water, is a classical example of mode water (Hanawa & Talley, 2001), featuring a thick vertical layer characterized by near-homogeneous properties including temperature, salinity, and oxygen concentration.It has a typical thickness of 200-300 m, located at 300 m depth (Deremble & Dewar, 2013;Worthington, 1958).It is formed during winter, when surface buoyancy loss leads to convection events that deepen the mixed layer in the Sargasso Sea, and through cross-frontal mixing in the southern flanks of the Gulf Stream (Davis et al., 2013;Joyce et al., 2013).Spring stratification caps off NASTMW again, causing it to act as an interannual buffer of wintertime atmospheric anomalies of temperature and carbon (Bates et al., 2002).Gyre circulation and eddy-induced advection allow NASTMW to spread horizontally southwards, causing it to occupy an area much larger than its formation location (Gary et al., 2014).This makes NASTMW also a key regulator of temperature (Sugimoto et al., 2017), organic carbon (Sugimoto et al., 2017), and nutrients (Palter et al., 2005) in the interior of the subtropical gyre.Subsequent destruction of NASTMW occurs primarily through vertical mixing at the top of the layer, but also through diapycnal mixing and alongisopycnal stirring (Billheimer & Talley, 2016).
Current understanding of the role of NASTMW in oceanic carbon uptake is either based on sparse observations (Bates, 2012;Bates et al., 2002;Billheimer et al., 2021) or is inferred from insights into physical mechanisms such as its formation, ventilation and pathways (Davis et al., 2013;Gan et al., 2023;Gary et al., 2014;Kwon et al., 2015;Li et al., 2022).However, a process-based view of how dissolved inorganic carbon (DIC) is transported along pathways from the ocean surface through NASTMW into the ocean interior is lacking.
We investigate how DIC concentrations change along pathways of NASTMW during its formation, persistence, ventilation, and physical export to higher-density water masses to better understand which processes alter carbon concentrations along this conduit between the atmosphere and the ocean interior.To do so, we trace virtual parcels of water along pathways into, out of, and within NASTMW using a coupled physical-biogeochemical, eddypermitting ocean model.Along the pathways of these flow-following Lagrangian parcels, we disentangle the influence of different vertical mixing and biogeochemical processes on the local DIC concentration.Specifically, we split biogeochemical processes into soft-tissue and carbonate components.We then quantify the timescales and magnitudes of DIC depletion and enrichment regimes, defined between local minima and maxima in DIC anomaly time series.Rather than only looking at the bulk change in DIC concentrations along each pathway, we also consider how these changes are distributed between processes and pathways, as well as in time and between a range of timescales.This allows us to better understand the complexity by which vertical mixing and biogeochemical processes affect the DIC content of NASTMW at different moments and timescales.
We focus on timescales of the order of years and less, as observations and modeling studies show that most NASTMW parcels have residence times shorter than a year (Fratantoni et al., 2013;Gary et al., 2014).We consider parcels that subduct into NASTMW, ventilating parcels, persisting NASTMW parcels, and parcels that are exported due to increases in density.Parcels in this last class are relevant candidates for longer sequestration on timescales of years to decades (the timescale of the gyre interior; Levine et al., 2011).

Ocean Model Data
To compute Lagrangian parcel trajectories and along-trajectory DIC changes, we use gridded ocean physics and biogeochemistry output data from a global hindcast ocean model at 1/4°horizontal resolution, comprised of the FREEGLORYS2V4 physics and FREEBIORYS2V4 biogeochemistry products developed by Mercator Ocean International (MOi) for the Copernicus Marine Service.These coupled products have an eddy-permitting resolution, resolving part of the mesoscale eddy regime, which plays a role in mode water formation (Davis et al., 2013;Xu et al., 2014).Ideally, we would use a model set that resolves the full mesoscale, to better represent NASTMW (Gan et al., 2023), or even parts of the submesoscale, as these play important roles in biogeochemical cycles (Lévy et al., 2024).However, an eddy-permitting resolution accommodates the high computational and storage demands from the physical-biogeochemical run (similar to Atkins et al., 2022), and can provide us with a mesoscale process-level understanding.The horizontal resolution of 1/4°is representative of that used in state-ofthe-art earth system models (Hewitt et al., 2020), and thus also bears relevance to their dynamics.We use versions of both data products on their native Arakawa C-grid, allowing for more precise Lagrangian trajectory computations (Delandmeter & van Sebille, 2019).
Ocean physics obeys conservation of mass, momentum, and biogeochemical budgets and the hindcast ocean model does not include data assimilation.We use a time series between 1995 and 2017, excluding earlier spin-up years.The length of the time series allows us to take into account interannual variability.Since the model is not constrained by observations after its initialization, it may exhibit drift from observed conditions over time.Thus, we use the model in the context of process understanding rather than precise reproduction of observational data.
Ocean physics in FREEGLORYS2V4 are simulated with NEMO version 3.1 (Madec et al., 2013) with the ORCA025 configuration, having a 22 km horizontal resolution at Cape Hatteras and 75 vertical levels (Bernard et al., 2006).Vertical mixing is parameterized using an adaptation of the turbulent closure model by Blanke and Delecluse (1993).Physics are initialized from the EN4 data product (Good et al., 2013) and atmospheric forcings come from 3-hourly ERA-interim reanalysis products from ECMWF (Dee et al., 2011).FREEGLORYS2V4 has an assimilated counterpart, GLORYS2V4, which is extensively described in Garric and Parent (2017).A comparison of the model with observations is found in Text S1 in Supporting Information S1.
Biogeochemistry in FREEBIORYS2V4 is modeled using the intermediate complexity PISCES-v2 model (Aumont et al., 2015).PISCES simulates the carbon cycle, carbonate chemistry, main nutrients (P, N, Fe, and Si), and the lower trophic levels of marine ecosystems (phytoplankton, microzooplankton, and mesozooplankton) using 24 prognostic variables in total.These tracers are advected and vertically mixed using the hydrodynamics from FREEGLORYS2V4, without horizontal diffusive mixing.Nitrate, phosphate, oxygen, and silicate are initialized using data from the World Ocean Atlas (National Oceanographic Data Center (U.S.) Ocean Climate Laboratory et al., 2002) and DIC and alkalinity are initialized with the GLODAP climatology (Key et al., 2004).The model includes atmospheric deposition and riverine input of Fe, Si, N and P, as well as Fe input from sediment.Although PISCES-v2 imposes a stoichiometric ratio of C:N:P = 122:16:1, cycles of phosphorus and nitrogen are different due to nitrogen fixation, denitrification, and differences in atmospheric deposition of N and P. Atmospheric pCO 2 is prescribed at the air-sea interface, computed from monthly global CO 2 mole fractions (Lan et al., 2023).A biogeochemical model verification at the global scale is found Perruche et al. (2019).
Near the Bermuda Atlantic Time-series Study (BATS) site (32°10′N, 64°10′E) FREEBIORYS2V4 exhibits a trend of increasing salinity-normalized DIC at 10 m depth, by +0.2 μmol/L/year (Figure S12 in Supporting Information S1).This is roughly a factor 5 smaller than the observed salinity-normalized DIC trend of +1.08 ± 0.05 μmol/kg/year (Bates et al., 2012).Therefore, the input of the biogeochemical model does not accurately represent the observed decadal increase of upper ocean DIC concentrations due to climate change in the Sargasso Sea.We thus only consider the two decades of model data in the context of interannual variability and focus on a process-based understanding of DIC changes along NASTMW pathways.Due to a bias in the model upper ocean salinity trend, the DIC trend at the surface of the Sargasso Sea is mostly between 0 and 1 μmol/L/year (Text S1 in Supporting Information S1), which is small compared to the high background DIC concentrations of ∼2055 μmol L 1 .

Definition of NASTMW
NASTMW was first identified by Worthington (1958) as a water mass with a uniform temperature around 18°C, giving it the moniker Eighteen Degree Water.It is most commonly classified using a temperature range centered around 18°C, typically 17-19°C (Forget et al., 2011;Kwon & Riser, 2004;Maze et al., 2009), with an added stratification constraint that delineates the vertical homogeneity of mode water (Klein & Hogg, 1996;Kwon & Riser, 2004).Alternatively, it is defined through a potential density range with a threshold to delimit low potential vorticity (PV) (Billheimer & Talley, 2016;Talley & Raymer, 1982).
We choose the following constraints for marking Lagrangian parcels as part of NASTMW.
1. Temperature at the location of the parcels is bounded between 17 and 20.5°C.The upper bound is higher than the typical 19°C, since the temperature stratification is slightly stronger in the NASTMW region in the model than in observations, due to model biases (see Text S1 in Supporting Information S1), similar to the model study of Gary et al. (2014).2. The local temperature stratification ∂T/∂z is smaller than 0.01°C m 1 .This stratification limit is weaker than the constraint of ∂T/∂z < 0.006°C m 1 stratification of Kwon and Riser (2004), but the same as in Gary et al. (2014).3. Parcels reside in NASTMW layers of at least 50 m thickness, to exclude thin mixed layers 4. Parcels reside in a contiguous volume of NASTMW of at least 1 × 10 11 m 3 .While this is only on the order of 10 4 times the typical winter NASTMW volume, it excludes many small volumes that are shed off from the main NASTMW volume.5. We only consider NASTMW west of 35°W, to exclude Madeira Mode Water, which has similar properties to NASTMW (Siedler et al., 1987).
Constraints 1-3 are similar to those used in Gary et al. (2014) and Kwon et al. (2015), who also investigate Lagrangian pathways of NASTMW in an ocean model, except that they use a slightly lower temperature upper bound of 20°C.A sensitivity analysis of the NASTMW constraints is found in Text S2 in Supporting Information S1.With the constraints used, we find an average yearly maximum volume of 9.0 × 10 14 m 3 , close to the NASTMW volume of 9.1 × 10 14 m 3 found by Joyce (2012) based on observations (Text S2 in Supporting Information S1).Although the volume of NASTMW exhibits strong seasonal and interannual variability, we find a strong decrease in NASTMW volume starting in 2010, which agrees with the observed decrease in mode water formation found by Stevens et al. (2020).
Figure 1 shows March and September snapshots of modeled NASTMW thicknesses defined using the above criteria.The imprint of mesoscale eddies on the NASTMW structure can be clearly observed (Fratantoni et al., 2013;Gary et al., 2014).Due to model biases, the modeled NASTMW has its core located farther eastward with respect to observations (Text S1 in Supporting Information S1), but the modeled NASTMW is used here to gain a process-level understanding of how DIC concentrations change along NASTMW pathways.

Initialization and Simulation of Lagrangian Parcels
To compute DIC changes along Lagrangian pathways, we simulate the movement of virtual Lagrangian parcels using the Parcels Lagrangian framework (Delandmeter & van Sebille, 2019) (version 2.4.1).These Lagrangian parcels have no defined size and behave like point particles that are advected using ocean model output velocities, while tracer concentrations are equal to those of the ambient water.Along the Lagrangian pathways, we sample molar concentrations of DIC and its precomputed vertical mixing flux, as well as alkalinity, nitrate, phosphate, and their mixing fluxes (Section 2.5).We also sample temperature, salinity, mixed and mixing layer depth, and NASTMW criteria (stratification and contiguity criteria) as these help us distinguish the different NASTMW pathways.
We homogeneously initialize parcels in all parts of NASTMW yearly on 1 September between 1995 and 2015, using the criteria of Section 2.2.Parcels are horizontally spaced apart 0.25°in the zonal and meridional directions, matching the horizontal model resolution.Vertically, parcels are released at fixed depths in NASTMW at 30 m intervals starting at 100 m and reaching down to 460 m (the minimum and maximum summer NASTMW depths in the model).
Parcels are advected forward and backward in time: forward simulations serve to investigate ventilating parcels, persisting NASTMW parcels, and parcels that are exported due to increases in water density.Backward simulations allow us to investigate which parcels have subducted from the mixed layer since the previous summer.Simulations use time steps of Δt = 90 min.For maximum velocities of the order ∼1 m s 1 and a nominal grid resolution of 20 km, this is well below the limit of Δt = Δx/U = 6 hr during which a parcel may travel distances at the grid scale.Parcels are simulated for 3 years, although for most of the analysis in this study, we use only the first year of integration data.Parcel locations and biogeochemical concentrations are saved at daily intervals.In total, we simulate a total of 2 × 861,164 trajectories (2 × 20, 504 ± 6, 487 per year, depending on the NASTMW volume), with the factor 2 indicating forward and backward trajectories.The parcel spacing and temporal output are chosen to balance statistical accuracy with the large computational and storage demands from sampling many biogeochemical fields and identifying individual DIC enrichment and depletion regimes (see Section 2.6).
Parcel trajectories are computed without adding any stochastic displacements that simulate vertical mixing or subgrid-scale isoneutral dispersion (Reijnders et al., 2022).Instead, the parcels represent the (grid-scale) mean flow.We can then see how subgrid-scale vertical mixing fluxes, which are sampled at parcel locations, influence carbon concentrations at the larger grid-scale.

NASTMW Pathway Definition
Parcels may enter and exit NASTMW at infinite points in space and time.We delineate four specific pathways into, within, and out of NASTMW that are typical to its life cycle.Rather than accounting for all possible NASTMW pathways, our study focuses on DIC transformations along these characteristic pathways, acknowledging their potential overlap and non-exhaustiveness.To determine which pathway a parcel follows, we use the daily output snapshots of the parcel trajectory computation.We restrict ourselves to pathways starting on 1 September, as this date coincides with a low NASTMW volume and ensures that the NASTMW maximum in March occurs midway through our simulations.Figure 2 shows the four pathways used in this study.They are defined as follows.
1. Subduction: Parcels were in the mixed layer on 1 September in the previous year and end up in NASTMW on the following 1 September.Here, the mixed layer in NEMO is defined as the layer where the temperature is within 0.2°C of the temperature at 10 m depth.We use the mixed layer rather than the mixing layer (Brainerd & Gregg, 1995) since it is sufficient that a parcel has recently been mixed in summer.Relatively few parcels are expected to subduct on timescales of a year.However, using a longer time scale allows parcels to travel larger horizontal distances, thus widening the domain of parcel origin, possibly far beyond the NASTMW formation region.We opt for specifically investigating short subduction time scales of 1 year in order to keep the parcel origin close to the NASTMW region.Text S3 in Supporting Information S1 discusses results using longer subduction timescales of two and 3 years.2. Persistence: Parcels persist in NASTMW throughout the year until next September.If a parcel temporarily exits the NASTMW-for instance, for a duration shorter than the interval between two consecutive daily snapshots-and then re-enters, it is still considered to have persisted, provided it is in NASTMW at all daily snapshots.3. Ventilation: Parcels from September NASTMW at one point reach the mixing layer, defined by the turbocline depth, and are present in NASTMW again next September.The turbocline depth in NEMO is computed by a transition in vertical mixing regimes, where the vertical eddy diffusivity drops below a predefined threshold (Madec et al., 2013).Air-sea heat and carbon fluxes act in the model's upper layer and propagate by vertical mixing throughout the mixing layer.Note that the mixing layer can partially overlap with NASTMW.A portion of persisting NASTMW parcels may thus ventilate as well, meaning parcels may be double-counted.4. Export: Parcels leave NASTMW and acquire a potential density that is higher than their last value within NASTMW (σ > σ NASTMW ).Leaving NASTMW can thus also be the result of not meeting the stratification criterion anymore when NASTMW is destroyed.We are particularly interested in parcels that persistently maintain their higher densities for a full year, to exclude parcels that only densify temporarily.Since not all parcels will densify immediately when the simulation starts, we here require that parcels have been densified out of NASTMW for at least a continuous year, 2 years after their initialization.We view these parcels as candidates for sequestration since they represent previous NASTMW parcels that are transformed and exported to higher-density waters.When parcels leave NASTMW, their densities may undergo slightly negative fluctuations.We relax the criterion slightly to allow for this: σ > σ NASTMW Δσ, with Δσ = 0.01 kg m 3 .σ is computed using TEOS-10 (McDougall & Barker, 2011).We discuss variations of Δσ in Text S4 in Supporting Information S1.
Note that we focus on the total, time-integrated, change in DIC along a Lagrangian pathway, indicated by ΔDIC, and on timescale distributions of DIC concentration changes along these pathways.The Lagrangian pathways of NASTMW parcels have already been extensively described from a physical perspective by Gary et al. (2014) and Kwon et al. (2015).
As mentioned, the above four pathways are not exhaustive.For example, between consecutive summers a parcel may temporarily leave NASTMW for a few days without reaching the mixing layer.Such a parcel would not fall into any of the above categories.Other examples are parcels that densify but do not remain denser than their NASTMW exit density for a full year or parcels that subduct over timescales longer than 1 year.However, the clear definitions of the four pathways in this study helps create a process-based understanding of carbon fluctuations throughout the life cycle of NASTMW.

Disentangling DIC Changes
The total change in DIC concentrations along a Lagrangian trajectory can be decomposed as follows: Here, the left-hand side is the Lagrangian or total derivative of DIC along the pathway, D[DIC]  Dt .
[DIC] in this study is expressed as a molar concentration of DIC, with units μmol L 1 , so that seawater density changes do not affect DIC concentrations in a control volume.The Lagrangian derivative evolves due to vertical mixing into and out of the water parcel, along-trajectory biogeochemical sources and sinks (indicated by a Lagrangian derivative), airsea exchange (only in the surface layer), and residual terms.
The vertical mixing fluxes are computed from the model output vertical diffusivity coefficient k z and vertical gradients in DIC: We compute these as daily Eulerian fields and sample them along Lagrangian pathways.Because the model includes a horizontal diffusive mixing parameterization only for temperature and salinity, but not for biogeochemical tracers, we do not include a horizontal mixing term.
Not all PISCES state variables are stored as output in FREEBIORYS2V4: only chlorophyll, nitrate, phosphate, silicate, DIC, total alkalinity, and dissolved oxygen are available.Therefore, we need to approximate the biogeochemical sources and sinks from the available variables.We follow Sarmiento and Gruber (2006) and compute the biogeochemical term as (3) Here, the first term on the right-hand side corresponds to soft-tissue production and remineralization.These values are estimated from changes in phosphate concentrations, with r C:P being the stoichiometric ratio of carbon to phosphorus in PISCES.The second term estimates changes in DIC due to calcite formation and dissolution from changes in total alkalinity, correcting for changes in total alkalinity due to changes in nitrate originating from soft-tissue processes (Brewer et al., 1975).The right-hand side of Equation 3 again consists of Lagrangian derivatives.Note that vertical mixing can also increase or decrease concentrations of phosphate, alkalinity, and nitrate along a Lagrangian trajectory.Rather than a biogeochemical effect, this is a physical effect on the DIC concentrations, as already captured in the mixing term in Equation 1.The Lagrangian derivatives of the tracers in Equation 3 are therefore calculated by subtracting their precomputed diffusive mixing fluxes from the total alongtrajectory change in tracer concentrations, such that biogeochemical effects are isolated.For example: with similar equations for total alkalinity and nitrate.
Explicit air-sea exchange in the ocean model occurs only in the uppermost layer (1 m).Changing DIC concentrations in this uppermost layer will influence the concentrations below by vertical diffusive mixing, as in Equation 2. In our Lagrangian simulations, parcels generally do not reach the uppermost layer: no parcels do so when simulated forward in time, and an annual average of 0.027% does so backward in time.This means that parcels do not experience explicit air-sea exchange.Instead, air-sea exchange only indirectly affects parcel DIC concentrations through strong diffusive mixing in the mixing layer.The air-sea exchange term from Equation 1 in our analysis thus becomes part of the vertical mixing term and is not treated separately.We also cannot differentiate between natural DIC and anthropogenic carbon (C ant ; Gruber et al., 1996) since it is not included as a separate tracer in FREEBIORYS2V4.
The residual term captures changes in DIC concentrations that cannot be accounted for by the mixing and biogeochemical terms in Equation 1.We compute it by subtracting the biogeochemical and vertical mixing terms in Equation 1 from the total DIC derivative along the trajectory.The residual contains the discrepancies between the biogeochemical DIC term computed from Equation 3 and the actual biogeochemical contribution as it occurs in PISCES.For the very few parcels that reach the upper model level, missing air-sea exchange will be captured by the residual.It also accounts for unconstrained numerical mixing of DIC and other biogeochemical state variables.Atmospheric and riverine deposition of phosphorus and nitrogen are neglected when computing their biogeochemical changes, thus also leaving an imprint on the residual.We cannot isolate the effects of local freshening or evaporation along Lagrangian parcel trajectories, because these terms cannot be constrained: usually, these can be estimated from salinity, but salinity is explicitly horizontally mixed, unlike biogeochemical tracers.Thus, we do not normalize DIC by local salinity in our analysis because horizontal mixing of salinity would cause a drift in the budget over time.In addition, we are explicitly interested in the effect of vertical mixing on DIC concentrations.Evaporation-and precipitation-related freshening directly impact nutrient concentrations only in the model surface layer.Below the first meter, the effect of mixing of fresher and more saline waters on DIC concentrations is part of the mixing term (Equation 2).
With the implicit treatment of air-sea exchange in the vertical mixing flux, and with the division of biogeochemical contributions into a soft-tissue term and a carbonate term, the total DIC derivative is then decomposed as: Journal of Geophysical Research: Oceans

10.1029/2023JC020814
To obtain an illustrative comparison of the contribution of each term to D[DIC] Dt , we compute the sum of the magnitudes of each component in the right-hand side Equation 5 for each trajectory.Figure 3a shows the percentage by which each component contributes to the sum of their magnitudes, computed from forward-in-time 1year trajectories from all years (30-day segments of disentangled time series in the year 2000 are found in Figure S35 in Supporting Information S1).Note that the net biogeochemical contribution may be lower than the sum of its constituents, as these may have opposite signs.We aim to discern the significance of each term in determining the ΔDIC across entire trajectories.To achieve this, we exclude time steps that cumulatively account for less than 5% of the sum of magnitudes, thereby focusing on the time steps for which the total rate of change predominantly influences ΔDIC.This method effectively screens out instances with minimal D [DIC]   Dt values, since these are susceptible to disproportionate impacts from numerical inaccuracies, thus distorting the residual's impact.The residual accounts for about 19% of the sum of magnitudes of each term, which per timestep is higher than the contribution of the mixing and carbonate terms.However, our main analysis will show that the residual will fluctuate in sign between timesteps and trajectories, such that its overall influence is limited.Figure 3a shows that biogeochemical terms are the dominant contributor to the total rate of change at each time step, with soft-tissue processes in turn being the main constituent.As expected, mixing fluxes become increasingly important in the mixing layer, where the contribution is twice as high as in the mixed layer and 13 times higher than below the mixing layer.Carbonate terms are roughly six times smaller than soft-tissue terms.
Because the residual term is composed of multiple unconstrained constituents (see above), we computed correlations between the residual term and other terms.This can shed light on whether the residual is systematically linked to one of the other processes.Figure 3b shows that the residual is uncorrelated with the other constituents, meaning none of these are systematically over-or underestimated.There is a weak correlation between the total DIC rate of change and the residual, which, as discussed earlier, could be related to some unconstrained process that is missing in the budget.We note that the soft-tissue contribution has a very high correlation with the total biogeochemical term, meaning that it almost always dominates this term.The mixing term more often than not is of the opposite sign as the biogeochemical terms, indicating net counteraction.
For each water parcel trajectory, we use the disentangled rates of change to reconstruct time series of DIC anomalies with respect to the initial concentration for soft-tissue processes, carbonate processes, vertical mixing, and residual processes.

Identifying Enrichment and Depletion Regimes and Timescales
One of our aims is to investigate the timescales and strengths at which DIC concentrations are depleted or enriched along NASTMW pathways.We examine both the cumulative ΔDIC along the pathway and intermittent DIC enrichment and depletion regimes affecting DIC concentrations.
We have opted for a straightforward approach to define DIC enrichment and depletion regimes.Specifically, these regimes are defined by the intervals in the time series between local minima, which indicate the start of enrichment, and local maxima, which indicate the start of depletion.Thus, the duration of regimes corresponds to the time intervals between these local minima and maxima, whereas the regime's magnitude is the change in DIC concentration during these intervals.In our analyses, we do not truncate regimes that have their start date before the pathway's end date: we analyze trajectories from 1 September till 1 September in the subsequent year (or 2 years, in case of export), but if a regime starts before this end date, we still include its entire timescale in our analysis.
To reduce the impact of minor fluctuations occurring over short periods (a few days or less), we apply a centered moving average to the time series after its disentanglement into constituent terms.Although this approach smooths the series, it does not completely eliminate short-time variability.Instead, it emphasizes significant changes in DIC concentration, minimizing the influence of brief, minor fluctuations at time scales shorter than the window size.Therefore, the window size partially sets the scale for which regimes are deemed significant.
The primary advantage of this methodology lies in its simplicity, offering a clear lens to assess the main dynamics of DIC variations over time.Given our model's nominal resolution of 1/4°, it does not resolve the submesoscale, which is associated with nutrient transport and biogeochemical structuring at timescales of the order of days (Lévy et al., 2012).Instead, we will apply a window length of 10 days, which is still much shorter than the lifetimes of mesoscale eddies (months) by which nutrients are supplied (McGillicuddy et al., 1998), and instead is of the same order as typical remineralization timescales (Siegel et al., 1999).This allows us to resolve processes on timescales of the order of a week and higher.Additionally, Text S5 in Supporting Information S1 repeats our analysis without any smoothing, and with smoothing using window sizes of 6 and 20 days.Especially when no smoothing is applied, the bulk of the regime lengths are shorter than 10 days, while their magnitudes are also much smaller.This illustrates the need for smoothing to shed light on processes at longer timescales.

Results: DIC Enrichment and Depletion Along NASTMW Pathways
Table 1 summarizes the mean ΔDIC for each pathway, averaged over all years, and the mean fraction of parcels that follow each pathway.This shows that parcels undergo the largest change in DIC while subducting from the mixed layer over the course of a year, while at the same time, only 1.2% of parcels make up this pathway.Forward in time, the most prominent pathway is persistence, accounting for 25.9% of the trajectories, followed by export and ventilation.Since the ΔDIC does not follow a normal distribution, we do not include statistics of variance in Table 1.Instead, these distributions are provided and discussed directly in the following subsections.
The above pathways are non-exhaustive and non-exclusive (Section 2.4), but the trajectories that are not accounted for can follow a myriad of paths, hindering their analysis.For comparison, Text S6 in Supporting Information S1 gives a brief overview of the total ΔDIC and timescale distributions for all trajectories emanating from NASTMW in the forward and backward simulations.In summary, for all trajectories starting in the mode water, ΔDIC = 5.4 μmol L 1 over the following year, while for all trajectories that ended up in the mode water, ΔDIC = 4.3 μmol L 1 over the previous year.
Each pathway is discussed individually in the following subsections.We plot the distribution of total ΔDIC and show the relative contribution of biogeochemical, mixing, and residual processes for different total ΔDIC magnitudes.The integrated monthly strength of each process is also discussed.Lastly, we examine distributions of enrichment and depletion events spread across timescales and processes.

The Subduction Pathway
As can be seen from Figure 4, subduction of parcels into NASTMW has a large impact on DIC concentrations in the water parcels: the mean increase is ∼100 μmol L 1 , though with a large spread for different parcels (Figure 4a).The dominant contribution comes from vertical mixing, which acts chiefly when the parcel is in the mixing layer.The mixing contribution grows from September until December, after which it decreases to nearzero around May, where it remains steady for the rest of the year (Figure 4f).This increased mixing coincides with increased downwelling (not shown): upon initialization, parcels experience downwelling at a mean rate of 0.1 m per day in September, increasing steadily to a maximum of approximately 1 m per day in March, after which downwelling velocities reduce again to 0.1 m per day in May.The strength of the vertical mixing flux is then proportional to the downwelling rate and the vertical gradient at the depth of the parcel (Figure S9 in Supporting Information S1).Most parcels exit the mixed layer in April when the mixing layer shoals again and mixing decreases drastically.In the winter months, when the mixing layer deepens and entrains nutrients, there is a small negative soft-tissue carbon contribution from primary productivity (Figure 4e).Although the spring bloom is visible as a minimum in the mean yearly fifth percentile, subducting NASTMW parcels experience a mean increase in DIC in April, as they move below the mixing layer again.Remineralization continues over the course of the following months, as subducting NASTMW is rich in semi-labile dissolved organic carbon that can be remineralized (Carlson et al., 1994;Krémeur et al., 2009).We investigate the timescales and magnitudes associated with the enrichment and depletion regimes in Figure 5. Vertical mixing is the dominant contributor across regime timescales (Figure 5a).Only at timescales of 10 days or less, the net ΔDIC is slightly negative.DIC depletion regimes at timescales of 30 days or less primarily occur around March (not shown), during peak primary production.The regime distribution has a long positive tail with around 43% of DIC changes associated with timescales of more than 100 days, mostly associated with vertical mixing.When regime detection is applied specifically to the mixing-related DIC anomaly time series, almost 90% of the ΔDIC is associated with regimes with these long timescales (Figure S41 in Supporting Information S1).This shows that vertical mixing steadily increases DIC concentrations as parcels subduct into NASTMW, adjusting to the ambient vertical DIC distribution, with larger DIC concentrations at depth that can supply the parcel with DIC from below (Figure S9 in Supporting Information S1).
While the soft-tissue, carbonate, and residual terms have a relatively minor contribution to the total ΔDIC (Figure 4c), their contributions are of similar order to the total ΔDIC of the other pathways (shown later).The residual and carbonate processes here exhibit the largest contribution (and spread) when parcels are in the mixing layer (Figures 36a and S36b in Supporting Information S1).
A yearly average of 1.2% of backtracked parcels originate from the mixed layer on the previous September 1st.
We also investigated subduction occurring over two and 3 years.As this allows parcels to subduct over longer timescales, more parcels meet this criterion (3.7% and 5.8% respectively).This is discussed in Text S3 in Supporting Information S1.In summary, these parcels on average experience a total ΔDIC of similar order (110 μmol L 1 in both cases; see Figures S13 and S15 in Supporting Information S1).While vertical mixing is still the dominant contributor in both cases, soft-tissue processes progressively make up a higher share of the mean total ΔDIC: for subduction over 1 year, these make up 9% of ΔDIC, with this contribution increasing to 20% and 30% when subduction occurs over 2 and 3 years, respectively.These soft-tissue contributions are an important factor in shaping the vertical distribution of DIC (Sarmiento & Gruber, 2006).The longer a parcel takes to subduct, the more time soft-tissue remineralization processes have to directly increase the parcel's DIC, whereas during quick subduction the parcel will instead adapt its DIC to ambient conditions through mixing.

The Persistence Pathway
Figure 6 shows the total transformation of DIC concentrations within parcels that persist in NASTMW throughout a full year.On average, this accounts for 25.9% of all parcels and thus agrees with the model study of Gary et al. (2014), where 74% of NASTMW parcels exit the mode water within a year.Note, however, the large interannual range of parcels that comprise this pathway (Figure 6a).The minimum of 2.9% is associated with parcels initialized in September 2013, where the following year 2014 marked a strong decline in the modeled Outliers are defined as regimes with magnitudes that fall outside the whisker ranges.
NASTMW volume toward its minimum in the summer of 2014 (Figure S14 in Supporting Information S1).
Interannual variability in NASTMW formation and volume is commonly observed (Billheimer & Talley, 2013;Stevens et al., 2020).Generally, we find that the percentage of parcels that persist in NASTMW is correlated with the volume in the next year, with a Pearson-R of 0.88 (p < 0.001); fewer parcels can persist in NASTMW if its volume shrinks from 1 year to the next.
Figure 6b shows that positive contributions are dominated by soft-tissue remineralization, which has a slightly positive monthly mean contribution year-round (Figure 6e).Vertical mixing leaves a distinctly negative imprint on the ΔDIC of persisting NASTMW parcels, meaning that it depletes parcels of carbon.This occurs specifically in winter (Figure 6f), when the mixing layer deepens, and causes 19% of parcels to have a negative ΔDIC.Because winter mixing is a primary driver of NASTMW formation, some parcels that persist in NASTMW may in fact reside in well-mixed newly formed NASTMW.Vertical mixing can then act to deplete DIC from these parcels as it is supplied to the euphotic zone.We find that the vertical displacement of a parcel is a predictor for the total ΔDIC (Pearson-R of 0.51, p < 0.001): parcels that move deeper, are more likely to have increased DIC concentrations.This can be due to a smaller likelihood of being temporarily entrained in NASTMW regions that are in contact with the mixing layer.
While the net residual term is smaller than the soft-tissue and mixing terms, the carbonate term is smaller than these residual terms, so we neglect it in our discussion for this pathway.Both the carbonate and residual terms show no clear yearly cycle (Figures S40c and S40d in Supporting Information S1).
From Figure 7 we see that net depletion is associated with timescales of 30 days and less.For timescales between 10 and 30 days, about half of the depletion is attributed to vertical mixing.This is largely associated with the winter convection.The contribution of photosynthesis, a soft-tissue process, has its mode at the 10-20 days timescale.The associated timescales of the order of weeks correspond to observations of the spring bloom in the Sargasso Sea (Nelson et al., 2004).
Soft-tissue DIC enrichment, associated with remineralization, has its mode at the 20-30 days timescale, but its tail extends over longer timescales than DIC depletion, with a contribution of almost half the total net ΔDIC at timescales longer than 100 days (see the timescales associated with the soft-tissue DIC anomaly, Figure S42 in Supporting Information S1). Figure 7b shows that for timescales less than 20 days, the mode ΔDIC of each individual regime is close to zero at timescales less than 20 days and gradually increases to 6.2 μmol L 1 for timescales longer than 100 days.The net negative ΔDIC at timescales less than 30 days then suggests that strongly negative "outliers" are responsible for a net decrease at this timescale.For regimes with durations of around a month, we find that these outliers are largely concentrated around March, coinciding with the spring bloom.The spring bloom can thus be linked to strong anomalous DIC depletion for parcels residing in NASTMW.

The Ventilation Pathway
On average, parcels that ventilate by temporarily reaching the mixing layer undergo a negligible net ΔDIC (Figure 8).The mean total ΔDIC of 0.7 μmol L 1 is smaller than the interannual standard deviation.On aggregate, this means that the ventilation of existing NASTMW parcels does not influence the NASTMW carbon content appreciably.However, the ΔDIC distribution for individual trajectories in Figure 8a shows a wide spread between positive and negative values.We find that the total ΔDIC along this pathway is correlated with the net downward displacement in the water column with a Pearson-R of 0.60 (p < 0.001).Thus, over one yearly ventilation cycle, the net deepening of a parcel is a predictor of its increase in DIC.
Figures 8b and 8c show a strong counteraction of DIC enrichment from soft-tissue remineralization and a negative contribution for vertical mixing.Figure 8f shows how winter mixing is responsible for the decrease in DIC, as the parcel exchanges its DIC with the mixing layer, supplying nutrients for primary production in the euphotic zone, as well as equilibrating with the upper layer in which air-sea fluxes allow for atmospheric gas exchange.Although the mean monthly soft-tissue term never becomes negative, the mean 5th percentile has a minimum around February and March.This coincides with a modeled maximum in net primary production of phytoplankton, associated with the spring bloom.Only some parcels experience this negative soft-tissue contribution directly, as not all parcels can reach the euphotic zone where primary production occurs.Instead, many parcels are linked to the spring bloom indirectly, supplying it with nutrients from the deeper parts of the mixing layer.Following this, the mean soft-tissue term has a slight positive maximum in April (Figure 8e), as the mixing layer shoals and moves above the parcel, allowing organic carbon to remineralize.The net DIC term remains positive over the following months (Figures 8d and 8e).Carbonate processes have a small positive yearly contribution, also peaking in April, after most parcels have left the mixing layer (Figure S22e in Supporting Information S1).Note that residual processes have no net effect on the ΔDIC for this pathway since they effectively cancel each other out (Figure 9a and Figure S40f in Supporting Information S1).Mixing is the main contributor to the net depletion of DIC at timescales shorter than 40 days (Figure 9a).The distribution of DIC enrichment and depletion regimes for ventilating parcels is somewhat similar to that of persisting NASTMW parcels, albeit with a larger contribution from mixing at the short timescales.Figure 9b shows that the interquartile range and whiskers are symmetric with the median around 0 for timescales less than 20 days, with the median becoming positive at longer timescales.However, since the normalized ΔDIC for regimes at timescales less than 30 days is negative, this net negative contribution must be due to strongly depleting outlier regimes, associated with vigorous mixing.In Figure S43 in Supporting Information S1, we identify regime timescales and magnitudes based on the DIC anomaly due to mixing processes.Interestingly, we find that about   40% of the net contribution of mixing processes has regime timescales of more than 100 days.When assuming a regime-based view of the total DIC anomaly, the relatively steady DIC depletion due to mixing during wintertime ventilation can be temporarily counteracted by local soft-tissue remineralization, such that mixing is not able to cause the total DIC anomaly to persistently decrease for such long timescales.

The Export Pathway
NASTMW parcels that get exported to denser surroundings mostly undergo a net positive ΔDIC (Figure 10), with a yearly mean of 9.9 μmol L 1 .We recall that here, parcels are integrated for a total of 2 years, during which their potential density remains, for at least the entire second year, higher than their potential density upon exiting NASTMW (Section 2.4).
The distribution of ΔDIC is asymmetric, with a longer tail in the positive direction (Figure 10a).Soft-tissue processes make up the bulk of the ΔDIC for trajectories in the positive tail (Figure 10b).Note that for small and negative ΔDIC, the residual becomes more prominent, indicating that the DIC budget becomes less well constrained by our disentanglement method (Section 2.5).Mixing fluxes again are concentrated in the winter months of the first year.Due to the way we select exported parcels, some may be temporarily entrained into the mixing layer in the first year of integration.However, outside of the first winter months, the mean DIC changes are solely governed by soft-tissue processes.Note that after the first year, these are effectively zero, albeit with the 5th and 95th percentile ranging between values of around ±2 μmol L 1 (Figure 10e), while the residual has a range of ±1 μmol L 1 (Figures S40g and S40h in Supporting Information S1).Since these ranges are of the same order, while the mean is close to zero, we conclude that after the NASTMW parcels remain exported, they undergo no clear net DIC depletion or enrichment in the second year; the change in ΔDIC occurs before.
Looking at enrichment and depletion regimes and timescales (Figure 11), we see that at short timescales of less than a month, positive and negative ΔDIC regimes nearly balance one another.The net ΔDIC of regimes with longer timescales becomes positive, dominated by soft-tissue remineralization.Mixing mainly acts to deplete exporting NASTMW parcels of their DIC at timescales of around a month and less, while, photosynthesis and remineralization balance each other out at timescales up to 2 weeks.At short timescales, the enrichment and depletion magnitudes show a symmetric distribution, also in terms of outliers (Figure 11b), meaning that at these timescales, there is a balanced counteraction.Exported parcels remain in waters of higher density than the local NASTMW and thus are situated at deeper depths.Compared to the persistence pathway (Figure 7a), it can therefore seem counter-intuitive that exported parcels experience more photosynthesis-dominated regimes at timescales of a few weeks, even if these regimes are balanced by remineralization (Figure 11a).However, note that in the thinner flanks of NASTMW (Figure 1b), parcels are less likely to persist if the NASTMW shrinks between years, such that they get exported instead.On average, exported parcels in such regions may then be located at shallower depths and experience more photosynthesis than parcels that in such years persist in the central core of NASTMW, which extends to deeper depths.
While the number of parcels that we consider exported depends on the potential density threshold, Δσ, these parcels qualitatively exhibit similar total ΔDIC, regime magnitudes, and timescales (Text S4 in Supporting Information S1).

Summary and Discussion
We adopted a Lagrangian frame of reference to study how DIC concentrations are altered along different NASTMW pathways in an eddy-permitting model that reproduces NASTMW for process study purposes.Of the four pathways considered, parcels experience the largest increase in DIC during subduction from the mixed layer into NASTMW, with an order of magnitude of ∼100 μmol L 1 .However, only 1.2% of parcels subduct over the course of a year.This percentage increases to 3.7% and 5.8% if we allow parcels to subduct over longer timescales of 2 and 3 years respectively, while the increase in DIC remains similar at ∼110 μmol L 1 .This implies that parcels generally take longer to move from the summer mixed layer to summer NASTMW, while the total enrichment along the subduction path is roughly unaffected.The persistence pathway makes up the largest fraction at 25.9%, which is in agreement with the findings of Gary et al. (2014).The subduction, ventilation, and export pathways had not been previously explored, and no previous study had yet investigated DIC concentrations along mode water pathways.Kwon et al. (2015) explore the timescales between the re-outcropping of NASTMW parcels, which somewhat resembles our ventilation pathway.However, the main differences are that our study requires parcels to be directly in the mixing layer instead of in an outcropping NASTMW column, and to reside in NASTMW for two consecutive summers, which complicates direct comparisons.While our results do not apply directly to other subtropical mode waters, they may provide initial qualitative hints of their DIC dynamics until further research is available.
Generally, a downward displacement of parcels is associated with net DIC enrichment.This is in line with the positive downward gradient of observed DIC in the water column, set by the remineralization of organic material.It then makes sense that the net ΔDIC for each of the four pathways is positive (Table 1), as the Subtropical Gyre is associated with net sinking (Berglund et al., 2022;Spall, 1992).Subducting parcels adjust to the vertical DIC gradient mostly through vertical mixing, but as longer subduction times are allowed (Text S3 in Supporting Information S1), the contribution from mixing shrinks in favor of direct remineralization within the parcel.Following the spring bloom, downward-moving persisting NASTMW parcels and exported NASTMW parcels also experience enrichment by remineralization directly.Ventilating parcels on the whole experience only a small appreciable net DIC enrichment during 1 year, which is of a similar order of magnitude as the surface ocean model DIC trend (Figure S12 in Supporting Information S1).This could indicate that ventilating trajectories propagate DIC increases at the surface toward the NASTMW interior.Unfortunately, model biases (Text S1 in Supporting Information S1) prevent us from directly comparing these values to the observed Eulerian DIC trend at the BATSsite (1.5 μmol kg 1 yr 1 ; Bates, 2012).These ventilating parcels follow a cycle where they first supply DIC to the euphotic zone by vertical mixing in late winter, followed by remineralization of organic material in the following months.The net sinking distance over a ventilation cycle is also correlated with enrichment.
The timescales at which photosynthesis depletes parcels of DIC are of the order of weeks, corresponding to observations of the spring bloom in the Sargasso Sea (Nelson et al., 2004).This depletion is most prominent in persisting, ventilating, and exported parcels.For these pathways, mixing also has a depleting component at similar timescales, which is associated with an upward supply of DIC toward the surface.Mixing-associated depletion peaks in March when the mixed layer depth is at its maximum, re-establishing a link between NASTMW and the surface.This highlights the importance of NASTMW as a reservoir that can store and resupply carbon to the ocean surface (Bates et al., 2002), with relevance to ocean biology and ocean-atmosphere carbon exchange.
While the four pathways chosen for this study are non-exclusive and non-exhaustive, they still provide a good picture of carbon cycling in NASTMW.This is corroborated by the fact that when all trajectories are analyzed together (Text S6 in Supporting Information S1), they qualitatively resemble the persistence pathway and aspects of the ventilation and export pathways.Backward-in-time trajectories also resemble the persistence pathway, with the minor contribution from the subduction pathway superimposed on this.We think that the trajectories that are unaccounted for will often resemble the persistence pathway, because the aggregate of all trajectories qualitatively resembles the persistence pathway strongly, while it only accounts for about a quarter of all pathways.
These "unaccounted trajectories" can for instance be parcels leaving NASTMW only temporarily while not being ventilated or exported, or parcels that remain in the direct vicinity of NASTMW where carbon dynamics may be somewhat similar.
In the context of carbon sequestration in the Sargasso Sea, our findings show that once NASTMW parcels are exported, they experience only small changes in DIC concentrations, which average out each other for at least 1 year.Some of these parcels may re-outcrop during strong convection in the following winters (Kwon et al., 2015), which re-liberates their carbon (Bates et al., 2002).However, exported parcels may also further sink and travel out of the NASTMW ventilation region to the Subpolar Gyre over timescales >10 years (Gary et al., 2014).Berglund et al. (2022) show that transport between the Subtropical Gyre and Subpolar Gyre occurs in a downward-spiraling fashion over timescales of several decades.The long-term fate of carbon sequestered through the export of NASTMW parcels provides a direction for future research.We find that before such sequestration, DIC concentrations are set primarily due to mixing and to a lesser extent due to remineralization during a parcel's journey from the mixed layer to NASTMW.
While this study did not consider NASTMW trends in the face of global (ocean) warming, we know that NASTMW is influenced by climate change.Observations show that ventilation and formation of NASTMW have decreased due to surface warming (Stevens et al., 2020).Under a strong warming scenario (SSP5-8.5),CMIP6 models report a shallowing of the winter mixed layer in the NASTMW formation region (Fox-Kemper et al., 2021), which would further decrease the ventilation of NASTMW parcels.This would decrease NASTMW formation and thus volume, decreasing the amount of heat and carbon that it can buffer, in the latter case allowing less carbon to potentially sequester through export of NASTMW parcels.We did not investigate these effects due to model biases (Text S1 in Supporting Information S1) and because of the limited temporal extent in the face of high interannual and decadal NASTMW variability (Deremble & Dewar, 2013;Kwon & Riser, 2004;Stevens et al., 2020).However, a future similar study using trajectories in a state-of-the-art Earth system model could shed light on this.
Of the four processes considered, soft-tissue processes and vertical mixing dominate DIC concentration changes.Carbonate processes play a much smaller role at these NASTMW depths and sub-annual timescales (Williams & Follows, 2011), sometimes with a smaller contribution than the residual.We note that here we estimated the carbonate term from alkalinity, while in reality it is also affected by nitrification.Still, we could not detect a systematic bias in the carbonate term (Section 2.5).Future studies that are interested in how carbonate processes influence DIC concentrations along Lagrangian trajectories would benefit from using a biogeochemical model where all the required state variables are stored for constraining the carbonate term.Exact calculation of the carbonate term can for example, become important when studying Lagrangian carbon dynamics in the deep ocean where carbonate dissolution plays a more prominent role.
Time series of carbon along Lagrangian trajectories have previously been studied by Brady et al. (2021) in the context of upwelling in the Southern Ocean.Like Cetina-Heredia et al. ( 2018) did for nitrate, Brady et al. (2021) characterized DIC timescales by the Lagrangian decorrelation timescale.However, decorrelation metrics assume that Lagrangian biogeochemical time series are stationary, while in reality biogeochemical depletion and enrichment are highly dependent on the spatial (vertical and horizontal) and temporal location of the water parcel.For example, the DIC time series of a parcel can be influenced by the entrainment into an eddy that experiences high primary productivity, or by a quick decrease in vertical mixing after subduction through the thermocline.These processes are highly non-linear and non-stationary.This is why we introduced a straightforward approach where we define enrichment and depletion regimes between local minima and maxima in smoothed Lagrangian DIC anomaly time series.Also, unlike methods based on spectral analysis, our method reveals regime timescales while staying agnostic about any periodicity, which may not be present.
Model data constraints prevent us from investigating the full mesoscale spectrum or submesoscale processes and variability, as higher resolution model data is not available for large regions over the span of decades.To deal with model biases in our eddy-permitting set-up, we tuned our NASTMW definition criteria for a realistic model volume and thickness, but, for example, our winter NASTMW core is too shallow (Text S1 in Supporting Information S1), possibly leading to an underestimation of the export pathway, which would be influenced by larger vertical excursions of NASTMW depths between summer and winter.Gan et al. (2023) show that an eddy-rich model would more faithfully reproduce the observed NASTMW spread and volume.Moreover, the unresolved submesoscale dynamics have large implications for biogeochemistry, for example, by creating fronts that provide short-lived nutrient pulses of just a few days (Lévy et al., 2012;Mahadevan, 2016).The findings in this study instead are chiefly related to mesoscale ocean dynamics, with a resolution similar to physical ocean components of state-of-the-art earth system models (Hewitt et al., 2020).While a similar study at a submesoscale resolution is at present computationally unfeasible, more process-based submesoscale Lagrangian studies, such as by Freilich et al. (2022) on phytoplankton growth rates, provide first steps in understanding Lagrangian carbon dynamics at smaller scales.
Our findings highlight that individual parcels undergo DIC enrichment and depletion regimes over a range of timescales and magnitudes, due to a complex interplay of vertical mixing and biogeochemical processes.We find that on short timescales of the order of weeks, enrichment and depletion often oppose one another, which for the subduction, persistence, and export pathways leads to a relatively small net ΔDIC at these timescales, which is why the largest net changes in DIC occur on longer timescales.Bulk Eulerian studies average out this underlying complexity of enrichment and depletion unfolding over different timescales.Our approach can thus complement Eulerian approaches when investigating the carbon cycle.This also applies beyond the context of NASTMW or mode waters entirely; our Lagrangian process and timescale decomposition can be applied to any study aimed at better understanding carbon dynamics along water pathways.For example, while Eulerian approaches are the dominant method to study marine carbon dioxide removal through ocean alkalinity enhancement (Fennel et al., 2023), our methods may be used to study how alkalinity enhancement interventions influence DIC at different timescales along pathways of water downstream from intervention sites.Similarly, one could disentangle influences on pCO 2 , phytoplankton, or even ecosystems if more model state variables are available.

Figure 2 .
Figure 2. Sketch of the four Lagrangian pathways in, within, and out of NASTMW as covered in this study in latitude-depth space.NASTMW is indicated by the dark blue blob.Actual NASTMW boundaries and outcropping locations also exhibit longitudinal, seasonal, and interannual variation.

Figure 3 .
Figure 3. (a) Percentage at which each term in Equation5contributes to the sum of their magnitudes, averaged over each time step.Percentages are computed using forward trajectories in all years, selecting only the time steps for which the total rate of change is responsible for at least 95% of the sum of all total magnitudes.We also examine the total biogeochemical rate of change, composed of soft-tissue and carbonate terms, which may have opposite signs.(b) Correlations between each of the terms, including the total derivativeD[DIC]  Dt .

Figure 4 .
Figure 4. Transformation of DIC concentrations along pathways of parcels that subduct and reach NASTMW.(a) Distribution of total ΔDIC per trajectory for all initialization years 1995-2015.Error bars indicate standard deviation for each bin per year.The average number of trajectories of this pathway is indicated as a percentage of all simulated trajectories per year, with min-max ranges indicated in brackets.(b) Relative contribution of each process to the average ΔDIC, for different ΔDIC strengths.The sum of contributions is always 1, meaning that contributions greater than 1 are balanced by contributions of the opposite sign.The relative contribution is computed only for bins with their edges between the 1st and 99th percentile of ΔDIC.(c) Mean of yearly average ΔDIC of trajectories, with the standard deviation across years in brackets.(d-f) Trajectory statistics per month, averaged across the different years for the total DIC rate of change (d), softtissue processes (e), and mixing fluxes (f).The black lines show the monthly integrated DIC change, averaged over trajectories and the orange shading the 5th and 95th trajectories percentiles.These trajectory statistics are averaged over the different years.The blue shading show the inter-annual range of the mean monthly integrated DIC changes (black line).Carbonate and residual terms are much smaller and are shown in Figures S40a and S40b in Supporting Information S1.

Figure 5 .
Figure 5. ΔDIC contribution of each timescale for parcels that subduct and reach NASTMW.(a) Relative ΔDIC of regimes of each timescale.This quantity is computed by summing the magnitudes of each positive and negative regime for all trajectories across years and then normalizing by the sum of ΔDIC of each whole trajectory."Net" shows the positive minus negative normalized ΔDIC.Because the distribution has a long tail, regimes longer than 100 days are grouped together.(b) Boxplot of ΔDIC magnitudes for each regime for each timescale.Maxima and minima of outliers are indicated by triangles.The number of positive and negative outliers is indicated as a percentage of the total number of regimes, which is indicated above.The boxplot follows the classical definition: whiskers are defined as Q 1 1.5*IQR and Q 3 + 1.5*IQR.Outliers are defined as regimes with magnitudes that fall outside the whisker ranges.

Figure 6 .
Figure 6.Similar to Figure 4, but for parcels that persistently remain in NASTMW.(a) Distribution of total ΔDIC per trajectory, averaged over initialization years 1995-2015.(b) Relative contribution of each process to the average ΔDIC, for different ΔDIC strengths.Note that when the average ΔDIC is negative, positive contributions to ΔDIC (e.g., soft-tissue remineralization) have a negative relative contribution.(c) Mean of yearly average ΔDIC of all trajectories.(d-f) Trajectory statistics per month, averaged across the different years for the total DIC rate of change (d), soft-tissue processes (e), and mixing fluxes (f).Carbonate and residual terms are much smaller and are shown in Figures S40c and S40d in Supporting Information S1.

Figure 7 .
Figure 7. ΔDIC contribution of each timescale for persisting NASTMW parcels.(a) Relative ΔDIC of regimes of each timescale.(b) Boxplot of magnitudes of each regime for each timescale.

Figure 8 .
Figure 8. Similar to Figure 4, but for ventilating NASTMW parcels.(a) Distribution of total ΔDIC per trajectory, averaged over initialization years 1995-2015.(b) Relative contribution of each process to the average ΔDIC, for different ΔDIC strengths.(c) Mean of yearly average ΔDIC of all trajectories.(d-f) Trajectory statistics per month, averaged across the different years for the total DIC rate of change (d), soft-tissue processes (e), and mixing fluxes (f).Carbonate and residual terms are much smaller and are shown in Figures S40e and S40f in Supporting Information S1.

Figure 9 .
Figure 9. ΔDIC contribution of each timescale for ventilating NASTMW parcels.(a) Relative ΔDIC of regimes of each timescale.(b) Boxplot of magnitudes of each regime for each timescale.

Figure 10 .
Figure 10.Similar to Figure 4, but for exported NASTMW parcels.(a) Distribution of total ΔDIC per trajectory, averaged over initialization years 1995-2015.(b) Relative contribution of each process to the average ΔDIC, for different ΔDIC strengths.(c) Mean of yearly average ΔDIC of all trajectories.(d-f) Trajectory statistics per month, averaged across the different years for the total DIC rate of change (d), soft-tissue processes (e), and mixing fluxes (f).Carbonate and residual terms are much smaller and are shown in Figures S40g and S40 h in Supporting Information S1.

Figure 11 .
Figure 11.ΔDIC contribution of each timescale for exported NASTMW parcels.(a) Relative ΔDIC of regimes of each timescale.(b) Boxplot of magnitudes of each regime for each timescale.

Table 1
Yearly Mean ΔDIC for Each Pathway in the Hindcast Model, Including Mean Occurrences