The formation and ventilation of an oxygen minimum zone in a simple model for latitudinally alternating zonal jets

An advection–diffusion model coupled to a simple dynamical ocean model is used to illustrate the formation and ventilation of an oxygen minimum zone. The advection–diffusion model carries a passive tracer with a source at the western boundary and a Newtonian damping term to mimic oxygen consumption. The dynamical model is a non-linear one-and-a-half-layer reduced-gravity model. The latter is forced by an annually oscillating mass flux confined to the near-equatorial band that, in turn, leads to the generation of mesoscale eddies and latitudinally alternating zonal jets at higher latitudes. The model uses North Atlantic geometry and develops a tracer minimum zone remarkably similar in location to the observed oxygen minimum zone in the eastern tropical North Atlantic (ETNA). This is despite the absence of wind forcing for a subtropical gyre and the shadow zone predicted by the ventilated thermocline theory. Although the model is forced only at the annual period, the model nevertheless exhibits decadal and multidecadal variability in its spun-up state. The associated trends are comparable to observed trends in oxygen within the ETNA oxygen minimum zone. Notable exceptions are the multi-decadal decrease in oxygen in the lower-oxygen minimum zone and the sharp decrease in oxygen in the upper-oxygen minimum zone between 2006 and 2013.

1 Introduction

As the ocean warms due to anthropogenic forcing, oxygen levels in the ocean are expected to decrease, putting pressure on the ocean ecosystem and fisheries (Stramma et al., 2008; Keeling et al., 2010; Schmidtko et al., 2017; Breitburg et al., 2018; Oschlies et al., 2018). In particular, this deoxygenation process is expected to lead to an expansion of the oxygen minimum zones (OMZs). These are large naturally occurring volumes of low-oxygen waters situated in the subsurface of the subtropical oceans (e.g. Karstensen et al., 2008). Yet, the fate of the OMZs depends on the delicate interplay between biological oxygen consumption and physical oxygen supply. One critical aspect is therefore to understand how oxygen is supplied to OMZs and how this ventilation mechanism varies in time. Figure 1a shows the mean oxygen distribution (black contours) at roughly 500 m depth in the tropical Atlantic Ocean. Prominent are the regions of low-oxygen concentration in the eastern basin separated by a band of higher concentration centred on the Equator and the much higher levels of concentration found near the western boundary. The off-equatorial regions of low concentration at the eastern boundary are commonly associated with the shadow zones predicted by the theory for the ventilated thermocline (Luyten et al., 1983; Brandt et al., 2015) and are regions with sluggish ocean circulation that are weakly ventilated. Brandt et al. (2015) discuss the eastern tropical North Atlantic (ETNA) OMZ, north of the Equator, in some detail. These authors note the important role played by latitudinally alternating zonal jets for supplying oxygen to this region. These jets can be seen in Fig. 1a and appear as the alternating bands of eastward (red) and westward (blue) zonal velocity in a multi-year average. The importance of these jets has also been demonstrated by Calil (2023), who found a much better representation of the Atlantic OMZs in a high-resolution model simulation that supports latitudinally alternating zonal jets compared to a lower-resolution simulation that does not. Furthermore, Delpech et al. (2020) have noted the importance of these jets for shaping the structure of tracer fields in the tropical Pacific. Brandt et al. (2015) document a multi-decadal downward trend in ETNA oxygen concentration that could be a consequence of anthropogenically driven ocean deoxygenation (Stramma et al., 2008, 2010; Keeling et al., 2010; Schmidtko et al., 2017). Brandt et al. (2015) note, however, that there are also stronger trends, both upward and downward, on the decadal timescale (see Fig. 1b and c taken from Fig. 19 in Brandt et al., 2015). It has been suggested that these trends are associated with changes in the zonal jets and, in particular, the strength of these jets (Brandt et al., 2010; Hahn et al., 2017).

