# The ocean flows downhill near the seafloor and recirculates upward above

René Schubert, Jonathan Gula, Esther Capó, Pierre Damien, M. Jeroen

Molemaker, Clément Vic, James C Mcwilliams

► To cite this version:

René Schubert, Jonathan Gula, Esther Capó, Pierre Damien, M. Jeroen Molemaker, et al. The ocean flows downhill near the seafloor and recirculates upward above. Nature Communications, 2025, 16 (1), pp.5873. 10.1038/s41467-025-61027-2. hal-05196317

HAL Id: hal-05196317

https://hal.univ-brest.fr/hal-05196317v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The ocean flows downhill near the seafloor and recirculates upward above

Received: 15 November 2024

Accepted: 10 June 2025

Published online: 01 July 2025

<figure><img src="image_3.png" /></figure>

Check for updates

René Schubert Ⓐ, Jonathan Gula Ⓑ<fnref n="1" />,²,³, Esther Capó², Pierre Damien², M. Jeroen Molemaker², Clément Vic<fnref n="1" /> & James C. McWilliams²

The ocean's circulation redistributes heat, salt, biota, dissolved gases, microplastics, and sediments on Earth. The abyssal ocean, in the lowest 1000 m above the seafloor, moves on average with the deeper seafloor to its left in the Northern Hemisphere and to its right in the Southern Hemisphere. This finding has received little attention and its consequences for the abyssal vertical circulation have remained largely unexplored. Here, we show, using current-meter measurements and numerical simulations, that the interior flow, $O(100 \text{ m}) - O(1000 \text{ m})$ above the seafloor, is deflected within the bottom boundary layer, the lowest $O(10 \text{ m})$, into a widespread downhill flow. This flow intensifies with the steepness of the seafloor. We further reveal that typical local changes in seafloor steepness lead to a shallow divergence and a deep convergence of this downhill flow. These are connected by an overlying upward recirculation forming closed overturning cells that extend on average over the lowest 1000 m of the ocean. Our study improves the understanding of the oceanic abyssal circulation and the climate-relevant overturning. Future research should focus on quantifying the transports of heat, particles, and dissolved chemicals associated with these abyssal slope overturning cells.

The vertical circulation of the ocean is of great importance for the Earth system, as it is associated with a huge transport of heat and dissolved gases that are exchanged with the atmosphere at the sea surface<fnref n="1" />,². Since the beginning of the industrial era, anthropogenic gas emissions have caused a decrease in the outgoing long-wave radiation, resulting in an excess heat in the Earth system. Most of this excess is absorbed by the ocean and transported into the interior³,⁴. One contributor to this transport is the global overturning circulation, which is initialized by surface cooling at high latitudes, leading to increasing densities that induce a net downwelling of water. In the Labrador, Greenland, Nordic, Ross and Weddell Seas, water sinks into the abyss and flows equatorward into the deep basins of the global ocean⁵. Tides and winds provide the bulk of the energy that fuels the diapycnal mixing required to upwell the abyssal waters back to the surface⁶⁻⁸. Measured mixing rates are orders of magnitude stronger above steep and rough topography than above flat topography⁹⁻¹¹. Further, the mixing rates above rough topography are observed to increase with

depth in the lowest few hundreds of meters above the bottom, particularly where strong flows reach the bottom. In the absence of compensating horizontal changes in mixing rates, this implies diapycnal downwelling¹²⁻¹⁴. Furthermore, the abyssal mixing reduces the stratification, which is maintained by restratifying processes, including submesoscale baroclinic mixed-layer instabilities in the bottom boundary layer¹⁵,¹⁶. The fact that mixing vanishes at the solid seafloor, implies that the mixing rate decreases with depth, resulting in diapycnal upwelling right above the bottom¹². This near-bottom diapycnal upwelling has been simulated¹⁷ and observed with a dye-release experiment¹⁸ in thin canyons.

In this study, we investigate the Eulerian ocean circulation (in horizontal and vertical space), which is nonlinearly related to the diapycnal circulation (across density surfaces). We present a systematic analysis of the orientation of the time-mean oceanic flow with respect to the underlying sloping seafloor. We thereby focus on the abyss, the lowest 1000 m of the ocean, where oceanic flows have been shown to

<footnote n="1">¹Univ Brest, CNRS, IRD, Ifremer, Laboratoire d'Océanographie Physique et Spatiale (LOPS), IUEM, Plouzané, France. ²UCLA, Los Angeles, CA, USA. ³Institut Universitaire de France (IUF), Paris, France. ✉e-mail: rene.schubert@univ-brest.fr</footnote>
be important for the transports of biota, chemicals, sediments, and microplastics¹⁹⁻²². In the presence of sloping seafloors, idealized mixing-driven bottom boundary layer models show that abyssal mixing induces a time-mean flow with near-bottom uphill and interior downward across-isobath components²³⁻²⁶. Across-isobath is the direction perpendicular to the isobaths, the isolines of the bathymetry. In thin fracture zones and rift valleys of mid-ocean ridges, near-bottom current-meter measurements also showed horizontal time-mean currents with across-isobath components directed toward the shallower seafloor, implying uphill flows¹⁸,²⁷⁻²⁹. However, above the continental slope and above seamounts, near-bottom current-meter measurements showed across-isobath flows towards the deeper seafloor, implying downhill near-bottom flows¹⁹,³⁰,³¹. Further, high-resolution regional simulations of the western Mediterranean Sea also show mainly downhill near-bottom flows³². Here, we systematically investigate a global collection of current-meter measurements and model simulations that cover large parts of the global ocean to identify the general orientation of the abyssal circulation.

A scalar parameter that measures the alignment of an arbitrary flow with the isobaths below is the topostrophy³³⁻³⁵

$$ \tau = \operatorname{sign}(\theta) \left( \bar{u} \frac{\partial H}{\partial y} - \bar{v} \frac{\partial H}{\partial x} \right), \quad (1) $$

where sign(θ) is the sign of the latitude θ, overlines indicate time averages, u is the zonal (west-east) velocity component in direction x (positive eastward), v is the meridional (south-north) velocity component in direction y (positive northward), and H is the seafloor depth. Positive-topostrophic flows have the deeper seafloor to the left (right) and negative-topostrophic flows have the deeper seafloor to the right (left) on the Northern (Southern) hemisphere (Fig. 1a). As the amplitude of τ increases, the flow more closely follows the isobaths. A flow with τ=0 flows in the across-isobath direction. The along-isobath component of the time-mean flow in idealized mixing-driven bottom boundary layer models has the deeper seafloor to the right (left) on the Northern (Southern) hemisphere²⁴,²⁵. In contrast to the negative-topostrophic prediction of the mixing-driven models, the present and previous studies based on current-meter measurements, floats, and ocean general circulation model simulations show predominantly positive topostrophy in the oceanic interior³³,³⁴,³⁶⁻³⁸. This is consistent with the barotropic theory of mean-flow driven by eddies, time-varying geostrophic flows, above spatially variable seafloor topography³⁹⁻⁴². Ocean general circulation model simulations with higher resolution, improved parameterization of unresolved effects, and energy conserving momentum advection schemes are associated with increased topostrophy. This means that the flow is more widespread positive topostrophic, and/or is more (less) aligned with the bathymetry in positive- (negative-)topostrophic direction, and/or has a higher (lower) speed in positive- (negative-)topostrophic direction⁴³⁻⁴⁵. Furthermore, current-meter measurements and simulations show larger τ with increasing depth⁴⁶,⁴⁷. If eddies mainly drive a positive topostrophic flow and mixing a negative topostrophic flow, the predominantly positive-topostrophic abyssal current-meter measurements indicate that the effects of eddies overcomes the effect of mixing on the abyssal time-mean flow in most parts of the global ocean. Imposing a positive-topostrophic along-isobath flow in a one-dimensional boundary layer model further leads to a downhill instead of an uphill near-bottom flow³²,⁴⁷.

In this study, we first confirm on the basis of current-meter measurements and numerical ocean model simulations of the Atlantic, the Pacific, and the Western Mediterranean Sea (Med. Sea) that the interior flow is widespread positive-topostrophic (Section “Positive-topostrophic time-mean interior flow”). The simulations are integrated with kilometer-scale grid-spacing in the horizontal and vertical grids that follow the seafloor and that are refined near the seafloor. These vertical

