Observations of diapycnal upwelling within a sloping submarine canyon

BLT Recipes experiment

The data presented here were collected during expedition DY132 of RRS Discovery⁵¹ between 19 June 2021 and 29 July 2021, supported by the Natural Environment Research Council (NERC) and National Science Foundation (NSF)-funded Boundary Layer Turbulence and Abyssal Recipes project. The experiment consisted of a dye release together with surveys and moorings measuring hydrographic data, velocity, and shear and temperature microstructure.

Dye release

A fluorescent dye (fluorescein C₂₀H₁₀O₅Na₂) was used for this experiment. The mixture used for the release consisted of 64 l of isopropyl, 6 l of seawater and 149 l of 40% fluorescein liquid. The mixture was chosen to be denser than the local water at the desired depth to ensure that it was released properly and did not rise without mixing. The dye-release system (InkBot; University of Exeter) comprises a 219-l drum equipped with a Seabird SBE 911plus CTD, an AQUAtracka III fluorometer (Chelsea Technologies Group), an altimeter and an OCEANO 2500S Universal acoustic release. For deployment, the InkBot was attached to the CTD sea cable and lowered to 10 m above the bottom using the ship winch. Once in position, the drum was emptied by simultaneously flipping the drum and opening the drum lid using the acoustic release. The drum was upside down for 15 min, while readings from the onboard CTD enabled the winch operator to ensure that the InkBot remained at the release temperature. The dye release began at 03:01 am on 1 July 2021. The BLT dye release is novel in its depth and deep-ocean environment, where common dye sampling techniques such as aerial multi-spectral imagery could not be used. Instead, we performed rapid-repeat fluorometer measurements using the FCTD and operated the ship in a pattern that dynamically followed the dye.

FCTD

The FCTD system⁵² is a tethered profiler with a Seabird SBE49 CTD, a dual needle micro-conductivity probe (built in-house at Scripps Institution of Oceanography by the Multiscale Ocean Dynamics group) and, for this experiment, a Turner C-FLUOR fluorometer and an altimeter. The instrument was raised and lowered with a direct drive electric winch at vertical speeds of approximately 3 m s⁻¹ while the ship was steaming at a speed of 0.5–1 knots. Data from both up- and down-casts were used. The goal of the FCTD survey was to sample the dye release as many times as possible before dye concentrations became too low to be detected. Transects were focused along the canyon axis to transit from one end of the dye patch to the other before turning around and repeating. In addition, two cross-canyon transects and a 12-hour time-series station were completed to better understand the full extent of the dye patch. During a transect, the FCTD profiled over a roughly 100-m vertical range and to within less than 10 m of the seafloor. Data from the FCTD were monitored and used interactively to guide the survey.

The micro-conductivity probe was used to estimate the dissipation rate of temperature variance (χ). Conductivity gradients can be used as temperature dominates the conductivity variance⁵³.

MP1 and MP2 moorings

The MP1 mooring was deployed on the canyon axis at 2,028-m water depth at 54° 14.312′ N, 11° 56.923′ W. The mooring was 600-m tall and consisted of a downwards-looking RDI Longranger 75-kHz acoustic Doppler current profiler (ADCP) mounted in the top float and a McLane moored profiler outfitted with a Seabird SBE52 CTD, a Falmouth Scientific acoustic current meter and an epsilometer turbulence package (Multiscale Ocean Dynamics group, Scripps Institution of Oceanography)⁵⁴ to measure χ from temperature gradients, and a Turner C-FLUOR fluorometer. The profiler was deployed on 28 June 2021, and profiled continuously for 7 days, collecting one profile every 30 min. The mooring was redeployed (as MP2) on 7 July 2021, without the fluorometer and epsilometer at 54° 10.938′ N, 11° 50.572′ W at a depth of 1,676 m. This deployment lasted until 6 October 2021.

MAVS1, TCHAIN and MAVS2 moorings

The MAVS1 and MAVS2 moorings were deployed for approximately 3 months from 6 July 2021 to 7 October 2021. The MAVS moorings were each 300-m tall and consisted of 8 modular acoustic velocity sensors (MAVS; Woods Hole Oceanographic Institute), 80 RBR Solo or Seabird SBE56 thermistors, a Seabird SBE37 CTD and an RDI Longranger 75-kHz ADCP mounted on the top float. MAVS1 was deployed at 54° 11.849′ N, 11° 51.719′ W at 1,612-m depth and MAVS2 was deployed at 54° 10.938′ N, 11° 50.572′ W at 1,466-m depth. The TCHAIN mooring (Hans van Haren, Royal Netherlands Institute for Sea Research), deployed at 54° 11.413′ N, 11° 51.137′ W at 1,529 m, consisted of a 150-m thermistor chain with 102 pre-attached thermistors and a 75-kHz RDI ADCP mounted on the top float. TCHAIN was deployed from 6 July 2021 to 11 August 2022. Detailed analysis of the long-term moorings will be discussed elsewhere.

