## Abstract

We present a simplified theory using reduced-gravity equations for North Atlantic Deep Water (NADW) and its variation driven by high-latitude deep-water formation. The theory approximates layer thickness on the eastern boundary with domain-averaged layer thickness and, in tandem with a mass conservation argument, retains fundamental physics for cross-equatorial flows on interannual and longer forcing time scales. Layer thickness anomalies are driven by a time-dependent northern boundary condition that imposes a southward volume flux representative of a variable source of NADW and damped by diapycnal mixing throughout the basin. Moreover, an outflowing southern boundary condition imposes a southward volume flux that generally differs from the volume flux at the northern boundary, giving rise to temporal storage of NADW within the Atlantic basin. Closed form analytic solutions for the amplitude and phase are provided when the variable source of NADW is sinusoidal. We provide a nondimensional analysis that demonstrates that solution behavior is primarily controlled by two parameters that characterize the meridional extent of the southern basin and the width of the basin relative to the equatorial deformation radius. Similar scaling applied to the time-lagged equations of Johnson and Marshall provides a clear connection to their results. Numerical simulations of reduced-gravity equations agree with analytic predictions in linear, turbulent, and diabatic regimes. The theory introduces a simple analytic framework for studying idealized buoyancy- and wind-driven cross-equatorial flows on interannual and longer time scales.

## 1. Introduction

One of the most important functions of the oceans in the global climate system is the meridional transport of heat and freshwater (Buckley and Marshall 2016). These transports are required to connect regions of net heat or freshwater gain with regions of net heat or freshwater loss, which are often located at great distances. One clear and climatically important example of this is the surface heat flux in the Atlantic Ocean (Talley 2003). The net heat loss in the subpolar North Atlantic and Nordic Seas is balanced by a net northward heat flux across the equator of approximately W (Trenberth and Caron 2001). This northward flux represents a significant fraction of the total northward heat flux in the coupled ocean–atmosphere system. In addition to heat, the ocean also transports freshwater and climatologically important passive tracers, such as and oxygen, into the Northern Hemisphere (Kawase and Sarmiento 1986).

Much of this northward heat transport at low latitudes is carried by the meridional overturning circulation (MOC; Hall and Bryden 1982). Northward flow of warm, salty water is balanced by southward flow of colder, fresher water at middepths (e.g., Schmitz 1995). These middepth waters are isolated from interaction with the atmosphere until they outcrop in the Southern Ocean, feeding upwelling-driven by winds over the Antarctic Circumpolar Current (Toggweiler and Samuels 1998; Wolfe and Cessi 2011). In the limit of weak diapycnal mixing in the ocean interior, the source waters formed by deep convection at the high latitudes of the North Atlantic are transported, nearly adiabatically, to these outcrops in the Southern Ocean. Strong diapycnal mixing has been measured near and above rough topography (Polzin et al. 1997), and there is of course strong mixing and water mass modification in the surface mixed layer, but mixing at the middepths of North Atlantic Deep Water (NADW) is generally much weaker.

The MOC has received much attention, both as a feature of the ocean circulation and as a component of the coupled climate system. Here, we focus on the dynamics of the middepth southward flow connecting the high-latitude North Atlantic with the high-latitude Southern Ocean. The basic mechanisms by which the lower limb of the Atlantic meridional overturning circulation (AMOC) is established have been explored in several previous studies. Kawase (1987) described the adjustment process of the deep ocean to a high-latitude mass source as a combination of coastal and equatorial Kelvin waves rapidly communicating the pressure signal to the eastern boundary, where it is then slowly propagated into the basin interior by baroclinic Rossby waves. In his model and theory, these Rossby waves were damped by a diapycnal mass flux. In a pair of papers, Johnson and Marshall (Johnson and Marshall 2002a,b) develop a theory for the time-dependent MOC response to variable water mass formation rates or Ekman transport driven at high latitudes. Their principal finding, as it relates to the present study, is that the equator acts to dampen the meridional propagation of transport anomalies in the hemisphere opposite the forcing hemisphere. They found that the degree of damping depends strongly on the forcing frequency. For parameters typical of the North Atlantic, variability in deep-water formation at periods less than roughly decadal remain largely within the North Atlantic. Analogously, variability in Ekman transport in the Southern Ocean at similar frequencies does not penetrate into the North Atlantic basin. This basic mechanism was applied to multibasin, global domains by Johnson and Marshall (2004) and to heat content and sea level change by Zhai et al. (2011).

The present study builds on these previous results by considering the influences of eddies and diapycnal mixing using a simplified formulation that yields closed-form analytic solutions for the spatially averaged layer thickness . The theory developed below details a new perspective on the propagation of NADW anomalies that is consistent with, and a parallel and simplified description of, previous theoretical work by Johnson and Marshall (2002b). Nondimensionalization is used to identify two controlling parameters and to provide a simple interpretation of the system behavior. We show that the same two parameters also appear in nondimensional forms of the Johnson and Marshall (2002a) governing equations for the eastern boundary layer thickness , providing a connection between the two approaches. The simplified formulation is allowed by exploiting the approximation which, in tandem with a mass conservation argument, retains the fundamental physics involved in cross-equatorial flows on interannual and longer time scales.

## 2. Governing equations

The 2.5-layer reduced-gravity equations are used to develop the theory presented below. The upper and bottom layers are assumed motionless, and the active middle layer represents the NADW responsible for transporting flow southward as part of the lower limb of the AMOC. While this is a significant idealization of the three-dimensional AMOC, we have performed simulations with a northward-flowing upper layer (with and without wind stresses) and a southward-flowing bottom layer. In this case, each layer yields similar behavior for layer thickness and volume transports as a single active middepth layer. This motivates a detailed study of the simplified configuration illustrated in Fig. 1, in which only the middle layer is in motion.