Figure 1 (a) Contours show the mean oxygen concentration (µmol kg⁻¹) in the tropical Atlantic at 500 m depth (close to the deep-oxygen minimum) as obtained from the World Ocean Atlas 2018 (Garcia et al., 2019). Colour shading shows the 1993–2016 time mean zonal velocity from the GLORYS12v1 reanalysis product averaged between the 454 and 541 m depth levels (CMEMS, 2023). (b, c) Oxygen anomalies for the region 9–15° N, 20–26° W (orange box in panel a) for 150–300 m (upper OMZ; panel b) and 350–700 m (lower OMZ; panel c). Grey circles represent all available data, whiskers show the interquartile range of data within each year, and the black squares are annual medians. Trends are calculated using annual medians weighted by the square root of available data within each year for the periods 1900–2013 (solid red line) and 2006–2013 (solid blue line). The dashed lines mark the standard errors in the trends. Panels (b) and (c) are taken from Fig. 19 in Brandt et al. (2015).

Latitudinally alternating zonal jets, as seen in Fig. 1a, are an ubiquitous feature of the ocean (Maximenko et al., 2005, 2008). Their identification using Argo float data is shown in Fig. 3 of Ascani et al. (2015), and their central role in the tropical ocean circulation is illustrated in Fig. 1 of Ménesguen et al. (2019), who also discuss possible mechanisms for supporting these jets. One of these mechanisms involves the net eddy-induced westward mass flux by mesoscale eddies and argues that the zonal jets are required to provide a net eastward mass flux in order to conserve mass (Marshall et al., 2013).

Here, we work with the simplest possible dynamical model that can generate mesoscale eddies and that, in turn, supports zonal jets, with the aim to shed light on both the formation of the ETNA OMZ and the decadal-to-interdecadal variations in oxygen concentration within the OMZ. The model is a non-linear one-and-a-half-layer reduced-gravity model and is very similar to that used by Qiu et al. (2013) to illustrate the generation of zonal jets in the tropical Pacific Ocean through Rossby wave resonant triad interactions. While the model includes the realistic coastlines of the North Atlantic, we note that very similar results are obtained using rectangular basin geometry, as shown by Köhn (2018), and that the effect of varying ocean bathymetry is excluded. The dynamical model is coupled to a simple advection–diffusion model for a passive tracer that includes a tracer source near the western boundary and a simple Newtonian relaxation term to mimic oxygen consumption. By coupling to the advection–diffusion model, we illustrate how the latitudinally alternating zonal jets can advect oxygen from the western boundary into the interior of the ocean, with their effect diminishing eastwards and facilitating the formation of an OMZ near the eastern boundary. Due to their inherent non-linearity originating in the mesoscale eddy field, these jets can exhibit their own internal variability, even on decadal timescales, analogous to what is known about the atmospheric circulation (James and James, 1989). There is, therefore, the potential to influence oxygen concentration in the OMZs on similar timescales.

The structure of the paper is as follows. Section 2 describes the model set-up and forcing, Sect. 3 presents the results, and Sect. 4 provides a summary and discussion.

2 Methods

2.1 The model

The governing equations for the dynamic model are those for a non-linear one-and-a-half-layer reduced-gravity model, namely

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial u}{\partial t}}-qhv=-{\displaystyle \frac{g^{\prime }}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial E}{\partial \mathit{\lambda }}}+M_u,\text{(2)}& {\displaystyle \frac{\partial v}{\partial t}}+qhu=-{\displaystyle \frac{g^{\prime }}{a}}{\displaystyle \frac{\partial E}{\partial \mathit{\theta }}}+M_v,\text{(3)}& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\theta }}}\left[{\displaystyle \frac{\partial \left(hu\right)}{\partial \mathit{\lambda }}}+{\displaystyle \frac{\partial \left(hv\cos \mathit{\theta }\right)}{\partial \mathit{\theta }}}\right]=F_{\mathrm{h}}.\end{array}$$

Here (λ,θ) are longitude and latitude, respectively; (u,v) are the eastward and northward velocities, respectively; h is the active layer thickness; $E=\frac{\mathrm{1}}{\mathrm{2}}(u^{\mathrm{2}}+v^{\mathrm{2}})+g^{\prime }h$ is the Bernoulli function, where g^′ is the reduced gravity; $q=(f+\mathit{\zeta })/h$ is the potential vorticity, where $\mathit{\zeta }=\frac{\mathrm{1}}{a\cos \mathit{\theta }}(\frac{\partial v}{\partial \mathit{\lambda }}-\frac{\partial u\cos \mathit{\theta }}{\partial \mathit{\theta }})$ is the relative vorticity; f=2Ωsin θ is the Coriolis parameter; $\mathrm{\Omega }=\mathrm{7.29}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹ is the rotation rate of the Earth; and a is the radius of the Earth. (M_u,M_v) are the lateral mixing of momentum terms. These are given, following Shchepetkin and O'Brien (1996), by