Fluorometer calibration

The Turner C-FLUOR fluorometer calibration uses the linear relationship C = (V − a)b where C is the concentration in ppb, V is the measured voltage, a is an offset and b is the calibration coefficient. The offset voltage was lower than the factory-provided value for both fluorometers used. For MP1, we used profiles before the dye release, and for the FCTD we used the upper 500 m of down-casts from the surface to measure the mean background voltage when zero dye was present. We used the factory calibration coefficient for both fluorometers. This gave calibrations of C_FCTD = (V_FCTD − 0.0140) × 31.1831 and C_MP = (V_MP − 0.0139) × 31.3309. The minimum detectable level of the fluorometers was determined as three standard deviations above the mean background level. This gave minimum detectable concentrations of 0.06 ppb for the fluorometer on the FCTD and 0.006 ppb for the fluorometer on MP1. The factor of 10 difference between detection limits is probably owing to individual sensor differences and an electrically noisier channel on the FCTD. Observed dye concentrations at MP1 were up to three times the detection limit. At the end of the FCTD survey, the levels were down to twice the detection limit.

Sampling uncertainty

Owing to high flow speeds—up to 0.4 m s⁻¹ along the canyon—the dye was advected rapidly along the canyon, requiring the dye survey to focus on two-dimensional transects through the patch. A source of uncertainty in our results comes from the resulting under-sampling of the dye patch. The two cross-canyon transects completed during the survey (Extended Data Fig. 5h,i) show that lateral distribution of the dye varied in time or along the length of the patch. During the first cross-canyon transect, the dye was spread fairly uniformly across the canyon width; however, in the second transect, the dye was banked on the northeast side of the canyon. Thus, along-canyon transects, which focused on the canyon axis, may at times have sampled through only the edge of the patch.

With only the velocity measurements from MP1 (down-canyon of the dye for most of the experiment), it is difficult to estimate where the dye patch might have been in the cross-canyon direction. Along-canyon transects may have sampled only the front or back edge or a small section of the patch at times, particularly during later transects when the dye patch was large. On the basis of dye measurements from MP1 (Fig. 2), the entire patch was about 4-km long 20 hours after the release. Later, along-canyon transects were all shorter than this (Extended Data Fig. 8). It is possible that MP1, located in the centre of the canyon, may have missed a large concentrated patch of dye if it were banked against one wall.

To investigate the effect of under-sampling on estimates of the centre of mass, we subsampled the cross-canyon transects (Extended Data Fig. 5h,i), the longest along-canyon transects (Extended Data Fig. 5b,c) and the time series (Extended Data Fig. 5l). Subsampling involved systematically selecting ten consecutive profiles along a transect. This approach ensured that the subsamples represented a fraction of the dye patch. The centre of mass was then estimated for each subsample. The results are shown in Extended Data Fig. 9. The effect of subsampling varies depending on the transect. The size of these standard deviations across different samples (Extended Data Fig. 9) varies from transect to transect in a similar way to the standard deviations of the dye-weighted density of each transect (bars in Extended Data Fig. 7). This is likely because the length and temperature range of the transects impact the overall standard deviation as well. As we cannot get an accurate estimate of the error due to sampling from all of the transects, we will use the standard deviation of the dye-weighted average to estimate the error.

This method may be particularly insufficient for the shorter transects. Looking at transects e, f and g (Fig. 3e–g), the temperatures measured are all warmer than 3.85 °C (green temperature contour) and, particularly in transects f and g, dye concentrations were relatively weak. However, during transect h, which was in the cross-canyon direction, there is more concentrated dye colder than 3.85 °C. This indicates a concentrated, colder part of the dye patch not measured in the three previous transects, and the true centres of mass may be colder than was measured. There is no way for us to account for this error accurately.

Another artefact of our sampling was due to the phase of the tide. For example, during the second transect (Extended Data Fig. 5b), we sampled the dye in a down-canyon direction as the dye itself was moving down-canyon, causing the patch to appear spread over a larger extent than it was. This may be another source of error in our estimates.

