1 Introduction

Equatorial deep jets (EDJs) were first observed in the Indian Ocean by Luyten and Swallow (1976) and later in the Atlantic and Pacific by Leetmaa and Spain (1981) and Eriksen (1982), respectively. The EDJs are a major feature of the equatorial current system and consist of a series of vertically stacked zonal jets that alternate in direction with depth. The Atlantic EDJs have a vertical length scale of around 300–700 m and their period is estimated to be 5 ± 1 years by Johnson and Zhang (2003), later refined to 4.5 years by Brandt et al. (2011) with an uncertainty of about 150 days estimated by Claus et al. (2016) for the, then available, 10 year time series. The observed downward phase propagation of the EDJs implies, according to linear theory, an upward energy propagation, a feature confirmed in a nonlinear model study by Matthießen et al. (2015) (see their Fig. 4). Several authors have stressed the impact of the EDJs on the local ocean and climate system. For example, Gouriou et al. (2001) indicated that some North Atlantic Deep Water is transported eastward along the equator by the EDJs and Brandt et al. (2012) identified the EDJs as a possible pathway for oxygen to the tropical oxygen minimum zones. Indeed, Dietze and Loeptien (2013) argued that the persistent problem of too little oxygen in the eastern equatorial oceans in Earth System Models can be attributed to the failure of models to correctly represent the equatorial current system, including the EDJs. Getzlaff and Dietze (2013) supported this argument by using an enhanced zonal diffusivity along the equator, which led to improved distributions of temperature, salinity, oxygen and nutrients in the equatorial region in their model. Brandt et al. (2011) showed the indirect impact of the EDJs on local climate through variations in the sea surface temperature (SST) and Matthießen et al. (2015) further highlighted the North Equatorial Countercurrent (NECC) where there appears to be a surface impact of the EDJs affecting the meandering of the current. The lack of EDJs in Earth System Models noted above is not limited to these particular models. Even high resolution, realistic ocean models have difficulty representing EDJs (Eden and Dengler 2008; Ascani et al. 2015). Nevertheless, idealised model set-ups do support EDJs (e.g. D’Orgeville et al. 2007; Ménesguen et al. 2009; Ascani et al. 2015; Matthießen et al. 2015) and hopefully analysis of these model runs can shed light on why realistic models are unsuccessful in simulating EDJs.

Studies on the theoretical background of the EDJs have led to their description in terms of resonant basin modes (Cane and Moore 1981; Yang and Liu 2003) with characteristic period \(T=\frac {4L}{c_{n}}\), where L is the basin width and c n the gravity wave speed associated with the dominant vertical normal mode (the time taken for a Kelvin wave and the reflected gravest long Rossby wave to cross the basin). The analogy with basin modes was first noted by Ascani et al. (2006) in an idealised model experiment but has since become generally accepted as a model for the observed EDJs (D’Orgeville et al. 2007; Greatbatch et al. 2012; Claus et al. 2014; Youngs and Johnson 2015; Claus et al. 2016). Indeed, basin modes appear to be important for describing not only the EDJs in the equatorial Atlantic Ocean but also the semi-annual and annual cycles (Brandt et al. 2016). Ascani et al. (2015) also noted that in their nonlinear model the EDJs are well represented by the gravest basin mode for a high baroclinic vertical normal mode. Furthermore, Matthießen et al. (2015) showed that in their idealised model experiments, most of the energy associated with the zonal velocity along the equator has both a time scale and vertical structure associated with resonant basin modes.