$$\begin{array}{}\text{(4)}& M_u=-{\displaystyle \frac{\mathrm{1}}{h}}\left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}{P}_{\mathit{\lambda }\mathit{\lambda }}+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}{P}_{\mathit{\lambda }\mathit{\theta }}-{\displaystyle \frac{\mathrm{2}\tan \mathit{\theta }}{a}}{P}_{\mathit{\lambda }\mathit{\theta }}\right],\end{array}$$

and

$$\begin{array}{}\text{(5)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }{M}_v=\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{h}}\left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}{P}_{\mathit{\lambda }\mathit{\theta }}+{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}{P}_{\mathit{\theta }\mathit{\theta }}-{\displaystyle \frac{\mathrm{2}\tan \mathit{\theta }}{a}}(P_{\mathit{\lambda }\mathit{\lambda }}-P_{\mathit{\theta }\mathit{\theta }})\right],\end{array}\end{array}$$

where

$$\begin{array}{}\text{(6)}& {\displaystyle }{P}_{\mathit{\lambda }\mathit{\lambda }}=-A_{\mathrm{h}}h\left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}u-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}v-{\displaystyle \frac{\tan \mathit{\theta }}{a}}v\right],\text{(7)}& {\displaystyle }{P}_{\mathit{\lambda }\mathit{\theta }}=P_{\mathit{\theta }\mathit{\lambda }}=-A_{\mathrm{h}}h\left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}v-{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}u+{\displaystyle \frac{\tan \mathit{\theta }}{a}}u\right],\end{array}$$

and $P_{\mathit{\theta }\mathit{\theta }}=-P_{\mathit{\lambda }\mathit{\lambda }}$, with A_h being the horizontal eddy viscosity. F_h is the model forcing described below.

The model domain spans the latitude/longitude range 50° N, 20° S and 65° W, 15° E; uses the realistic coastlines of the Atlantic Ocean at the western and eastern boundaries; and has a horizontal resolution of 0.1°. The northern and southern boundaries are closed, as is the Strait of Gibraltar, and a no-slip boundary condition is used with a horizontal eddy viscosity of A_h=100 m² s⁻¹. The model is set up with an undisturbed layer depth of H₀=500 m and a reduced gravity of $g^{\prime }=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{2}}$ m s⁻², giving a gravity wave speed $c=\sqrt{g^{\prime }{H}_{\mathrm{0}}}=\mathrm{2.7}$ m s⁻¹ that is within the range of first baroclinic mode gravity wave speeds for the tropical Atlantic (Chelton et al., 1998). The forcing corresponds to a mass flux in and out of the dynamically active layer that is confined near the Equator and oscillates at annual period $\frac{\mathrm{2}\mathit{\pi }}{\mathit{\omega }}$ and is given by