Typically, long-term chemical tracer release surveys are designed to provide information on the three-dimensional spread of the tracer. Observations are often objectively mapped to higher-resolution grids⁵⁵. These maps provide an estimate of the fraction of tracer found. Doing such an inventory is difficult, given the two-dimensional nature of our sampling pattern. However, as this study focuses on the first rather than the second moment of the tracer, the results are less impacted by outliers.

Calculations

Estimating adiabatic versus diapycnal upwelling

Consider the situation in our canyon, which is equivalent to that in many past buoyancy and mass budget calculations^3,56: a volume bounded below by a sloping seafloor, on the sides by the canyon walls and the top by a neutral surface γ_o, which is a mean distance H above the bottom and has an area A_z out to the location of the in-flow. A lateral flow u_in(t) is incident beneath γ_o through a cross-sectional canyon area A_x. Volume conservation for an incompressible fluid requires one of two things to happen: γ_o can rise adiabatically at a rate w_adia, which requires no mixing, or fluid can exit the volume via a divergent turbulent buoyancy flux, J_b, which produces a turbulent diapycnal velocity \({w}^{* }=\frac{1}{{N}^{2}}\frac{\partial {J}_{b}}{\partial z}\). Hence, at all times, the following balance must hold:

Assuming u_in, w_adia and w* are constant in space, and the canyon has a rectangular cross-section with constant width, equation (2) can be simplified to

where the ratio A_x/A_z is approximately equal to tan α, where α is the slope of the bathymetry. On very long timescales, the observed constancy of the ocean’s stratification requires w_adia = 0, giving a balance between u_in and w*.

We verify that this balance holds at both MP1 and MP2 on the approximately 3-day timescales of the dye observations. The in-flow velocity, u_in, was estimated as the along-canyon velocity below an isotherm. The 3.7 °C isotherm and 4.2 °C isotherms were chosen for MP1 and MP2, respectively, as they were on average about 100 m above the bottom and therefore within the region where near-boundary mixing occurs. Given the displacement of these isotherms (η), the adiabatic velocity is then w_adia = dη/dt. Velocities and temperatures were low-pass filtered at a period of 48 hours with a fourth-order Butterworth filter to remove the diurnal and semi-diurnal tides.

Extended Data Fig. 3 shows that both u_in and w_adia vary in time but that the dominant balance is between u_in and w*. The average of w_adia at MP2 is approximately zero, and u_in is within one standard deviation of the mean of w_adia just 4% of the time. These results bolster our interpretation that in-flow is balanced by diapycnal transport on subtidal timescales as in longer-term tracer releases such as the Brazil Basin. By contrast, adiabatic motions on tidal timescales are, in fact, key for effecting the observed water-mass transformation and exchange with the interior, as argued in the main text.

Centre of mass

For a dye or tracer of concentration C, the tracer-weighted average operator is defined as

Here we use the first moment of the tracer-weighted average of the potential density anomaly (\(\overline{{\sigma }_{\theta }}\)), which represents the centre of mass of the dye in density space, to describe the location of the dye patch^40,41. Ideally, the integral is performed on the full extent of the three-dimensional dye patch. In practice, however, we are limited to the spatial information of the survey. The integrals were estimated as sums in the vertical and along-transect directions for each transect through the dye patch.

Upwelling rate

Given the first density moment, the dye-weighted diapycnal velocity is given by \({w}_{{\rm{dye}}}=-\frac{{\partial }_{t}\overline{{\sigma }_{\theta }}}{\left|\overline{\nabla {\sigma }_{\theta }}\right|}\) (refs. ^40,41). As mixing acts on dye gradients as well as density gradients, the dye-weighted diapycnal velocity yields twice the dye-weighted density velocity, \({\partial }_{t}\overline{{\sigma }_{\theta }}=2\overline{\dot{{\sigma }_{\theta }}}\), where \(\,{\rm{D}}{\sigma }_{\theta }/{\rm{D}}t=\dot{{\sigma }_{\theta }}\) is the material derivative of the potential density anomaly⁴⁰. A total of 12 transects were completed over 3 days before dye concentrations became too low to detect. Below, we describe four methods for estimating the diapycnal upwelling rate, the results of which are shown in Table 1.