Following a presentation of the reduced-gravity framework, a description of the theory for the domain-averaged layer thickness is given, an outline of analytic solutions is provided and their salient features discussed.

### a. Reduced-gravity equations

The reduced-gravity equations are placed on the equatorial beta plane in a simple rectangular domain free from wind stress and topography. In dimensional form, the 2.5-layer reduced gravity equations with motionless upper and bottom layers simplify to

Subscripts denoting layers are dropped with the understanding that the dynamics being considered are those belonging to layer 2 in Fig. 1. The horizontal coordinates are along zonal and meridional directions, respectively, and is the unit vector normal in the vertical. The horizontal extent of the domain is taken to be , with the equator at . Horizontal velocity is denoted as , and layer thickness is denoted as *h*. The reduced gravity is , is the gradient of planetary vorticity, and is the constant eddy viscosity coefficient.

If the flow is nonadiabatic, then the right-hand side of (2) is included where is the imposed spatially uniform thickness relaxation time scale to the prescribed target layer thickness . Previous studies have implemented diapycnal mixing using the functional form , where is the diapycnal diffusivity, which acts to increase layer thickness (Allison et al. 2011). However, the form of diapycnal mixing used by Kawase (1987) and used in this study is approximately the same as the parameterization for appropriate choices of *γ* and . In fact,

and the appropriate choices for *γ* and are and . This is valid for ; therefore, may be chosen to be any thickness such that , where is the equilibrium thickness. In the analysis given below we provide analytic predictions for , allowing diapycnal mixing in the form of a relaxation as a substitution for the form involving the diffusivity , that is,

### b. Domain-averaged continuity equation

This section presents the theoretical model for layer thickness of NADW in the form of a first-order, nonlinear, ordinary differential equation. Spatially averaging (2), and noting that zonal velocities vanish at western and eastern boundaries, gives

where and are the distances from the equator to the northern and southern boundaries. We have introduced the averaging operator

where *F* is some field variable, and *A* is the domain area. Equation (4) relates the time rate of change of domain-averaged layer thickness to the zonally integrated volume transports at and and the diapycnal mass flux. The volume transports at these latitudes are defined to be

Since we aim to model NADW, flow at these boundaries is such that , that is, flow is southward and is driven by remote downwelling in the north and upwelling in the south. In particular, upwelling and downwelling processes are presumed to exist in regions beyond the meridional extent of the domain. In addition, volume fluxes at the northern and southern boundaries are imposed in zonal geostrophic balance where inertial accelerations are omitted and zonal velocities are set to zero, implying that meridional gradients of layer thickness vanish.

Volume flux into the domain through the northern boundary at is imposed and provides the forcing method by which the flow is put into motion. This imposed forcing is related to, and is an idealized simplification of, downwelling at high latitudes. This applied volumetric flux is prescribed as a time-dependent condition that varies sinusoidally with amplitude , frequency , and mean , that is,

Note that is the forcing period over which fluctuations of the AMOC occur. The case when and will also be considered and, as our results will show, the system response in the limit will tend toward the response observed when . For the southern boundary, flow is fluxed southward from the domain in zonal geostrophic balance

where is the planetary vorticity at the southern boundary.

We now make note of a key assumption that greatly simplifies the equation for and allows for closed form analytic solutions. This assumption takes the layer thickness on the eastern boundary as representative of the basin-averaged layer thickness, . To justify this and to show when the relation may be valid, we assume further that mixing is weak and the ocean interior is in geostrophic balance. Therefore, if no normal flow boundary conditions are imposed, then it is expected that layer thickness will be *y* independent along the eastern boundary, that is, . Since we are concerned with interannual and longer time scales, we assume that the combination of coastal and Kelvin waves are fast compared to the forcing period . In this sense, pressure signals originating from the high-latitude forcing are *immediately* propagated south along the western boundary, east along the equatorial waveguide, and finally poleward along the eastern boundary (Kawase 1987; Johnson and Marshall 2002b). Marshall and Johnson (2013) distinguish between Kelvin and Rossby waves as the responsible coastal mechanism for the adjustment of the MOC. These boundary waves have speeds given by , where is the gravity wave speed, is the deformation radius, and is the Munk layer thickness. For the numerical simulations provided below, when is the equatorial deformation radius and when is deformation radius at the highest latitudes in the numerical domain. Therefore, pressure signals propagate quickly about the basin. If we expect the approximation to hold, then sufficient time must be allowed for the interior thickness to respond, via westward-propagating Rossby waves, to the pressure signal along the eastern boundary so that is indeed representative of . Otherwise, layer thickness in the interior may differ significantly from , and there is no a priori justification for the relation .

To quantify the practicality of the relation , consider, for , the layer thickness differential

where is the layer thickness evaluated at a distance equal to the width of the western boundary layer from the western boundary. If can be shown to be small relative to the signal of interest, then the relation may be justified. To this end, we note that for a geostrophic interior the time rate of change of layer thickness may be related to zonal gradients through

where is the Rossby wave speed (Johnson and Marshall 2002b). Suppose *η* represents the amplitude for variations of , *T* is the time scale for variability of , and , then we may use (11) to write (10) as

Therefore, ≪ 1 if ≪ 1, and the assumption that basin-averaged layer thickness is representative of the layer thickness on the eastern boundary is valid in the limit that the time scale *T* for variability of is long compared to the time scale for a Rossby wave to propagate across the basin. For the forcing given by (8) and in the limit of rapidly propagating Kelvin waves, this condition requires the forcing period to be long compared to the Rossby wave basin-crossing time, that is, ≪ 1. However, note that ≫ 1 is not necessarily true for high-frequency forcing since the amplitude for the variation of , as determined by reduced-gravity dynamics, may be small.