There are different theories for the generation of EDJs. D’Orgeville et al. (2007) showed, in an idealised study, that destabilisation of Yanai waves (also called Mixed Rossby Gravity waves) generated near the western boundary and that have eastward group velocity into the interior of the ocean, can drive a high baroclinic mode structure following the argument of Hua et al. (2008). Eden and Dengler (2008) also investigated the generation of vertically stacked jets in a high resolution ocean model where the energy source was thought to be fluctuations of the Deep Western Boundary Current in the model. Ascani et al. (2015) described a mechanism whereby Yanai waves shed from tropical instability waves interact to support the EDJs. In particular, a high vertical mode intraseasonal Yanai wave, with the same vertical scale as the EDJs, interacts with a low vertical mode intraseasonal Yanai wave to maintain the deep jets through the meridional advection term in the momentum equation. This mechanism implies that deep jets are supported by the intraseasonal waves over a considerable depth range. In fact Claus et al. (2016) have diagnosed the forcing for the deep jets from the mooring data at 23∘W on the equator and find that the implied power input that maintains the deep jets takes place over the whole depth range occupied by the jets. Concerning the vertical phase propagation of the jets, we note that in the particular case studied by D’Orgeville et al. (2007), the vertical phase propagation is not in agreement with observations, although it should be noted that their model runs are shorter than those considered by Ascani et al. (2015) and Matthießen et al. (2015) and possibly do not show fully developed EDJ cycles. In D’Orgeville et al. (2007), the phase propagation is such as to take energy away from the depths where the imposed destabilizing Yanai waves have relatively large amplitude to depths where the destabilizing Yanai waves have relatively low amplitude (see Fig. 6 in their paper and the discussion thereon). Furthermore, the vertical symmetry of their experiment (the experiment is continually forced by a second baroclinic mode, intraseasonal Yanai wave, uses vertically uniform stratification and free slip/no normal flux boundary conditions at the top and bottom boundaries) suggests that this would not change in a longer model run. We also note a distinction between a generation mechanism, that sets the scale of the deep jets, like the theories of D’Orgeville et al. (2007) or Hua et al. (2008), and a sustaining mechanism, that can maintain the deep jets in a quasi-equilibrium, as in Ascani et al. (2015) or Muench and Kunze (2000). The two mechanisms are thought to operate together and are not in any way contradictory to each other. Indeed, the mechanism described by Ascani et al. (2015) somewhat resembles, and may simply be another way to interpret, the resonant triad mechanism described in Appendix A of Hua et al. (2008).

In this paper, we revisit the numerical model setup used by Ascani et al. (2015) and Matthießen et al. (2015) and focus on the spin-up process through which the EDJs emerge in the model. The model setup is described in Section 2 and the model analysis is presented in Section 3. Section 4 provides a summary and discussion.

2 Model setup

We are using the MITgcm (Marshall et al. 1997) in five different configurations (see Table 1). The time stepping uses the Adams-Bashforth scheme and the advective terms use a second order centred (non-dissipative) scheme (see Marshall et al. (1997) for the details). The basic configuration (EDJ_Ref) is the same as the wide basin case (Model A) in Matthießen et al. (2015) and differs from Ascani et al. (2015) only in the model used (they used the POP model, we use the MITgcm) and the number of vertical levels (they used 100, we use 200). Indeed, the model results obtained from EDJ_Ref do not differ significantly from those of Ascani et al. (2015) during the overlap period.

Table 1 Model configurations used in this study
Full size table

The model domain is a rectangular basin, extending 72∘ in longitude and 40∘ in latitude, centred on the equator, with a flat bottom and vertical walls on all sides and a horizontal resolution of \(\frac {1}{4}^{\circ } \times \frac {1}{4}^{\circ }\). There are no sponge layers used in our model configurations. The model has 200 levels in the vertical, with a finer scale at the top and a more coarse resolution at the bottom and a total depth of 5000 m. Vertical mixing is parameterised following Pacanowski and Philander (1981) with background vertical diffusivity and viscosity of 10− 5 m2 s−1 and biharmonic mixing is used in the horizontal for both tracers and momentum with coefficient − 2 × 1010 m4 s−1. An additional sensitivity experiment is described in which the Pacanowski and Philander (1981) is replaced by the KPP mixing scheme (Large et al. 1994; see Table 1 for a complete list of the experiments). Only potential temperature is used to set the density stratification (the salinity is uniform and constant throughout the integrations) and a linear equation of state is used. The initial profile, in experiments initialised with a state of rest, is obtained by area averaging potential temperature at each depth from the World Ocean Atlas (Levitus et al. 2013) within the region of the Atlantic Ocean corresponding to the model domain. The surface temperature in all model runs is relaxed to its initial condition with a timescale of 30 days. The wind stress used to drive the model is time independent and is the annual mean zonally-averaged meridional and zonal wind stress computed from the NCEP reanalysis (Kalnay et al. 1996).

The reference experiment, EDJ_Ref, is initialised from a state of rest, as in Ascani et al. (2015) and Matthießen et al. (2015). We also include an experiment, EDJ_Mean, initialised with the density field averaged over the last 30 years of EDJ_Ref, without spatial averaging (but with the velocity set to zero, as before), in order to investigate the sensitivity of the emerging EDJs in the model to the choice of the initial state. A third configuration, EDJ_Low, is the same as EDJ_Ref except that the wind stress is scaled by a factor of 10− 3 in order to reduce the importance of the nonlinear terms. It should be noted that this model run also uses the Pacanowski and Philander (1981) mixing scheme and, as such, since the velocities are very weak, the vertical mixing for both momentum and potential temperature is given by the background value of 10− 5 m2 s−1 and hence is very weak. A fourth configuration EDJ_Const is the same EDJ_Ref but using a uniform stratification for the initial state with N = 0.003 s−1 corresponding to a potential temperature of 27 ∘C at the surface and 4 ∘C at the bottom of the basin (see Table 1). The advantage of using a uniform stratification is that the vertical structure functions for the vertical baroclinic modes can be computed analytically allowing us to make an analytical reconstruction of the spin-up. A fifth configuration, EDJ_KPP, is the same as EDJ_Ref except that KPP mixing scheme is used.