$$\begin{array}{}\text{(8)}& F_{\mathrm{h}}(\mathit{\lambda },\mathit{\theta })=\left\{\begin{array}{ll}-A_{\mathrm{0}}\sin \left(\mathit{\omega }t\right)\exp [-\mathit{\beta }{y}^{\mathrm{2}}/(\mathrm{2}c\left)\right]& \\ & \mathrm{\text{if 15°\hspace{0.17em}S }}\le \mathit{\theta }\le \text{15°\hspace{0.17em}N}\\ \mathrm{0}& \mathrm{\text{otherwise}},\end{array}\right.\end{array}$$

where $A_{\mathrm{0}}=\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, corresponding to ≈2.6  m d⁻¹. y is the meridional distance from the Equator, and β is the meridional gradient of f. The e-folding distance around the Equator, $\sqrt{\mathrm{2}c/\mathit{\beta }}$, is 4.6° latitude corresponding to the equatorial radius of deformation. This choice of forcing is a simple and efficient way to induce an annual cycle of layer thickness and zonal velocity along the Equator (in reality associated with the annual cycle of wind stress) and is not intended to be realistic.

The dynamical model is coupled to an advection–diffusion model for which the governing equation is

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Ch}{\partial t}}& {\displaystyle }+\mathrm{\nabla }\cdot \left(h\mathbf{u}C\right)=\\ {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \left({\mathit{\kappa }}_{\mathrm{h}}h\mathrm{\nabla }C\right)-JCh-h\mathit{\gamma }(C-C_{\mathrm{0}})+C{F}_{\mathrm{h}},\end{array}\end{array}$$

where C is the tracer concentration (mimicking oxygen), κ_h is the tracer diffusivity (set equal to, for simplicity, the eddy viscosity A_h), and ∇⋅(…) is the divergence operator in the spherical (latitude/longitude) coordinates. The tracer source is given by the hγ(C−C₀) term, where

$$\begin{array}{}\text{(10)}& C_{\mathrm{0}}(\mathit{\lambda },\mathit{\theta })=\left\{\begin{array}{ll}\mathrm{1}& \text{if 65°\hspace{0.17em}W }\le \mathit{\lambda }\le \text{55°\hspace{0.17em}W},\\ \mathrm{0}& \text{otherwise}\end{array}\right.\end{array}$$

implying a strong source within 10° of the western boundary (because of the use of realistic coastlines, the source is confined to the region north of 10° N; see Fig. 2d). To prevent a sharp edge in the tracer concentration at the step in the C₀ field, the restoring factor γ is set up to be spatially dependent and is given by

$$\begin{array}{}\text{(11)}& \mathit{\gamma }(\mathit{\lambda },\mathit{\theta })=\left\{\begin{array}{ll}{\mathit{\gamma }}_{\mathrm{0}}/\mathrm{2}(\mathrm{1}-\tanh [\frac{\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{0}}}{{\mathit{\lambda }}_{\mathrm{s}}}\left]\right)& \\ & \mathrm{\text{if 65°\hspace{0.17em}W }}\le \mathit{\lambda }\le \text{55°\hspace{0.17em}W}\\ \mathrm{0},& \mathrm{\text{otherwise}},\end{array}\right.\end{array}$$

where λ₀=63°, λ_s=1°, and ${\mathit{\gamma }}_{\mathrm{0}}=\mathrm{1}/\mathrm{8.35}$ d⁻¹, indicating a strong restoring within a few degrees of the western boundary. Placing the source near the western boundary is consistent with both the observed oxygen distribution (Fig. 1a) and the notion that the most recently ventilated water is associated with the western boundary current system (see, for example, Fratantoni and Richardson, 2006, and Kirchner et al., 2009).

Figure 2 Spin-up showing (a–c) the layer thickness, h, and (d–f) the tracer concentration, C, after 6, 60, and 600 months, respectively.

Tracer consumption is modelled using a simple Newtonian relaxation term, the strength of which is determined by the tracer consumption rate J. Here J is set to $J=\mathrm{6.5}\times {\mathrm{10}}^{-\mathrm{10}}$ s⁻¹, is associated with a relaxation timescale of roughly 50 years, and is, as a result, much longer than the timescale associated with the source term. This value for the consumption rate is half that diagnosed by fitting a one-dimensional advection–diffusion model to observations in the Pacific by Van Geen et al. (2006) and used by Brandt et al. (2010, 2012) in studies of the oxygen field in the equatorial Atlantic over the depth range 300–500 m. Using a smaller value ensures that the tracer is able to cross the entire model basin without being damped out and is also consistent with the notion of a lower apparent oxygen utilisation rate at greater depth (Karstensen et al., 2008), where the model is likely to be most applicable. It should be noted that the use of a Newtonian relaxation term with a constant spatially uniform consumption rate is a simplification that is sufficient for our purposes in the simple model set-up used here but, nevertheless, a crude representation of the complex biogeochemical processes that determine oxygen consumption in reality. It follows that care should be taken when interpreting the tracer modelled here as a proxy for oxygen.

The equations are discretised on an Arakawa C-grid (Arakawa, 1972) using an Adams–Bashforth third-order time-stepping scheme. The details can be found in Köhn (2018). The model is initialised with a state of rest with the undisturbed layer depth of H₀=500 m and zero tracer concentration C and is then run for 400 years with a time step of 500 s.

2.2 Thickness-weighted averaging

The most natural way to average the tracer Eq. (9) is to use thickness-weighted averaging (e.g. De Szoeke and Bennett, 1993; Greatbatch and McDougall, 2003), where, for any variable α,

$$\begin{array}{}\text{(12)}& \widehat{\mathit{\alpha }}={\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\alpha }h}}{\stackrel{\mathrm{‾}}{h}}},\end{array}$$