The slope of a weighted linear regression is used to estimate the evolution of the dye’s centre of mass over time (\({\partial }_{t}\overline{{\sigma }_{\theta }}\); Extended Data Fig. 7). Weights chosen for the regression were \({\mathcal{W}}=\frac{n{s}_{i}^{-2}}{{\sum }_{i=1}^{n}{s}_{i}^{-2}}\) where s_i are standard deviations of dye-weighted density for each transect i = 1, …, n and n = 12, such that transects with large standard deviations are weighted low. Linear regression (solid blue line in Extended Data Fig. 7), with R² = 0.852, yields \({\partial }_{t}\overline{{\sigma }_{\theta }}=-0.0341\pm 0.0038\,{\rm{kg}}\,{{\rm{m}}}^{-3}\,{{\rm{d}}}^{-1}\) where the standard error of the slope of the fit gives the error. We approximate the density gradient with its vertical component (\(\overline{{\partial }_{z}{\sigma }_{\theta }}\)). Here the horizontal gradient of the density is approximately ten times smaller than the vertical component and, therefore, has a negligible impact on the vertical velocity. We calculate the dye-averaged vertical density gradient for each transect and then use the weighted average over all the transects. The weights are equivalent to \({\mathcal{W}}\) but are a function of the standard deviations of the dye-weighted density gradient. The density gradient used for the calculation of w_dye is then \(\overline{{\partial }_{z}{\sigma }_{\theta }}=-1.4\times 1{0}^{-4}\pm 0.3\times 1{0}^{-4}\,{\rm{kg}}\,{{\rm{m}}}^{-4}\). For comparison, using all the FCTD data, without weighting by the concentrations, the average density gradient is −2.0 × 10⁻⁴ ± 6.4 × 10⁻⁴ kg m⁻⁴. The weighted linear regression gives an upwelling velocity of 250 ± 75 m d⁻¹, assuming that errors associated with the time rate of change of the centre of mass and the centre of mass of the vertical gradient are independent. Temporal and spatial density gradients in this location are not independent owing to the strong influence of the tide on the stratification; however, quantifying this effect on the centre of mass is difficult.

The linear regression can also be done without weighting (Extended Data Fig. 7, solid red line). In this case, \({\partial }_{t}\overline{{\sigma }_{\theta }}=-0.0207\pm 0.0055\,{\rm{kg}}\,{{\rm{m}}}^{-3}\,{{\rm{d}}}^{-1}\) and \(\overline{{\partial }_{z}{\sigma }_{\theta }}=-2.0\times 1{0}^{-4}\pm 4.7\times 1{0}^{-5}\,{\rm{kg}}\,{{\rm{m}}}^{-4}\), yielding an upwelling velocity of 101 ± 50 m d⁻¹. For this fit, R² = 0.587.

An alternative method for estimating the upwelling rate is to calculate pairwise estimates between transects²². Given 12 individual estimates of the centre of mass, there are 66 estimates of the upwelling rate between pairs of transects. The time evolution of the dye-weighted density (\({\partial }_{t}\overline{{\sigma }_{\theta }}\)) is estimated by finite differencing between each observation, where the observation time is chosen to be the average time for the transect. That is, \({({\partial }_{t}\overline{{\sigma }_{\theta }})}_{i,j}=\frac{{(\overline{{\sigma }_{\theta }})}_{j}-{(\overline{{\sigma }_{\theta }})}_{i}}{{t}_{j}-{t}_{i}}\). The average over all pairwise estimates is \({\partial }_{t}\overline{{\sigma }_{\theta }}=-0.0164\pm 0.0052\,{\rm{kg}}\,{{\rm{m}}}^{-3}\,{{\rm{d}}}^{-1}\) where the error is the standard error over all pairwise estimates. The spatial density gradient (\(| \nabla \overline{{\sigma }_{\theta }}| \)) is taken to be the average of the pair of tracer-weighted vertical density gradient estimates. The overall average for the vertical density gradient is the same as for the unweighted linear regression case. Taking the average of all 66 estimates of the upwelling rate gives 125 ± 31 m d⁻¹ where the error is the standard error on the mean.

A final, fourth estimate of the upwelling rate can be calculated from the change in depth of the centre of mass during the experiment; we consider depth changes relative to the bottom to avoid spurious upwelling due to tidal aliasing. Using the dye-weighted average of the height above bottom, the height of the centre of mass during the first transect was 28 m above the bottom. Similarly, for the final observation during the time series, the average height of the centre of mass was 151 m above the bottom. Using the same time convention as the pairwise method, the time difference between observations was 46 hours, giving an upwelling rate of 64 m d⁻¹. The change in density between the first and last observations was −0.0312 kg m⁻³ d⁻¹.

All four of these methods are imperfect: the weighted linear regression uses weights that may not be representative of the actual sampling error; neither the unweighted linear regression nor the pairwise methods account for the sampling error at all; and the change in depth estimate does not account for the diapycnal component. However, all these methods give a significant positive velocity of \({\mathcal{O}}(100\,{\rm{m}}\,{{\rm{d}}}^{-1})\) and together provide confidence in our assertion that we observed diapycnal upwelling.