Effects of ²³⁸U variability and physical transport on water column ²³⁴Th downward fluxes in the coastal upwelling system off Peru

The eastern boundary region of the southeastern Pacific Ocean hosts one of the world's most dynamic and productive upwelling systems with an associated oxygen minimum zone (OMZ). The variability in downward export fluxes in this region, with strongly varying surface productivity, upwelling intensities and water column oxygen content, is however poorly understood. Thorium-234 (²³⁴Th) is a powerful tracer to study the dynamics of export fluxes of carbon and other elements, yet intense advection and diffusion in nearshore environments impact the assessment of depth-integrated ²³⁴Th fluxes when not properly evaluated. Here we use vessel-mounted acoustic Doppler current profiler (VmADCP) current velocities, satellite wind speed and in situ microstructure measurements to determine the magnitude of advective and diffusive fluxes over the entire ²³⁴Th flux budget at 25 stations from 11 to 16^∘ S in the Peruvian OMZ. Contrary to findings along the GEOTRACES P16 eastern section, our results showed that weak surface wind speed during our cruises induced low upwelling rates and minimal upwelled ²³⁴Th fluxes, whereas vertical diffusive ²³⁴Th fluxes were important only at a few shallow shelf stations. Horizontal advective and diffusive ²³⁴Th fluxes were negligible because of small alongshore ²³⁴Th gradients. Our data indicated a poor correlation between seawater ²³⁸U activity and salinity. Assuming a linear relationship between the two would lead to significant underestimations of the total ²³⁴Th flux by up to 40 % in our study. Proper evaluation of both physical transport and variability in ²³⁸U activity is thus crucial in coastal ²³⁴Th flux studies. Finally, we showed large temporal variations on ²³⁴Th residence times across the Peruvian upwelling zone and cautioned future carbon export studies to take these temporal variabilities into consideration while evaluating carbon export efficiency.

1 Introduction

Isotopes of thorium (Th) are widely used as tracers for particle cycling in the oceans (Waples et al., 2006). In particular, ²³⁴Th has been extensively used to trace particle dynamics and export fluxes in the upper ocean and to quantify the marine budgets of important macro- and micronutrients such as carbon (C), nitrogen (N), phosphorus (P) and iron (Fe) (e.g., Bhat et al., 1968; Buesseler et al., 1992; Coale and Bruland, 1987; Lee et al., 1998; Le Moigne et al., 2013; Cochran and Masqué, 2003; Van Der Loeff et al., 2006; Black et al., 2019). ²³⁴Th has a relatively short half-life (${\mathit{\tau }}_{\mathrm{1}/\mathrm{2}}=\mathrm{24.1}$ d) that allows studies of biological and physical processes occurring on timescales of days to weeks. Unlike its radioactive parent uranium-238 (²³⁸U, ${\mathit{\tau }}_{\mathrm{1}/\mathrm{2}}=\mathrm{4.47}$ Ga) that is soluble in seawater, ²³⁴Th is highly particle reactive with a particle–water partition coefficient of 10³ to 10⁸ (Santschi et al., 2006, and references therein) and is thus strongly scavenged by particles (Bhat et al., 1968). Generally, a deficit of ²³⁴Th relative to ²³⁸U is observed in the surface ocean and reflects net removal of ²³⁴Th due to particle sinking, whereas secular equilibrium between ²³⁴Th and ²³⁸U is observed for intermediate and deep waters. Integrating this surface ²³⁴Th deficit with depth yields the sinking flux of ²³⁴Th and, if elemental : ²³⁴Th ratios are known, the sinking flux of elements such as C, N, P, Si and trace metals (e.g., Bhat et al., 1968; Buesseler et al., 1998, 1992, 2006; Coale and Bruland, 1987; Weinstein and Moran, 2005; Owens et al., 2015; Black et al., 2019; Puigcorbé et al., 2020).

Various ²³⁴Th models have been put forward to study adsorption–desorption, aggregation and export, but single box models that assume negligible ²³⁴Th fluxes due to physical transport are commonly used to calculate oceanic ²³⁴Th-derived particle fluxes (see detailed review by Savoye et al., 2006). This assumption is typically appropriate in open-ocean settings where ²³⁴Th fluxes due to advection and diffusion are small relative to the downward fluxes of ²³⁴Th associated with particle sinking. However, in upwelling regions such as the equatorial Pacific and coastal systems, advective and diffusive ²³⁴Th fluxes may become increasingly important (e.g., Bacon et al., 1996; Buesseler et al., 1998, 1995; Dunne and Murray, 1999). For example, in the equatorial Pacific, strong upwelling post El Niño could account for ∼50 % of the total ²³⁴Th fluxes (Bacon et al., 1996; Buesseler et al., 1995). Ignoring the upwelling term could thus lead to an underestimation of ²³⁴Th fluxes by a factor of 2. Conversely, horizontal diffusion carrying recently upwelled, ²³⁴Th-replete waters has been shown to balance the upwelled ²³⁴Th fluxes in the central equatorial Pacific (Dunne and Murray, 1999). To the contrary, advective and diffusive ²³⁴Th fluxes were minimal off the Crozet Islands in the Southern Ocean due to limited horizontal ²³⁴Th gradients, long residence time of water masses, and low upwelling rates and diffusivities (Morris et al., 2007).

The dynamic nature of coastal processes requires that physical terms be included in ²³⁴Th flux calculation whenever possible. Accurate measurements of current velocities and diffusivities are however challenging, and thus direct observations of the effects of physical processes on ²³⁴Th distributions in coastal regions are scarce. Limited studies have incorporated advection and diffusion in the nearshore zones of the Arabian Sea (Buesseler et al., 1998), Gulf of Maine (Gustafsson et al., 1998; Benitez-Nelson et al., 2000), South China Sea (Cai et al., 2008) and Peruvian oxygen minimum zone (OMZ) (Black et al., 2018). In the Arabian Sea, coastal upwelling during the southwest monsoon season could account for over 50 % of the total ²³⁴Th flux (Buesseler et al., 1998). Horizontal advection has been shown to be substantial in the inner Casco Bay of the Gulf of Maine (Gustafsson et al., 1998), whereas offshore advection and diffusion are only important in late summer (Benitez-Nelson et al., 2000). Therefore, the importance of physical processes on the ²³⁴Th flux estimate is highly dependent on the seasonal and spatial variability of the current velocities, diffusivities and ²³⁴Th gradients. In terms of the Peruvian OMZ, Black et al. (2018) showed that coastal upwelling accounts for >50 % of total ²³⁴Th fluxes at 12^∘ S; however, how upwelling ²³⁴Th fluxes vary seasonally and spatially in this region is unclear.

Another uncertainty in ²³⁴Th flux calculations in such regions stems from variations on dissolved ²³⁸U activities. Generally speaking, U behaves conservatively under open-ocean oxic conditions and is linearly correlated with salinity (Chen et al., 1986; Ku et al., 1977; Owens et al., 2011). However, numerous studies have shown that such a correlation breaks down in various marine environments including the tropical Atlantic (Owens et al., 2011), Mediterranean Sea (Schmidt and Reyss, 1991) and Arabian Sea (Rengarajan et al., 2003). Although it is generally accepted that deviations from the linear ²³⁸U–salinity correlation will lead to differences in the final calculated ²³⁴Th fluxes, there is currently little knowledge on how significant these differences could be.

