Origins of material in the ocean are commonly identified by tracing Lagrangian particle trajectories backward-in-time in two or three dimensions. While this is mathematically consistent, numerical computations are hampered by numerical round-off and truncation errors, leading to discrepancies between forward- and backward-in-time trajectories. The chaotic nature of ocean flows amplifies these errors. We identify an additional issue with Lagrangian backtracking, related to the reversal of stability with regards to velocity convergence and divergence. Trajectories near convergent regions are stable to numerical errors when calculated forward-in-time but become unstable backward-in-time. The timescales at which trajectories reside in convergent zones are thus underestimated in backward-in-time computations, meaning convergent regions (downwelling zones) become underrepresented and divergent zones (upwelling zones, river mouths) overrepresented as trajectory sources. Using mesoscale experiments representing common set-ups, we show that already for timescales of less than half a year, this leads to systematic biases in the regions identified as particle origins. These biases can extend over distances of thousands of kilometers. While this stability bias is linked to divergence, it is not only limited to 2D trajectories in 3D flows, as we discuss how inappropriate treatment of surface boundary conditions in 3D Lagrangian studies can also introduce an effective non-zero divergence. These findings have consequences for source-attribution modeling, for example in the context of water mass tracing, ecology, and pollution studies. Backtracking is typically applied to material that has accumulated in convergent zones, where the stability bias is especially relevant, which further impedes source attribution.

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. While the mean ΔDIC for parcels that persist within NASTMW in one year is relatively small at +6 µmol/L, this masks underlying dynamics: individual parcels undergo interspersed DIC depletion and enrichment, spanning several timescales and magnitudes. The strongest ΔDIC is during subduction of water parcels (+101 µmol/L  in one year), followed by transport out of NASTMW due to increases in density in water parcels (+10 µmol/L). 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.

To capture the effects of mesoscale turbulent eddies, coarse-resolution Eulerian ocean models resort to tracer diffusion parameterizations. Likewise, the effect of eddy dispersion needs to be parameterized when computing Lagrangian pathways using coarse flow fields. Dispersion in Lagrangian simulations is traditionally parameterized by random walks, equivalent to diffusion in Eulerian models. Beyond random walks, there is a hierarchy of stochastic parameterizations, where stochastic perturbations are added to Lagrangian particle velocities, accelerations, or hyper-accelerations. These parameterizations are referred to as the 1st, 2nd and 3rd order ‘Markov models’ (Markov-N), respectively. Most previous studies investigate these parameterizations in two-dimensional setups, often restricted to the ocean surface. On the other hand, the few studies that investigated Lagrangian dispersion parameterizations in three dimensions, where dispersion is largely restricted to neutrally buoyant surfaces, have focused only on random walk (Markov-0) dispersion. Here, we present a three-dimensional isoneutral formulation of the Markov-1 model. We also implement an anisotropic, shear-dependent formulation of random walk dispersion, originally formulated as a Eulerian diffusion parameterization. Random walk dispersion and Markov-1 are compared using an idealized setup as well as more realistic coarse and coarsened (50 km) ocean model output. While random walk dispersion and Markov-1 produce similar particle distributions over time when using our ocean model output, Markov-1 yields Lagrangian trajectories that better resemble trajectories from eddy-resolving simulations. Markov-1 also yields a smaller spurious dianeutral flux.