3 Results

Regarding the first experiment EDJ_Ref, Fig. 1 shows the zonal velocity in the middle of the basin directly on the equator as a function of depth and time. The Equatorial Undercurrent (EUC) is visible as a strong eastward current near the surface; the slow downward drift is due to the long time span of the model run and the slow diffusion of the thermocline. Below the EUC are the EDJs; that is, the system of zonal jets that vertically alternate direction with depth and time. Our interest here is mostly in the initial spin-up period during the first 50 years of the model run during which energy appears to first propagate downward and then be reflected upward from the ocean bottom in the model. The first episode of the spin-up, from the beginning of the model run to year 12, shows mostly upward phase propagation, while the second episode, from year 12 to year 50, is dominated by downward phase propagation. After year 50, the model run remains dominated by downward phase propagation, although the amplitude is subject to some decadal modulation. The downward phase propagation in the model is consistent with the observations of Brandt et al. (2011) from their mooring data at 23∘W on the equator.

Fig. 1
figure 1

Hovmoeller diagram showing zonal velocity in m s−1 as a function of time and depth at a location near the centre of the basin on the equator in experiment EDJ_Ref. The black dashed line marks year 12 of the model run

Full size image

Figure 2a, b shows snap-shot vertical sections across the equator of zonal velocity from EDJ_Ref at different phases of the deep jet cycle. The stacked zonal jets within 2∘ either side of the equator are visible and are broadly consistent with what we know from observations (compare with Fig. 1b in Eden and Dengler 2008). The most obvious difference is that the deep jets in the model extend over the whole depth range whereas in the observations they are confined above the depth reached by the Mid-Atlantic Ridge (about 2000 m). Figure 2c shows a snap-shot of the zonal velocity along the equator from the model. This shows the strong zonal coherence of the jets and also the slight tendency to tilt down towards the east, consistent with Fig. 6 in the reconstruction of Claus et al. (2016) and the analysis of Youngs and Johnson (2015) based on data from Argo floats.

Fig. 2
figure 2

Snapshots of zonal velocity from EDJ_Ref in m s−1. Panels a and b show meridional sections across the equator in the middle of the basin. Panel a shows a section during year 17 of the model run and panel b during year 27. Panel c shows a section along the equator during year 27 of the model run

Full size image

As in the observations, the EDJs in EDJ_Ref are mostly composed of basin modes (Fig. 3). To produce Fig. 3, the normal mode decomposition, as described in Gill (1982), has been applied to

$$ u(z,t)-u_{mean}(z,t)=\sum\limits_{n} a_{n}(t) \hat{p}_{n}(z), $$
(1)

where the \(\hat {p}_{n}(z)\) are the vertical structure functions calculated from the initial density field by solving the Sturm-Liouville problem

$$ \frac{d}{dz}\left[\frac{1}{N^{2}}\frac{d \hat{p}_{n}}{dz}\right] + \frac{1}{{c_{n}}^{2}}\hat{p}_{n} = 0 $$
(2)

with boundary conditions \(\frac {1}{N^{2}}\frac {d\hat {p}_{n}}{dz} = 0\) at the surface z = 0 and at the bottom z = −H appropriate for the baroclinic modes (Gill 1982). The vertical structure functions are normalised so that \(\int \limits _{-H}^{0} \hat {p}_{n}(z)\hat {p}_{m}(z)dz=\delta _{nm}H\), where δ n m = 1 when n = m and is zero otherwise. u(z,t) is the zonal velocity shown in Fig. 1 and u m e a n (z,t) is a 30-year running mean of this zonal velocity and is used to remove the signal of the EUC (before year 15, the average at time t is taken from 0 to 2*t ). As is clear from Fig. 3, most of the energy sits along the line corresponding to the gravest resonant basin mode for each baroclinic normal mode, with the energy peaks between time scales of 3.5 and 5 years corresponding to the EDJs. Note that the whole model run excluding the first 15 years was used to compute Fig. 3.

Fig. 3
figure 3

Fourier spectra of the zonal velocity in m s−1 as a function of vertical mode number at the location on the equator used to produce Fig. 1. The blue line is the period of the gravest basin mode, T = 4L/c n , where c n is the gravity wave speed associated with the nth vertical normal mode

Full size image