In this study, we report vertical profiles of ²³⁴Th and ²³⁸U along four transects perpendicular to the coastline of Peru (i.e., shelf–offshore transects). We evaluate the ²³⁸U–salinity correlation in low-oxygen waters and how deviations from this correlation impact final ²³⁴Th flux estimates. We also assess the spatial and temporal importance of advection and diffusion on ²³⁴Th flux estimates.

2 Sampling and methods

2.1 Seawater sampling and analysis

Seawater samples were collected at 25 stations along four shelf–offshore transects between 11 and 16^∘ S in the Peruvian OMZ during two cruises, M136 and M138, on board the R/V Meteor (Fig. 1). Cruise M136 took place in austral autumn (11 April to 3 May 2017) along two main transects at 12 and 14^∘ S (Dengler and Sommer, 2017). Two stations from M136 (stations 458 and 495) were reoccupied within a week (repeat stations 508 and 516, respectively) to evaluate the steady-state assumption in the ²³⁴Th flux calculation. The surface sample of the repeat station 508 (reoccupied 4.5 d after station 458) was missing so only results from repeat stations 495 and 516 (occupation interval 1.5 d) were compared and discussed in terms of the non-steady-state model (Sect. 3.3). ²³⁴Th sampling during cruise M138 was carried out in austral winter (1 June to 4 July 2017) and focused on four shelf–offshore transects at 11, 12, 14 and 16^∘ S.

Figure 1 Maps showing (a) locations of each station from M136 (white squares) and M138 (grey circles) and (b) monthly-averaged current field in the top 15 m from 16 April to 15 May 2017 derived from altimetry measurements (http://marine.copernicus.eu/, last access: 8 June 2020; product ID: MULTIOBS_GLO-PHY_REP_015_004). Color boxes in (a) schematically divide the four shelf–offshore transects. Map (a) was created with Ocean Data View (Schlitzer, 2014). The white box in (b) highlights our study area.

At each station, a stainless-steel rosette with Niskin bottles (Ocean Test Equipment^®) was deployed for sampling of total ²³⁴Th in unfiltered seawater and dissolved ²³⁸U (0.2 µm pore size, AcroPak^® polycarbonate membrane). High-vertical-resolution sampling was performed in the upper 200 m where most of the biological activity occurs; additional depths were sampled down to 600 m, or 50 m above the seafloor. Deep seawater at 1000, 1500, and 2000 m was sampled at three stations to determine the absolute β counting efficiency. Salinity, temperature, oxygen concentrations and fluorescence data (Table S1 in the Supplement) were derived from the sensors (Seabird Electronics^® 9plus system) mounted on the CTD frame (Krahmann, 2018; Lüdke et al., 2020).

Sample collection and subsequent chemical processing and analysis for total ²³⁴Th followed protocols by Pike et al. (2005) and the SCOR working group RiO5 cookbook (https://cmer.whoi.edu/, last access: 27 March 2019). Briefly, a ²³⁰Th yield tracer (1 dpm) was added to each sample (4 L) before Th was extracted with MnO₂ precipitates. Precipitates were filtered onto 25 mm quartz microfiber filters (Whatman^® QMA, 2.2 µm nominal pore size) and dried overnight at 50 ^∘C, after which they were counted at sea on a Risø^® low-level beta GM multicounter until uncertainty was below 3 % and again 6 months later at the home laboratory for background ²³⁴Th activities. After the second beta counting, filters were digested in an 8 M HNO₃ ∕ 10 % H₂O₂ solution (Carl Roth^®, trace metal grade). A total of 10 dpm of ²²⁹Th was added to each sample at the beginning of digestion to achieve a 1 : 1 atom ratio between ²²⁹Th and ²³⁰Th. Digested samples were diluted in a 2.5 % HNO₃ ∕ 0.01 % HF mixture, and ²²⁹Th ∕²³⁰Th ratios were measured using an inductively coupled plasma mass spectrometer (ICP-MS) (ThermoFisher^® Element XR) to determine the chemistry yield and final ²³⁴Th activities. The average yield was calculated to be 97 % ± 6 % (n=247). For a subset of samples (marked in Table S1) whose analysis failed during initial ICP-MS measurement, anion chromatography (Bio-Rad^® AG1x8, 100–200 mesh, Poly-Prep columns) was performed to remove Mn from the sample matrix before another ICP-MS analysis. This subset of samples also included three samples (marked in Table S1) whose initial ICP-MS measurement was successful to test whether anion chromatography affects final ICP-MS results. Identical ²²⁹Th ∕²³⁰Th ratios were measured for samples with and without column chromatography (see Table S1 footnotes for details).

Each ²³⁸U sample was acidified to pH∼1.6 at sea and transported home for analysis. Samples of dissolved ²³⁸U were diluted 20 times in 1 N HNO₃ at the home laboratory and spiked with an appropriate amount of ²³⁶U spike to achieve ²³⁶U:²³⁸U ∼ 1 : 1. Ratios of ²³⁶U:²³⁸U were analyzed by ICP-MS (ThermoFisher Element XR), and activities of ²³⁸U were calculated using isotope dilution. Seawater certified reference materials (CRMs), CASS-6 and NASS-7, and the International Association for the Physical Sciences of the Oceans (IAPSO) standard seawater were analyzed routinely for uranium concentrations.

2.2 Flux calculation

Assuming a one-box model, the temporal change of ²³⁴Th activities is balanced by production from ²³⁸U, radioactive decay of ²³⁴Th, removal of ²³⁴Th onto sinking particles, and transport into or out of the box by advection and diffusion (Bhat et al., 1968; Savoye et al., 2006; and references therein):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial A_{\mathrm{Th}}}{\partial t}}=\mathit{\lambda }(A_{\mathrm{U}}-A_{\mathrm{Th}})-P+V,\end{array}$$

where A_U and A_Th are respectively the activities of dissolved ²³⁸U and total ²³⁴Th, λ is the decay constant of ²³⁴Th, P is the net removal flux of ²³⁴Th, and V is the sum of advective and diffusive fluxes. It is recommended that the time interval between station occupations should be >2 weeks in order to adequately capture the temporal variability of the mean spatial gradients rather than small local changes (Resplandy et al., 2012). The solution of Eq. (1) (Savoye et al., 2006) is

$$\begin{array}{}\text{(2)}& P=\mathit{\lambda }\left[{\displaystyle \frac{A_{\mathrm{U}}\left(\mathrm{1}-e^{-\mathit{\lambda }\mathrm{\Delta }t}\right)+A_{\mathrm{Th}\mathrm{1}}\cdot e^{-\mathit{\lambda }\mathrm{\Delta }t}-A_{\mathrm{Th}\mathrm{2}}}{\mathrm{1}-e^{-\mathit{\lambda }\mathrm{\Delta }t}}}\right],\end{array}$$