grids allow for a better representation of the seafloor topography and associated effects on the oceanic flows compared to other vertical grids used in ocean models. Second, we show that the ocean flows on average downhill near the seafloor and recirculates upward above (section “Consequences of interior positive topostrophy for the near-bottom flow”). Third, we show that the along-isobath interior flow and the downhill flow right above the bottom are stronger above larger topographic slopes. Typically, the seafloor becomes steeper from the maximum elevation of the seafloor to about mid-depth and less steep below. Both facts combined imply that there is near-bottom divergence, and thus downwelling into the near-bottom layer, in the former regions and near-bottom convergence, and thus upwelling into the interior, in the latter regions. The so established abyssal slope overturning cells are closed by an upward recirculation of the water, which can reach far up into the water column (section “Overturning above sloping seafloor”).

## Results

### Positive-topostrophic time-mean interior flow

The interior topostrophy is positive at most abyssal current-meter measurement locations and in most of the simulated domains (Fig. 1a). This means that the abyssal flow mainly has the deeper seafloor to the left (right) on the Northern (Southern) hemisphere. Topostrophy is larger above steeper seafloor slopes and reaches its largest values at seafloor hills, trenches, and continental slopes (see section “Numerical Simulations” for a map of the seafloor slope). Topostrophy is reduced on the shelf and between 5°S and 5°N, where the planetary vorticity goes to zero. Regions between 5°S and 5°N and where the seafloor is shallower than 500 m are excluded from the computation of domain averages throughout the study. This excludes about 18/19/28% of the ocean area in the Atlantic/Pacific/Med. Sea simulation. A vertical section at 53°N in the Atlantic shows that topostrophy is positive in the lower 2000 m of the water column for many locations and shows only weak vertical structure (Fig. 1b). The vertical section cuts through some of the rare observations available on the Mid-Atlantic Ridge.

Current-meter observations and simulations are in good agreement in the interior (Fig. 1a), at the vertical section (Fig. 1b), and at most other locations (Supplementary Fig. 7). Furthermore, domain averages of topostrophy as a function of height above the bottom and seafloor depth, for regions where the topostrophy is positive or negative, also show a good agreement with the rare observations (Fig. 1c and d). Here and in the following, for a comparison as a function of height above the bottom, the simulated quantities are interpolated onto the same heights above the bottom every 20 m in the lowest 1000 m and every 100 m between 1000 and 3000 m above the bottom. If the aiming height lies below the lowest simulated value, this lowest value is used for the horizontal velocity components. The vertical velocity is linearly extrapolated below the lowest data point using a value of zero at the seafloor. The amplitude of the domain-averaged topostrophy increases slightly with depth down to about 100 m above the bottom (Fig. 1c). Domain averages as a function of the seafloor depth show that the interior topostrophy weakens with increasing seafloor depth (Fig. 1d). This is true for seafloor depths of up to 5500 m. One explanation for this could be that the strength of the interior eddy field decreases on average with depth. In regions with deeper seafloor, this may lead to a weaker abyssal time-mean flow resulting in smaller abyssal topostrophy. In regions deeper than 5500 m, which are mainly the deep-sea trenches, the interior topostrophy increases however with increasing seafloor depth. About 70-90% of the measurements and of the simulated domains show positive topostrophy in the deep interior (Fig. 1e). With increasing distance from the bottom, the percentage decreases to about 50% at 3000 m above the bottom for the Pacific, Atlantic, and observations, and to less than 10% for the Med. Sea. The latter is due to the dominant cyclonic near-surface circulation in the Western Med. Sea⁴⁸.
<figure><img src="image_4.png" /></figure>

**Fig. 1** | Topostrophy $\tau = \text{sign}(\theta)(\bar{u}\partial_H/\partial y - \bar{v}\partial_H/\partial x)$, computed from three simulations (shading) and current-meter observations (circles). In (a), the interior topostrophy is shown at 300 m above the bottom for the simulations and from current-meter measurements taken between 100 and 1000 m above the bottom. If there are several current-meter measurements at one location in this layer, the deepest value is shown. In (b), a vertical section of topostrophy is shown at 53°N from the simulation of the Atlantic and from current-meter measurements in the latitudinal band 53°N ± 0.5°. The inlays right of (a, b) illustrate the along-isobath component of positive and negative topostrophic flows on the Northern Hemisphere (NH) and the Southern Hemisphere (SH). In (c), the domain-averaged topostrophy is shown as a function of height above the bottom. In (d), the domain-averaged topostrophy

300 m above the bottom is shown as a function of seafloor depth for the simulations. For the observations, averages from interior measurements between 100 m and 2000 m above the bottom are shown. Domain averages are shown separately for all positive-topostrophic and all negative-topostrophic locations. For the simulations, the volume-weighted means are shown. e shows the percentage of the positive-topostrophic area or measurements as a function of height above the bottom for the same regions. For (c-e), regions shallower than 500 m and between 5°N and 5°S are ignored. For (c, d), the distributions of the measurements are shown with gray shading in the background. Hatching indicates where no measurements are associated with the respective values. In (c, d), thin lines show the weighted 10% (90%) quantiles for negative- (positive-) topostrophic volumes in the simulations.

## Consequences of interior positive topostrophy for the near-bottom flow

For a uniform, geostrophic interior flow in a homogeneous fluid above a flat bottom, the bottom friction leads to a cyclonic deflection of the interior flow⁴⁹. This means that the direction of the flow rotates counter-clockwise (clockwise) with decreasing distance to the bottom on the Northern (Southern) hemisphere. In this case, the maximum rotation angle is 45° at the bottom⁵⁰ and the transport-weighted mean angle in the Ekman layer (the layer, where the flow is deflected due to bottom friction) is 11°. However, stratification and sloping seafloors

affect the Ekman layer structure²⁴,³²,⁴⁷. A previous collection of observations show 5–100 m thick Ekman layers that are associated with cyclonic veering angles of 10° to 20°⁵¹. In Fig. 2a the distribution and the average of the angle between the isobaths and the horizontal flow, $\alpha$, from the current-meter measurements are shown as a function of height above the bottom. $\alpha$ is positive (negative), if the flow is directed towards the shallower (deeper) seafloor and thus near the bottom indicative for uphill (downhill) flow. $\alpha$ is associated with a relatively large spread between about –45° and 45°. For heights above the bottom below 80 m, positive angles are very rare, while negative values
<figure><img src="image_5.png" /><figcaption>**Fig. 2** | Consequences of interior positive topostrophy for the near-bottom flow. **a** Shows the angle between the isobaths and the time-mean horizontal flow averaged over all measurements and the model domains in 20 m height above the bottom bins. The angle is defined as positive if the flow is directed towards shallower regions and negative if the flow is directed towards deeper seafloor depths. The inlay to the right of (a) shows the number of measurements in each 20 m height above the bottom bin. **b** shows the time-mean vertical velocity, $\bar{w}$, as a function of height above the bottom. **c** shows the 30 km</figcaption></figure>

bin-averaged $\bar{w}$ in the bottom grid cells of the simulations (shading) and estimated from current-meter measurements in the lowest 60 m of the water column with relation (2) (circles). **d** shows the distributions of the near-bottom vertical velocities. For (a, b), the distributions of the measurements are shown in the background with gray shading. Measurements and simulated data in the equatorial region between 5°S and 5°N, and in shallow regions with seafloor depths less than 500 m are ignored for (a, b). In (a, b), thin lines show the 10% and 90% quantiles that are area-weighted for the simulations.
dominate, and have larger amplitudes. The angle averaged over all measurements in 20 m bins decreases from zero (the along-isobath direction) between 80 m and 100 m to about -4° below 20 m. In summary, the near-bottom flow is mainly directed towards the deeper seafloor indicating a widespread near-bottom downhill flow.