Figure 4 focuses on the spin-up period and compares the zonal velocity at the equator in the fully nonlinear run, EDJ_Ref, against that in the low amplitude run EDJ_Low for which the amplitude of the wind forcing was reduced by a factor of 10− 3. When looking at Fig. 4, it should be noted that the model output from EDJ_Low has been scaled up by a factor 103 so that the amplitude in the two plots can be compared directly. Also, whereas the reference case (top panel) uses a standard linear scale, the low amplitude case (lower panel) is plotted using a logarithmic scale in order to include both the intense (when the velocity is scaled) currents near the surface and the weak zonal flow (associated with basin modes) below. A feature of the two plots is their similarity in shape, despite the much weaker role for the non-linear terms in the low amplitude case. We shall return to this point later when we attempt to reconstruct the spin-up shown here using linear superpositions of basin modes. At the same time, the difference in amplitude seen in the plots below the top few hundred metres clearly shows the role played by the nonlinear terms in amplifying the model response in EDJ_Ref compared to EDJ_Low. This is consistent with the central role played by nonlinearity in the formation of the EDJs in the theory of D’Orgeville et al. (2007), Hua et al. (2008) and Ascani et al. (2015). We also note that the analyses presented in the studies of Ascani et al. (2015) and Matthießen et al. (2015) were confined to the second episode of the spin-up (years 12–50), characterised by downward phase propagation and did not consider the episodes of downward phase propagation after the spin-up (after year 50) that can be seen in Fig. 1.

Fig. 4
figure 4

Hovmoeller diagram showing zonal velocity in m s−1 as a function of time and depth at the same location on the equator as in Fig. 1. a The first 50 years from EDJ_Ref and b the first 50 years from EDJ_Low, where the velocity is multiplied by a factor of 103 (the inverse factor of what was used to reduce the wind stress used to drive this experiment). The velocity in b is plotted using a logarithmic scale to highlight the deep jet structure when the amplitude is low

Full size image

Concerning the relatively strong currents near the surface in EDJ_Low, we note that these correspond to the EUC in this model run but rather than being a uni-directional current, the model EUC here consists of several bands of currents that alternative in direction, eastward and westward, with a very short vertical scale. Both the relatively large amplitude and the unusual vertical structure are a consequence of the very weak vertical mixing in this experiment (recall that the Pacanowski and Philander (1981) mixing scheme was used and that because the forcing has been scaled down by 103, the velocities used for the calculation of the Richardson Number are very weak, leaving only the background eddy viscosity and diffusivity as the vertical mixing). Nevertheless, in the actual model run itself, the maximum amplitude of the EUC is only 0.03 m s−1 and the EUC is always stable with no tropical instability waves in this experiment. Likewise, as can be seen from Fig. 4b, the actual model velocities in EDJ_Low below the EUC are tiny, of order 0.0001 m s−1 or less. Figure 5 compares the zonal velocity in both model runs vertically averaged over the region below the EUC (noting that, as in Fig. 4b, the velocity in EDJ_Low has been scaled up by a factor of 103 in order to enable comparison with EDJ_Ref). From this, it is clear that the scaled velocities in EDJ_Low are relatively higher during the first 10 years than in EDJ_Ref and that at later times, relatively much higher zonal velocities are found in EDJ_Ref associated with the deep jets. The decline in the zonal velocity in the low amplitude case at depth is associated with the adjustment of the low amplitude model, as in a linear dynamical adjustment, towards Sverdrup balance, as described by McCreary (1981), with the EUC being here strongly confined near the surface. In particular, because of the weak dissipation in this model run, the weak velocities at depth are not a consequence of dissipation but rather of the linear dynamical adjustment. Concerning the higher (scaled) velocities in EDJ_Low compared to EDJ_Ref in the early years of the model runs, it should be noted that in EDJ_Ref the vertical mixing is much higher than in EDJ_Low, implying much more dissipation.

Fig. 5
figure 5

Mean absolute zonal velocity below 200 m in EDJ_Low (multiplied by 1000, blue line) and in EDJ_Ref (green line) in m s−1

Full size image

Figure 6a shows a mode diagram for the zonal velocity shown in Fig. 1 and given by Eq. 1. The mode diagram shows that there is a shift of energy from modes 1–5 to modes 13–16 over the first 12 years, followed by a slower shift towards higher modes with energy reaching modes 18–20 by year 60. Thereafter, a somewhat similar process can be observed with intermittent bursts of high vertical mode energy propagating towards higher modes over a timescale of 5–10 years, albeit with weaker amplitude than the initial process. The same modal analysis applied to the low amplitude configuration EDJ_Low is shown in Fig. 6b. In this case, the energy level falls off after only 3 years, quite different from the fully nonlinear case. It should also be noted that the barotropic mode (mode 0) contains very little energy in both experiments and is of no interest here.