and

$$\begin{array}{}\text{(13)}& \mathit{\alpha }=\widehat{\mathit{\alpha }}+{\mathit{\alpha }}^{\prime \prime }.\end{array}$$

Here an overbar refers to a long time average, and we note that $\stackrel{\mathrm{‾}}{h{\mathit{\alpha }}^{\prime \prime }}=\mathrm{0}$. Averaging Eq. (9) then gives

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \widehat{C}\stackrel{\mathrm{‾}}{h}}{\partial t}}=& {\displaystyle }-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h\mathbf{u}C}\right)+\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{\kappa }}_{\mathrm{h}}h\mathrm{\nabla }C}\right)\\ {\displaystyle }& {\displaystyle }-J\widehat{C}\stackrel{\mathrm{‾}}{h}-\mathit{\gamma }\stackrel{\mathrm{‾}}{h}(\widehat{C}-C_{\mathrm{0}})+\stackrel{\mathrm{‾}}{C{F}_{\mathrm{h}}}.\end{array}\end{array}$$

The first term on the right-hand side of Eq. (14) can be understood by noting that the thickness-weighted mean velocity, $\widehat{\mathbf{u}}$, is the primary velocity variable after averaging (Greatbatch and McDougall, 2003) and can be decomposed into two parts

$$\begin{array}{}\text{(15)}& \widehat{\mathbf{u}}=\stackrel{\mathrm{‾}}{\mathbf{u}}+{\mathbf{u}}^{*},\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathbf{u}}$ is the Eulerian mean velocity, and

$$\begin{array}{}\text{(16)}& {\mathbf{u}}^{*}={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathbf{u}}^{\prime }{h}^{\prime }}}{\stackrel{\mathrm{‾}}{h}}}\end{array}$$

is the eddy-induced velocity, and here $h=\stackrel{\mathrm{‾}}{h}+h^{\prime }$. The eddy-induced velocity is analogous to the Stokes drift associated with surface gravity waves (Stokes, 1847) and is an additional velocity that needs to be added to the Eulerian mean velocity to account, after averaging, for the Lagrangian advection of tracer (see Marshall et al., 2013, for a further discussion). The first term on the right-hand side of Eq. (14) can now be decomposed as

$$\begin{array}{}\text{(17)}& {\displaystyle }\begin{array}{rl}{\displaystyle }\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h\mathbf{u}C}\right)& {\displaystyle }=\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h}\widehat{\mathbf{u}}\widehat{C}\right)+\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h{\mathbf{u}}^{\prime \prime }{C}^{\prime \prime }}\right)\\ {\displaystyle }& {\displaystyle }=\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h}\stackrel{\mathrm{‾}}{\mathbf{u}}\widehat{C}\right)+\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h}{\mathbf{u}}^{*}\widehat{C}\right)+\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h{\mathbf{u}}^{\prime \prime }{C}^{\prime \prime }}\right).\end{array}\end{array}$$

The three terms on the far right-hand side correspond to the divergence of the advective flux associated with the Eulerian mean velocity, the eddy-induced velocity, and the transient eddies (the eddy-mixing term), respectively.

3 Results