where Δt is the time interval between repeat occupations of a station, and A_Th1 and A_Th2 are respectively total ²³⁴Th activities during the first and second occupation. At times when repeat sampling is not possible within an adequate cruise timeframe, steady-state conditions are generally assumed, i.e., $\frac{\partial A_{\mathrm{Th}}}{\partial t}=\mathrm{0}$. In this case, Eq. (1) is simplified into

$$\begin{array}{}\text{(3)}& P=\underset{\mathrm{0}}{\overset{z}{\int }}\mathit{\lambda }(A_{\mathrm{U}}-A_{\mathrm{Th}})\mathrm{d}z+V.\end{array}$$

The vertical flux of ²³⁴Th, P (dpm m⁻² d⁻¹), is integrated to the depth of interest. Earlier studies generally used arbitrarily fixed depths (e.g., the base of mixed layer or ML, and 100 m) for ²³⁴Th and particulate organic carbon (POC) flux estimates (e.g., Bacon et al., 1996; Buesseler et al., 1992). Recent studies emphasized the need to normalize POC flux to the depth of the euphotic zone (EZ), which separates the particle production layer in the surface from the flux attenuation layer below (Black et al., 2018; Buesseler and Boyd, 2009; Rosengard et al., 2015). In the open ocean, the depth of EZ is generally similar to ML depth. The PAR (photosynthetically active radiation) sensor was not available during both of our cruises, so it was not possible to identify the base of the EZ. For the purpose of this study, the slight difference of the exact depth chosen (ML vs. EZ) was of little relevance to the significance of physical processes and ²³⁸U variability. Due to sampling logistics, however, we did not sample at the base of the ML but 5–20 m below the ML. This depth corresponded closely to the EZ depth used in Black et al. (2018) in the same study area during austral spring 2013. For the purpose of comparison with earlier studies which reported ²³⁴Th fluxes at 100 m, we also calculated ²³⁴Th fluxes at 100 m in this study.

2.3 Quantification of the physical fluxes

The physical term V in Eq. (2) is expressed as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}V& =\underset{\mathrm{0}}{\overset{z}{\int }}\left(w{\displaystyle \frac{\partial \mathrm{Th}}{\partial z}}-u{\displaystyle \frac{\partial \mathrm{Th}}{\partial x}}-v{\displaystyle \frac{\partial \mathrm{Th}}{\partial y}}\right)\mathrm{d}z\\ & +\underset{\mathrm{0}}{\overset{z}{\int }}\left(K_x{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{Th}}{\partial x^{\mathrm{2}}}}+K_y{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{Th}}{\partial y^{\mathrm{2}}}}-K_z{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathrm{Th}}{\partial z^{\mathrm{2}}}}\right)\mathrm{d}z,\end{array}\end{array}$$

where w is the vertical (i.e., upwelling) velocity (m s⁻¹); u and v are respectively the zonal and meridional current velocities (m s⁻¹); and K_x, K_y, and K_z represent eddy diffusivities (m² s⁻¹) in zonal, meridional and vertical directions, respectively. $\frac{\partial \mathrm{Th}}{\partial z}$, $\frac{\partial \mathrm{Th}}{\partial x}$ and $\frac{\partial \mathrm{Th}}{\partial y}$ are vertical and horizontal ²³⁴Th gradients (dpm L⁻¹ m⁻¹); and $\frac{{\partial }^{\mathrm{2}}\mathrm{Th}}{\partial x^{\mathrm{2}}}$, $\frac{{\partial }^{\mathrm{2}}\mathrm{Th}}{\partial y^{\mathrm{2}}}$ and $\frac{{\partial }^{\mathrm{2}}\mathrm{Th}}{\partial z^{\mathrm{2}}}$ are respectively the second derivative of ²³⁴Th (dpm L⁻¹ m⁻²) on the zonal, meridional and vertical directions.

2.3.1 Estimation of upwelling velocities

In the Mauritanian and Peruvian coastal upwelling regions, there is strong evidence that upwelling velocities in the mixed layer derived from satellite scatterometer winds and Ekman divergence (Gill, 1982) agree well with those from helium isotope disequilibrium (Steinfeldt et al., 2015). The parameterization by Gill (1982) considers the baroclinic response of winds blowing parallel to a coastline in a two-layer ocean. Vertical velocity (w) at the interface yields

$$\begin{array}{}\text{(5)}& w={\displaystyle \frac{\mathit{\tau }}{\mathit{\rho }fa}}{e}^{-x/a},\end{array}$$

where τ is the wind stress (kg m⁻¹ s⁻²) parallel to the coast line, ρ the water density (1023 kg m⁻³), f the Coriolis parameter (s⁻¹) as a function of latitude, a the first baroclinic Rossby radius (km) and X the distance (km) to the coast.

Upwelling velocities were calculated at stations within 60 nautical miles of the coast, where upwelling is the most significant (Steinfeldt et al., 2015). We used a=15 km for all stations based on the results reported by Steinfeldt et al. (2015) for the same study area. The magnitude of monthly wind stress was estimated from the monthly wind velocities (Smith, 1988):

$$\begin{array}{}\text{(6)}& \mathit{\tau }={\mathit{\rho }}_{\mathrm{air}}{C}_{\mathrm{D}}{U}^{\mathrm{2}},\end{array}$$

where ρ_air is the air density above the sea surface (1.225 kg m⁻³), C_D the drag coefficient (10⁻³ for wind speed <6 m s⁻¹) and U the wind speed.

Monthly wind speed (m s⁻¹) fields from the MetOp-A ASCAT scatterometer sensor with a spatial resolution of 0.25^∘ (Bentamy and Croize-Fillon, 2010) were retrieved from the Centre de Recherche et d'Exploitation Satellitaire (CERSAT), at IFREMER, Plouzané (France) (data version numbers L3-MWF-GLO-20170903175636-01.0 and L3-MWF-GLO-20170903194638-01.0). We assumed a linear decrease of w from the base of the mixed layer toward both the ocean surface and 240 m depth (bottom depth of our shallowest station). Upwelling rates at any depth between 0 and 240 m at individual stations could thus be determined once w was estimated. Following Rapp et al. (2019), an error of 50 % was assigned to estimated upwelling velocities to account for uncertainties associated with the spatial structure and temporal variability of the wind field, as well as the satellite wind product near the coast.

2.3.2 Estimation of upper-ocean velocities

During both cruises a phased-array vessel-mounted acoustic Doppler current profiler (VmADCP; 75 kHz Ocean Surveyor, Teledyne RD Instruments) continuously measured zonal and meridional velocities in the upper 700 m of the water column (Lüdke et al., 2020). Postprocessing of the velocity data included water track calibration and bottom editing. After calibration, the remaining uncertainty of hourly averages of horizontal velocities is smaller than 3 cm s⁻¹ (e.g., Fischer et al., 2003). For the horizontal advective flux calculation (Eq. 3), velocities collected within a 10 km radian at inshore stations (stations 353, 428, 458, 475, 508, 904 and 907) and within a 50 km radian at offshore stations (Lüdke et al., 2020) were averaged. Data collected at the same positions within 5 d due to station repeats were also included in the velocity average. As representative for the near-surface flow, we extracted the velocity data from the top 30 m for M136 stations and top 50 m for M138 stations (defined as the “top layer” hereafter); these depths correspond to 5–20 m below the base of the ML during each cruise.

