1 Introduction

Accurate prediction of the future of the Earth’s climate relies on our understanding of ocean circulation because of its role in governing air-sea exchange and storage of heat and carbon. Transient ocean tracers can provide unique information about ocean circulation as they are sensitive to various processes integrated over a range of timescales.

Radiocarbon (14C) in dissolved inorganic carbon is a tracer with both natural and anthropogenic influences. 14C is produced naturally in the atmosphere by the interaction of nitrogen atoms and cosmic rays and it decays with a 5700-year half-life. Fossil fuel emissions beginning from the Industrial Revolution reduce the ratio 14C /C (∆14C) in atmospheric CO2 by adding CO2 with no 14C. Then, nuclear weapon testing in the 1950s and 1960s abruptly increased 14C in the atmosphere (Rafter and Fergusson 1957; Nydal 1963).

Observations of bomb-produced 14C in the ocean from the 1970s-1990s have provided strong constraints on the rate of air-sea gas exchange (Broecker et al. 1985; Sweeney et al. 2007). After entering the ocean via gas exchange, the ∆14C distribution in the interior ocean is governed by advection–diffusion processes and radioactive decay. Radiocarbon has been used for many decades to understand ocean ventilation over a range of timescales corresponding to both natural and bomb radiocarbon, and to evaluate model representation of transport and water mass formation processes (Jain et al. 1995; England and Maier-Reimer 2001; Duffy et al. 1995; Guilderson et al. 2000; Gnanadesikan et al. 2004; Matsumoto and Key 2004; Grumet et al. 2005; Galbraith et al. 2011). Air-sea exchange was the dominant control on the rate of oceanic uptake of bomb-produced 14C until the 1990s, but then the dominant control shifted to shallow-to-deep transport and mixing, which also regulate the uptake of anthropogenic CO2 and are important for heat uptake (Graven et al. 2012; Khatiwala et al. 2013; Winton et al. 2013; DeVries et al. 2017).

Previous studies using radiocarbon to study recent ocean circulation have focused on the steady-state or average ocean circulation and have not examined variability in ocean circulation on decadal timescales. Tracer-based studies of recent ocean circulation variability have mainly relied on chlorofluorocarbons, which are entirely man-made and have been measured in the ocean since the 1970s. These studies have revealed variations in upper-ocean ventilation, which have impacts on ocean CO2 uptake, although attributing the variability to internal or external drivers is challenging (Waugh et al. 2013; DeVries et al. 2017; Lester et al. 2020).

Radiocarbon may provide unique information about ocean circulation on decadal timescales, however, the relationships between recent ocean circulation and 14C variability has not yet been explored by global ocean models. In this study, we simulate ∆14C in the ocean and examine how it is affected by circulation changes. We describe the implementation of radiocarbon tracer in the ocean model (Nucleus for European Modelling of the Ocean (NEMO)) and evaluate the performance of the model in simulating ∆14C by comparing it to observations. We then analyze the differences in ∆14C when the model is driven with fixed or time-varying atmospheric forcings and we suggest some potential mechanisms to explain the impact of variable forcing, although we do not perform a detailed analysis of the specific ocean circulation changes in the model. We use the simulations to explore the sensitivity of ∆14C to ocean circulation changes, and to assess whether ∆14C observations would likely be useful for data-driven detection of ocean circulation change.

2 Data and method

2.1 Ocean abiotic radiocarbon modelling

We used the global ocean circulation model Nucleus for European Modelling of the Ocean (NEMO) version 3.6 (Madec 2014), which is coupled to the Los Alamos sea-ice model (CICE, Rae et al. (2015)) with ∼ 1 degree resolution. Vertically, the full depth water column is divided into 75 levels with thickness varying from 1 m (top 10 m depth) to ∼ 200 m at depths below 3500 m. A detailed description of the model set-up and evaluation of its physical dynamics are given in Storkey et al. (2018), and briefly in Marzocchi et al. (2021). Therefore, here, we only provide a description of the new incorporation of 14C in the passive tracer module (Tracers in Ocean Paradigm, TOP) in NEMO.

Following the Ocean Model Intercomparison Project (Orr et al. 1999, 2017), radiocarbon is simulated using two prognostic tracers: dissolved inorganic carbon (DIC) and radiocarbon in DIC (DI14C). We use the “abiotic” formulation, which neglects any changes in dissolved inorganic carbon components due to biological activities and ignores mass-dependent isotopic fractionation during air-sea exchange (Toggweiler et al. 1989). ∆14C in the ocean is then calculated as 1000(DI14C/DIC − 1), with units in ‰. This quantity can be directly compared to the observational data reported as ∆14C since biological effects are presumed to be small (Joos et al. 1997).