The condition ≪ 1 is generally satisfied at low latitudes for interannual or longer period forcing because waves propagate quickly near the equator. For a 6000-km-wide basin with a typical gravity wave speed of 2 m s^{−1}, the basin-crossing time for a mode-1 baroclinic equatorial Rossby wave is about 100 days and is fast compared to the time scales of interest. From mid- to high latitudes Rossby waves are much slower, and so this assumption becomes more tenuous. To a rough approximation, the relative error for is given by , where is a nondimensional measure of the basin area with long crossing time, is the latitude at which , and *H* is the spatial and temporal mean layer thickness. Even in cases where the high-latitude crossing time is comparable to , that is, when , the condition ≪ 1 will be satisfied over most of the basin and ≪ 1. Moreover, as will be shown below, the amplitude *η* of variability in layer thickness is small for forcing periods less than the high-latitude basin crossing time scale; thus, ≪ 1 and is valid in parameter space for which the forced response is large. In fact, at high-frequency forcing, both full numerical solutions and the analytic solutions below predict a small-amplitude response; therefore, the errors incurred by assuming are small and are not fundamental to the overall predictions of the theory. Favorable comparisons are also made with the more complete solutions of Johnson and Marshall (2002a) in section 4.

Zonally integrating the product in (9) yields the following for the volume flux at the southern boundary

Here, the approximation has been made, and is a prescribed layer thickness at the southwestern corner of the domain, . Information propagating equatorward along the western boundary and originating from outside our model domain is responsible for setting , and therefore, it is appropriate to impose this as a boundary condition. Using the expressions for and given in (13) and (8), the domain-averaged continuity equation [(4)] becomes

This first-order differential equation becomes our governing equation for layer thickness of NADW. While we will be particularly focused on the case in which *γ* is sufficiently large such that the flow may be considered adiabatic, we will also consider the case for nonadiabatic flow. In each case, closed-form analytic solutions are obtained.

### c. Analytic solutions for

We briefly outline the analytic solutions that solve (14) for four different cases of interest, which we split into two categories: adiabatic and nonadiabatic flow. For each category we consider two possibilities for the volume flux : either is constant in time [with in (8)] or varies sinusoidally (with ). See the appendix for details.

#### 1) Adiabatic flow

In the case of adiabatic flow in the presence of constant volume flux through the northern boundary, with , solutions for domain-averaged layer thickness take the form

The equilibrium layer thickness is

and is approached over the intrinsic time scale

and is a constant determined by the initial condition . The time scale *τ* may be interpreted as the equilibrium volume of NADW divided by the transport into the domain from the north. This time scale was also found by Johnson and Marshall (2002b) for the equilibration of layer thickness on the eastern boundary.

If is allowed to vary sinusoidally, then solutions resemble the case when is constant; however, they contain sinusoidal variations that are phase-shifted from the imposed perturbation by the amount , so that variations are approximately in phase with for . In addition, the amplitude of variation is given by

These solutions take the form

where and the expressions for *H*, *τ*, and are the same as those given for the case of constant volume flux. The solution consists of the exponential transition to the equilibrium solution (first term), a phase-shifted sinusoidal oscillation (second term), and a transient associated with the initial condition (third term).

#### 2) Nonadiabatic flow

Solutions for above are valid in the limit ≫ ; however, solutions for take a similar form when nonadiabatic effects are considered, that is, when . If the volume flux varies sinusoidally, then nonadiabatic solutions become

where the solution form is similar to (19). Parameter may be interpreted similarly to *H* in (16): the equilibrium thickness obtained as . On the other hand, may be interpreted as an equilibrium thickness obtained if time is allowed to march backward, that is, . The distinct pair results from the quadratic form of (14) and corresponds to layer thicknesses that satisfy . Expressions for equilibrium layer thicknesses , perturbation thickness , and time scale are

### d. Analytic solutions for

Once the analytic form for is known, this can be used to determine the closed form solutions for by substitution of into (13). We note that this solution is fully nonlinear, which, in the adiabatic limit ≫ and for ≫ , is of the form

where we have taken . The terms omitted in taking the long-time limit ≫ simply describe the transient adjustment process as the system equilibrates from the initial condition, a process illustrated in the results provided in section 3a.

The first term for is the imposed mean volume transport while the second gives fluctuations about the mean with identical phase to that of ; however, using (18), the amplitude is seen to differ from the imposed amplitude at the northern boundary, through (8), by the scalar factor . Therefore, the amplitude of this fluctuation is easily seen to decay to zero for high-frequency forcing, that is, when ≫ 1, and obtains its maximum value in the limit of low-frequency forcing, that is, when ≪ 1. This analytic result is similar to the time-lagged expression for volume transport found by Johnson and Marshall (2002a). This acts as a low-pass filter to variability propagating out of the source hemisphere, a mechanism termed the “equatorial buffer” by Johnson and Marshall (2002b).

The remaining terms result from nonlinear projections of the imposed frequency onto the mean zero-frequency mode and the higher-frequency mode with phase . Each of these projections have amplitude that is smaller than the lower-order fluctuations [the second term in (24)] by the amount . The projection onto the zero-frequency mode acts to increase mean southward transport at the southern boundary when the forcing is sinusoidal, an effect that increases with increasing , or equivalently, increasing . Instances when , for , transport at the southern boundary is equal to the imposed mean transport, that is, . Instances when , for *n* odd, the magnitude of obtains its largest [] and smallest [] values.