2.3.3 Estimation of vertical and horizontal eddy diffusivities

While the strength of ocean turbulence determines the magnitude of diapycnal or vertical eddy diffusivities, the intensity of meso- and submesoscale eddies determines the magnitude of lateral eddy diffusivities. During the R/V Meteor cruise (M136) and the follow-up cruise (M137) in the same region, the strength of upper-ocean turbulence was measured using shear probes mounted to a microstructure profiler. The loosely tethered profiler was optimized to sink at a rate of 0.55 m s⁻¹ and equipped with three shear sensors; a fast-response temperature sensor; an acceleration sensor; two tilt sensors; and conductivity, temperature, depth sensors sampling with a lower response time. On transit between each CTD station 3 to 9 microstructure profiles were collected. Standard processing procedures were used to determine the dissipation rate of turbulent kinetic energy (ε) in the water column (see Schafstall et al., 2010, for a detailed description). Subsequently, turbulent vertical diffusivities K_z were determined from $K_z=\mathrm{\Gamma }\mathit{\epsilon }{N}^{-\mathrm{2}}$ (Osborn, 1980), where N is stratification and Γ is the mixing efficiency for which a value of 0.2 was used following Gregg et al. (2018). Stratification (buoyancy frequency) was calculated using CTD data retrieved from microstructure profilers and following the gsw_Nsquared function from the Gibbs Sea Water library (McDougall et al., 2009; Roquet et al., 2015). A running mean of 10 dbar was applied to avoid including unstable events due to turbulent overturns. The 95 % confidence intervals for averaged K_z values were determined from Gaussian error propagation following Schafstall et al. (2010).

Altogether, 189 microstructure profiles were collected during M136 (Thomsen and Lüdke, 2018) and 258 profiles during the follow-up cruise M137 (unpublished data; 6–29 May 2017). An average turbulent vertical diffusivity profile was calculated from all inshore (<500 m water depth) profiles and another one from all offshore (>500 m water depth) profiles (Fig. S1 in the Supplement). Microstructure profiles collected during cruise M138 were not available, but there were very small variations amongst the cruise average inshore and offshore microstructure profiles from M136 and M137 despite the drastic change in the intensities of the poleward Peru–Chile Undercurrent (Lüdke et al., 2020). It thus appears appropriate to apply these average vertical diffusivities also to stations during M138.

Horizontal eddy diffusivity could not be determined from data collected during the cruises. Surface eddy diffusivities in the North Atlantic OMZ were estimated to be on the order of a few 1000 m² s⁻¹, which decrease exponentially with depth (Hahn et al., 2014). A similar magnitude of eddy diffusivities was estimated for the eastern equatorial South Pacific based on surface drifter data and satellite altimetry (Abernathey and Marshall, 2013; Zhurbas and Oh, 2004). We thus consider an eddy diffusivity of 1000 m² s⁻¹ as a good approximation in this study for the evaluation of horizontal diffusive ²³⁴Th fluxes.

2.4 Residence time of ²³⁴Th

The residence time (τ_Th) of total ²³⁴Th represents a combination of the time required for the partition of dissolved ²³⁴Th onto particulate matter and that for particle removal. In a one-box model, the residence time of an element of interest can be estimated by determining the standing stock of this element and the rates of elemental input to the ocean or the rate of elemental removal from seawater to sediments (Bewers and Yeats, 1977; Zimmerman, 1976):

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_{\mathrm{Th}}={\displaystyle \frac{A_{\mathrm{Th}\left(\mathrm{mean}\right)}\cdot Z}{P}}.\end{array}$$

For the case of ²³⁴Th, A_Th(mean) is the averaged ²³⁴Th activities of the surface layer, Z is the depth of top layer and P is the removal flux of ²³⁴Th.

3 Results

3.1 Profiles of dissolved ²³⁸U, total ²³⁴Th, oxygen and fluorescence

The vertical profiles of ²³⁸U and ²³⁴Th activities are shown in Fig. 2 and tabulated in Table S1. Data from station 508 were reported in Fig. 2 and Table S1 but excluded in the Discussion section, because the surface sample at 5 m from this station was missing, which prevents any flux calculation. Also tabulated in Table S1 are temperature, salinity, and concentrations of oxygen and fluorescence obtained from the CTD sensors. Uranium concentrations of CRMs and the IAPSO standard seawater are reported in Table S2.

Figure 2 Profiles of ²³⁸U (black) and ²³⁴Th (orange squares – M136; orange circles – M138) along with concentrations of oxygen (grey) and fluorescence (green). Profiles are organized by cruises, transects, and distance to shore from left to right and top to bottom, indicated by east (E) to west (W) arrows. Error bars for both ²³⁸U and ²³⁴Th are indicated. Dashed red lines indicate the depth of the mixed layer. The start of the oxygen-deficient zone is where oxygen diminishes. Bottom depths are indicated for stations whose bottom depths are shallower than 600 m.

Activities of ²³⁸U showed small to negligible variations with depth, averaging 2.54±0.05 dpm L⁻¹ (or 3.28±0.07 ng g⁻¹, 1 SD, n=247) at all stations. The vertical distributions of ²³⁸U did not appear to be affected by water column oxygen concentrations or the extent of surface fluorescence maxima (Fig. 2). Average U concentrations of both CASS-6 (2.77±0.04 ng g⁻¹, 1 SD, n=5) and NASS-7 (2.86±0.05 ng g⁻¹, 1 SD, n=5) measured in this study agreed well with certified values (2.86±0.42 ng g⁻¹ and 2.81±0.16 ng g⁻¹, respectively). Average ²³⁸U concentration measured in our IAPSO standard seawater (OSIL batch P156) (3.24±0.06 ng g⁻¹, 1 SD, n=27) is slightly higher than that reported in Owens et al. (2011) (3.11±0.03 ng g⁻¹, 1 SD, n=10, OSIL P149) and may reflect slight differences in U concentrations between different OSIL batches.

Total ²³⁴Th activities varied from 0.63 to 2.89 dpm L⁻¹ (Fig. 2). All stations showed large ²³⁴Th deficits in surface waters with ²³⁴Th ∕²³⁸U ratios as low as 0.25 (Fig. 3). The extent of surface ²³⁴Th deficits did not vary as a function of depths of either mixed layer or the upper oxic–anoxic interface or as a function of the magnitude of surface fluorescence concentrations (Table 1, Fig. 2). ²³⁴Th at all stations generally reached equilibrium with ²³⁸U at depths between 30 and 250 m (Table 1). The equilibrium depths were slightly shallower toward the shelf at the 11, 12 and 16^∘ S transects (Fig. 3). At station 912, deficits of ²³⁴Th extended beyond 600 m depth (Fig. 2). The following stations (stations 428, 879, 898, 906, 907, 915, 919) displayed a secondary ²³⁴Th deficit below the equilibrium depth, indicative of ²³⁴Th removal processes. A small ²³⁴Th excess at depth was only observed for station 559 at 100 m. Ratios of ²³⁴Th ∕²³⁸U for deep samples at 1000, 1500, and 2000 m varied between 0.95 and 1.02 (1.00±0.04, 1 SD, n=11), suggesting that ²³⁴Th was at equilibrium with ²³⁸U at these depths.

