Sea ice loss drives a regime shift in Arctic Ocean nitrogen biogeochemistry | Communications Earth & Environment

Datasets used in this study

In this study we collate repeated sections across the Fram Strait (78.83–79.5°N) of concurrent temperature, salinity and nutrient (nitrate, phosphate and silicate) measurements collected between 1998 and 2023 from April to September (Table S3; Fig. S4). Note that silicate concentrations in 2008 (iAOOS) were not measured and those in 2018 (FS2018) were found to be unreliable due to instrument failure. Data were obtained from PANGAEA, NOAA, the British Oceanographic Data Centre (BODC) and directly from researchers involved in the Fram Strait Arctic Outflow Observatory, a long-term Fram Strait monitoring programme maintained by the Norwegian Polar Institute. The repeated sections allowed direct comparison of hydrographic and biogeochemical changes throughout the years. All nutrient data have been quality controlled, and their accuracy meets publication standards. Across the time series, nutrient measurements have been made by a range of different laboratories (Table S3) using standardised protocols and certified materials. The drop in average nitrate concentrations in the polar surface water does not correspond to a specific shift in methodologies. For example, the consistent low-nitrate concentrations detected from 2009 onwards have been measured by four different laboratories some of which were also responsible for the earlier measurements10. Nutrient samples were analysed with an accuracy of 0.05 μM, 0.01 μM and 0.18 μM for average nitrate, phosphate and silicate, respectively, across all cruises9,10,58,81,82.

Water mass classification

This study focuses on biogeochemical changes in the cold, fresh Arctic-sourced outflowing waters. These waters are typically classified as Polar Surface Waters (PSW) using density (σθ ≤ 27.7 kg m−3) and temperature (θ ≤ 0 °C) constraints83,84,85,86. To isolate it from local changes in AW intrusion and mixing, which can obscure upstream changes, this study uses a lower density constraint of σθ ≤ 26.9 kg m−3 but will continue to refer to this water mass as PSW for simplicity. Note that this difference limits PSW depth extent from 110 m to 70 m (Fig. S5), allowing for the purest form of PSW to be identified (Fig. S6). One exception is cruise JR17005; in 2018, Atlantic water intrusion was stronger than usual (Fig. 2d), so a lower density constraint was used (σθ ≤ 25.8 kg m−3). In addition, the top 10 m was excluded from all cruises to avoid the influence of sea ice surface melting and local biological uptake. To avoid Greenland meltwater influence, PSW westwards of 16°W were excluded.

Nutrient time series

The challenge of any time series biogeochemical study is in deriving annually representative values from seasonal and spatial variations (both in the vertical and horizontal space) associated with the signal. Here, we use calculated mean averages and range of values on a per cruise basis. PSW is identified by the physical properties described (temperature and density), and the data is grouped to calculate mean averages and the ranges as one standard deviation per cruise. In the time series figures, the mean averages are presented with markers and the range of values associated with a specific time frame of change (pre and post-regime shift) as shaded envelopes. These envelopes represent the mean plus one standard deviation of cruises in pre and post-regime shift (see Tables S1and S2). The standard errors are also estimated for all nutrient time series and are represented as error bars.

Change-point analysis

In order to statistically identify regime shifts in the biogeochemical properties of PSW and physical features of Arctic circulation from 1998 to 2023, we use the R package “Detection of Structural Changes in Climate and Environment Time Series” (EnvCPT87,88). Within the EnvCPT package, a number of predefined linear models are fit to a given time series. We include 8 statistical models in our analysis: a constant mean, a constant linear trend, a constant mean with first-order autocorrelation, a constant linear trend with first-order autocorrelation, multiple change points in the mean, multiple change points in the linear trend, multiple change points in the mean with first-order autocorrelation, and multiple change points in the linear trend with first-order autocorrelation. Subsequently, all model fits are ranked with the Akaike Information Criterion (AIC), which combines the maximum likelihood estimate for each statistical model. Thereby, the statistical model the lowest AIC score corresponds to the model that is best suited to describe the time series in question. We apply his methodology to the time series of (i) nitrate concentrations (Fig.2a), (ii) N:P ratios (Fig. 2b), (iii) Si:N ratios (Fig. 2c), (iv) stratification (Fig. 2d) and (v) total benthic denitrification signal at the time of PSW outflow from the shelves (Fig. 3c). For each of these quantities, the model including multiple change points in the linear trend is the best to describe the time series (Fig. S7).

Time series of stratification