Air-sea gas exchange is calculated with a quadratic wind speed parameterization of the piston velocity with a coefficient of 0.257 cm/h m2/s2 from Wanninkhof (2014), scaled by fractional ice cover. The alkalinity, for the estimation of the aqueous concentration of CO2, is computed as a function of salinity, and the concentration of phosphate and silicate are set constant (SiO2 = 7.7e−3 mol/m3 and PO4 = 5.1e−4 mol/m3).

2.2 Simulations

2.2.1 Spinup simulation

We ran a long spin-up simulation to allow the model physics and ocean carbon cycle to reach a steady state, even in the deep ocean. The boundary conditions for atmospheric CO2 and ∆14C were prescribed to 1850 values of 284 ppm and ∼ − 4‰ respectively (Meinshausen et al. 2017; Graven et al. 2017). ∆14C in the ocean was initialized at 0‰ and DIC was initialized from OMIP fields based on observations (http://omip-bgc.lsce.ipsl.fr/). The spin up was forced by repeating JRA-55-do neutral year (JRA-NY, Stewart et al. 2020) atmospheric data, which is a one year span (01-May-1990 to 30-April-1991) from JRA-55-do (version 1.3) data when major climate modes such as the North Atlantic Oscillation, Southern Oscillation and Southern Annular Mode were largely neutral. The model was integrated for 9000 years, when the drift in ∆14C was less than 0.001‰/year over more than 99% of the oceanic volume, satisfying the OMIP criterion for steady-state (Orr et al. 2017). Temperature and salinity are also spun up after 9000 years so they do not have any trend or inter-annual variation. The DIC, DI14C, temperature and salinity fields at the end of the 9000 years spin-up were then used as initial conditions to perform historical simulations.

2.2.2 Historical simulations

We ran two simulations with historical atmospheric CO2 and ∆14C boundary conditions. The CO2 concentrations for 1850 to 2014 are from (Meinshausen et al. 2017) and the data over 2015 to 2020 are from NOAA’s Global Monitoring Laboratory https://gml.noaa.gov/ccgg/trends/global.html. Atmospheric ∆14C is from Graven et al. (2017) for the period 1850 to 2015 and from the simulated transient ∆14C for SSP2-4.5 over 2016–2020 (Graven et al. 2020). Global mean values are used for the CO2 boundary condition while three zonal bands (90° S-30° S,30° S-30° N,30° N-90° N) are used for ∆14C. One simulation ("Hist") was forced with inter-annually varying JRA-55-do atmospheric reanalysis data (version 1.5, Tsujino et al. (2018)), which is available for the period 1958 through 2020 at 3-hour intervals. The JRA-55-do data is cycled over 1850 to 2020, where the 1850 atmospheric forcing corresponds to the year 1975. The other simulation ("Hist-NYF") was forced with the JRA-NY atmospheric data (‘NY’ stands for Neutral Year and ‘NYF’ means Neutral Year Forced) in every year. This includes seasonal variation in the forcing data but no inter-annual variability.

Then, we ran two simulations with fixed atmospheric CO2 and ∆14C boundary conditions, using 1850 values. Again, one simulation ("Fix") was forced with inter-annually varying JRA-55-do atmospheric reanalysis data and the other ("Fix-NYF") was forced with JRA-NY atmospheric data.

The "Hist" simulation represents the real world ocean while the "Fix-NYF" simulation is a control simulation continuing the spin-up simulation, and the "Hist-NYF" and "Fix" simulations allow us to examine the effect of changing circulation on total ∆14C and on natural ∆14C.

2.3 Ocean observations

We use hydrographic survey data for ∆14C from the Global Ocean Data Analysis Project (GLODAP v2.2020) data inventory (Key et al 2004; Olsen et al. 2020) to evaluate the model. We focus on two meridional sections, one in the Atlantic along roughly 30° W (A16) and the other in the Pacific Ocean along roughly 150° W (P16). We compare the simulated and observed changes in ∆14C over three decades (year and months of sampling can be found in Table 1). We bin the observations by 2° in latitude and discrete vertical depths over the years the 1990s, 2000s and 2010s available in Table 1 and linearly interpolated ∆14C for bins with no data. In the results, figures related to observational data show colour for each bin and have a blocky appearance, whereas figures related to only model output show smooth contours.

Table 1 Date of ∆14C repeat sections available in GLODAP and used for decadal comparisons
Full size table

3 Results

3.1 Simulated bomb 14C inventory and air-sea fluxes

First, we evaluate the globally-integrated oceanic inventories of bomb 14C, defined as the excess radiocarbon relative to the year 1950 (Broecker et al. 1985), in the NEMO model. The simulated inventories are consistent with the observation-based estimates of Naegler (2009) for the GEOSECS data from the 1970s and the WOCE data from the 1990s (Fig. 1a). The inventory stops increasing around 2012 and starts to decrease. The time at which the 14C inventory peaks depends on the efficiency of mixing processes and varies between models (Graven et al. 2012). For example, the CCSM model inventory showed an earlier peak around 2000, while the ECCO model inventory peaks a bit later in 2015–17 (Khatiwala et al. 2018). Comparing the Hist and Hist-NYF simulations, the Hist-NYF inventory is slightly higher for the period 1980–2020.

Fig. 1
figure 1

(a) Globally integrated oceanic radiocarbon inventory relative to the year 1950, (b) zonally integrated air-to-sea radiocarbon flux in 5 degree latitude bins for the Hist simulation, and (c) radiocarbon flux anomaly between Hist and Hist-NYF simulations

Full size image

Patterns of net flux of 14C into the ocean are similar to other models (Fig. 1b; Graven et al. 2012). There is a dramatic increase in 14C flux after the bomb testing in the 1950s and 1960s, with the largest fluxes in the Southern Ocean. The net air-sea flux changed its sign to negative in the low latitudes by the mid-1990s in response to ocean 14C uptake and decreasing atmospheric ∆14C that caused the ocean to have higher ∆14C than the atmosphere. The area of 14C outgassing extends progressively poleward with time. However, the Southern Ocean remains a strong sink for 14C.

Comparing the Hist and Hist-NYF simulations (Fig. 1c) shows that the largest effect of the inter-annual forcing is to decrease the influx of 14C to the Southern Ocean in the 1960s-70s. The effect of the inter-annual forcing is positive south of 65° S and between 25 and 45° latitudes in both hemispheres, and negative north of 45° N. To a large extent the positive or negative differences are sustained through the whole period, however, in the tropics and the Southern Ocean the fluxes switch between higher and lower fluxes with the inter-annual vs neutral year forcing. Sustained differences could suggest the JRA-NY forcing is significantly different from the average forcing over the JRA period 1958–2020 in those regions.

3.2 Modeled ∆14C vs observations

∆14C at the end of the long spin-up is ∼40‰ lower than the estimate of natural radiocarbon in GLODAP (Key et al 2004) at depths below the thermocline in most regions (Fig. 2). Part of the discrepancy can be attributed to a bias in the estimate of natural radiocarbon below the thermocline in GLODAP (Sweeney et al. 2007; DeVries and Holzer 2019). However, ∆14C in NEMO is still too low, indicating that the ventilation of the deep ocean in NEMO is too sluggish.

Fig. 2
figure 2

Basin-wide mean profiles of ∆14 C from our 9000 year spinup simulation and the data provided by OMIP-BGC, based on GLODAP

Full size image

Figures 3 and 4 show the vertical ∆14C profiles for the Hist simulation and the observations along the two meridional sections in the Atlantic and Pacific Oceans. Here, the model data is re-sampled to match the location and time of the observations (see Table 1) and binned in the same way as the observations (Sect. 2.3).

Fig. 3
figure 3

Observed (a–c) and simulated (d–f) ∆14C along the A16 section in the Atlantic Ocean in the three time periods (rows). No observations are available for the North Atlantic in the 2000s. (g–i) represents the bias, modeled—observed ∆14C. The units for all panels are ‰

Full size image
Fig. 4
figure 4

Same as in Fig. 3 but along the P16 section in the Pacific Ocean

Full size image

The observed ∆14C profiles in the Atlantic Ocean (Fig. 3a–c) show high ∆14C in the subtropical gyres down to 500m and in the shallow tropical ocean. In the North Atlantic, high ∆14C penetrates deeply, below 2000m, due to deep water formation while in the Southern Ocean ∆14C is low in the deep ocean up to the surface due to upwelling and limited air-sea equilibration. The profile is similar in the Pacific, but without the deep water formation and high ∆14C in the north (Fig. 4a–c). A decrease in ∆14C from the 1990s to 2010s is observed above 500m in both the Atlantic and the Pacific following the decrease in atmospheric ∆14C.

These features are qualitatively well captured in the model simulation (middle column d-f in Fig. 3 and 4) but there are biases in the model. The model biases are largely negative (>−50‰), especially in deep waters and in the Southern Ocean (right column e–f in Figs. 3 and 4). This negative bias in recent decades is similar to the bias in natural radiocarbon in the spin-up simulation (Fig. 2). Patches of positive bias are present, primarily in the tropics and high latitudes. The largest positive bias (>40‰) is found in the north Pacific at depths shallower than 500m. With negative biases at depth and positive biases in the shallow ocean, NEMO overestimates the vertical gradients in ∆14C in these regions.

We also compare the changes in ∆14C (∆∆14C) over decades in the simulations with the observations. Figure 5 shows the normalized change in ∆14C from the 1990s to the 2010s in observations and in the Hist simulation along A16 and P16. ∆∆14C is negative in shallow water (surface to ∼300 m) and positive in the thermocline depths. The negative values in the surface layer arise from downward transport and from outgassing of 14C while the positive values in the thermocline indicate 14C has been transported from the surface layers and accumulated (Graven et al. 2012). This pattern is well-captured in the model. However, the rate of change is underestimated. In the model, the decreases in shallow layers are too weak and the increases at depth are also too weak, by ∼5‰ per decade for the strongest changes. Underestimates of ∆∆14C may indicate that the model is too efficient at dispersing 14C, similar to the ECCO model (Graven et al. 2012). In contrast, the CCSM model overestimated the magnitude of changes for the 1990s to 2000s in Graven et al. (2012), indicating that CCSM did not transfer enough bomb 14C into deep layers, which may have resulted in part from errors in modeled density. Interestingly, this result seems to contrast with the comparison for natural radiocarbon (Fig. 2), which indicated NEMO was transporting 14C to the deep ocean too slowly, suggesting there may be different biases for decadal-scale vs millennial-scale transport in the model.

Fig. 5
figure 5

Changes in ∆14 C between 1990s and 2010s along A16 (a, b) and P16 (c, d). Observations are shown in the top row and the Hist simulations in the bottom row. Note that the values are normalised to ‰ per decade and the Hist values are re-sampled to match the exact location and time in GLODAP observations

Full size image

3.3 Effect of variable forcing on ∆14C changes along the sections

In this section, we discuss how the radiocarbon distribution is impacted by the change in ocean circulation or gas exchange arising from variable atmospheric forcing. Here we show the distribution along the transects in the model grid itself rather than re-gridding to match observations. We show the average over the years sampled on each section when sections were observed over multiple years (Table 1). To show the effect of variable forcing on total ∆14C changes, we plot the difference: ∆∆14C simulated in Hist minus ∆∆14C simulated in Hist-NYF (∆∆∆14C). So, the difference between Hist and Hist NYF ignores the effect of temporal changes in atmospheric ∆14C and describes how the distribution of ∆14C changes in the inter-annually varying ocean compared to the steady-state ocean. We also consider the effect on natural radiocarbon by plotting ∆∆14C simulated in Fix minus ∆∆14C simulated in Fix-NYF, which is included in the main text. In the supplement, we additionally plot the effect on bomb radiocarbon. We note that the differences in ∆∆14C between approximately decadal hydrographic cruises can reflect variations in forcing and circulation over a range of timescales: inter-annual, decadal and multi-decadal.

In the Pacific, along the P16 transect, the decadal changes in atmospheric forcing caused significant variability in ∆14C distribution in the upper 1000 m water depths (Fig. 6). Potential density also exhibits significant variability (Fig. 6e, f) in the upper 1000 m water depths associated with the changes in atmospheric forcings. Moreover, the differences in ∆∆14C are largely anti-correlated with changes in potential density (Fig. 6e, f). Since ∆14C typically decreases with depth and potential density increases with depth, much of the variation may arise from displacement of isopycnal surfaces.

Fig. 6
figure 6

The effect of variable forcing on changes in ∆14C over decades along P16 in the Pacific Ocean. Difference in ∆∆14C between Fix and Fix-NYF, representing the normalised changes in natural ∆14C (a, b). Difference in ∆∆14C between Hist and Hist-NYF, representing the normalised changes in total ∆14C (c, d). The changes in potential density in the Hist and Hist-NYF simulation (e, f)

Full size image

∆∆∆14C, the difference in the ∆14C change from one decade to the next in the JRA-interannually-forced simulation minus the JRA-NY-forced simulation, exhibits a dipole pattern in the Northern Hemisphere. From the 1990s to 2000s, ∆∆14C with variable forcing is 15‰ decade−1 higher around 20° N at 200-1000 m depths, while around 40° N it is >10‰ decade−1 lower (Fig. 6c). From the 2000s to 2010s, the sign has flipped, at 20° N ∆∆14C is >10‰ decade−1 lower and at 40° N >5‰ decade−1 higher (Fig. 6d). Potential density fields also exhibit similar feature (Fig. 6e, f). This indicates a change of circulation pattern between the two decades as it represents the change in ∆∆14C between dynamic and steady-state simulations and it is reflected in potential density. Change in subtropical gyre subduction related to changes in atmospheric forcings (Schneider et al. 1999; Qu and Chen 2009) could be causing the ∆∆14C and potential density changes. In the Southern Hemisphere, ∆∆∆14C is small for the 1990s-2000s. For the 2000s–2010s, ∆∆∆14C is stronger and there is a pattern of lower ∆∆14C south of 50° S and higher ∆∆14C north of 50° S. This pattern could reflect faster ventilation of Subantarctic Mode Water and faster upwelling of Circumpolar Deep Water, as found by Waugh et al. (2013) using CFCs but in the previous period (1990s-2000s). Both natural (Fig. 6a, b) and total radiocarbon (Fig. 6c, d) show similar patterns in response to changes in atmospheric forcing, showing that natural (Fig. 6a, b) and bomb ∆14C (Supplementary SFig. 1) are affected in similar ways.

The difference in the decadal rate of change of ∆14C due to circulation changes along the A16 transect in the Atlantic (Fig. 7a–d) is different from the Pacific transect. In the Atlantic, changes occur over the entire water column, while in the Pacific significant changes are visible only in the upper 2000m water depths. This is because strong transport in the North Atlantic even in the deep ocean occurs on decadal timescale, while Pacific deep ocean circulation happens on longer timescales. The simulation with the variable atmospheric forcing shows higher ∆∆14C in the North Atlantic north of 40° N, 500–2000m depth, with smaller and more variable differences in other regions. Again, the patterns are similar for natural ∆14C (Fig. 7a, b) and for total ∆14C (Fig. 7c, d), but the primary driver of the signal is bomb 14C.

Fig. 7
figure 7

Same as in Fig. 6 but along A16 in the Atlantic Ocean. The box in the figure points to the correlation ambiguity between ∆∆∆14C and ∆∆σ

Full size image

For the region north of 40° N, 500-2000m, shown in the box in Fig. 7, higher ∆14C is accompanied by both decreases (1990s-2000s) and increases (2000s–2010s) in potential density (Fig. 7e, f). Elsewhere, ∆∆∆14C and potential density changes are anti-correlated, as in the Pacific transect. We will explore the North Atlantic region in the box in Fig. 7 in more detail in the next section.

3.4 North Atlantic ∆14C

Figure 8a shows the time series of the difference between ∆14C in the Hist simulation and ∆14C in the Hist-NYF simulation, averaged in the box shown in Fig. 7 (50° –60° N, 500–1000 m). ∆∆14C (the difference between Hist and Hist-NYF) exhibits a decreasing tendency from 1958 to 1988 with decadal-scale oscillations that peak during 1972–74 and then during 1983–85. After 1988, ∆∆14C shows a monotonic increase without oscillations. It reached ∼−15‰ in 1988 and increased to greater than zero after 2000.

Fig. 8
figure 8

Time series of simulated variables in the North Atlantic, averaged over latitude 50-60° N, longitude 20–40° W, depth 500–1000m (a) anomaly in ∆14C and (b) anomaly in potential density (σ) and AMOC (latitude-depth mean); latitude–longitude box averaged (c) winter mean (January–February–March-JFM) mixed layer depth (MLD) (d) winter mean net heat flux (NHF) and e winter mean wind-stress (TAU)

Full size image

The pattern of ∆14C change is shown in more detail for the vertical column down to 3000m (50°–60° N) in Fig. 9. In the 1970s, high ∆14C is much more concentrated at the surface in Hist than in Hist-NYF. After 2000, ∆14C is higher in Hist than in Hist-NYF down to about 2000 m.

Fig. 9
figure 9

Plot of ∆14C vs time and depth averaged in the box 50–60° N and 20–40° W (north Atlantic) for the (a) Hist simulation and (b) Hist-NYF

Full size image

The ∆∆14C peaks in 1972–74 and 1983–85 align with deep winter mixed layers that transported more bomb 14C into deeper layers (Figs. 8c, 9). The lack of such events after 1988 may result from the weaker vertical gradient in ∆14C (Fig. 9), and/or from less variability in the mixed layer depth (Fig. 8c). The mixed layer depth in Hist-NYF in this region is actually much deeper than the average in Hist (Fig. 8c) which appears to have accentuated the 14C transport into the interior North Atlantic during the period 1960–90 for Hist-NYF as compared to Hist.

The longer-term pattern of decreasing ∆∆14C until 1988 and then increasing ∆∆14C appears to be more in phase with AMOC than with the mixed layer depth, potential density, or net heat flux (Fig. 8). There is a good correlation between ∆∆14C and the difference in AMOC over 1958–2020 (r = −0.53). A stronger AMOC could ventilate ∆14C quickly to and away from the location, possibly reducing the amount of radiocarbon there, whereas a weaker AMOC could result in higher ∆14C because of slow ventilation. In contrast, the longer term trend of mixed layer depth is to increase while ∆∆14C is decreasing, which is the opposite relationship as found for the individual years with high mixed layers that also had high ∆∆14C.

Thus, our hypothesis is that the variability in MLD acts to determine the ∆14C distribution at local and annual scales (when the vertical gradients are large) while AMOC variability impacts the ∆14C distribution on a multi-decadal timescale. The variability in winter MLD is largely modulated by net heat flux (NHF, Fig. 8d) and the wind-stress (TAU, Fig. 8E; Alexander et al. 2000; Pookkandy et al. 2016). Therefore, the short time-scale variability in rate of accumulation of ∆14C in the interior ocean is influenced by the atmospheric conditions at the time of air-sea exchange (Marzocchi et al. 2021).

3.5 Global map of difference in radiocarbon distribution due to circulation changes

Changes in ocean circulation resulting from varying atmospheric forcing cause changes in the ∆14C distribution in the surface and interior ocean. Figure 10a–d shows the mean difference in ∆14C between Hist and Hist-NYF 1991 to 2020 at the surface, 300m, 1000m, and 3000m depths. Figure 10e–h shows the standard deviation of the difference over 1991 to 2020.

Fig. 10
figure 10

At different depths, the mean difference in ∆14C (Hist—Hist-NYF) from 1991 to 2020 (a–d) and the standard deviation of the difference (e–h)

Full size image

The areas of the surface ocean with the largest differences and the largest standard deviations of the differences between Hist and Hist-NYF are upwelling areas such as the eastern Pacific, the subpolar gyres in the North Pacific and North Atlantic, and the Southern Ocean (Fig. 10). These areas are where vertical exchanges of deeper 14C-depleted waters can vary and affect the surface ∆14C values. On average, surface 14C was higher in the Northeast Atlantic, Southeast Pacific and north of the Ross Sea, and lower in the North Pacific and north of the Weddell Sea over 1991–2020 in Hist as compared to Hist-NYF. The differences are up to ±12‰. Differences at 300m are largely similar in pattern to differences at the surface, but with stronger magnitudes of up to ±20‰ or more.

In the North Pacific shallow subpolar waters (<300m), ∆14C in the Hist simulation is 20‰ lower than Hist-NYF (Fig. 10). Over 1991–2020, these shallow ocean regions are decreasing in ∆14C, so the negative value indicates a stronger decrease in Hist due to differences in circulation from Hist-NYF that could be related to stratification, gyre dynamics and/or upwelling.

At ∼1000m depth, the largest differences between Hist and Hist-NYF and variability in Hist are found in the Southwest Pacific and Southern Indian Oceans, with the mean difference showing a dipole structure in the North Atlantic that is opposite to the pattern at 300 m (Fig. 10). The subduction of radiocarbon rich water in the Southern Ocean as mode and intermediate water (SAMW and AAIW) or so called thermocline ventilation are the primary cause of high anthropogenic carbon sequestration in the Southern hemisphere midlatitude intermediate depths (Downes et al. 2009, 2017; Sallee et al. 2010, 2012; Jenkins et al. 2010). Moreover, Marzocchi et al. (2021) have presented evidence (using dye tracer experiments) that dye tracer in these locations is associated with a timescale of more than 20–25 years leading to accumulation over time.Hence, observing variation in radiocarbon in these regions could help to identify and understand ventilation changes in Southern Mode and Intermediate waters (Fig. 6; Waugh et al. 2013).

Radiocarbon variability is small at 3000m depths except for the deep western boundary current of the Atlantic, which shows lower ∆14C in Hist than Hist-NYF. This is consistent with our hypothesis of the effect of slower AMOC-related ventilation on 14C in Hist over 1990–2020 or previous decades, leading to more 14C at mid-depths of 500–1000m in the North Atlantic but less 14C further along the ventilation route in the deep western boundary current. Deeper mixed layers in the North Atlantic in Hist-NYF than in Hist may also have driven stronger transport of 14C into North Atlantic Deep Water.

3.6 Discussion and summary

The NEMO model simulations explored here show how ∆14C in the ocean responds to changes in atmospheric forcing. The simulations are integrated from 1850 to 2020 using JRA-55-do inter-annually varying and neutral-year atmospheric reanalysis data, where the initial boundary fields of ocean DIC, DI14C, dynamic and thermodynamic fields are from a 9000-year spin-up simulation. The model captures the large-scale ∆14C pattern and its change over time along the cruise tracks in the Pacific (P16) and Atlantic (A16) oceans. However, the model is biased low by more than 20‰ in the interior ocean and the magnitude of changes between the 1990s and 2010s is underestimated. The biases are somewhat different to the findings for ECCO and CCSM, where weak temporal changes were associated with high biases in deep ocean ∆14C (Graven et al. 2012), perhaps indicating that biases in ocean circulation can be different for different timescales.

Analyzing the ∆14C changes over decades for the variable and fixed circulation shows a strong large-scale association with changes in potential density, indicative of vertical or horizontal shifts in water distribution. Both natural (Fix-Fix_NY) and contemporary radiocarbon (Hist-Hist_NY) show similar response patterns to decadal changes in atmospheric forcing but a larger magnitude for Hist simulations due to the added nuclear bomb effect. Along the P16 transect, the rate of change in ∆14C due to circulation changes exhibits a dipole pattern, which flips sign over the decades. A dipole pattern in ∆14C was also reported by Galbraith et al. (2011), extending from surface to ∼300m depths, related to changes in gyre circulation and shows some correlation with the Pacific Decadal Oscillation (PDO). While we show the effects of variable forcing here and discuss some potential mechanisms, more detailed analysis is needed to confirm the specific drivers of the ∆14C distribution changes, including the effects of climatic oscillations and specific ocean transport and mixing processes.

The North Atlantic is well known for the convective overturning circulation (AMOC) and the NEMO model simulates ∼5Sv less than the RAPID array estimates at 26.5° N and 1000m (Moat et al. 2020). We hypothesize that the AMOC trend is an important driver of the circulation-driven ∆14C differences on multi-year timescales. Strong AMOC causes rapid ventilation of the absorbed 14C at the air-sea interface leading to quick dissipation from the North Atlantic intermediate depths (Goslar et al. 1995; Stocker and Wright 1996). Recent studies show the importance of AMOC for the accumulation of anthropogenic carbon and other tracers in the North Atlantic (Brown et al. 2021; Romanou et al. 2017). The NEMO simulations also show inter-annual variations in ∆14C related to mixed layer depth reaching 500-1000m depth.

Our study shows how ocean radiocarbon distribution depends on physical properties, such as mixed layer depth and ocean currents which are potentially sensitive to climate variability and change. This response varies among climate models and influences the transient response to ∆14C distributions, so that examining the response of different ocean models and using different atmospheric boundary conditions would help to identify robust patterns and characteristics. We show that key ocean regions for CO2 and heat uptake such as the Northern subpolar gyres and the Southern Mode and Intermediate Water appear to show the strongest signals of circulation change in ∆14C, and these signals may be large (≥10‰) relative to the measurement uncertainty (2–5‰). Therefore, observing and simulating the radiocarbon distribution in the ocean and how it is affected by various ocean dynamics can contribute to better understanding and prediction of circulation impacts on ocean biogeochemistry.