Fig. 6
figure 6

Mode diagram showing the amplitude of the projection of the zonal velocity, as given by Eq. 1, in m s−1 and as a function of time and vertical mode number at the same location as in Fig. 1 in EDJ_Ref (a) and EDJ_Low (b). Note the different time axis used in (b) compared to (a)

Full size image

Turning to Experiment EDJ_Const, Fig. 7a shows a Hovmoeller diagram of the zonal velocity similar to Fig. 1. The upward phase propagation at the beginning and the downward phase propagation later on is also visible here, as in EDJ_Ref, and evolves on a similar time scale to that seen in EDJ_Ref. A model run with uniform stratification has the advantage that the vertical structure functions, \(\hat {p}_{n}\), can be calculated analytically as

$$ \hat{p}_{n}= \sqrt{2}\cos{\frac{\pi n z}{H}} $$
(3)

with a corresponding gravity wave speed of \(c_{n}=\frac {N H}{n \pi }\) for n = 1,2,3,... (Gill 1982). Figure 8 shows the projection for the zonal velocity onto each vertical normal mode in this model run, showing a very similar behaviour to that seen in Fig. 6a, with energy transferring rapidly from modes 1–5 to mode 15 over the years 12–15 of the run.

Fig. 7
figure 7

Hovmoeller diagrams showing a zonal velocity in m s−1 as a function of time and depth at the same location as in Fig. 1 for EDJ_Const; b the analytic reconstruction of the zonal velocity using modes 15 and 16 and c the analytic reconstruction using modes 14, 15 and 16

Full size image
Fig. 8
figure 8

As Fig. 6 but for EDJ_Const

Full size image

We now show that it is possible to reconstruct the upward phase propagation at the beginning of the model run and the subsequent downward phase propagation later on using a linear superposition of basin modes. We noted from Fig. 3 that almost all the variability in EDJ_Ref is confined to basin modes and this is also true for the uniform stratification run being considered here (not shown). We now consider the linear superposition of two basin modes U 1 and U 2, given by

$$ U=U_{1} + U_{2}= \sin{(\omega_{1} t)}\cos{(n_{1} z)}+ \sin{(\omega_{2} t)} \cos{(n_{2} z)} $$
(4)

where \(z=\frac {\pi z^{*}}{H}\) (z ∗ being the vertical coordinate measured positive upwards with z ∗ = 0 at the surface), n is the vertical mode number, \(\omega =\frac {2 \pi c_{n}}{4L}\) is the corresponding basin mode frequency, and the amplitude is set to one for both modes for convenience (note that the choice of amplitude is not important, only that the amplitude of both modes be the same). This combination of two basin modes can be reformulated using standard trigonometric identities to give

$$\begin{array}{@{}rcl@{}} U &=& 2\cos(n_{1} z)\sin\left[\frac{(\omega_{1}+\omega_{2}) t}{2}\right]\cos\left[\frac{(\omega_{1}-\omega_{2}) t}{2}\right] \\ &&- 2\sin(\omega_{2} t)\sin\left[\frac{(n_{2}+n_{1}) z}{2}\right]\sin\left[\frac{(n_{2}-n_{1}) z}{2}\right] \end{array} $$
(5)

Figure 7b shows the resulting pattern of zonal velocity for the 15th and 16th vertical normal modes (that is, n 1 = 15 and n 2 = 16). The similarity to the evolution of the zonal velocity shown in Fig. 7a is clear, a similarity that is improved by adding the 14th mode but with only half the amplitude in Fig. 7c. From Eq. 5, it is clear that the low-frequency modulation in time comes from the \(\cos \left (\frac {(\omega _{1} - \omega _{2}) t}{2}\right )\) term with period \(T=\frac {4 \pi }{\omega _{1}-\omega _{2}}\) of about 25 years. This can be clearly seen near the surface and bottom where the second term is small because of the sine dependence on z. It follows that the time evolution of the zonal velocity shown in Fig. 7a can be largely reconstructed using just a few basin modes.

An interesting aspect of the reconstruction is that the phase of the basin modes used for the reconstruction is set by putting u = 0 at the start of the model integration. It therefore appears that these basin modes are excited at the start of the model integration but that energy is transferred into these basin modes during the first 10 years, or so, of the model run, as indicated by the mode diagram shown in Fig. 8. During the time energy transfers into these modes, they can be regarded as “pre-existing” basin modes, excited at the start of the model integration. This idea of pre-existing modes that get amplified through nonlinear interactions is supported by comparing the time evolution of the flow in EDJ_Ref and EDJ_Low in Fig. 4a, b. Note, in particular, the appearance of downward phase propagation in the low amplitude model run after year 10, similar to that in the reference run except for the enhanced amplitude in the reference run.