Table 1 ²³⁴Th fluxes due to production and decay, upwelling, and vertical diffusion below the mixed layer and at 100 m. Horizontal advective fluxes were not quantified at 100 m. Refer to text for details.

Figure 3 Shelf–offshore distributions of ²³⁴Th ∕²³⁸U along the four studied transects, as shown in Fig. 1, for M136 (a) and M138 (b). White dots denote station location.

3.2 Vertical and horizontal ²³⁴Th gradients

Discrete vertical ²³⁴Th gradients in each profile (or the curvature of the profile) were estimated by the difference in ²³⁴Th activities and that in sampling depths. As such, vertical ²³⁴Th gradients varied greatly amongst stations and were larger at shallow depths ranging from 0.003 to 0.085 dpm L⁻¹ m⁻¹ (median 0.013 dpm L⁻¹ m⁻¹). Vertical ²³⁴Th gradients were essentially negligible at and below equilibrium depths.

While calculation of the vertical ²³⁴Th gradient is straightforward, the same is hardly true for the determination of horizontal ²³⁴Th gradient. Mean ²³⁴Th activities in the top layer (see Sect. 2.3.2 for depth definition) of the water column are highly variable amongst stations (Table 3, Fig. 4) and likely reflect variations occurring at small temporal and spatial scales in the Peruvian OMZ. Quantification of the horizontal ²³⁴Th gradient between individual station may thus not be adequate to evaluate large-scale advection and eddy diffusion across the study area. Therefore, alongshore ²³⁴Th gradients on a larger spatial scale (1^∘ apart) were instead calculated by grouping stations into 1^∘ by 1^∘ grids and averaging ²³⁴Th activities of each grid for the top layer. Alongshore ²³⁴Th gradients in the top layer at nearshore stations for M138 are fairly consistent, ranging from $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{6}}$ to $\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{6}}$ dpm L⁻¹ m⁻¹, with a slightly stronger gradient in the north compared to the south. The net difference in alongshore ²³⁴Th gradient is merely $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}$ dpm L⁻¹ m⁻¹. A slightly smaller alongshore ²³⁴Th gradient of $\mathrm{4.8}\times {\mathrm{10}}^{-\mathrm{7}}$ dpm L⁻¹ m⁻¹ was observed for M136. The magnitude of the net difference in alongshore ²³⁴Th gradient for M136 cannot be adequately quantified, due to smaller spatial sampling coverage. Judging on the similarity in the spatial distributions of mean ²³⁴Th between cruises M136 and M138 (Fig. 4), it is reasonable to assume that the net difference in alongshore ²³⁴Th gradient remained similar during both cruises.

Figure 4 Distributions of averaged ²³⁴Th activities during M136 (a, top 30 m) and M138 (b, top 50 m).

3.3 Steady-state vs. non-steady-state models

The relative importance of ²³⁴Th fluxes due to advection and diffusion was assessed here assuming steady-state conditions, which assume negligible temporal ²³⁴Th variability. But how valid is this assumption in the Peruvian upwelling zone? Profiles of temperature and oxygen at repeat stations 458 and 508 showed that a lightly cooler and oxygen-depleted water mass dominated at the upper 50 m at station 508 (Fig. 5). However, an assessment of the ²³⁴Th fluxes at these two stations was not possible as the surface sample from station 508 was missing. Repeat stations 495 and 516 show substantial temporal variations in ²³⁴Th activities at each sampled depth in the top 200 m, while temperature and salinity profiles confirmed that similar water masses were sampled during both occupations (Fig. 5). Particularly, the surface ²³⁴Th deficit was more intense at station 495 (${}^{\mathrm{234}}\mathrm{Th}/{}^{\mathrm{238}}\mathrm{U}=\mathrm{0.44}$) compared to station 516 (${}^{\mathrm{234}}\mathrm{Th}/{}^{\mathrm{238}}\mathrm{U}=\mathrm{0.73}$). Correspondingly, ²³⁴Th fluxes decreased substantially from station 495 to station 516. At 100 m, the difference in ²³⁴Th fluxes between these two stations was ∼30 % (3200±90 dpm m⁻² d⁻¹ at station 495 and 2230±110 dpm m⁻² d⁻¹ at station 516). At 200 m where ²³⁴Th resumed equilibrium with ²³⁸U at both stations, the ²³⁴Th flux difference was ∼25 % (4510±220 dpm m⁻² d⁻¹ at station 495 and 3455±200 dpm m⁻² d⁻¹ at station 516). Taking the non-steady-state term in Eq. (1) into consideration (see details in Resplandy et al., 2012, and Savoye et al., 2006, for the derivation of flux formulation and error propagation) increased total ²³⁴Th at station 516 by 40 % to 3110±1870 dpm m⁻² d⁻¹ at 100 m (or 45 % to 5040±2290 dpm m⁻² d⁻¹ at 200 m), which is indistinguishable within error from fluxes at station 495. The large errors associated with the non-steady-state calculation due to the short duration between station occupations prevent a meaningful application of this model in the current study (also see discussion in Resplandy et al., 2012). As estimation of the physical fluxes is independent of the models chosen between steady and non-steady states, the following results and discussion sections regarding physical effects on the ²³⁴Th flux estimates are based on the steady-state model only.

Figure 5 Profiles of temperature (solid lines) and salinity (dashed lines) for repeated stations (a) 458 (purple) and 508 (yellow) and (d) 495 (blue) and 516 (orange). Panels (b) and (c) are respectively profiles for stations 458 and 508 of ²³⁸U (black), ²³⁴Th (color squares), and concentrations of oxygen (grey) and fluorescence (green). Panels (e) and (f) are respectively profiles for stations 495 and 516 of ²³⁸U (black), ²³⁴Th (color squares), and concentrations of oxygen (grey) and fluorescence (green).

3.4 Export fluxes of ²³⁴Th

Fluxes of ²³⁴Th due to radioactive production and decay (hereafter “production flux”), upwelling, and vertical diffusion were reported in Table 1 and Fig. 6 for both depths 5–20 m below the ML and at 100 m. The production fluxes of ²³⁴Th at 5–20 m below the ML ranged from 560 to 1880 dpm m⁻² d⁻¹, whereas at 100 m they were much higher at 850 to 3370 dpm m⁻² d⁻¹. There is no discernable trend regarding the production fluxes between the shelf and offshore stations, similar to those seen along the eastern GP16 transect (Black et al., 2018).