The strength of stratification was obtained by calculating the difference in potential density between 10 and 200 m depth (∆ρ)58. The thermal (∆ρT) and haline (∆ρS) contributions to this density difference, where ∆ρ ∼ ∆ρT + ∆ρS, were estimated as follows:

$$\Delta {\rho }_{T}=-\alpha \rho ({T}_{2}-{T}_{1})$$

(1)

$$\Delta {\rho }_{S}=-\beta \rho ({S}_{2}-{S}_{1})$$

(2)

where α and β are the thermal expansion and haline contraction coefficients, respectively, calculated at the average temperature, salinity, and pressure between 10 and 200 m. Subscripts 1 and 2 refer respectively to surface (10 m) and 200 m depth values of temperature (T) and salinity (S) and ρ is the average potential density between 10 and 200 m.

Denitrification rate calculations

Yearly averages of benthic denitrification rates in the Chukchi, East Siberian, Laptev and Kara Seas from 1998 to 2022 were estimated for two different net primary productivity (NPP) scenarios: (1) constant NPP over time, i.e., control experiment, and (2) steady increase in NPP between 1998 and 2008 followed by a shift to a different rate of steady increase between 2009 and 2022. The latter is based on satellite observations (ref. 5). Changes in primary productivity rates (in TgC yr−1) were used to estimate changes in benthic denitrification rates (in TgN yr−1), assuming that they are directly proportional. This proportionality has been widely observed31,89 and assumed in models21,90,91,92,93,94 to calculate denitrification rates. Sources of error in our benthic denitrification estimates include the propagated uncertainties of the literature-derived NPP estimates (ref. 24) and their rates of change (ref. 5) on which our BD calculations are based. Satellite-based NPP observations have an uncertainty of 35%29 and neglect the contribution of under-ice blooms, which have been observed in the Chukchi sea95,96. Consequently, satellite-based NPP, and hence BD estimates, on this shelf may be underestimated, although this effect is likely limited because the marginal ice zone (and not under-ice) blooms dominate NPP in the southern, shallower portion of the Chukchi shelf97 where BD exerts the strongest control. Also, the trend in NPP increase (from ref. 5), and by extension the BD trends derived from it, remains robust given that under-ice NPP has not increased over time97.

Constant denitrification rate scenario

Benthic denitrification rates in the “control” scenario – where primary productivity rates over the shelves do not change over time – were assumed to stay constant at values in ref. 21. Following this study, we calculated the average yearly denitrification rate over each shelf, x, using the rate-depth relationship shown in Eq. (3), where shelf depths for Chukchi, East Siberian, Laptev and Kara Seas are 80, 58, 48 and 131 m, respectively. This exponential relationship was derived from in-situ denitrification rate measurements during spring 2004 in the Chukchi Sea and is assumed to hold for the other shelves21. Denitrification rate measurements were based on the downward diffusive flux of nitrate from sediment cores98 and corrected using whole core incubations where the flux of nitrogen gas out of the sediments was measured21.

$$\bar{{{denitrification}\,{rate}}_{x}}=60.42\times {\left({depth}\right)}^{-0.944}$$

(3)

$$\bar{{{denitrification}\,{rate}}_{x}}=\bar{{{denitrification}\,{rate}}_{{Chukchi}}}{{{\rm{x}}}}\left(\frac{N{{PP}}_{x}}{{{NPP}}_{{Chukchi}}}\right)$$

(4)

The average denitrification rates were then scaled by the relative proportion of the primary production at each shelf x to the Chukchi Sea. This is expressed in Eq. (4), where NPPx is the primary productivity of shelf x and NPPChukchi is the primary productivity of the Chukchi Sea. Primary productivity rates for Chukchi, East Siberian, Laptev and Kara Seas in 2004 used were 18.4, 7.9 ± 2.0, 4.0 ± 1.0, 8.8 ± 2.3 Tg C yr−121,24. Resulting denitrification rates were 3.1, 1.8, 1.1 and 0.9 Tg N yr−1, respectively (Table S4). The direct proportionality between these primary productivity and denitrification rates in each shelf is carried forward in the second scenario of this study (see 2.2).

Steady denitrification rate increased from 1998 to 2008 and 2009 to 2019

Benthic denitrification rates in the scenario where primary productivity rates (NPP) in the shelves increase at the same steady rate between 1998 and 2008 and then shift to a different steady state between 2009 and 2018 were derived using reported NPP rates of change (in change year−1, see ref. 5) for the two periods (Table S5) and taking primary productivity rates observed in 2004 as fixed reference values (Table S4, see ref. 21). Note that as NPP rates of change were only measured until 2018, rates beyond this were obtained by linear extrapolation.