In the model, the annual period forcing along the Equator generates annual period equatorial Rossby and Kelvin waves that, in turn, excite annual period coastal Kelvin waves that propagate poleward along the eastern boundary of the model domain. These, in turn, radiate annual period Rossby waves into the interior. These Rossby waves destabilise by means of a resonant triad instability (see Qiu et al., 2013, where a detailed discussion can be found) and shed eddies that subsequently propagate westward towards the western boundary. Figure 2 shows snapshots at 6, 60, and 600 months into the model simulation. The top panel shows the layer depth, h, corresponding dynamically to the pressure, and the lower panels show the tracer concentration, C. After 6 months, the front excited by switching on the forcing is clearly visible in Fig. 2a, extending southwestwards towards the Equator from the eastern boundary. The much greater westward extension towards the Equator is because the westward propagation speed for (long) Rossby waves is given by $\mathit{\beta }{g}^{\prime }h/f^{\mathrm{2}}$ (Anderson and Killworth, 1979) and hence increases rapidly towards the Equator, where f=0. At later times (Fig. 2b and c), this front has destabilised and is shedding eddies that propagate westward and already, at 600 months, fill the model domain. Concurrent with these developments (lower panels), the tracer is confined entirely to the source region near the western boundary at 6 months (Fig. 2d). At 60 months (Fig. 2e), the eddies are already penetrating the western boundary region, leading to eastward-spreading filaments of relatively high tracer concentration into the interior and away from the western boundary. At 600 months (Fig. 2f), these filaments fill the basin to the west of the destabilised front and are clearly spreading the tracer into the interior of the model domain away from the western boundary. An interesting point is that while the eddies propagate westward, the tracer filaments extend eastward. As discussed by Marshall et al. (2013), the westward propagation of the eddies is associated with a net westward mass flux associated with the eddy-induced velocity (see Eq. 16 and the discussion). Marshall et al. (2013) note that in order to conserve mass, there must be a corresponding eastward mass flux, and it is this eastward mass flux that manifests itself through the eastward extending filaments. Another interesting point is the advection of tracer along the Rossby wave front at the location where this front starts breaking down into eddies. This can be clearly seen in Fig. 2f, where the eastward-spreading filaments turn northeastward near 20° N.

Figure 3 shows the mean fields averaged over the last 80 years of model integration (that is, years 321 to 400). The tracer distribution shown in Fig. 3b is not dissimilar to that shown in Fig. 1a and, intriguingly, indicates the presence of a tracer (oxygen) minimum zone east of the destabilising front noted in Fig. 2b and c. This region of the tracer minimum is remarkably similar in appearance to the observed OMZ shown in Fig. 1a and is almost exactly colocated with the observed OMZ. In particular, the concentration reaches its minimum southeast of the Cabo Verde islands, hugs the African coast, and increases again towards the Equator. In the model, its presence is due to the lack of penetration into this region of tracer filaments spreading from the west. The mean thickness-weighted zonal velocity $\widehat{\mathbf{u}}$ shown in Fig. 3c shows alternating bands of eastward and westward flow, corresponding to the latitudinally alternating zonal jets noted in the introduction. The corresponding zonal bands in the mean layer thickness reflect the meridionally varying pressure field required to geostrophically balance the zonal jets. These jets are strongest west of the destabilising Rossby front, as we expect, and less coherent at the latitude band of the tracer minimum on the eastern boundary. The mean meridional velocity shown in Fig. 3d is very different, is of much smaller magnitude over most of the model domain, and is important only when very close to the western boundary and in the western part of the equatorial region, most likely associated with the rectification of (short) Rossby waves that are excited on the western boundary and have eastward group velocity.

Figure 3 Mean fields of (a) layer thickness, h, (b) the tracer concentration, and (c) the zonal and (d) meridional velocity components over the last 80 years of model integration. The layer depth is averaged in a normal sense, while panels (b), (c), and (d) show thickness-weighted averages.