In all three simulations, 200 m above the bottom, the flow deviates on average by about 1.25° from the along-isobath direction and is directed towards shallower seafloors (Fig. 2a). Above, α decreases to zero at about 2000 m above the bottom in the Atlantic and the Pacific simulations. Below, the horizontal velocity vector rotates cyclonically towards the bottom, crossing the along-isobath direction at about 100 m in the Pacific and Atlantic simulations, and at about 30 m in the Med. Sea simulation. This suggests that, on average, the ocean also flows downhill near the seafloor in the simulations. Below 100 m, the average α from the observations agrees very well with that of the Atlantic simulation, where the bulk of the observations in the bottom 100 m come from. This gives evidence that the relevant near-bottom dynamics are reasonably well resolved by the simulation. The 10% and 90% quantiles from the simulations agree with the observations.

The upward-positive vertical velocity, w, shows strong near-bottom downwelling and weak interior upwelling in the simulations when averaged in time over the domain (Fig. 2b). Similar to α, the interior upwelling in the Atlantic and the Pacific decreases to zero at 2000 m above the bottom. The near-bottom zero crossings also occur at similar, though slightly greater, heights above the bottom as in the case of α. The vertical velocity can not be observed directly by standard moored acoustic Doppler current profilers, as the four beams are set-up to measure only the horizontal velocities and the resulting signal to noise ratio is not sufficiently large for meaningful vertical velocities⁵². The bottom boundary condition gives

$$ \bar{w} = - \bar{u}^h \cdot \nabla^h H, \quad (2) $$

where $\bar{u}^h = (\bar{u}, \bar{v})$, H is the seafloor depth, and $\nabla^h$ is the horizontal gradient operator. Although this relation is only true at the very bottom, it can be used to estimate $\bar{w}$ from $\bar{u}^h$ in the lowest 60 m with a median

error of less than 40%, which increases with distance from the bottom (see Supplementary section 2). This error is large enough that the amplitudes of the resulting vertical velocities must be interpreted with caution, but small enough to estimate the sign of the vertical velocity. The domain-averaged estimate of $\bar{w}$ from the observations clearly shows downwelling in the lowest 60 m of the water column (Fig. 2b). The amplitude of the estimate is about an order of magnitude larger than in the simulations, partly because the observations are biased towards energetic regions.

The simulated $\bar{w}$ is widespread downwards near the bottom (Fig. 2c), particularly in regions of strong interior positive topostrophy away from the equator (Fig 1a). For Fig. 2c, small-scale variations in the simulated $\bar{w}$ in the bottom cells have been removed by bin-averaging over -30 km × 30 km subdomains (5 × 5 grid cells in the Pacific, 10 × 10 grid cells in the Atlantic, 20 × 20 grid cells in the Med. Sea) to highlight the dominance of near-bottom downwelling. Consistently, most observations in the lowest 60 m of the water column also show downward estimated vertical velocities (circles in Fig 2c). However, there are also observations showing uphill near-bottom flows, especially over flat seafloor, where also the simulations show very weak and sometimes positive near-bottom vertical velocities. Distributions of near-bottom $\bar{w}$ show that near-bottom downhill flows are more widespread than uphill flows (Fig. 2d). This skewness is further enhanced for the 30 km bin-averages for which 73/69/85% of the area of the Atlantic/Pacific/Med. Sea, that is deeper than 500 m, is associated with near-bottom Eulerian downwelling.

## Overturning above sloping seafloor

The along-isobath component of interior positive-topostrophic flows is stronger above larger seafloor slopes. At 300 m above the bottom, especially for slopes greater than 1% which are present in most regions, the along-isobath velocity increases with increasing slope (Fig. 3a). For the Pacific simulation, this is also true for relatively small slopes below 1%. Observations between 100 and 2000 m follow the Atlantic simulation curve. Consistent with the fact that the strength of the idealized Ekman flow is a linear function of the interior flow⁵⁰, and consistent

<figure><img src="image_6.png" /><figcaption>**Fig. 3 | The impact of the seafloor steepness.** a shows the time- and domain-averaged interior along-isobath velocity at 300 m above the bottom from the simulations and averaged over the measurements taken between 100 and 2000 m above the bottom as a function of the seafloor slope. b shows the time-mean vertical velocity, w, in the bottom grid cells from the simulations and averaged over the measurements in the lowest 30 m of the water column. The inlays show the distribution of the slope under the measurements and for the simulations. The</figcaption></figure>

distributions of the measurements are shown in the background with gray shading. Hatched lines indicate that there are no measurements associated with these values. Thin lines show the weighted 10% and 90% quantiles for the simulations. Measurements and simulated data in the equatorial region between 5°S and 5°N, and in shallow regions with seafloor depths less than 500 m are ignored for the averages and quantiles.
a) a sketch of the abyssal overturning cells above sloping seafloor

mixing-driven flow

<figure><img src="image_7.png" /></figure>

negative-topostrophic
interior flow
NH: into the page
SH: out of the page

<figure><img src="image_8.png" /></figure>

following
Garrett (1990)

<figure><img src="image_9.png" /><figcaption>**Fig. 4 | Overturning above sloping seafloor.** In (a), a sketch of the abyssal overturning cells is shown. Abbreviations mean divergence (DIV), convergence (CONV), Northern Hemisphere (NH), and Southern Hemisphere (SH). In the right side of (a), the abyssal flow predicted by a bottom boundary layer model solely driven by mixing²⁴ is illustrated. The domain-averaged divergence in the Atlantic (b) and Pacific (c) simulations are shown as a function of cross-isobath seafloor curvature and mean depth within each curvature band (shading). Arrows show domain averages of vertical velocities and of the cross-isobath velocity component. Note</figcaption></figure>

that orange and black arrows are differently scaled. The respective scales are shown in the inlay of (b). Gray contours show the amplitude of the domain-averaged along-isobath velocity component in cm s⁻¹. Volumes between 5°S and 5°N, in shallow regions with seafloor depths less than 500 m, with cross-isobath seafloor slopes <0.01, and above continental slopes (where the bathymetry convolved with a 200 km top-hat kernel is associated with slopes greater than 0.006) are ignored for the domain averages.

with relation (2), the near-bottom downwelling is also stronger the greater the seafloor slope in both the simulations and the estimates from observations using relation (2) (Fig. 3b). Again, the vertical velocities and their increase with slope are larger for slopes steeper than 1%.

The slope of typical bathymetric features locally increases with depth in their shallower part and decreases in their deeper part. This, combined with the corresponding change in the downhill flow, implies shallow divergence and deep convergence (Fig. 4a). To visualize the domain-averaged divergence and circulation above the slopes in the Atlantic and the Pacific simulation, the domain averages are first computed as a function of height above the bottom and the cross-

isobath seafloor curvature

$$C = \frac{-\frac{\partial^2 H}{\partial a^2}}{\left(1 + \frac{\partial H^2}{\partial a^2}\right)^{\frac{3}{2}}}, \quad (3)$$

where *a* is the direction towards greater depth ($\frac{\partial H}{\partial a} > 0$). The curvature is negative, if the seafloor becomes steeper with depth ($\frac{\partial^2 H}{\partial a^2} > 0$) and positive, if the seafloor becomes flatter with depth ($\frac{\partial^2 H}{\partial a^2} < 0$). Second, the domain-averaged seafloor depth is computed as a function of *C*. If
this is done for the whole basins, the continental slope introduces a strong minimum in depth around $C=0$ as it is associated with large areas of moderate slopes that are partially located at much lower depth compared to most other regions. For visualization purposes, we thus exclude the continental slope regions from the domain averages shown in Fig. 4 by ignoring regions where the bathymetry, low-pass filtered with a 200 km top-hat kernel, is associated with slopes larger than 0.006. The overturning cells are however also present over the continental slopes. Furthermore, we ignore relatively flat regions with cross-isobath slopes less than 1%, where $\bar{w}$ is small, noisy and shows a weak increase with increasing slope (Fig. 3b). For the included domain, the average seafloor depths for the Atlantic and the Pacific simulation show an increase with increasing $C$ (Fig. 4). Note that $C$ and not the shape of the average seafloor depths is representative of the change in seafloor steepness. The former is deformed by the non-linearly distributed seafloor and its derivatives. Finally, the domain-averaged quantities are shown in Fig. 4b, $c$ at this average depth plus the respective height above the bottom.