## 3. Numerical experiments

We briefly outline the setup for the numerical simulations of the reduced-gravity equations [(1) and (2)] and the parameters used for their calculation. The numerical domain is an idealized rectangular basin with dimensions (approximately wide extending from S to N) and flat topography. The domain is centered about the equator such that . The reduced-gravity equations are discretized and solved on a C-grid using finite differences on a uniform grid with spatial resolution km. The numerical time-stepping scheme used is the explicit third-order Adams–Bashforth scheme with time step s (when is large) or s (when is small). At eastern and western boundaries, no-slip and no normal flow conditions are applied.

Numerical simulations are performed with two layers in anticipation of high Reynolds number calculations in which baroclinic instability may become important. However, first we consider flow in the laminar regime () where the upper layer is essentially passive. In this study we fix Sv (1 Sv = ), Sv, and with forcing period years.

Within 400 km of the northern boundary the meridional velocity is restored toward a hyperbolic secant profile that is geostrophically balanced. The layer thickness assumes a hyperbolic tangent decay from the modeled layer thickness at . At the imposed velocity is , where is the horizontal scale of the western boundary layer, and the amplitude is determined such that the desired southward volume flux through the northern boundary is obtained. Meridional velocity at the southern boundary is restored in a similar fashion. The profile for layer thickness assumes a hyperbolic tangent decay from the modeled layer thickness at on the southern boundary and zero at the western boundary (). Specifically, the meridional velocity is restored toward . This boundary condition allows all the transport that is approaching the southern boundary to exit the domain along the western boundary. We have found the results to be insensitive to the details of boundary layer width or restoring time scale. We set so the transport at the southern boundary is initially small and allowed to adjust from zero to a signal with mean transport *S*; however, the results are not qualitatively sensitive to this choice. The following results illustrate the validity of the analytic solutions for outlined in section 2c.

### a. Results

#### 1) Adiabatic flow

##### (i) Steady forcing

When the flow is adiabatic and the northern volume flux is steady, the domain-averaged layer thickness is given by (15). Figure 2a shows numerical solutions (solid curves) and analytic solutions (dotted curves) for . We note an incremental disagreement between these curves as decreases, corresponding to weakly stratified flows. This disagreement is due to a decreasing deformation radius for which, particularly at high latitudes, . Additional simulations (not shown) confirm that unresolved deformation radii contribute to the disagreement seen in Fig. 2a. Despite a declining resolution with decreasing , analytic solutions for provided in section 2c agree with numerical simulations of the reduced-gravity equations.

Figure 2b shows volume transports at the southern boundary and the imposed constant volume transport at the northern boundary (shown in gray). While we set , values seen in Fig. 2b slightly overestimate this due to numerical quadrature and interpolation used to diagnose transport. Moreover, southern transports equilibrate to the northern volume flux, balancing the incoming mass flux. Analytic and numerical solutions to are also seen to coincide (not shown). As decreases the volume flux at the southern boundary requires more time for equilibration. This trend, which goes like , can be related to the time required for Rossby waves at the southern boundary to travel across the basin from the eastern boundary to the western boundary. Since Rossby waves propagate with speed at midlatitudes and at speeds less than at the equator, weakly stratified adiabatic flows require longer periods than strongly stratified flows to equilibrate.

From Fig. 2b the flow displays evidence of early time instabilities, discernible as oscillations in the time series of for years. To aid in quantifying the degree of turbulence achieved by the flow, we define the Reynolds number

where, for these numerical solutions, . For years the layer thickness is small and . Eventually, these early time instabilities vanish as the layer thickness increases and the Re decreases below a critical value (Springer and Kawase 1993), where viscous effects dominate inertial accelerations.

The appearance of planetary vorticity and area *A* in the expressions for *H* and *τ* indicates that basin geometry plays a role in setting these quantities. Defining the area as

increases in and can be shown to impede equilibration. For example, , whereas *τ* scales linearly with and . That basin geometry can control equilibration is an indication of equatorial dampening effects found by Johnson and Marshall (2002a,b). We discuss these effects further in section 4 from a nondimensional perspective that illuminates the role of basin geometry.

##### (ii) Periodic forcing

If has the sinusoidal form given in (8), then domain-averaged layer thickness evolves according to (19). Similar to Fig. 2a, Fig. 3a provides a comparison of numerical solutions (solid curves) and analytic solutions (dotted curves), but Fig. 3a illustrates the effect of time-dependent volume flux at the northern boundary. Again, marginally resolved deformation radii at high latitudes in the numerical solutions are responsible for minor disagreement when the flow is weakly stratified. The analytic solutions in section 2c agree well with numerical solutions of the reduced-gravity equations.

Figure 3b shows the effect of sinusoidal volume transport at the northern boundary (shown in gray) on the transport at the southern boundary . Again, southern transports equilibrate to the northern volume flux (on average) and analytic and numerical solutions to are seen to coincide; however, high-frequency oscillations due to turbulence are not captured (not shown). The same leading-order features observed when persist when varies sinusoidally. Notably, the amplitude of and differ, indicating that southward meridional transport in the western boundary layer is lost to/gained from the basin interior. We return to the issue of loss of transport in section 4, where we connect this to basin geometry and propagation period of Rossby waves by extending the analysis of Johnson and Marshall (2002a,b).

The amplitude of variations for are given by the expression for in (18). For ≫ , the amplitude of variations scale like , so that high-frequency forcing has only a weak influence on domain-averaged thickness. When is small tends toward the maximum value . Equation (18) shows how the geometry of the domain also determines the magnitude of . For example, , but when and when . This reduced amplitude of for increasing reflects the equatorial damping discussed by Johnson and Marshall (2002a); however, the dependence on suggests a further dependence on the basin geometry. We return to this point, and the connection between our model and that of Johnson and Marshall (2002a), in section 4.