Figure 6 Bar charts of ²³⁴Th fluxes due to production and decay (blue), upwelling (orange) and vertical diffusion (grey) for the depths at 5–20 m below the ML (top) and 100 m below sea surface (bottom). Color boxes corresponds to individual transects in Fig. 1. Within each transect, stations from west (offshore) to east (nearshore) are listed from left to right. Error bars (1 SE) are indicated.

Alongshore winds were unusually weak off Peru preceding and during our sampling campaign as a result of the 2017 coastal El Niño (Echevin et al., 2018; Lüdke et al., 2020; Peng et al., 2019), which resulted in nominal upwelling in the water column. At nearshore stations, upwelling rates at the base of the ML varied between $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{7}}$ and $\mathrm{9.7}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹, whereas upwelling rates at offshore stations were on the order of 10⁻¹⁰ to 10⁻⁸ m s⁻¹ and essentially negligible. As a result, upwelled ²³⁴Th fluxes at 5–20 m below the ML were only significant at stations closest to shore; these stations were 428 (130 dpm m⁻² d⁻¹), 883–12 (80 dpm m⁻² d⁻¹) and 904–16 (280 dpm m⁻² d⁻¹), whose upwelled ²³⁴Th fluxes accounted for 10 %, 11 % and 25 % of the total ²³⁴Th fluxes respectively (Fig. 6). Upwelled ²³⁴Th fluxes at the rest of the stations accounted for less than 2 % of the total ²³⁴Th fluxes (6 % at stations 353 and 907–11) and were insignificant. At 100 m, both vertical ²³⁴Th gradients and upwelling rates were significantly smaller compared to shallower depths. As a result, upwelled ²³⁴Th fluxes were less than 70 dpm m⁻² d⁻¹, or less than 4 % of total ²³⁴Th fluxes.

Similarly, vertical diffusivities, shown as running mean over 20 m in Fig. S1, were an order of magnitude higher at shallow stations ($\mathrm{3.2}\times {\mathrm{10}}^{-\mathrm{4}}\pm \mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹; 1 SD, 27 to 100 m below sea surface) compared to those at deep stations ($\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{5}}\pm \mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹; 1 SD; 34–100 m below sea surface). Within the upper 27 to 33 m layer at offshore deep stations, vertical diffusivities decreased exponentially by an order of magnitude within a few meters; below this depth, vertical diffusivities remained relatively stable (Fig. S1). This is not surprising as wind-driven turbulence is most significant at the ocean surface (Buckingham et al., 2019). In this study, the sampling depths immediately below the ML were generally 30 and 60 m. A few high-vertical-diffusivity values around 30 m at deep stations were not likely representative for the 30–60 m water column layer. We thus opted to only apply vertical diffusivities below 33 m at deep stations. Relative standard errors (RSEs) associated with diffusivity estimates varied from 35 % to 55 %. Vertical diffusive ²³⁴Th fluxes at 5–20 m below the ML, determined using both vertical diffusivity and vertical ²³⁴Th gradient, varied greatly amongst stations. At shallow stations 428, 458 and 883–12, vertical diffusive ²³⁴Th fluxes made up 37 % (490 dpm m⁻² d⁻¹), 14 % (160 dpm m⁻² d⁻¹) and 21 % (160 dpm m⁻² d⁻¹) of total ²³⁴Th fluxes, respectively (Fig. 6). At the rest of the stations, vertical diffusive ²³⁴Th fluxes appeared to be insignificant, ranging between 1 % and 10 % in the total ²³⁴Th flux budget. At 100 m, vertical diffusive ²³⁴Th fluxes at station 428, 458 and 883–12 remained high at 390, 150 and 120 dpm m⁻² d⁻¹, respectively, whereas those at the rest of the stations accounted for <2 % of the total ²³⁴Th flux.

Horizontal advective and diffusive ²³⁴Th fluxes were both very small. Average alongshore current velocities (Lüdke et al., 2020) for the top layer varied from 0.06 to 0.34 m s⁻¹. At the periphery of a freshly formed anticyclonic eddy (station 915-1), alongshore current velocities could be as high as 0.53 m s⁻¹. Taking the mean alongshore velocity of 0.2 m s⁻¹ and the net difference in the alongshore ²³⁴Th gradient of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}$ dpm L⁻¹ m⁻¹, the resulting net horizontal advective ²³⁴Th flux at the top layer is ∼50 dpm m⁻² d⁻¹, a mere 3 %–9 % of the total ²³⁴Th fluxes.

Horizontal diffusive ²³⁴Th flux was estimated using an average eddy diffusivity of 1000 m² s⁻¹ (see methods Sect. 2.3.3) and the alongshore ²³⁴Th gradient. A maximum value of 10 dpm m⁻² d⁻¹ was calculated, which accounted for <1 % of total ²³⁴Th flux at all stations. Note that the horizontal advective and lateral diffusive fluxes presented here are a rough estimate and should only provide an idea of their order of magnitude. Due to the uncertainty inherent to the estimates, we refrain from adding these values to Table 1.

4 Discussion

4.1 Lack of linear ²³⁸U–salinity correlation in the Peruvian OMZ

The water column profiles of ²³⁸U in the Peruvian OMZ (Fig. 2) are similar to those seen in the open ocean (see compilations in Owens et al., 2011, and Van Der Loeff et al., 2006, and references therein). It thus appears that water column suboxic/anoxic conditions alone are not sufficient to remove U, in contrast to sedimentary U studies underlying low-oxygen waters where soluble U(VI) diffuses downward into subsurface sediments and is reduced to insoluble U(IV) (Anderson et al., 1989; Böning et al., 2004; Scholz et al., 2011). Our inference is in accord with water column ²³⁸U studies in intense OMZs in the eastern tropical North Pacific (Nameroff et al., 2002) and the Arabian Sea (Rengarajan et al., 2003), where ²³⁸U concentrations remain constant over the entire upper water column studied.

Dissolved ²³⁸U and salinity across the entire Peruvian OMZ displayed poor linear correlation regardless of seawater oxygen concentrations (Fig. 7a–b). The general consensus is that U behaves conservatively in oxic seawater in the open ocean, and early observations have shown that ²³⁸U activities can be calculated from salinity based on a simple linear correlation between the two (e.g., Chen et al., 1986; Ku et al., 1977). Compilations in Van Der Loeff et al. (2006) and Owens et al. (2011) further demonstrated that the majority of uranium data points in the global seawater dataset follow a linear correlation with seawater salinity. The ²³⁸U-salinity formulations from either Chen et al. (1986) or Owens et al. (2011) are thus generally appropriate for open-ocean conditions and have been widely used in ²³⁴Th flux studies. However, this linear ²³⁸U–salinity correlation breaks down in the Peruvian OMZ. Furthermore, the measured ²³⁸U activities in this study correlated poorly with those calculated from salinity using the Owens formulation regardless of water column oxygen concentrations (Table S2, Fig. 7c), with the former significantly higher than the projected values and differences up to 10 %. This evidence suggested that non-conservative processes have introduced significant amount of dissolved U into the water column.