A detailed analysis of the tracer budget for the same 80 years shows that the main balance is between the divergence of the total advective flux, given by the first term on the right-hand side of Eq. (14), and the consumption term (see Fig. 4a and b); the local time derivatives play a negligible role, consistent with the model being in a statistically steady state. This balance is particularly well established in the latitude band between the Equator and 30° N, which is where the tracer minimum is located. Splitting the divergence of the total advective flux into its separate parts shows that the dominant terms are associated with the Eulerian mean velocity, $\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h}\stackrel{\mathrm{‾}}{\mathbf{u}}\widehat{C}\right)$ (Fig. 4d), and the transient eddies, $\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{h{\mathbf{u}}^{\prime \prime }{C}^{\prime \prime }}\right)$ (Fig. 4f), with these two terms showing a strong tendency to cancel each other out. The term associated with the eddy-induced velocity (Fig. 4e) is of a much smaller magnitude. Nevertheless, as can be seen from Fig. 4a, the net effect of the advective tracer flux is to provide a source of the tracer. The importance of the Eulerian mean velocity and the transients is, perhaps, not surprising when looking at the instantaneous snapshot in Fig. 2f, with filaments extending eastwards from the western boundary region but also being advected meridionally by the transient eddies.

Figure 4 Tracer budget terms for the last 80 years of model integration. See the panel labelling for the individual panel content. The decomposed advective terms in panels (d), (e), and (f) are summed in panel (a) following Eq. (17).

Although the dynamical model is forced only at the annual period, because it is a non-linear model, it nevertheless exhibits variability on much longer (decadal to interdecadal) timescales (James and James, 1989). This low-frequency variability is reflected in the tracer concentration, as can be seen in Fig. 5b. Variations in tracer concentration of up to 5 % from the mean take place on decadal to interdecadal timescales. Köhn (2018) carried out a detailed analysis of the trends in the model and concluded that the observed trends in the ETNA OMZ, namely those discussed by Brandt et al. (2015) and Hahn et al. (2017), are largely within the range of trends in the model. The two major discrepancies are the pronounced oxygen decrease in the upper (150–300 m depth) OMZ during 2006–2013 (an ≈14 % decrease) and the multidecadal decrease in the lower (350–700 m depth) OMZ (an ≈14 % decrease, this time over 50 years) that can be seen in Fig. 1b and c, respectively. It should, nevertheless, be remembered that the model used here is very idealised and that results from much more complex models may lead to different conclusions. The important point is the presence of decadal to multidecadal trends in the tracer concentration exhibited by the statistically steady state of the model. These trends are driven by changes in the advective forcing of the tracer concentration, consistent with the suggestion of Brandt et al. (2010).

Figure 5 (a) Thickness-weighted average of tracer concentration along 23° W over the last 240 years of model integration and (b) the temporal evolution of the corresponding anomalies in tracer concentration, which is shown as a percentage of the mean concentration shown in panel (a). The anomalies are convolved with a 4-year (i.e. 49-month) kernel to remove the seasonal cycle.

4 Summary and discussion

Calil (2023) has shown that using a high-resolution ocean model to resolve the latitudinally alternating zonal jets leads to a much better representation of the eastern tropical North Atlantic (ETNA) oxygen minimum zone (OMZ) than when these jets are not resolved. Earlier, Brandt et al. (2010) had suggested that changes in these jets, in particular their strength, might be able to explain decadal trends in the oxygen concentration within the ETNA OMZ. Here, we have presented results from a simple dynamical model that supports mesoscale eddies and the corresponding latitudinally alternating zonal jets with a view to shedding light on both the formation of the ETNA OMZ and the subsequent decadal-to-interdecadal variability in oxygen concentration within the OMZ. The simple dynamical model, which uses the geometry of the North Atlantic Ocean, is coupled to an advection–diffusion model to mimic oxygen. The tracer carried by the advection–diffusion model has a source near the western boundary of the North Atlantic and a relaxation term to mimic consumption. The results show that filaments of relatively high tracer concentration are carried by the Eulerian mean jets and the eddies eastward away from the western boundary source towards the eastern boundary. The mean tracer distribution in statistically steady state is remarkably similar to that in the observed ocean with a minimum in tracer concentration in the same general location as the observed OMZ. This region of tracer minimum is situated to the east of the Rossby wave front which destabilises and is the source of the eddies in the model, exactly as described by Qiu et al. (2013).