### b. Nonadiabatic flow

While the focus of this study is on adiabatic flow, for completeness, we make a brief mention of the results when diapycnal mixing is a leading-order effect. Specifically, whereas previous adiabatic considerations presumed ≫ , we now consider . This regime of nonadiabatic effects is comparable to flows studied by Kawase (1987), where Rossby waves are damped by diapycnal mass fluxes. For these numerical simulations, we set the spatially uniform relaxation time scale years and the relaxation layer thickness m.

#### 1) Steady forcing

Figure 4a compares curves from numerical simulations of the reduced-gravity equations (solid curves) and analytic solutions (dotted curves). As before, good agreement between these curves is evident. From (21) it can be shown that for ≪ 1, recovering the adiabatic result, while for ≫ 1. An important parameter is the ratio of the damping time scale *γ* to the time scale for Rossby waves to propagate across the basin (Kawase 1987). For the values of in Fig. 4 this ratio spans the interval , and Rossby waves are damped by diapycnal mass fluxes. It is clear how the equilibrium time scale depends on .

Figure 4b provides curves for volume transport at the northern boundary (shown in gray) and at the southern boundary . If the domain-averaged layer thickness equilibrates to a value different from , then a diapycnal mass flux into or out of the layer is observed, the sign being determined by whether is greater or less than . For example, when , Fig. 2a shows that in the absence of diapycnal mixing. With spatially uniform diapycnal mixing, Fig. 4a shows the curve for with tends toward and Fig. 4b shows that , implying the loss of mass. Similarly, when , tends toward and differs from by the amount of diapycnal mass fluxes.

#### 2) Periodic forcing

When sinusoidal transport is imposed in the presence of diapycnal mass fluxes, much of the above behavior is observed. However, in this case , and this can allow to fluctuate about . This behavior can be seen in Fig. 5. When , the basin-averaged layer thickness is such that for , and the flow experiences a diapycnal mass loss. Volume transport via diapycnal fluxes for other values of oscillate in time between loss (when ) and gain (when ); however, time-mean diapycnal mass fluxes for simulations shown in Fig. 5 are precisely those observed for simulations shown in Fig. 4, where .

Since in this case, (18) may be written

where, similar to (21), the ratio indicates the significance of diapycnal mass fluxes. Here, *τ* denotes the intrinsic time scale for equilibration in the adiabatic limit given by (17). When ≪ 1, diapycnal fluxes are weak and (27) approaches its adiabatic value given by (18). In the limit ≫ 1, diapycnal mass fluxes become significant and cause to decrease monotonically to zero. Therefore, if is allowed to vary from small to large, basin-averaged layer thickness anomalies are maximal in the adiabatic limit and increasingly negligible as diapycnal mass fluxes act to dampen these anomalies when is large. Since volume flux at the southern boundary is a function of [see (13)], the amplitude of will decrease when ≫ 1, requiring a loss of meridional volume transport into the domain interior. The adiabatic analog of the issue of loss of meridional volume transport and amplitude of is discussed further in section 4.

### c. Turbulent flow

To this point, all comparisons between numerical solutions of the reduced-gravity equations and the analytic solutions of basin-averaged layer thickness have been made in a laminar regime for which viscous effects dominate inertial accelerations. Since Earth’s oceans are replete with nonlinear behavior, we consider turbulent flow and how well analytic solutions for compare in this case. Here, we set and consider . This value of lateral viscosity permits instabilities to develop and thus extends the resolved physics beyond the geostrophic approximation used in the theory. However, we note that submesoscale processes and higher-order vertical modes are still neglected, and the Reynolds number defined in (25) depends on the eddy viscosity coefficient , which is not well constrained theoretically or observationally.

Figure 6 shows an instance of the volume transport when . The Reynolds number is an order of magnitude larger than the critical value . Eddies dominate the flow along the western boundary south of the equator, as expected from stability theory (Edwards and Pedlosky 1998), similar to findings by Goes et al. (2009) and consistent with observations in the Brazil basin (Dengler et al. 2004). A comparison of analytic and numerical versions of yields a difference during the initial establishment of the circulation where appears to increase more rapidly than analytic solutions predict. This is shown by the dashed curve for in Fig. 2a. This difference is due to an artificial mass source in regions where the layer thickness is prohibited from becoming negative. Early in these turbulent calculations, where initially, eddying motion causes fluctuations in layer thickness that, if permitted, result in . As increases and instances of no longer occur, agreement between analytic and numerical values of is observed.

A stability analysis by Isachsen et al. (2007) reports that baroclinic Rossby waves are unstable to barotropic perturbations resulting in increased westward phase velocities. A decomposition of layer thickness and velocities into time-mean and time-fluctuating (eddy) components shows eddy fluxes of layer thickness are negative (westward). These fluxes are strongest during the initial spinup when is small and the Reynolds number is large. This suggests that stirring motion due to time-dependent eddies acts to expedite the interior response to transport signals beyond what is achievable with baroclinic Rossby waves alone. Simulations of steady () and time-dependent () flows with and initial layer thickness were computed. The time-dependent flow shows an increased variation of compared to that of the steady flow, however, this increased variation is accompanied by a decrease in variation of the transport . This is opposite to the behavior for given in (24), where an increase in corresponds to an increase in the variation of . Instead, the effect of eddies on projects variations onto fast time scales, diminishing the variability on slow time scales and requiring the observed increase in variation of .