Figure 7 Cross plots of measured ²³⁸U activities vs. salinity for M136 (a) and M138 (b), showing poor linear relationship between ²³⁸U and salinity. Panel (c) shows a direct comparison between measured and salinity-based ²³⁸U to further highlight the large difference between the two. The solid blue line indicates the 1 : 1 ratio between measured and projected ²³⁸U. Dashed blue lines indicate the ± errors reported in Owens et al. (2011). Error bars for measured ²³⁸U activities are smaller than symbols.

It is likely that this poor ²³⁸U–salinity correlation in the water column is not a unique feature off the coast of Peru. Poor correlations between dissolved ²³⁸U and salinity have been previously observed in open-ocean settings such as the Arabian Sea (Rengarajan et al., 2003) and the Pacific Ocean (Ku et al., 1977), as well as shelf–estuary systems such as the Amazon shelf (McKee et al., 1987; Swarzenski et al., 2004). It is possible that the narrow range of salinity within any single ocean basin precludes a meaningful ²³⁸U–salinity correlation (Ku et al., 1977; Owens et al., 2011). For the Peruvian shelf system, two possible scenarios may further explain the lack of linear ²³⁸U–salinity correlation in the water column. Firstly, authigenic U within the sediments may be remobilized under El Niño–Southern Oscillation (ENSO)-related oxygenation events. In reducing pore water, U reduction and removal from pore water is usually seen within the Fe reduction zone (Barnes and Cochran, 1990; Barnes and Cochran, 1991; Scholz et al., 2011). As such, a downward diffusive flux of U across the water–sediment interface is expected in reducing sedimentary environment. However, pore water and bottom water geochemistry measurements during two previous cruises (M77-1 and M77-2) along an 11^∘ S transect off Peru showed large diffusive fluxes of U out of the Peruvian shelf sediments despite the fact that both Fe reduction and U reduction took place in the top centimeters of sediments (Scholz et al., 2011). It was suggested that a minute increase in bottom water oxygen concentration induced by El Niño events would be sufficient in shifting the U(VI)/U(IV) boundary by a few centimeters and remobilize authigenic U (Scholz et al., 2011). Preceding and during our sampling campaign, a coastal El Niño event, with coastal precipitation as strong as the 1997–1998 El Niño event, had developed rapidly and unexpectedly in January and disappeared by May 2017 during cruise M136 (Echevin et al., 2018; Garreaud, 2018; Peng et al., 2019). This strong coastal El Niño event could induce an oxygenation event large enough to remobilize authigenic U along the Peruvian shelf. Secondly, resuspension of bottom sediments and subsequent desorption of U from ferric oxyhydroxides could affect the ²³⁸U-salinity relationship, similar to that seen on the Amazon shelf at salinity above 10 (McKee et al., 1987) and in laboratory experiments (Barnes and Cochran, 1993). Fe reduction and release from the Peruvian shelf sediments (Noffke et al., 2012; Scholz et al., 2014) could release additional U to overlying waters. The magnitude of such, however, has not been quantified.

The consequence of the notable difference between measured ²³⁸U in this study and salinity-based ²³⁸U to ²³⁴Th flux according to Eq. (2) is neither linear nor straightforward, because the vertical gradients of both ²³⁸U and ²³⁴Th strongly affects the impacts of ²³⁸U variations on ²³⁴Th fluxes. In this study, ²³⁴Th fluxes at 100 m derived from salinity-based ²³⁸U lead to significant underestimation of ²³⁴Th fluxes by an average of 20 % and as high as 40 % (Table 2). These differences in ²³⁴Th fluxes will have direct consequences for ²³⁴Th-derived elemental fluxes such as C, N, P and trace metals. It is thus important to note that U concentrations in coastal systems are highly sensitive to bottom water oxygen concentrations and redox-related U addition, the variability of which is expected to intensify with future climate change (Shepherd et al., 2017). Relatively minor variations in dissolved ²³⁸U could account for substantial overestimation/underestimation of the depth-integrated ²³⁴Th fluxes. We thus encourage future ²³⁴Th flux studies in such environments to include seawater ²³⁸U analysis.

Table 2 Comparison of ²³⁴Th fluxes at 100 m calculated with measured ²³⁸U activities and those with salinity-based ²³⁸U.

^* For the purpose of evaluating flux difference due to uranium variability, we report here ²³⁴Th fluxes only due to radioactive production and decay.

4.2 Dynamic advective and diffusive ²³⁴Th fluxes

The significance of advection and diffusion in the total ²³⁴Th flux budget highly depends on the upwelling rate, current velocity, vertical diffusivity, and ²³⁴Th gradient on the horizontal and vertical directions. Our results demonstrated that physical processes off Peru during and after the 2017 coastal El Niño have very limited impact on the downward fluxes of ²³⁴Th (Fig. 6).

Our findings are in reasonable agreement with those from the GEOTRACES GP16 eastern section along 12^∘ S from Peru to Tahiti, in which Black et al. (2018) quantified both horizontal and vertical advective ²³⁴Th fluxes. Horizontal advective fluxes for the upper 30 m water column estimated during GP16 were ∼180 dpm m⁻² d⁻¹ for all nearshore and offshore stations, similar in magnitude to those estimated in our study (∼50 dpm m⁻² d⁻¹). Upwelling fluxes along GP16 eastern section were suggested to account for 50 % to 80 % of total ²³⁴Th fluxes at the base of the euphotic zone (Black et al., 2018), a depth similar to or slightly deeper than ML depths in the current study where upwelling fluxes accounted for less than 25 % of total ²³⁴Th fluxes. Total ²³⁴Th fluxes along the GP16 eastern section, ranging from 4000 to 5000 dpm m⁻² d⁻¹ at the base of the euphotic zone, were much higher than those in our study (560 to 1900 dpm m⁻² d⁻¹ 5–20 m below the ML). This difference could be related to the period of sampling (austral autumn and winter 2017 in our study vs. austral spring 2013 for the GP16 section). We note that the estimated vertical mixing rates based on ⁷Be isotope at the base of the euphotic zone along the GP16 section (Kadko, 2017) were at least an order of magnitude higher than the upwelling rates at the base of the ML at nearby stations in our study. This difference could stem from different methods used to estimate upwelling rates at different timescales and may also reflect the dynamic upwelling system off Peru in which upwelling rates vary greatly seasonally and interannually. During cruises M136 and M138, upwelling favorable easterly winds off Peru were weak, resulting in negligible coastal upwelling. Coastal upwelling in the same general area was also suggested to be negligible in austral summer 2013 during cruise M92 due to nominal surface wind stress (Thomsen et al., 2016). Results from studies conducted in the same year (October to December 2013, Kadko, 2017; December 2012, Steinfeldt et al., 2015; January 2013, Thomsen et al., 2016) indicate that seasonal upwelling rates vary drastically in the Peruvian upwelling zone. The seasonal dynamics of coastal upwelling off Peru are similar to those seen in the Arabian Sea, where large upwelled ²³⁴Th fluxes only occurred during the middle-to-late southwest monsoon at stations close to shore (Buesseler et al., 1998). Our findings lend further support to earlier studies that advection and diffusion are seasonally important for ²³⁴Th fluxes in regions with high upwelling velocities and diffusivities such as the equatorial Pacific (Bacon et al., 1996; Buesseler et al., 1995; Dunne and Murray, 1999) and coastal sites such as the Arabian Sea (Buesseler et al., 1998) and offshore Peru (Black et al., 2018; this study).