As expected, in the lowest 100 m, where the flow is downhill (Fig. 2b), it is associated with a shallow divergence where $C<0$ and a deep convergence where $C>0$. This is consistent with the increasing (decreasing) strength of the interior flow (contours in Fig. 4b and c). The near-bottom divergence leads to downwelling and convergence in the interior. The near-bottom convergence leads to upwelling and divergence in the interior. In the interior, upwelling and downwelling extend about 1000 m into the water column and are connected by an upward recirculation (Fig. 4b and c). Near the bottom, the thickness of the downwelling layer decreases with increasing mean depth, due to the decrease in divergence and the subsequent increase in convergence of the downhill flow. The domain-averaged abyssal slope overturning cell exists in both the Atlantic and the Pacific simulation in a similar way (Fig. 4b and c). However, in the Pacific, there is an anomaly around $C = -1 \pm 0.5$, where the upward recirculation branch is very weak and the along-isobath speed shows a local maximum. The reasons of these deviation remain unclear. We recomputed the domain averaged vertical velocity as a function of height above the bottom and sea-floor curvature but included regions of the continental slopes and multiplied it by the respective areas. The maximum of the so derived vertical transport integrated over the deepest 200 m, at heights above the bottom where it is directed downward, provides a rough estimate of the total transport associated with the abyssal slope overturning cells: 314 Sv in the Atlantic and 589 Sv in the Pacific.

## Discussion

The here revealed time-mean abyssal overturning cells above sloping seafloor are widespread in the global ocean. If the driver of the overturning cells, the time-mean interior flow, is mainly a consequence of the interactions of mesoscale eddies, it is expected to follow mainly the meso- and larger-scale bathymetry and not the smaller bathymetric features. Thus, the here discovered abyssal overturning cells are part of the abyssal flow averaged over mesoscale horizontal distances. Locally, the abyssal flow can differ strongly from the overturning cell structure. Furthermore, there are regions of sloping seafloor where the overturning cells might not develop, for example, if the mesoscale interior eddies do not extend down to the bottom boundary layer or if topographic blocking suppresses its near-bottom deflection. This is for example the case in canyons of submesoscale width in the region of the Mid-Atlantic Ridge$^{28}$. In such areas, mixing is expected to induce near-bottom uphill flow and interior negative topostrophy$^{24}$ and therefore abyssal overturning cells of opposite direction$^{25}$. In regions of strong mixing and weak interior flow, the mixing-induced overturning might overcome the interior-induced overturning. Both, topographic blocking and strong mixing, occur in thin canyons with strong tidal activity such as the ones around the Rockall Trough$^{18}$. Although our simulations of the Atlantic with 3-km grid spacing do not

resolve such canyons well, they show time-mean upward flow and reduced topostrophy near the bottom there (see section "Numerical Simulations"). More moored measurements accompanied by mesoscale dye-release experiments in regions away from narrow canyons and strong mixing are necessary to improve the ability to validate ocean simulations and to disentangle the roles of mixing and interior eddies on the abyssal circulation. Also more recent acoustic Doppler current profiler with a fifth beam, dedicated to measure the vertical velocity directly$^{52}$, will lead to further insight in this area.

The major difference between the widely held mixing-driven view presented in Garrett (1990) (Fig. 4a) and our findings lies in the following hypothesis: Garrett (1990) assumes that the ocean is at rest above a layer whose dynamics are driven by mixing processes$^{24}$. These processes, in turn, drive a negative-topostrophic flow that extends out of the bottom mixed layer. On the contrary, our study is based on the observational evidence that the ocean is not at rest above the seafloor but that there is a mean flow that is dominantly positive-topostrophic and that drives the specific dynamics in the bottom boundary layer through frictional processes. In essence, in the vertical plane, Garrett (1990) proposes a "bottom-up" vision$^{24}$ whereas we propose a "top-down" vision. Our study focuses on the Eulerian flow and its orientation to the sloping bathymetry and not on diapycnal velocities and the advection of buoyancy. The implications of positive-topostrophic interior and near-bottom downhill flow for the abyssal buoyancy balance as well as the associated diapycnal velocities and buoyancy advection are investigated in a companion study$^{32}$.

Regions with very strong negative-topostrophic boundary currents, such as the Gulf Stream south of Cape Hatteras, the Brazil Current, or the West Australian Current, show negative interior topostrophy both in observations and simulations (Fig. 1a). This indicates that the gyre-scale forced deep extensions of the boundary currents overcome the effect of the positive-topostrophy inducing interior dynamics in these regions. Furthermore, in the Pacific simulation, regions with steep slopes are associated with negative interior topostrophy. The reason for this deviation remains unclear. A sensitivity test in the form of a parallel simulation of the Atlantic, using a doubled horizontal grid spacing, which is in the horizontal almost identical to that of the Pacific run, shows no negative topostrophy above steep slopes (see section "Numerical Simulations").

In the tropics, where the planetary vorticity goes to zero, the interior topostrophy is either weak or of spatially varying sign (Fig. 1a). The near-bottom vertical velocities are mainly upward (Fig. 2c) if a $K$-profile parameterization (kpp) is used for the mixing in the simulation, where $K$ is a mixing coefficient. But the near-bottom $\bar{w}$ is found to be downward if a $K-\epsilon$ parameterization is used, where $\epsilon$ is the turbulent kinetic energy (section "Numerical Simulations"). This points to a strong mixing-dependence of the abyssal tropical circulation. Away from the tropics, we expect that the abyssal overturning cells occur in most parts of the global ocean. This includes also the regions that have not been directly addressed in the present study due to the limited computational power, such as the Indian Ocean, the Southern Ocean, and the Arctic Ocean. The reduction in mesoscale eddy size at high latitudes may result in the abyssal time-mean flow aligning with smaller-scale bathymetric features compared to lower latitudes. This would be particularly relevant to the Arctic Ocean.

Differences between the abyssal circulation in the three investigated domains of the Atlantic, Pacific, and the Med. Sea can be attributed to differences in the model configurations. With decreasing horizontal grid-spacing (6 km in the Pacific, 3 km in the Atlantic, and 1.5 km in the Med. Sea), the topostrophy (Fig. 1c) as well as the near-bottom downward vertical velocity (Fig. 3b) increases due to a better representation of eddy dynamics and bathymetry. The refinement of the vertical resolution near the bottom increases from the Atlantic over the Med. Sea to the Pacific simulation and improves the representation of the bottom boundary layer. The height above the bottom,
at which the upward recirculation changes on average to the down-
ward near-bottom flow decreases accordingly (Fig. 2b). The effects of
tides, mixing parameterization, grid spacing, as well as bottom
roughness on the abyssal circulation in an ocean model are investi-
gated in more detail with a series of sensitivity simulations and are
presented in the “methods” section. The results consistently show that,
on average, the ocean flows downhill near the bottom and recirculates
upward above. The effect of tides on near-bottom vertical velocity and
topostrophy is complex and depends on the region and the mixing
parameterization. Using the kpp mixing scheme is associated with
mainly small effects of tides, while using the K-ε mixing scheme shows
larger differences of both signs. Refining the models grid-spacing or
using a finer bathymetry leads to larger topostrophy and stronger
abyssal overturning cells.

The here presented overturning cells potentially affect the
transport of near-bottom materials such as sediments, larvae, litter,
and microplastics as they are transported by the moving water, as well
as the vertical exchange of heat, salt, and dissolved gases, which are of
importance for climate and ecosystems. In the future, this importance
needs to be quantified and respective effects need to be embedded in
climate and Earth system simulations either through resolving them
explicitly or through a parametrization of their effects.

Methods

Quantities from current-meter measurements