The equation of a line representing the rate of change in primary productivity (NPP) was found by solving a system of six linear Eqs. (5)-(10). The unknowns in the equation are: NPP in 2008 (NPP2008), NPP in 1998 (NPP1998), gradient (m) and intercept (c) of the line. The resulting equation would have the form y = mx + c, where y is the NPP rate for a given year (x).

$${{{{\rm{NPP}}}}}_{1998}+({\Delta {{{\rm{NPP}}}}}_{1998-2008}( \% )/100* {{{{\rm{NPP}}}}}_{1998})={{{{\rm{P}}}}}_{2008}$$

(5)

$${{{{\rm{NPP}}}}}_{1998}

(6)

$${{{{\rm{NPP}}}}}_{2008} > {{{{\rm{NPP}}}}}_{2004{{\_{{{\rm{C}}}\& D}}}}$$

(7)

$${{{\rm{m}}}}=({{{{\rm{NPP}}}}}_{2008}-{{{{\rm{NPP}}}}}_{1998})/\left(2008-1998\right)$$

(8)

$${{{\rm{c}}}}={{{{\rm{NPP}}}}}_{1998}-{{{\rm{m}}}}* 1998$$

(9)

$${{{\rm{m}}}}* 2004+{{{\rm{c}}}}={{{{\rm{NPP}}}}}_{2004{{\_{{{\rm{C}}}\& D}}}.}$$

(10)

Equation (5) assumes that the NPP in 2008 (NPP2008) is the same as NPP1998 plus its percentage change between 1998 and 2008. The change in NPP between 1998 and 2008 (∆NPP1998–2008 in %) is not directly available from the literature. Instead, it was derived by ∆NPP1998–2008 (%) = ((NPP2008 – NPP1998)/ NPP1998) *100 where NPP2008 = NPP1998 + (\(\Delta\)NPP1998–2008(years) * 10 yrs) and \(\Delta\)NPP1998–2008(years) is reported in the literature (in ref. 5; see Table S5). Equation (6) assumes that NPP in 1998 (NPP1998) is less than the NPP measured in 2004 (NPP2004_C&D) (in ref. 21; see Table S4). Equation (7) assumes that NPP in 2008 (NPP2008) is more than the NPP measured in 2004 (NPP2004_C&D) (in ref. 21; see Table S4). Equation (8) states that the gradient (m) of the line is defined as the differences in NPP from 1998 to 2008. Equation (9) states that the intercept of the line (c) is defined as NPP1998 – m * 1998. Equation (10) states that the NPP rate in 2004 should be the same as that in ref. 21 (see Table S4).

Once NPP2008 is solved for, the NPP rate in remaining years is found by adding the reported 2009–2018 rate of change per year in ref. 5 (see Table S5).

The estimated NPP for a given shelf (x) and year (t) \(\left(N{{PP}}_{x,t}\right)\) is then used to calculate the average benthic denitrification rate \((\bar{{{denitrification}\,{rate}}_{x,{t}}})\) by scaling the benthic denitrification rate at that shelf in 2004 (from ref. 21; see section 2.1) by the relative proportion of the primary production at shelf x at the time of interest \(\left(N{{PP}}_{x,t}\right)\) to the the primary production at shelf x in 2004 (\({{NPP}}_{x,\,2004}\); ref. 24), as expressed in Eq. (11) below.

$$\bar{{{denitrification}\,{rate}}_{x,\,t}}=\bar{{{denitrification}\,{rate}}_{x,\,2004}}{{{\rm{x}}}}\left(\frac{N{{PP}}_{x,t}}{{{NPP}}_{x,\,2004}}\right)$$

(11)

Monte Carlo uncertainty analysis of BD estimates

To assess whether uncertainties in NPP could affect the inferred temporal trends in benthic denitrification, we performed a Monte Carlo uncertainty propagation focused on the Chukchi Sea shelf, where the post-2008 increase in BD is most pronounced. Specifically, we use the Monte Carlo uncertainty propagation to address: (1) whether plausible uncertainties in NPP could alter the sign of the inferred post-2008 increase in BD in the Chukchi Sea, and (2) whether the post-2008 BD trend represents a strengthening relative to the pre-2008 trend. Annual NPP values were randomly varied by ± 35% using a Gaussian error distribution, consistent with reported uncertainties in satellite-derived NPP estimates relative to in situ measurements29. The NPP time series was re-anchored to the observed 2004 value (18.4 Tg C yr⁻¹)21,24, and BD was recalculated for each realisation assuming the same linear NPP–BD relationship used throughout this study.