It seems surprising that the model has memory of its initial condition at least out to year 30. It certainly retains its memory well beyond the time (10–15 years based on Fig. 4b) after which the low amplitude spin-up has faded away and even while the EDJs in the model are being sustained by nonlinearity in EDJ_Ref. It is also clear that the reconstruction captures the upward phase propagation (corresponding to downward energy propagation for linear waves) in the initial phase of the spin-up and the later downward phase propagation that is characteristic of the EDJs themselves. The reconstruction is even able to reproduce the exact time (near 25 years) where the upward phase propagation changes to downward phase propagation. Note that the introduction of the 14th mode (n = 14) does not change this transition time by much, but narrows the wave beam, which in turn leads to a more realistic reconstruction.

We have also reconstructed the low amplitude solution shown in Fig. 4b. Figure 9 shows the mean amplitude of the projection of the zonal velocity onto each mode, averaged over the first 10 years of the model integration. The maximum amplitude is in mode 10 and the reconstruction uses all modes from 8 to 12 with the amplitude taken from the figure (again with u = 0 at t = 0). The solution can be seen in Fig. 10. Obviously, the reconstruction does not reproduce the reduction in amplitude seen in Fig. 4b but it does capture the transition time for the change from predominantly upward phase propagation to downward phase propagation. We note too that even the low amplitude case exhibits downward phase propagation after year 10, and that Fig. 10 also serves as a reconstruction of the spin-up in the nonlinear experiment, EDJ_Ref. This is not surprising given the energy present around mode 10 in the first 20 years of EDJ_Ref (see Fig. 6a). However, the presence of energy around mode 15 between years 10 and 60 has the effect of extending the period of downward phase propagation during the spin-up compared to the low amplitude case, as well as reducing the vertical scale of the EDJs.

Fig. 9
figure 9

Amplitude of the projection of the zonal velocity onto each mode averaged over the first 10 years of EDJ_Low at the same location as used to produce Fig. 6b. The modes between the green and red lines are used for the reconstruction

Full size image
Fig. 10
figure 10

Hovmoeller diagram showing the analytic reconstruction of the zonal velocity in m s−1 for EDJ_Low using mode 8 to 12 (as seen in Fig. 9)

Full size image

Comparing the reconstruction shown in Fig. 10 with EDJ_Ref at later times (see Fig. 4a), it is clear that whereas the reconstruction shows alternating phases of upward and downward phase propagation, EDJ_Ref is dominated by downward phase propagation. The other feature of EDJ_Ref that is not captured by the reconstruction is the modulation of the episodes of downward propagation, in particular the reduction in energy around mode 15 that can be seen in Fig. 6a between years 60 and 80 and its reappearance again after year 80. Indeed, the impression from Fig. 1 is that the cycle that began around year 20 is in the process of repeating itself, starting around year 80, and perhaps also again after year 150.

To test the dependence of the solution on the initial state, we have also run the same model set-up as in EDJ_Ref but with the initial temperature and salinity fields replaced by the average over the last 30 years of EDJ_Ref. As can be seen from Fig. 11, and comparing with Fig. 4a, the EDJs in this case are a lot weaker than in EDJ_Ref with only very weak evidence for downward phase propagation, this being confined to the upper 1000 m or so. It follows that the initial state matters and has a clear influence on the whole of the subsequent model evolution (at least for the 50 years shown here).

Fig. 11
figure 11

Hovmoeller diagram showing zonal velocity in m s−1 as a function of time and depth at the same location on the equator as in Fig. 4a but for EDJ_Mean

Full size image

Finally in this section, we have repeated the first experiment EDJ_Ref but this time with the Pacanowski and Philander (1981) mixing scheme replaced by the KPP scheme (Large et al. 1994). The Hovmoeller diagram (Fig. 12a) for the zonal velocity at the same location as used for Fig. 1, looks quite similar to that shown in Fig. 1, but there are some differences in detail. These are most clearly revealed by looking at the projection of the zonal velocity onto the vertical modes, following Eq. 1, shown in Fig 12b. Comparing with Fig. 6a, it is clear that following the initial spin-up, the energy is more obviously being transferred to higher modes than in experiment EDJ_Ref, reaching even mode 25, and that this process tends to repeat itself in time, corresponding to the decadal modulation of the EDJs that can be seen in Fig. 12a. The reason for this difference in behaviour is not clear but obviously has to do with the use of a different vertical mixing scheme. One possibility is that the higher modes are more damped when using the Pacanowski and Philander (1981) scheme, preventing the energy from spreading to these modes in experiment EDJ_Ref.