For this study, we analyze time-mean horizontal velocity components derived from current-meter measurements (https://www2.whoi.edu/site/aomip/data/validation-data/currents-data/) and seafloor depth gradients derived from the SRTM15+^53 bathymetry product (ftp://ftp.nodc.noaa.gov/nodc/archive/arc0090/0150537/1.1/data/0-data/topo15.grd). The current-meter data is a compilation based on three previous data collections^36,54,55. For this study, we only use the time averages of measurement series that are longer than 3 months. The

<figure><img src="image_10.png" /><figcaption>**Fig. 5 | Measurement distribution.** The number of measurements as a function of height above the bottom (20 m bins) in the simulated domains of the Atlantic, the Pacific and the Mediterranean Sea. Measurements between 5°S and 5°N and in regions shallower than 500 m are not counted.</figcaption></figure>

bathymetric gradients are identified as the nearest value to those of
the SRTM15+ bathymetry, previously low-passed with a 10 km top-hat
convolution kernel (see also section “Numerical Simulations”). For
most analysis, regions between 5°S and 5°N and in regions shallower
than 500 m are excluded. The number of measurements in all other
regions in the model domains are compared in Fig. 5 as a function of
height above the bottom. There are not many measurements in the
Med. Sea. In particular near the bottom but also for many other height
above the bottom bins, more measurements are available in the
Atlantic than in the Pacific. Measurements from five moorings near the
North pole (θ > 89°N), that are all associated with heights above the
bottom of more than 1400 m, have been excluded from all analysis to
simplify the convolution procedure.

Numerical simulations

The simulations of the Atlantic, Pacific, and Med. Sea were selected because they use terrain-following vertical sigma levels, a relatively large ocean domain, and relatively small horizontal- and vertical grid spacings. These simulations were designed to realistically simulate the interior flow and its interaction with the bottom boundary layer over sloping terrain. The Atlantic and the Med. Sea simulations are performed with the Coastal and Regional Ocean Community model (CROCO, 10.5281/zenodo.7415055), which is based on the Regional Ocean Modeling System (ROMS^56). The Pacific run is performed with the UCLA version of ROMS^57,58. The simulations have been integrated on an Arakawa C-grid^59. The grids are orthogonal based on an oblique Mercator projection and have almost uniform spacings in both horizontal dimensions. The horizontal grid spacing varies from 2.574/4.230/1.496 to 3.500/6.250/1.500 km for the grid of the Atlantic/Pacific/Med. Sea simulations. In the vertical, 100/100/60 sigma-levels where used with different level distributions (Fig. 6a-c). The vertical grid of the Pacific simulation is associated with more vertical levels near the bottom than the Atlantic simulation.

The bathymetry is taken from the SRTM30+ dataset⁶⁰. Slopes are limited to 0.2 dH/H, where dH is the difference in seafloor depth from one grid-cell to the other, to reduce pressure gradient errors⁶¹. The resulting seafloor slope and depth are shown for the three simulations in Fig. 7.

The effect of bottom friction is parameterized with a logarithmic/logarithmic/quadratic law using a roughness length of 0.01/0.02/0.001 m. No explicit lateral diffusivities and viscosities are used. Mixing is parameterized with kpp, which includes a Richardson number based parameterization in the interior⁶². Hourly atmospheric forcing datasets from the Climate Forecast System Reanalysis with a horizontal resolution of 0.2" have been used in the Atlantic and the Med. Sea simulation and from ERAS (fifth generation European Center for Medium-Range Weather Forecasts reanalysis for the global climate and weather⁶³) with a horizontal resolution of 0.25" in the Pacific simulation. Wind stresses are estimated using bulk formulae⁶⁴, including a

<figure><img src="image_11.png" /><figcaption>**Fig. 6 | The distribution of vertical levels for the simulations.** The distribution of vertical levels for the simulations of the Pacific (a), the Mediterranean Sea (b), the Atlantic with 3 km horizontal grid-spacing and 100 vertical levels (c), the Atlantic with 6 km horizontal grid-spacing and 50 vertical levels (d), the Atlantic nest in the</figcaption></figure>

Southern Hemisphere (SH) (e), the Atlantic nest in the Northern Hemisphere (NH) for the investigation of the impact of the resolution (f), and the Atlantic nest in the NH for the investigation of the impact of bottom roughness (g). The number in brackets gives the number of vertical levels.
<figure><img src="image_12.png" /><figcaption>**Fig. 7 | Bathymetries.** Seafloor slope (a) and depth (b) from the Atlantic, Pacific, and Mediterranean Sea simulations. The boundaries of the domains are highlighted with black lines. Note that in the Pacific (Atlantic) simulation, the Atlantic (Pacific) part of the domain is not simulated.</figcaption></figure>

current feedback on surface wind stress⁶⁵. Initial and boundary conditions were interpolated from Simple Ocean Data Assimilation⁶⁶ for the Atlantic simulation and from MERCATOR-OCEAN global solution GLORYS (12v1) reanalysis⁶⁷ for the Pacific and the Med. Sea. In the Pacific (Atlantic) simulation, tidal motions of the ten (eight) dominant astronomical components are forced through the combination of a surface tidal geopotential force that includes self attraction and loading, and barotropic currents tidal fluctuations at the lateral boundaries. They are derived from TPXO9 for the Pacific simulation and from TPXO7.2 and GOT99.2b (ntrs.nasa.gov/citations/19990089548) for the Atlantic simulation⁶⁸. The Med. Sea simulation is not forced with tides.

The Pacific simulation was initialized in September 2011 and integrated until the end of 2015. The results above are shown for the time-mean over the last 3 years 2013–2015. The Med. Sea simulation was initialized on January 1st 2010 and ran for 4 years. The results above are also shown for the time-mean over the last 3 years 2011–2013. The Atlantic simulation was initialized in January 2004 and run for 3 years without tides. On model day July 1st 2005, a parallel simulation with tides has been initialized from the restarts of the simulation without tides and has been also run until the end of 2006. For the analysis above, the time-means for the year 2006 from the run with tides are presented. The comparison of bottom topostrophy for the year 2006 in the Atlantic runs with and without tides shows hardly any difference (Fig. 8a–c). However, the amplitude of the bottom vertical velocity is found to be slightly enhanced with tides (Fig. 8d–f).

Several sensitivities to changes in the model configuration are investigated by comparing additional parallel simulations. First, a parallel Atlantic simulation was initialized in January 2004 without tides using a $K-\epsilon$ mixing parameterization with a Canuto A stability function instead of kpp for the parameterization of mixing. This run has been integrated until the end of 2012. On March 25th 2007, a parallel run with tides has been initialized from the run without tides and integrated until the end of 2011. The comparison of the bottom topostrophy with and without tides for the period 2009–2011 shows larger differences compared to the kpp-runs (Fig. 9a–c). Tides increase the near-bottom topostropy in most regions. However, along the

crests of the Mid-Atlantic Ridge and a few parts of the continental slopes, tides decrease the near-bottom topostrophy. Upward bottom vertical velocities averaged over the same period occur more often with tides, in particular above flat topography (Fig. 9d–f). Above the Mid-Atlantic Ridge, bottom vertical velocities are reduced. On the other side, there are also regions with enhanced downhill flow, in particular above the continental slope in regions where the topostrophy is enhanced with tides.

For regions away from the tropics (excluding 5°S–5°N) that are deeper than 500 m, the domain-averaged topostrophies of the runs with kpp (with and without tides) and the run with $K-\epsilon$ without tides are similar (Supplementary Fig. 1a). The run with tides and $K-\epsilon$ shows on average enhanced topostrophies in the lowest 2000 m compared to the other runs. The vertical velocity profile as a function of height above the bottom is similar for all four runs when averaged over the same region (Supplementary Fig. 1b). The fraction of positive topostrophy as a function of height above the bottom is also similar for all four runs (Supplementary Fig. 1c).

In all four simulations, the location of a recent dye release experiment in a canyon of the Rockall trough west of Ireland¹⁸ (black dots in Figs. 8 and 9) is located in an area of reduced positive or negative topostrophy (inlays in Figs. 8 and 9). Although the canyon is not well resolved with a 3 km horizontal grid-spacing, in all four simulations the canyon is associated with regions of upward near-bottom vertical velocities as observed with the dye release experiment and indicated by current-meter measurements¹⁸. The upward near-bottom vertical velocity is strongly enhanced when tides are simulated. These results indicate that when the mixing effect on the abyssal circulation locally dominates that of the interior eddies, the abyssal positive topostrophy and near-bottom downwelling are reversed to negative topostrophy and near-bottom upwelling. Strong mixing and a topographic blocking that suppresses the near-bottom deflection of the interior-flow may contribute to the near-bottom uphill flow, observed and simulated in narrow canyons²⁸. However, away from narrow canyons, usually the effect of interior eddies dominates the effect of mixing as the interior topostrophy is found to be widespread positive and the near-bottom vertical velocity is usually downhill at
<figure><img src="image_13.png" /><figcaption>**Fig. 8** | Tidal impact using *K*-profile mixing parameterization (kpp) where *K* is a mixing coefficient. Topostrophy (a–c) and vertical velocity (d–e) are shown in the bottom grid-cells averaged over model year 2006 in an Atlantic-wide simulation without tides (a, d) and a parallel simulation with tides (b, e), as well as their difference (c, f). Mixing has been parameterized with kpp in both runs. The inlays</figcaption></figure>

show a zoom into the Rockall trough canyon. The location of a dye-release experiment¹⁸ is marked with a black dot in the inlays and the large maps. Contour lines in the inlays show the model bathymetry with a contour interval of 500 m with the deeper seafloor in the North.

least when horizontally averaged over the mesoscale range (see also Fig. 10). The abyssal overturning cells, computed in the same way as for Fig. 4b, show only very small differences for the four runs (Supplementary Fig. 2). Note, however, that in the simulation with tides and *K*-ε the domain-averaged overturning cell is weaker than in the other runs. A point-wise comparison of the simulation-based topostrophies at the measurement locations with those derived from the current-meter measurements shows that there are only small differences between the experiments at the measurement locations (Supplementary Fig. 3). Runs with *K*-ε mixing parameterization show slightly smaller RMSE than those with kpp. Runs without tides show slightly smaller RMSE than those with tides. The latter is unexpected as the simulation of tides should increase the realism of the simulations. This, the small RMSE differences, and the sparseness of the observations indicate that a conclusion on whether the kpp or the *K*-ε mixing parameterization leads to more realistic results cannot be drawn from this database.

Second, to investigate the effect of resolution, a parallel Atlantic simulation without tides and with *K*-ε mixing scheme has been initialized in January 2004 and has been integrated until the end of 2012 using a doubled grid-spacing of about 6 km and 50 vertical levels (Fig. 6d). The lower resolution run shows for the period 2008–2012 weaker time-mean topostrophy at 300 m above the bottom almost in the whole domain with amplitudes closer to those of the Pacific simulation, which uses also a grid-spacing of about 6 km (Fig. 1 and Supplementary Fig. 4). This is expected, as the higher horizontal

resolution is associated with a better representation of the interior eddies that are thought to induce the interior mean flow, as well as a better representation of the bathymetric slopes.

Additionally, two nested simulations are initialized from the Atlantic runs with *K*-ε mixing scheme. For the nested simulations, SRTM15+ instead of SRTM30+ is used for the bathymetry. The nested runs use the kpp mixing scheme. For the time-means of the nested simulations presented in the following, the first 10 days have been excluded. The first nested simulation was initialized from the basin-wide Atlantic run without tides on February 23rd 2008 and integrated for 111 days with a horizontal grid-spacing of about 1 km for a domain above the Mid-Atlantic Ridge in the Southern Hemisphere (thin violet box in Supplementary Fig. 4). For this simulation without tides, 300 vertical levels have been used with a distribution shown in Fig. 6e. The second nested simulation was initialized from the basin-wide Atlantic run with tides on August 26th 2008 and integrated for 230 days with a horizontal grid-spacing of about 800 m for a domain above the Mid-Atlantic Ridge in the Northern Hemisphere (thin violet box in Supplementary Fig. 4). For this simulation with tides, 200 vertical levels have been used with a distribution shown in Fig. 6f. Furthermore, to test the effect of small-scale bathymetry and pressure gradient effects, for this simulation, the bathymetry has been smoothed using a two-dimensional Gaussian kernel of a width that is twelve times the horizontal grid-spacing in contrast to the standard of four times the horizontal grid-spacing. To further investigate the effect of resolution,
<figure><img src="image_14.png" /><figcaption>**Fig. 9** | Tidal impact using *K-ε* mixing parameterization where *K* is a mixing coefficient and *ε* is the turbulent kinetic energy. The Figure is similar to Fig. 8 but shows averages over the model period 2009–2011 from two parallel simulations that use the *K-ε* mixing parameterization instead of the *K*-profile mixing</figcaption></figure>

parameterization. Topostrophy (a–c) and vertical velocity (d–e) are shown in the bottom grid-cells for an Atlantic simulation without tides (a, d) and a parallel simulation with tides (b, e), as well as their difference (c, f).

area-averages for subdomains in the nested regions are compared in Supplementary Fig. 5. All simulations show time-mean abyssal positive topostrophies and near-bottom downwelling. The higher the model resolution, the larger is the average positive topostrophy, the interior upward flow, and the near-bottom downward flow. However, the average vertical velocity in the 3 km Atlantic-wide run is very close to the one in the 800 m and the 1 km nested simulations. The sign change from interior upwelling to near-bottom downwelling occurs for the 3 km run at the same height above the bottom as for the nested runs. This and the comparable amplitudes indicate that the near-bottom vertical circulation is on average well represented in the 3 km run although the vertical resolution is relatively coarse (Fig. 6c). The nested run in the Northern Hemisphere with very high near-bottom vertical resolution shows dominantly positive topostrophy and downward flow in the bottom grid-cells (Fig. 10b, f). On average, the bottom vertical grid-spacing is only 1.7 m in this run.

Third, to investigate the effect of bottom roughness, we perform two additional nested simulations above the Mid-Atlantic Ridge in the Northern Hemisphere (thin violet box in Supplementary Fig. 4) that are run in parallel to the second nested simulation described above. Both are integrated for 50 days and use a different distribution of vertical levels compared to the second nested simulation (Fig. 6g). The third nested simulation is else identical to the second. The fourth nested simulation is identical to the third, besides that the bathymetry is smoothed only with the standard four times the horizontal grid-

spacing Gaussian kernel. Both smoothing procedures ensure that $dH/H < 0.2$, where $dH$ is the difference in seafloor depth from one grid-cell to the other, which is the classic criterion to keep pressure gradient errors within acceptable levels⁶¹. For simplicity, the second nested simulation is referred to as smoothed (long run) in the following, the third to as smoothed, and the fourth to as unsmoothed. For the analyzed time-mean we again ignore the first 10 days of the simulations. Smoothing the bathymetry leads to less strong positive topostrophy, reduced amplitudes of vertical velocities, and to a larger-scale structure of both fields (Fig. 10). Although there are plenty of regions with upward flow in the bottom grid-cells, on average, the near-bottom vertical velocity is strongly negative. Note that the short length of the experiments of only 50 days and in consequence the very short averaging period of 40 days may not be long enough so that the positive-topostrophic time-mean flow is fully established. A comparison of the smoothed run to the smoothed (long run), where the near-bottom vertical resolution is additionally enhanced, shows less regions with negative topostrophy (Fig. 10b and c) and less strong near-bottom upwelling in the respective regions (Fig. 10e–g).

For the computation of topostrophy from the current-meter measurements, seafloor gradients from the SRTM15+ bathymetry with a horizontal grid-spacing of 1/240' have been used. It is found that the number of measurement-derived negative topostrophies is reduced and the number of positive topostrophies enhanced, if the SRTM15+ bathymetry is smoothed with a two-dimensional 10 km top-hat kernel
<figure><img src="image_15.png" /></figure>

**Fig. 10** | **The impact of bottom roughness.** Time-mean topostrophy (a–d) and vertical velocity (e–h) from three 800 m horizontal resolution and 200 vertical level nested simulations above the Mid-Atlantic Ridge in the Northern Hemisphere (thin violet box in Supplementary Fig. 4). Domain averages are shown as a function of height above the bottom in (a, e); values in the bottom grid-cells in the other subplots. A 230 day run with smoothed bathymetry and very high near-bottom

vertical resolution (Fig. 6f) smoothed (long run) is compared to two parallel 50 day
runs with less vertical levels near the bottom (Fig. 6g) one with smoothed the
other one without additional bathymetry smoothing (unsmoothed). From each
simulation, the first 10 days are not included in the averaging period. All three
simulations simulated tides and use K-profile mixing parameterization, where K is a
mixing coefficient.

before the computation of the topostrophy (Supplementary Fig. 6). This indicates that the time-mean positive-topostrophic interior flow follows a slightly smoothed bathymetry. On the one hand side, this could be due to the fact that the time-mean flow is mainly driven by mesoscale interior eddies. On the other hand side, it might be that the flows inertia does not allow it to be sensitive to the seafloor roughness under a specific threshold scale. Smoothing with a 20 km kernel instead of a 10 km kernel leads to the same effect and additionally to a large percentage of values near zero which points to a compression of the distribution that maybe induced by less extreme amplitudes due to the stronger smoothing. A comparison to the distribution of the topostrophy in the Atlantic at the measurement locations shows that the latter agrees very well with the one of the observations when the SRTM15+ bathymetry is smoothed with a 10 km kernel (Supplementary Fig. 6a).

Probability density functions for the comparison of the topos-
trophy at the grid points closest to the measurement locations show a
reasonable agreement between observations and simulations when
the SRTM15+ bathymetry is slightly smoothed before the computation
of topostrophy (Supplementary Fig. 7). The root-mean-square errors
(RMSE) drop from 1.11/0.93/4.29 to 0.46/0.46/0.41 for the Atlantic/
Pacific/Med. Sea simulations if the SRTM15+ bathymetry is smoothed
with a 10 km top-hat kernel. Here, the results from the Atlantic simu-
lation with tides and kpp mixing parameterization is shown and thus
from the simulation for which the results are presented in the main
part of the study. For the Atlantic and the Pacific simulation, the RMSE
slightly drops further, if a 20 km kernel is used instead of the 10 km
kernel. For all three simulations, the amplitudes of the simulated
topostrophies are on average slightly larger than those from the
observations.

Data availability

The data to reproduce the main figures of this study can be accessed at
https://doi.org/10.24433/CO.7868381.v1. The ocean model output
data is too large for an online upload. Contours showing coastlines

have been created with www.naturlearthdata.com/downloads/110m-
physical-vectors/.

Code availability

The code that has been applied to observational and numerical model
data can be accessed here: https://doi.org/10.24433/CO.7868381.v1.

References

1. Yu, L. & Weller, R. A. Objectively analyzed air-sea heat fluxes for the global ice-free oceans (1981–2005). *Bull. Am. Meteorol. Soc.* **88**, 527–540 (2007).

2. Wanninkhof, R., Asher, W. E., Ho, D. T., Sweeney, C. & McGillis, W. R. Advances in quantifying air-sea gas exchange and environmental forcing. *Annu. Rev. Mar. Sci.* **1**, 213–244 (2009).

3. Liu, W., Xie, S.-P. & Lu, J. Tracking ocean heat uptake during the surface warming hiatus. *Nat. Commun.* **7**, 10926 (2016).

4. Zanna, L., Khatiwala, S., Gregory, J. M., Ison, J. & Heimbach, P. Global reconstruction of historical ocean heat storage and transport. *Proc. Natl. Acad. Sci. USA* **116**, 1126–1131 (2019).

5. Marshall, J. & Speer, K. Closure of the meridional overturning circulation through southern ocean upwelling. *Nat. Geosci.* **5**, 171–180 (2012).

6. Munk, W. & Wunsch, C. Abyssal recipes II: energetics of tidal and wind mixing. *Deep Sea Res. Part I Oceanogr. Res. Pap.* **45**, 1977–2010 (1998).

7. De Lavergne, C., Madec, G., Le Sommer, J., Nurser, A. G. & Garabato, A. C. N. On the consumption of Antarctic Bottom Water in the abyssal ocean. *J. Phys. Oceanogr.* **46**, 635–661 (2016).

8. de Lavergne, C., Groeskamp, S., Zika, J. & Johnson, H. L. The role of mixing in the large-scale ocean circulation. in *Ocean Mixing* 35–63 (Elsevier, 2022).

9. Polzin, K., Toole, J., Ledwell, J. & Schmitt, R. Spatial variability of turbulent mixing in the abyssal ocean. *Science* **276**, 93–96 (1997).

10. Laurent, L. S. et al. Turbulence and diapycnal mixing in Drake Passage. *J. Phys. Oceanogr.* **42**, 2143–2152 (2012).
11. Waterhouse, A. F. et al. Global patterns of diapycnal mixing from measurements of the turbulent dissipation rate. *J. Phys. Oceanogr.* **44**, 1854–1872 (2014).

12. Ferrari, R., Mashayek, A., McDougall, T. J., Nikurashin, M. & Campin, J.-M. Turning ocean mixing upside down. *J. Phys. Oceanogr.* **46**, 2239–2261 (2016).

13. McDougall, T. J. & Ferrari, R. Abyssal upwelling and downwelling driven by near-boundary mixing. *J. Phys. Oceanogr.* **47**, 261–283 (2017).

14. Holmes, R. & McDougall, T. J. Diapycnal transport near a sloping bottom boundary. *J. Phys. Oceanogr.* **50**, 3253–3266 (2020).

15. Callies, J. Restrafitation of abyssal mixing layers by submesoscale baroclinic eddies. *J. Phys. Oceanogr.* **48**, 1995–2010 (2018).

16. Wenegrat, J. O., Callies, J. & Thomas, L. N. Submesoscale baroclinic instability in the bottom boundary layer. *J. Phys. Oceanogr.* **48**, 2571–2592 (2018).

17. Drake, H. F. et al. Dynamics of eddying abyssal mixing layers over sloping rough topography. *J. Phys. Oceanogr.* **52**, 3199–3219 (2022).

18. Wynne-Cattanach, B. L. et al. Observations of diapycnal upwelling within a sloping submarine canyon. *Nature* **630**, 884–890 (2024).

19. Simpson, J. H. & McCandliss, R. R. "The Ekman Drain": a conduit to the deep ocean for shelf material. *Ocean Dyn.* **63**, 1063–1072 (2013).

20. Harris, P. T. Shelf and deep-sea sedimentary environments and physical benthic disturbance regimes: a review and synthesis. *Mar. Geol.* **353**, 169–184 (2014).

21. Vic, C., Gula, J., Roullet, G. & Pradillon, F. Dispersion of deep-sea hydrothermal vent effluents and larvae by submesoscale and tidal currents. *Deep Sea Res. Part I: Oceanogr. Res. Pap.* **133**, 1–18 (2018).

22. Peng, G., Bellerby, R., Zhang, F., Sun, X. & Li, D. The ocean's ultimate trashcan: Hadal trenches as major depositories for plastic pollution. *Water Res.* **168**, 115121 (2020).

23. Wunsch, C. On oceanic boundary mixing. In *Deep Sea Research and Oceanographic Abstracts* Vol. 17, 293–301 (Elsevier, 1970).

24. Garrett, C. The role of secondary circulation in boundary mixing. *J. Geophys. Res.: Oceans* **95**, 3181–3188 (1990).

25. Ruan, X. & Callies, J. Mixing-driven mean flows and submesoscale eddies over mid-ocean ridge flanks and fracture zone canyons. *J. Phys. Oceanogr.* **50**, 175–195 (2020).

26. Peterson, H. G. & Callies, J. Rapid spinup and spindown of flow along slopes. *J. Phys. Oceanogr.* **52**, 579–596 (2022).

27. Clément, L., Thurnherr, A. M. & Laurent, L. C. S. Turbulent mixing in a deep fracture zone on the Mid-Atlantic Ridge. *J. Phys. Oceanogr.* **47**, 1873–1896 (2017).

28. Lahaye, N. et al. Deep currents in the rift valley of the north mid-Atlantic ridge. *Front. Mar. Sci.* **6**, 597 (2019).

29. Thurnherr, A. M., Clément, L., Laurent, L. S., Ferrari, R. & Ijichi, T. Transformation and upwelling of bottom water in fracture zone valleys. *J. Phys. Oceanogr.* **50**, 715–726 (2020).

30. Xie, X. et al. Observation of near-inertial waves in the bottom boundary layer of an abyssal seamount. *J. Phys. Oceanogr.* **53**, 635–645 (2023).

31. Wei, B. et al. Rapid down-slope transport of fresh dissolved organic matter to the deep ocean in the Eastern North Atlantic. *Geophys. Res. Lett.* **51**, e2024GL110349 (2024).

32. Capó, E. et al. Abyssal slope currents. *J. Phys. Oceanogr.* **54**, 2373–2392 (2024).

33. Holloway, G. A shelf wave/topographic pump drives mean coastal circulation. *Ocean Model.* **68**, 12–15 (1986).

34. Holloway, G. Representing topographic stress for large-scale ocean models. *J. Phys. Oceanogr.* **22**, 1033–1046 (1992).

35. Holloway, G. et al. Water properties and circulation in Arctic Ocean models. *J. Geophys. Res. Oceans* **112**, Nr. C4, (2007).

36. Holloway, G. Observing global ocean topostrophy. *J. Geophys. Res. Oceans* **113**, Nr. C7, (2008).

37. Voet, G., Quadfasel, D., Mork, K. A. & Søiland, H. The mid-depth circulation of the Nordic Seas derived from profiling float observations. *Tellus A Dyn. Meteorol. Oceanogr.* **62**, 516–529 (2010).

38. Bai, Y., Wang, Y. & Stewart, A. L. Does topographic form stress impede prograde ocean currents? *J. Phys. Oceanogr.* **51**, 2617–2638 (2021).

39. Bretherton, F. P. & Haidvogel, D. B. Two-dimensional turbulence above topography. *J. Fluid Mech.* **78**, 129–154 (1976).

40. Salmon, R., Holloway, G. & Hendershott, M. C. The equilibrium statistical mechanics of simple quasi-geostrophic models. *J. Fluid Mech.* **75**, 691–703 (1976).

41. Herring, J. R. On the statistical theory of two-dimensional topographic turbulence. *J. Atmos. Sci.* **34**, 1731–1750 (1977).

42. LaCasce, J., Palóczy, A. & Trodahl, M. Vortices over bathymetry. *J. Fluid Mech.* **979**, A32 (2024).

43. Merryfield, W. J. & Scott, R. B. Bathymetric influence on mean currents in two high-resolution near-global ocean models. *Ocean Model.* **16**, 76–94 (2007).

44. Penduff, T. et al. Influence of numerical schemes on current-topography interactions in 1/4 global ocean simulations. *Ocean Sci.* **3**, 509–524 (2007).

45. Maltrud, M. & Holloway, G. Implementing biharmonic neptune in a global eddying ocean model. *Ocean Model.* **21**, 22–34 (2008).

46. Tagliabue, A. et al. Mechanisms driving the dispersal of hydrothermal iron from the northern Mid Atlantic Ridge. *Geophys. Res. Lett.* **49**, e2022GL100615 (2022).

47. MacCready, P. & Rhines, P. B. Slippery bottom boundary layers on a slope. *J. Phys. Oceanogr.* **23**, 5–22 (1993).

48. Pinardi, N. & Masetti, E. Variability of the large scale general circulation of the Mediterranean Sea from observations and modelling: a review. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **158**, 153–173 (2000).

49. Ekman, V. W. On the influence of the earth's rotation on ocean-currents. *Ark. Mat. Astron. Och Fys.* **2**, 1–53 (1905).

50. Cushman-Roisin, B. & Beckers, J.-M. *Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects* (Academic Press, 2011).

51. Perlin, A. et al. Organization of stratification, turbulence, and veering in bottom Ekman layers. *J. Geophys. Res. Oceans* **112**, Nr. C5, (2007).

52. Comby, C. et al. Measuring vertical velocities with ADCPs in low-energy ocean. *J. Atmos. Ocean. Technol.* **39**, 1669–1684 (2022).

53. Tozer, B. et al. Global bathymetry and topography at 15 arc sec: Srtm15+ *Earth Space Sci.* **6**, 1847–1864 (2019).

54. Gregory, D. *Ocean Data Inventory (ODI): A Database of Ocean Current, Temperature and Salinity Time Series for the Northwest Atlantic* (Fisheries & Oceans Canada, Science, Canadian Science Advisory Secretariat, 2004).

55. Scott, R. B. et al. Total kinetic energy in four global eddying ocean circulation models and over 5000 current meter records. *Ocean Model.* **32**, 157–169 (2010).

56. Shchepetkin, A. & McWilliams, J. C. The Regional Oceanic Modeling System (ROMS): a split-explicit, free-surface, topography-following- coordinate ocean model. *Oceanogr. Meteorol.* **9**, 347–404 (2005).

57. Shchepetkin, A. F. & McWilliams, J. C. A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate. *J. Geophys. Res. Oceans* **108**, Nr. C3, (2003).

58. Molemaker, M. J., Dollery, D., Damien, D. & McWilliams, J. UCLA ROMS https://doi.org/10.5281/zenodo.8256776 (2023).

59. Arakawa, A. & Lamb, V. R. Computational design of the basic dynamical processes of the UCLA general circulation model. *Gen. Circ. Models Atmos.* **17**, 173–265 (1977).

60. Becker, J. et al. Global bathymetry and elevation data at 30 arc seconds resolution: SRTM30_PLUS. *Mar. Geod.* **32**, 355–371 (2009).
61. Auclair, F., Marsaleix, P. & Estournel, C. Sigma coordinate pressure gradient errors: evaluation and reduction by an inverse method. *J. Atmos. Ocean. Technol.* **17**, 1348–1367 (2000).

62. Large, W. G., McWilliams, J. C. & Doney, S. C. Oceanic vertical mixing: a review and a model with a nonlocal boundary layer parameterization. *Rev. Geophys.* **32**, 363–403 (1994).

63. Skamarock, W. C. et al. A description of the advanced research WRF version 3. *NCAR Tech. note* **475**, 113 (2008).

64. Large, W. B. Surface fluxes for practitioners of global ocean data assimilation. *Ocean Weather Forecasting: An Integrated View of Oceanography* 229–270 (Springer, 2006).

65. Renault, L. et al. Modulation of wind work by oceanic current interaction with the atmosphere. *J. Phys. Oceanogr.* **46**, 1685–1704 (2016).

66. Carton, J. & Giese, B. A reanalysis of ocean climate using Simple Ocean Data Assimilation (SODA). *Mon. Weather Rev.* **136**, 2999–3017 (2008).

67. Lellouche, J.-M. et al. Recent updates to the Copernicus Marine Service global ocean monitoring and forecasting real-time 1/12 high-resolution system. *Ocean Sci.* **14**, 1093–1126 (2018).

68. Egbert, G. D. & Erofeeva, S. Y. Efficient inverse modeling of barotropic ocean tides. *J. Atmos. Ocean. Technol.* **19**, 183–204 (2002).

## Acknowledgements

R.S. and J.G. gratefully acknowledge support from the French National Agency for Research (ANR) through the project DEEPER (ANR-19-CE01-0002-01). J.G. further gratefully acknowledges support from ISblue "Interdisciplinary graduate school for the blue planet" (ANR-17- EURE-0015). The Atlantic simulations were performed using high performance computing resources from GENCI-TGCC (Grants 2024-A0170112051), and from DATARMOR of "Pôle de Calcul Intensif pour la Mer" at Ifremer Brest France. Support for this research was provided by the Office of Naval Research (ONR grants #N00014-18-1-2599 P00001 and #N00014-23-1-2226) and the National Science Foundation (NSF grant #2241754). Computing support for the Mediterranean Sea and Pacific simulations was provided by the NSF ACCESS program.

## Author contributions

R.S., J.G., E.C., P.D., M.J.M., and J.C.M. designed the study. J.G., E.C., and P.D. performed the numerical model simulations. R.S. developed and performed all other analysis, produced the figures and wrote the main

part of the text. All authors, including C.V., contributed to the discussion of the results and to the text.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-61027-2.

**Correspondence** and requests for materials should be addressed to René. Schubert.

**Peer review information** *Nature Communications* thanks the anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025