## 4. Nondimensional analysis

In this section we consider the amplitude of the AMOC variability in a nondimensional framework using the analytic solutions presented in section 2c. This provides a general understanding of what controls variability and allows us to connect more directly the results from our analytic model, specifically (18) above, and the delay equations of Johnson and Marshall (2002a,b).

The dimensional solution for the adiabatic flow subject to sinusoidal forcing is now nondimensionalized using the basin width as the characteristic length, average layer thickness *H* as the characteristic height, the mode-1 baroclinic equatorial Rossby wave phase speed as the characteristic velocity, and as the characteristic time. The choice for *T* describes the time it takes the fastest available (long wave) equatorially trapped Rossby wave to cross the basin (Gill 1982). We note the choice for *U* and *T* are only true for domains with sufficiently large zonal extent, that is, for ; however, for domain widths similar to that of the Atlantic, the longest available Rossby waves travel at roughly the speed of *U*.

The solution for the variability of the domain-averaged thickness [(18)] is nondimensionalized using the above scaling, resulting in the nondimensional domain-averaged thickness variability

where is the nondimensional forcing period. Several nondimensional numbers arise: , , . If we define the equatorial deformation radius , then . Defined this way, *λ* and *μ* can be interpreted as nondimensional aspect ratios for the southern basin and equatorial band, respectively. If these parameters are to be representative of NADW, then , , and .

From this form, behavior can be seen to change at the transition from high-frequency forcing to low-frequency forcing, which occurs near , within the factor of . This is the period at which a Rossby wave can cross the basin at the southern latitude of the domain. However, the solution also depends on the northern extent, or overall area of the layer, which gives the scale factor of . For high-frequency forcing , while for low-frequency forcing . The amplitude of is shown in Fig. 7a as a function of *λ* and *P*. We fix and , which are roughly representative of the Atlantic. We find a similar behavior for variations in *μ* and *P* with in Fig. 7b. The closed form solution provides a simple interpretation of the transition from a weak response in the basin-averaged layer thickness (and outflow at the southern boundary) to a strong response.

We now make the connection between our closed form solutions for the average layer thickness and the delay equations derived by Johnson and Marshall (2002b) for the layer thickness on the eastern boundary and the meridional gradient of the transport streamfunction. Using the above scaling, the nondimensional equation for the layer thickness on the eastern boundary [Eq. (14) of Johnson and Marshall (2002b)] is written as

where *λ* and *μ* are as previously defined, is the nondimensional transport into the domain from the north, and *c* is the nondimensional baroclinic Rossby wave speed. Away from the equator, the dimensional Rossby wave speed is given by , while near the equator the wave speed is limited by the mode-1 baroclinic equatorial Rossby wave speed . The nondimensional wave speed *c* is thus defined as

The solution for may be written as

where, for simplicity, we have taken , although the general results are not sensitive to this choice. There are two integrated quantities, defined as

The approximation for is valid as long as , that is, if the meridional extent of the basin is such that . In that case, half of the contribution to comes from the equatorial waveguide and half comes from off the equator.

The characteristic scales , *U*, and *T* are used to derive a nondimensional equation for the change in meridional transport as a function of latitude [Eq. (13) of Johnson and Marshall (2002b)],

We present results as a function of *λ*, *μ*, and *P* to provide a broad perspective on the system response. The forcing is given by a mean inflow from the north, which scales as , and a periodic oscillation of amplitude , written in nondimensional form as

The nondimensional amplitude (scaled by ) of the layer thickness anomaly on the eastern boundary is shown in Fig. 8a as a function of *λ* and *P*. At higher-frequency forcing, the amplitude of the variability on the eastern boundary is small. This is consistent with the calculations reported by Johnson and Marshall (2002a). Small variations in correspond to small variations in the outflowing dimensional transport at the southern boundary since . Figure 8a shows that the amplitude of increases with increasing *P*, but the values of *P* that cause the increase in depend on *λ*. The transition is well predicted by the nondimensional closed form solution [(28)] described in section 2. Larger values of *λ* have longer forcing periods *P* that cause to increase. This dependence is most easily understood by consideration of the flow exchange between the western boundary layer and the interior, now calculated from (33).

First, recall the discussion in section 2b regarding the establishment of the circulation via the initial propagation of boundary waves and subsequent westward propagation of Rossby waves. Note that the initial boundary waves establish meridional pressure gradients about the equatorial waveguide that are later weakened by Rossby waves with propagation speeds that decrease with latitude (Kawase 1987; Johnson and Marshall 2002b). When there is little flow exchange between the western boundary layer and the basin interior, most of the transport exits the domain to the south, as discussed by Johnson and Marshall (2002a,b). This requires that the layer thickness on the eastern boundary have large-amplitude fluctuations in order to support a nondimensional meridional transport at , allowing flow to exit the domain. The total loss of transport from the western boundary, integrated from the southern boundary to the northern boundary, is shown, in nondimensional form, in Fig. 8b as a function of *λ* and *P*. At high-frequency forcing essentially all of the transport anomaly in the western boundary layer is fluxed into the basin interior. This is consistent with the weak signal seen in Fig. 8a for . Since little variability reaches the southern boundary, and differ markedly. This is evident from a comparison of (34) for and the nondimensionalized form of which, using (24) and neglecting nonlinear effects, can be written as