Fig. 12
figure 12

a Hovmoeller diagram showing zonal velocity in m s−1 as a function of time and depth at the same location as in Fig. 1 but using the KPP mixing scheme. b As Fig. 6a except for the experiment using the KPP mixing scheme

Full size image

4 Discussion

Let us begin by returning to the uniform stratification case, so that the vertical modes can be written as cosine functions. We first note that if two neighbouring vertical modes are oscillating in time at the same mean frequency (not necessarily the basin mode frequency), ω, and the same amplitude, then whether this results in the dominance of upward or downard phase propagation depends on the phase difference between the two oscillations in time. To see this, we note that

$$\begin{array}{@{}rcl@{}} U &\,=\,& \sin(\omega t)\cos(mz) + \sin(\omega t+\theta)\cos(nz) \end{array} $$
(6)
$$\begin{array}{@{}rcl@{}} &\,=\,& 2\sin\left[\omega t+\frac{\theta}{2}+\frac{(n+m)}{2}z\right]\cos\left[\frac{(m-n)}{2}z+\frac{\theta}{2}\right] \\ &&\!+ 2\sin\left[\omega t\,+\,\frac{\theta}{2}-\frac{(n\,+\,m)}{2}z\right]\cos\left[\frac{(n\,-\,m)}{2}z+\frac{\theta}{2}\right], \end{array} $$
(7)

where z is again given by \(z=\frac {\pi z^{*}}{H}\), as can be shown using standard trigonometric identities. Note that the amplitude of the downward propagating wave is given by \(\cos \left (\frac {m-n}{2}z+\frac {\theta }{2}\right )\), that of the upward propagating wave by \(\cos \left (\frac {n-m}{2}z+\frac {\theta }{2}\right )\), and the relative importance of upward and downward phase propagation then depends on the phase difference 𝜃.

We next note that if the two neighbouring modes are not oscillating with the same frequency, i.e.

$$ U=U_{1} + U_{2}= \sin{(\omega_{1} t)}\cos{(n_{1} z)}+ \sin{(\omega_{2} t)} \cos{(n_{2} z)} $$
(8)

then the effect is the same as having a time-dependent phase difference between two modes oscillating at the same frequency, as shown by

$$\begin{array}{@{}rcl@{}} &&{\kern5pt}\sin(\omega_{1}t)\cos(n_{1} z) + \sin(\omega_{2}t)\cos(n_{2} z)\\&&{\kern10pt} = \sin(\omega_{1}t)\cos({n_{1}} z) + \sin(\omega_{1}t+\Phi(t))\cos({n_{2}} z) \end{array} $$
(9)

where Φ(t) = (ω 2 − ω 1)t. It follows that when two neighbouring basin modes interact, as in Eq. 8, there must result alternating phases of upward and downward phase propagation. Consequently, if the energy always stays constant in the modes, there will be a cyclic stable system were phases of upward and downward phase propagation are following on one another. This cyclic stable system does not exist in EDJ_Ref, where downward phase propagation dominates after the initial spin-up (see Fig. 1). This implies that the energy cannot stay constant in each mode in EDJ_Ref, consistent with Fig. 6a.

Returning to the low amplitude experiment, EDJ_Low, we have seen that the amplitude of the solution decreases rapidly after year 10 (see Fig. 4b), which is not seen in EDJ_Ref (Fig. 1). This means that there is an energy enhancement mechanism that comes from the nonlinearities that sustain the deep jets in EDJ_Ref. Such a mechanism has been discussed by Ascani et al. (2015) and it involves transfer of energy from intraseasonal Yanai waves to the EDJs. The intraseasonal Yanai waves are, in turn, shed by tropical instability waves that develop near the surface of the model ocean as a consequence of the imposed wind forcing (see Ascani et al. 2015 for the details). Furthermore, as we noted earlier, since the reconstructions have their phase set at the beginning of the model run, it is as if the effect of the nonlinearity is to transfer energy into pre-existing basin modes that are generated as part of the spin-up (note that basin modes are excited even as part of a linear spin-up following the sudden application of the wind stress to a resting ocean, as in our model experiments). As we have seen, the dominance of downward phase propagation in EDJ_Ref after the initial spin-up means that the energy does not remain constant in the modes. Two questions arise. First, what is the fate of the energy that disappears around year 60 in Fig. 6a and secondly what happens when the energy reappears around mode 10, as happens in Fig. 6a around year 90 (year 70 in Fig. 12a). An analysis by Aiki et al. (2017) showed that for a single baroclinic mode, the main region for dissipation is the extratropics where extratropical Rossby waves are being radiated from the eastern boundary as part of the basin mode solution, a factor that may well play a role here. Regarding the second question, one possibility is that in addition to the enhancement mechanism for pre-existing basin modes, there is also a generation mechanism at work that re-establishes the energy in the higher modes. A candidate for the later generation mechanism is the mechanism of Hua et al. (2008), in which the vertical scale of the EDJs is set either by the destabilization of intraseasonal Yanai waves or by resonant triad interactions involving intraseasonal Yanai waves. All these topics require further investigation and are topics for future research.