Given that satellite-derived NPP errors may be similar from year to year or vary independently between years, we tested the two limiting cases: fully correlated errors applied uniformly across all years (representing systematic bias, i.e., algorithmic bias) and independent year-to-year errors (representing interannual variability, e.g., from cloud cover, light/ice conditions). For each case, 20,000 Monte Carlo realisations were generated, and BD trends were re-estimated for the pre-2008 (1998–2008) and post-2008 (2009–2018) periods. All simulations were performed using a fixed random number initialisation to ensure reproducibility. This analysis was designed to test the robustness of inferred BD trends to plausible NPP uncertainties rather than to provide precise probabilistic error bounds on absolute BD rates.

Design of the Lagrangian experiments

Lagrangian experiments were computed to analyse pathways of water mass transport by tracking virtual “particles” through time and space. Specifically, a backward approach was used to (1) trace where the Polar Surface Waters outflowing in the western Fram Strait come from, and to (2) evaluate how contributions from different source regions – particularly Siberian shelves – have changed over time. For this purpose, the open-source and open-access PARCELS (“Probably A Really Computationally Efficient Lagrangian Simulator”, v2.499) particle-tracking software package was used. This framework for computing Lagrangian particle trajectories can be found at http://www.oceanparcels.org. It works by advecting “virtual particles” along ocean currents by physical dynamical laws. We chose to use ocean velocities from the NEMO general ocean circulation model (Nucleus for European Modelling of the Ocean37). It allowed us to derive multi-annual trajectories, and sample other physical properties (i.e. sea-ice concentration, salinity, temperature) along the pathway, with a daily frequency. In particular, we used the version of NEMO corresponding to the Global Ocean Reanalysis and Simulations 2 version 4 (GLORYS2V4) with a 0.25° × 0.25° resolution in the Arctic Ocean. The dataset covers the time period from 1993 to 2022 and is a re-analysed product constrained by satellite and in-situ observations to minimise and correct for model bias. NEMO has been extensively validated throughout the Arctic and consistently reproduces sea ice coverage100, mixed layer depth101 and stratification102 observations. The velocity field dataset (“Global Ocean Ensemble Physics Reanalysis” product ID GLOBAL_MULTIYEAR_PHY_ENS_001_031) was compiled from the Copernicus Marine Environment Monitoring Service (http://marine.copernicus.eu) maintained by Mercator Ocean International, which develops, quality checks and calibrates the product following standard procedures established by the modelling community.

For each Lagrangian experiment, a total of 729 particles uniformly distributed along the Fram Strait (78 to 80 °N and −20 to 0° E) were initialised. This release region corresponds to the western Fram Strait and is designed to enclose the sampling locations where PSW was observed by in-situ hydrographical data throughout the time period of interest (Fig. S4). Particles were tracked for 10 years with positions recorded every day. The choice of a 10-year trajectory assumes that within this time frame, the particle will have completed its trajectory from its point of origin (or the outer bounds of Arctic Ocean) to the Fram Strait, thereby encapsulating its entire path evolution since formation.

“Releases” of particles took place every 31st of July between 2003 and 2022 with the same initial 729 particle grid used for each model run. The time of release was chosen for consistency with available time series of observations. In the first model run, particles were released in Fram Strait on 31st July 2003, and their daily locations in the preceding ten years (back to 1993) were outputted. In the last model run, particles were released on 31st July 2022, and source contributions were estimated from daily 2017 to 2022 trajectories (details on this in the “post-processing” section below).

Particles were initialised at 50 m depth and then advected via ocean currents in two dimensions (zonal and meridional directions). The model depth of 50 m was chosen as: (1) it is below the mixed layer across basins and seasonal variability (∼5–30 m in summer and ∼25 to 50 m in winter101), (2) it is within the Arctic upper halocline (which generally extends up to 80 m depth) and captures the transpolar drift core depth of 50 m33 and (3) PSW core depth is also at around 50 m according to this study’s Fram Strait time series of observations (Fig. S5). The bathymetry from IBCAO (International Bathymetric Chart of the Arctic Ocean103), as well as modelled temperature and salinity were sampled along trajectories.

Post-processing of particle trajectories