As the period increases and increases, less of the variability is lost to the interior and the southward flow is able to exit at the southern boundary. The transition takes place roughly in accord with the increase in (Fig. 8a). The period at which this transition takes place is given by the time it takes a baroclinic Rossby wave to cross the basin at midlatitudes. For the southern latitude of the model domain and from (28), this is given by , indicated in Fig. 8 by the white line. At forcing periods longer than this, midlatitude Rossby waves have sufficient time to propagate pressure anomalies from the eastern boundary to the western boundary before the forcing changes. This weakens the meridional pressure gradient on the offshore side of the deep western boundary layer, consequently reducing the zonal flow and results in more of the transport anomaly being transmitted to the south.

For very high frequency forcing, the amplitude of increases and slightly less transport is lost from the boundary layer into the interior. This trend becomes more pronounced for smaller values of *λ* (not shown). This is a result of a resonance between the forcing frequency and the time it takes a Rossby wave to cross the equator. At , the signal on the eastern boundary and the western boundary within the equatorial waveguide are the same, meaning that the meridional flow into the interior is weak and the transport remains along the western boundary. The transition to this resonant state begins at , below which the amplitude of begins to increase. For typical North Atlantic parameters, this occurs at a period of about 90 days, and so is not really relevant for AMOC variability and begins to violate the approximation and the assumption that Kelvin waves are fast compared to the forcing time scale. This regime of high-frequency forcing has been explored by Moore (1968) and Marshall and Johnson (2013), where layer thickness is no longer *y* independent and is shown to scale as .

Johnson and Marshall (2002b) found that the flow exchange between the western boundary layer and the interior was concentrated on the equator for high-frequency forcing but became more broadly distributed in latitude at lower frequencies. The degree of equatorial trapping is indicated here by the latitude , at which 75% of the total mass transport is lost from the western boundary layer in the Northern Hemisphere occurs. This is shown in Fig. 8c. Consistent with Johnson and Marshall (2002b), the mass loss is concentrated near the equator for high-frequency forcing and becomes more uniform with latitude ( increases) for low-frequency forcing. It is largely independent of (*λ*). However, somewhat unexpectedly, at very low frequency forcing the exchange again becomes trapped near the equator. This is explained by consideration of (33). For forcing periods much longer than the basin-crossing time scale, the layer thickness on the western boundary is close to that on the eastern boundary at all latitudes. However, the zonal transport is proportional to the change in layer thickness times the wave speed *c*, which is much larger near the equator than at higher latitudes, thus refocusing the exchange near the equator.

An example of the meridional structure of the flow exchange between the boundary layer and the interior is shown in Fig. 9 for three forcing frequencies, . For the exchange is concentrated near the equatorial waveguide and shows the influence of the phase shift just off the equator due to slower high-latitude Rossby waves. The transport amplitude decreases at higher latitudes because . At , the exchange is now broadly distributed over a wider range of latitudes such that . In both these cases the direction of flow exchange between the boundary layer and interior changes sign between the equator and midlatitudes due to latitude dependence for Rossby waves to cross the basin. At , the exchange is all of the same sign with latitude but decays away from the equator. The same qualitative behavior can be seen in the calculations of computed from numerical solutions to the reduced-gravity equations in Fig. 10 where , , and .

A similar analysis of (29) and (33) has been carried out by fixing and varying *μ* between 10 and , as summarized in Fig. 11. In general, increasing *μ* decreases the variability of , with the transition from weak to strong variability occurring for periods longer than (white line). The equatorial trapping at high frequencies is found to be stronger for larger *μ*, although the reemergence of trapping at low frequencies is only weakly dependent on *μ*.

These results from the delay equations of Johnson and Marshall (2002a) are in broad agreement with the closed form solution [(28)] provided here for the average layer thickness. This demonstrates the applicability of the layer-average approach to understanding the time-dependent response and validates the simple interpretation of the results derived from the closed form solutions.

## 5. Conclusions

We have presented a study of southward flowing North Atlantic Deep Water (NADW) at middepths. This study focuses on NADW dynamics in a framework that treats water masses above and below the NADW as motionless. While such a framework is a great simplification of the vertical structure of the meridional overturning circulation in the real ocean, numerical simulations of two-layer reduced-gravity equations with (and without) wind stresses and with buoyantly-driven upper and lower layers display the same behavior as the theory and model configuration presented here, so we feel that this is a useful approximation. The present study details a simplified and parallel description of the propagation of NADW anomalies to that provided by the work of Johnson and Marshall (2002a). The theory presented here permits closed-form solutions that are consistent with the numerical solutions to the delay equation provided by Johnson and Marshall (2002a) and numerical simulations of reduced-gravity equations for laminar, diabatic, and moderately turbulent regimes.

The theory for domain-averaged layer thickness of North Atlantic Deep Water derived from the continuity equation is valid in the limit of slow time-scale variability of the layer thickness on the eastern boundary to the relatively fast time scale for a Rossby wave to propagate across the basin, permitting the approximation . This allows closed form analytic solutions, which are provided when the southward volume flux at the northern boundary is constant or sinusoidal and both northern and southern boundary conditions are geostrophically balanced. However, numerical simulations and analytic solutions with high-frequency forcing predict a small-amplitude response for layer thickness; therefore, any errors incurred in this high-frequency limit by assuming do not diminish the predictions of the theory. To allow a detailed examination of the effect on the circulation, a southward transport signal with a fixed frequency is imposed at the northern boundary and is solely responsible for driving fluid motions. The imposed forcing is meant to mimic waters that downwell at high latitudes through surface buoyancy forcing or entrainment into descending overflows. Although the periodic forcing provides convenient analytic solutions, in the linear limit any time series of transport variability could be reproduced by the theory through a Fourier series.