5 Summary and Conclusions

Results from several idealised model experiments have been described. The experiments show the development of vertically alternating zonal jets along the equator that are analogous to the equatorial deep jets (EDJs) that have been observed, for example in the equatorial Atlantic (Brandt et al. 2011; Claus et al. 2016). We noted that the energy in the model runs associated with the zonal velocity along the equator is mostly confined to frequencies associated with resonant basin modes (Cane and Moore 1981; Yang and Liu 2003), a feature of the solution we have been able to exploit in order to describe the emergence of the EDJs in the model solutions. In particular, the initial spin-up, associated with first upward phase propagation and later downward phase propagation, can be understood as the superposition of two or three neighbouring basin modes that have their phase set at the start of the model run. It is as if, during the first phase of downward phase propagation, the nonlinearity in the model has the effect of putting energy into these pre-existing basin modes, the mechanism for which has been described by Ascani et al. (2015). It should be noted that the EDJs analysed by Ascani et al. (2015), and also Matthießen et al. (2015), belong to this first period of downward phase propagation arising from the spin-up, although there is no reason to suppose that the mechanism supporting the EDJs during this time is different from that applying later. We have also seen that the model solution is sensitive to the choice of initial conditions; in this model set-up, it does seem to be important to have basin modes excited at the start of the model run that can be subsequently amplified.

We noted that a prediction of the reconstructions we have presented is the appearance of alternating phases of upward and downward phase propagation in the future. The fact this does not occur (downward phase propagation dominates) means that the energy cannot remain constant in each vertical mode. Instead, as we have seen, the energy in the active vertical modes (typically between modes 15 and 20) fades out after the initial spin-up phase and then shows a tendency to reappear at later times. Associated with this decline and re-emergence, the EDJs themselves undergo modulation on decadal time scales. One possibility is that the re-emergence is associated with the destabilisation mechanism of Hua et al. (2008), although this remains speculation at this time.

The preference for downward phase propagation (implying upward energy propagation) is intriguing. We noted in the introduction that in the model run of D’Orgeville et al. (2007), energy tends to be taken from depths where the imposed intraseasonal wave has the largest amplitude (implying the largest energy input) to depths where the amplitude of the imposed intraseasonal wave is zero (implying no input), suggesting that the direction of energy propagation is determined by the vertical distribution of the energy sources and sinks, with energy, generally, propagating towards the sinks. In the experiments reported here, the energy propagation is upward (see Fig. 4 in Matthießen et al. 2015) and it seems likely that interaction between the deep jets and the Equatorial Undercurrent in the model provides a significant energy sink.

It should be noted that the energy source for the deep jets comes from intraseasonal waves that are shed by tropical instability waves in the surface layers (see Ascani et al. 2015) and that the forcing acts over a considerable depth range (see Claus et al. 2016 for a diagnosis of the forcing from observations and note that the forcing acts over the whole depth range occupied by the jets, something we have also diagnosed from the model and will be discussed in detail elsewhere). Note, in particular, that Ascani et al. (2015) identified the nonlinear interaction of two intraseasonal Yanai waves as being important, one with relatively large vertical scale and the other with the vertical scale of the EDJs themselves. The existence of the latter suggests that the intraseasonal Yanai waves and the EDJs are interacting with each other and it is possible that this mechanism at least partly explains why it is pre-existing basin modes that get energised during the spin-up phase. It should also be noted that Muench and Kunze (1999, 2000) have suggested that internal waves play a role in maintaining the EDJs, another issue requiring further investigation.

Finally, we note that the analysis carried out by Claus et al. (2016) identified a single frequency as being sufficient to describe the EDJs as observed on the equator at 23∘ W. We noted in Section 4, that two neighbouring modes oscillating at a single frequency can lead to either upward or downward phase propagation, depending on their phase difference, and that the analysis of Claus et al. (2016) effectively sets this phase difference (the observed EDJs have downward phase propagation). It is also the case that the available time series for the EDJs in observations is still short and there is no evidence yet of any low-frequency modulation of the EDJs in observations. Much longer time series (several decades at least) of observations are required to identify such modulation. Likewise, identifying such modulation is undoubtedly a requirement if more than one dominant frequency for the EDJs is to be extracted from the observations.