4.3 Residence time of ²³⁴Th in the Peruvian OMZ

The residence time calculated using Eq. (6) was based on a simplified one-dimension (1D) model of Zimmerman (1976). This 1D steady-state model is obviously an oversimplification of a multi-dimensional process; it however provides a good first-order estimate for understanding the highly dynamic nature of the ²³⁴Th residence time. It also provides a reasonable value that can be directly compared to values estimated in earlier ²³⁴Th flux studies that did not consider the physical processes. Furthermore, we showed in the Discussion (Sect. 4.2) that physical processes, namely upwelling and vertical diffusion, are only important at a few shelf stations. We thus consider this simple 1D model robust in estimating the residence time of total ²³⁴Th.

In this study, residence time of total ²³⁴Th in the top layer varied from 20 d at shallow stations to 95 d at deep stations (mean $\mathit{\tau }=\mathrm{51}\pm \mathrm{23}$ d, 1 SD, n=24; Table 3). These values were similar to those estimated within the California Current (Coale and Bruland, 1985) and the residence times of particulate organic carbon (POC) and nitrogen (PON) (Murray et al., 1989) but were much longer than predicted in nearshore shelf waters where residence times of total ²³⁴Th were on the order of a few days (Kaufman et al., 1981; Kim et al., 1999; and references therein). The longer residence times estimated in our study could reflect a combination of weak surface ²³⁴Th deficits (²³⁴Th=0.63 to 1.82 dpm L⁻¹) (Fig. 3) and low export fluxes (800 to 2000 dpm m⁻² d⁻¹, Fig. 7). Nearshore seawater samples during GP16 (Black et al., 2018) featured similar surface ²³⁴Th deficits (²³⁴Th = 0.63 to 1.33 dpm L⁻¹) but much higher downward ²³⁴Th fluxes (4000 to 5000 dpm m⁻² d⁻¹) as a result of strong upwelling, implying that residence times of total ²³⁴Th in the Peruvian OMZ during GP16 occupation would be 3–6 times shorter. Indeed, a quick reassessment of the GP16 data predicted a shorter residence time of total ²³⁴Th of 5–23 d within the euphotic zone of the coastal Peruvian OMZ.

Table 3 Residence time of total ²³⁴Th in the top layers of Peruvian OMZ.

^* Here “the top layer” refers to the top 30 m during M136 and top 50 m during M138.

These temporal variations on the residence times of total ²³⁴Th have important implications for the estimation of POC fluxes and quantification of carbon export efficiency. Firstly, seasonal changes in Th residence times reflect variations in particle removal over different integrated timescales. For example, POC produced in surface waters during GP16 (austral spring 2013) (Black et al., 2018) would have been exported out of the euphotic zone 3–6 times faster than during austral autumn 2017 (this study). Secondly, to properly evaluate carbon export efficiency, surface net primary production (NPP) should be averaged over a timescale similar to the residence time of total ²³⁴Th during station occupation. Applying a 16 d averaged NPP for the export efficiency estimate (Black et al., 2018; Henson et al., 2011) would likely not be appropriate in the current study in which total ²³⁴Th fluxes integrated timescales of several weeks. ²³⁴Th residence times should thus be properly quantified in coastal studies before deriving export efficiencies over varying NPP integration timescales.

5 Conclusions and implications for coastal ²³⁴Th flux studies

Advection and diffusion are important in coastal and upwelling regions with respect to ²³⁴Th export fluxes (Bacon et al., 1996; Buesseler et al., 1995, 1998; Dunne and Murray, 1999). Our findings show that their significance is subject to the seasonal variability of the current and upwelling velocities, diffusivities, and ²³⁴Th gradients and should be evaluated on a case-by-case basis. Advective fluxes are perhaps the most straightforward to estimate as current velocities can be obtained routinely from shipboard ADCP measurements and upwelling rates calculated from satellite wind stress (Steinfeldt et al., 2015; Bacon et al., 1996). Horizontal and vertical velocities derived from general ocean circulation models also provide a good first-order estimate for advective ²³⁴Th fluxes; this approach has been successfully demonstrated in a few studies (Buesseler et al., 1995, 1998). In addition, the anthropogenic SF₆ tracer and radium isotopes, widely used to quantify nutrient and Fe fluxes (Charette et al., 2007; Law et al., 2001), as well as ⁷Be isotope (Kadko, 2017), could be used independently to constrain horizontal and vertical exchange rates of ²³⁴Th (Morris et al., 2007; Charette et al., 2007; Buesseler et al., 2005). When in situ microstructure measurements are available (this study), vertical diffusivity can be directly calculated to estimate the vertical diffusive ²³⁴Th fluxes. Yet, microstructure analysis is not a routine measurement on oceanographic cruises. Earlier studies in the equatorial Pacific and the Gulf of Maine have shown that general ocean circulation models and a simple assumption on dissipation coefficients could provide a robust estimate on vertical and horizontal diffusivities (Benitez-Nelson et al., 2000; Gustafsson et al., 1998; Charette et al., 2001). Therefore, the calculation of physical fluxes is possible, though challenging, and ²³⁴Th fluxes due to physical processes should be carefully considered when conducting research in a coastal and upwelling systems.

A striking finding in this study is that the assumption of a linear ²³⁸U–salinity correlation could lead to one of the largest errors in ²³⁴Th flux estimates. In our study, using the salinity-based ²³⁸U activities resulted in significant underestimation of total ²³⁴Th fluxes by as much as 40 %. Because the translation of ²³⁸U activities to ²³⁴Th fluxes is not linear, larger differences between measured and salinity-based ²³⁸U do not necessarily contribute to greater overestimation or underestimation of ²³⁴Th fluxes. For example, a moderate difference of 3 %–6 % in ²³⁸U throughout the upper 100 m at station 898 leads to a 40 % difference in the final ²³⁴Th flux, while a 5 %–9 % difference in ²³⁸U at station 906 only resulted in a 16 % ²³⁴Th flux difference (Tables 2, S2). We would thus stress the importance of ²³⁸U measurements in future ²³⁴Th flux studies, particularly in coastal and shelf regions.

Finally, our study showed that the residence times of total ²³⁴Th in the Peruvian nearshore waters varied seasonally. Tropical OMZs are important hotspots for carbon sequestration from the atmosphere and enhanced sedimentary carbon preservation (Arthur et al., 1998; Suess et al., 1987). These OMZs are projected to intensify as a result of future climate change (Keeling and Garcia, 2002; Schmidtko et al., 2017; Stramma et al., 2008). Future studies should take into consideration the large temporal variations of the residence times of total ²³⁴Th in order to properly evaluate how carbon biogeochemical cycles and carbon export efficiency in these OMZs will respond to continuing ocean deoxygenation.