The choice of northern and southern boundary conditions in our theoretical model permits closed analytic solutions; however, more general boundary conditions may be used to incorporate forcing due, for example, to Ekman transports that may still yield simple analytic solutions. In fact, the inclusion of a southward Ekman transport for the southern boundary condition is a trivial modification to the differential equation governing and its solutions. At the same time these boundary conditions place limitations on the nature of the solutions for since their analytic forms determine the tractability of the domain-averaged continuity equation. If the choice for these boundary conditions yield a tractable differential equation for , then analytic solutions amenable to analysis may be obtained. We claim that the closed-form solutions for are significant to the extent that their expressions are simple, compact, and readily analyzed.

Similar to findings by Johnson and Marshall (2002b), basin geometry and stratification are found to set the intrinsic time scale *τ* and equilibrium layer thickness *H*. Moreover, when forcing at the northern boundary is sinusoidal, analytic solutions yield a closed form expression for the amplitude of domain-averaged layer thickness anomalies. A nondimensional analysis generalizes the dynamic behavior observed and relates the theory developed here to that of Johnson and Marshall (2002a,b). The nondimensional parameters related to basin geometry are the aspect ratios , , and , while *P* describes the nondimensional forcing period. These parameters are used to understand the system response to a range of forcing frequencies and aspect ratios, including values relevant to NADW. Particularly, loss (gain) of meridional transport to (from) the basin interior about the equator is explained within this nondimensional framework. Roughly speaking, if forcing periods are less than the time it takes a Rossby wave to cross the basin at the latitude of the southern boundary, the variability in average layer thickness is weak while for longer forcing periods the variability is large.

## Acknowledgments

This research was supported by the National Science Foundation under Grant OCE-1634468. The authors thank David Marshall and Helen Johnson for their helpful conversations. We extend thanks to an anonymous reviewer whose helpful comments improved the presentation and clarity of this work.

### APPENDIX

#### Analytic Solutions to the Domain-Averaged Reduced-Gravity Continuity Equation

Here we outline the derivation of the analytic solutions to the domain-averaged continuity equation and provide the associated intrinsic time scales and equilibrium thicknesses. To begin, note that variability of volume transport drives variability of layer thickness. We assume variability of volume transport into the domain be such that resulting layer thickness variability is small compared to the equilibrium thickness. This assumption can be expressed by introducing the dimensional parameter

where is the amplitude of the imposed volume transport variability, *A* is the area of the domain, *τ* is some time scale to be determined, and *H* is the asymptotic () equilibrium thickness. The parameter *ε* represents the characteristic scale for layer thickness perturbations due to the transport variability . We pose a solution to (14), using *ε* as the expansion parameter, of the form

Equation (A2) is dimensional where ; however, the second term in the sum represents layer thickness perturbations that scale as *ε*, and thus is nondimensional and . Substitution of (A2) in (14) yields

where, for the moment, we allow , solely to illustrate that the method of solution is similar when ≫ . At leading order, , the governing equation is

and describes nonadiabatic flow in the presence of steady forcing. If ≫ the damping term becomes a higher-order effect and the flow is adiabatic to leading order. To first order, , the governing equation is

This captures the effects of time-dependent volume transport at the northern boundary, the effect of outflow at the southern boundary, and the nonadiabatic effects.

While the order-by-order solution below continues to assume , adiabatic solutions can be recovered by setting . More precisely, adiabatic flow requires the time scale for imposed thickness relaxation, that is, *γ*, be sufficiently slow such that nonadiabatic effects are negligible. To enforce this in the context of adiabatic flow, we assume ≫ , where *τ* is the intrinsic time scale of the flow. Otherwise, adiabatic effects may be considered by allowing . This adiabatic regime is most comparable to the weak damping case studied by Kawase (1987); however, when ≫ damping via nonadiabatic effects is sufficiently weak such that Rossby waves propagate virtually unhindered.

##### a. : Nonadiabatic flow, steady forcing

The leading-order equation [(A4)] is separable with the right-hand side that can be written as a quadratic, that is,

with real constant coefficients

When (which is the case for a southern boundary residing south of the equator), then , , and it follows that and (A6) may be written as

The values and are the real and distinct roots of the quadratic in (A6). Integration of (A7) along with the initial condition yields the solution

where , and *τ* is the intrinsic time scale for the saturation of layer thickness to its equilibrium value . The roots correspond to layer thicknesses for which and obtain in the limits , respectively. The time scale *τ* and the equilibrium layer thicknesses are given by

##### b. : Nonadiabatic flow, periodic forcing

In the presence of imposed variability of volume transport into the domain, that is, small but nonzero *ε*, the solution to (A5) provides the dynamic response. While is now known, the solution for may be simplified by Taylor-expanding at large *t* in (A5) rendering the following linear equation

with initial condition . The first term on the right-hand side results from the applied volumetric flux prescribed as a time-dependent condition that varies sinusoidally with amplitude , frequency . Recalling that *ε* is the characteristic scale for layer thickness perturbations, the order-one variable is dimensionless and has solution

where and and *τ* is determined by (A9). The domain-averaged layer thickness for nonadiabatic flow with imposed periodic volume transport has solution

Note that in the limit , decreases, and solutions with time-dependent forcing resemble those with steady forcing.

##### c. Adiabatic flow, steady or periodic forcing

The adiabatic solutions follow as a special case of the nonadiabatic flow. The flow may be considered adiabatic when ≫ , which may be accomplished by setting the formally order-one quantity to zero. Doing so leaves the coefficient *a* unchanged, while and . This case has solution of the form given in (A8), however, with modified time scale and equilibrium thicknesses given by

In the case of periodic forcing, the solution has these same values for time scale and equilibrium thickness, but the solution is of the form given in (A13).

## REFERENCES

*Atmosphere–Ocean Dynamics.*Academic Press, 662 pp.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).