The Arctic Ocean was split into segments—as shown in Fig. 1 —to aid identification of Polar Surface Waters and quantification of shelf contributions to PSW via estimation of residence time along each of the Siberian shelves. This allowed deriving a benthic denitrification signal per particle upon outflow from the shelves. Such tasks were completed during post-processing of the Lagrangian model output, ultimately aiming to assess the role of changing circulation and benthic denitrification rates in driving the regime shift in PSW nitrate concentrations and nutrient stoichiometry.

PSW identification

In order to exclusively consider PSW (i.e., southwards flowing, Arctic-derived waters), modelled temperatures sampled along trajectories were used. Particles with temperatures of over 0 °C at time of arrival in Fram Strait were discarded as Atlantic Waters and masked. The percentage of PSW particles (out of all particles released) did not change over time (~81%) – as expected given that particles at 50 m should not be affected by Atlantification – thus supporting the validity of conditions above.

Siberian shelf contributions to PSW

PSW particles that were sourced in the Siberian shelves (Chukchi, East Siberian, Laptev and Kara seas) are referred to as Siberian shelf particles and satisfied the following condition: at any point during its trajectory, the particle crossed the shelf environment (depth less than 100 m) within the shaded shelf bounds.

To estimate residence time of Siberian shelf particles in each of the shelves, the following additional geographical limits were used (based on ref. 4). Note that residence time is defined as days spent by particles within the bounds of shelf of interest.

Chukchi Sea: South of 81.8 °N, between 155.5 °W and 180.0 °W.

East Siberian Sea: South of 81.8 °N, between 180 °W and 148.0 °E.

Laptev Sea: South of 81.8 °N, between 103.6 °E and 148.0 °E

Kara Sea: South of 81.8 °N, between 103.6 °E and 66.3 °E.

It is worth noting that if a particle reached Fram Strait from its shelf region of formation in less than 10 years (maximum length of the advection scheme), the particle would seem stationary in the shelf environment rather than “disappearing”. To correct for this Lagrangian model artifact, the time during which the particle was stationary was excluded from shelf residence time. Moreover, for consistency, particles were only tracked for a maximum of three years on the shelves. The same was applied when deriving benthic denitrification signal (see 3.2.3). Three years was selected as this is the average total time spent on the shelves throughout the time series for ten-year particle tracking (Fig. 3a).

The portion of total shelf residence time that each particle spent in each shelf (expressed as a “Shelf contribution to PSW” as a percentage) was then calculated for each model run as a measure of changes in particle pathways.

Modelled benthic denitrification signal

Only Siberian shelf particles carried a benthic denitrification signal, as they are exposed to strong nitrate losses as a result of denitrification. Given that benthic denitrification rates are different on each of the shelves (i.e., Chukchi, East Siberian, Laptev and Kara Sea shelves), a total denitrification signal upon outflow from the shelves is estimated in Tg per year to explore the influence of changes in contribution from different shelves from circulation changes on PSW N loss. For this, a BD signal is assigned to each shelf-influenced particle based on particle residence time on each of the shelves and estimated denitrification rates on the shelves during time of shelf transit. For instance, a particle that spent its whole shelf residence time in the Kara shelf would carry the benthic denitrification signal of this region (i.e., the benthic denitrification rate of this shelf at the year of transit, as shown in Fig. 2e). In most cases, particles shared their shelf residence time across different shelves. Thus, each shelf signal is estimated separately, see Eq. (12–15), based on BD rates on each shelf of interest at time of transit and adding them according to their relative contributions to PSW, using total shelf residence time (ca. 3 years) as a normalisation factor, as shown in Eq. (16).

In the equations below, RT corresponds to residence time, RTTotal is the total time spent on the Siberian shelves and DR is (benthic) denitrification rate.

$${{{\rm{Chukchi\; shelf}}}} ({{{\rm{CS}}}}) {{{\rm{BD\; signal}}}}=\left({{{{\rm{RT}}}}}_{{{{\rm{CS}}}}}/{{{{\rm{RT}}}}}_{{{{\rm{Total}}}}}\right)* {{{{\rm{DR}}}}}_{{{{\rm{CS}}}}}$$

(12)

$${{{\rm{East\; Siberian}}}}({{{\rm{ES}}}}){{{\rm{BD\; signal}}}}=\left({{{{\rm{RT}}}}}_{{{{\rm{ES}}}}}/\,{{{{\rm{RT}}}}}_{{{{\rm{Total}}}}}\right)* {{{{\rm{DR}}}}}_{{{{\rm{ES}}}}}$$