We noted in the Introduction that the OMZs in the ocean have been associated with the shadow zones that are predicted by the ventilated thermocline theory of Luyten et al. (1983) for the wind-driven circulation of the subtropical gyre, regions that are not ventilated in the theory and are thought to be characterised by a near-stagnant ocean. We note that in the theory, the ventilated thermocline is driven by the downward Ekman pumping due to the wind stress curl over the subtropical gyre. Nevertheless, the simple model we have described here excludes this downward Ekman pumping and is forced only near the Equator, with no forcing in the latitudinal band of the tracer minimum shown in Fig. 3b yet exhibiting an OMZ remarkably similar to the observed OMZ in the tropical North Atlantic. Indeed, the mechanism we have demonstrated offers an alternative explanation for the existence of OMZs that is independent of the ventilated thermocline theory. It is an interesting point that the ETNA OMZ we have focused on is located in a region where the mean wind stress curl has the opposite sign to that associated with the subtropical gyre (see, for example, Leetmaa and Bunker, 1978, or the more recent estimate by Risien and Chelton, 2008). It follows that the ventilated thermocline theory, which relies on the equatorward Sverdrup transport associated with a subtropical gyre, is unlikely to be applicable to the ETNA OMZ. However, this is not to say that wind forcing has no role to play in explaining the observed ETNA OMZ. In particular, Peña-Izquierdo et al. (2015), while noting the importance of zonal jets for ventilating the ETNA OMZ, note also that the upper OMZ is ventilated predominantly by water from the South Atlantic via the North Equatorial Countercurrent (NECC), a wind-driven current, and the North Equatorial Undercurrent (NEUC). Likewise, the lower OMZ was found to be ventilated predominantly by water originating in the North Atlantic subtropical gyre by means of deep zonal jets (Brandt et al., 2015) analogous to those simulated by our model. It follows that both the mechanism we have identified here and the wind-driven circulation play a role in ventilating the ETNA OMZ in reality but that the shadow zone predicted by the ventilated thermocline theory is not required to explain the ETNA OMZ. It also seems likely that the latitudinally alternating zonal jets associated with the mesoscale eddy field, such as we have discussed here, are most applicable for ventilating the lower OMZ, below roughly 250 m depth, with the direct effect of the wind-driven circulation being more important for ventilating the upper, roughly 250 m, water column. We also note that although the ventilated thermocline theory of Luyten et al. (1983) may not be relevant for explaining the ETNA OMZ, this does not mean that it does not play a role in the dynamics of other OMZs, in particular those in the Pacific Ocean.

We have shown that although the dynamical model has forcing only on the annual timescale, the modelled tracer concentration exhibits pronounced decadal and multidecadal variability. The analysis by Köhn (2018) suggests that the associated decadal and interdecadal trends offer a plausible explanation for the observed trends in the ETNA OMZ (for example, those documented by Hahn et al., 2017). Exceptions are the slow multidecadal decrease in oxygen in the lower OMZ (350–700 m depth) and the rapid decrease in the upper OMZ (150–300 m depth) observed between 2006 and 2013 (see Fig. 1b and c taken from Brandt et al., 2015). However, it should be remembered that the model we have used is very simplified and that further studies using more comprehensive models are needed. Nevertheless, we note that more comprehensive models generally underestimate rates of multidecadal oxygen decline, as well as interannual to decadal variability (Oschlies et al., 2018). One possible reason for this, supported by the work of Calil (2023), may well be an inability of the models to properly resolve the latitudinally alternating zonal jets.

Finally, we note that Chang et al. (2008) have argued that the amplitude of the annual cycle in the equatorial Atlantic can be expected to reduce as the strength of the Atlantic Meridional Overturning Circulation reduces under anthropogenic forcing. While acknowledging that the forcing we use to drive our model is very idealised, this suggests that a future climate scenario might be mimicked by reducing the amplitude of the forcing applied to the model. Qiu et al. (2013) have shown that in this case, the front, where the westward propagating Rossby waves destabilise, moves westward. Such a westward displacement would, in turn, reduce the ventilation of the OMZ in our model, expanding the OMZ westward. It is possible, therefore, that the ventilation of the OMZs in the ocean by the latitudinally alternating jets will be reduced in the future, contributing to the deoxygenation of the ocean that is expected as a result of anthropogenic climate change (Stramma et al., 2008, 2010; Keeling et al., 2010; Schmidtko et al., 2017).