(13)

$${{{\rm{Laptev\; Sea}}}}({{{\rm{LS}}}}){{{\rm{BD\; signal}}}}=\left({{{{\rm{RT}}}}}_{{{{\rm{LS}}}}}/{{{{\rm{RT}}}}}_{{{{\rm{Total}}}}}\right)* {{{{\rm{DR}}}}}_{{{{\rm{LS}}}}}$$

(14)

$${{{\rm{Kara\; Sea}}}}({{{\rm{KS}}}}){{{\rm{BD\; signal}}}}=\left({{{{\rm{RT}}}}}_{{{{\rm{KS}}}}}/{{{{\rm{RT}}}}}_{{{{\rm{Total}}}}}\right)* {{{{\rm{DR}}}}}_{{{{\rm{KS}}}}.}$$

(15)

$${{{\rm{Total}}}}\; {{{\rm{benthic}}}}\; {{{\rm{denitrification}}}}\; {{{\rm{signal}}}}= {{{\rm{BD}}}}\; {{{\rm{signal}}}}\; {{{\rm{in}}}}\; {{{\rm{CS}}}}+{{{\rm{BD}}}}\; {{{\rm{signal}}}}\; {{{\rm{in}}}}\; {{{\rm{ES}}}} \\ +{{{\rm{BD}}}}\; {{{\rm{signal}}}}\; {{{\rm{in}}}}\; {{{\rm{LS}}}}+{{{\rm{BD}}}}\; {{{\rm{signal}}}}\; {{{\rm{in}}}}\; {{{\rm{KS}}}}$$

(16)

The rate of benthic denitrification (DR) used varied according to the modelling scenario considered, namely (1) constant denitrification, and (2) steady increase from 1998–2008 and shift to a different steady increase from 2009 to 2022. In scenario 2, a temporal averaging method was employed to account for varying denitrification rates across different years. This was necessary as particles often travelled through a given shelf across years. BD signal was calculated as the weighted average of the yearly denitrification rates, where the weights corresponded to the proportional duration of residence at the shelf of interest in each year. Moreover, each shelf contribution to the BD signal was also calculated as a percentage to distinguish whether regions driving change varied over time. It is worth noting that no NPP rate estimates are available before 19985.

As a benthic denitrification signal was derived for each Siberian shelf particle, a mean total BD signal and standard deviation (σ) could be estimated by considering the total amount of Siberian shelf particles for each model run. The mean anomaly ± σ was plotted at the start of the model run time and encapsulates particle trajectory in the ten years since then.

Time series of model velocities

To better untangle the intrinsic role of NEMO large-scale circulation in explaining the shift in BD signal in Lagrangian model results, speed anomaly and speed linear trends maps were constructed to compare pre and post shift variability. This was possible using the meridional and zonal velocity daily averages from NEMO (u and v, respectively) and converting the u and v to mean current speed |U| for the given period of interest using the equation \(\sqrt{{u}^{2}+{v}^{2}}\). The |U| anomaly was then constructed by subtracting mean current speed |U| for the two periods. Linear trends were estimated using linear least-squares regression.

Influence of changes in shelf contribution

To investigate the relative importance of circulation changes on the benthic dentification signal carried by PSW, a “control” Lagrangian experiment was implemented. This experiment consisted in deriving a BD signal (following method in 3.2.3) assuming constant benthic denitrification rates on the shelves. Specifically, the BD rates reported in 2004 (from ref. 21; see Table S4) are assumed for the whole time series. Its resultant contribution to the total BD signal (which takes into account changes in BD on the shelves) can be derived as a percentage following \(\frac{{BD\; signal\; derived\; using\; constant\; BD\; rates}}{{BD\; signal\; derived\; using\; changing\; BD\; rates}}{{{\rm{x}}}}100\). The resulting percentage was used to construct the two shades of the envelope shown in Fig. 3c. The two sets of BD signals used to derive the percentages, namely one derived using constant BD rates and the other using changing BD rates, are shown in Fig. S8 (in grey and red, respectively).

More importantly, to estimate how much of the increase in the total BD signal can be attributed exclusively to changing contribution of shelf waters, linear regressions were derived for estimated BD signals in the two described scenarios (Fig. S8). Linear regressions are constructed for the period after 2004, when the BD signal from both scenarios diverge. By comparing their respective gradients, the influence of changes in shelf contribution due to faster surface currents to the change in BD signal can be estimated (i.e., 0.05/0.15 ×100 = 33%; Fig. S8).