## Abstract

Mesoscale eddies shape the Beaufort Gyre response to Ekman pumping, but their transient dynamics are poorly understood. Climate models commonly use the Gent–McWilliams (GM) parameterization, taking the eddy streamfunction to be proportional to an isopycnal slope *s* and an eddy diffusivity *K*. This local-in-time parameterization leads to exponential equilibration of currents. Here, an idealized, eddy-resolving Beaufort Gyre model is used to demonstrate that carries a finite memory of past ocean states, violating a key GM assumption. As a consequence, an equilibrating gyre follows a spiral sink trajectory implying the existence of a damped mode of variability—the eddy memory (EM) mode. The EM mode manifests during the spinup as a 15% overshoot in isopycnal slope (2000 km^{3} freshwater content overshoot) and cannot be explained by the GM parameterization. An improved parameterization is developed, such that is proportional to an effective isopycnal slope , carrying a finite memory *γ* of past slopes. Introducing eddy memory explains the model results and brings to light an oscillation with a period ≈ 50 yr, where the eddy diffusion time scale *T*_{E} ~ 10 yr and *γ* ≈ 6 yr are diagnosed from the eddy-resolving model. The EM mode increases the Ekman-driven gyre variance by *γ*/*T*_{E} ≈ 50% ± 15%, a fraction that stays relatively constant despite both time scales decreasing with increased mean forcing. This study suggests that the EM mode is a general property of rotating turbulent flows and highlights the need for better observational constraints on transient eddy field characteristics.

## 1. Introduction

The Beaufort Gyre, a major anticyclonic circulation feature in the Arctic Ocean, hosts a substantial fraction of the overall Arctic freshening (Haine et al. 2015). The large-scale gyre circulation has been directly linked to its freshwater content (FWC) via the process of Ekman pumping that converges relatively fresh surface waters and deepens the halocline (e.g., Proshutinsky et al. 2002). The Ekman pumping arises due to transient anticyclonic winds that cause significant gyre variability on interannual and longer time scales. However, observations indicate that the halocline depth (roughly equivalent to the FWC) varies, but does not always mimic, the variability in the strength of the anticyclonic wind stress (e.g., Proshutinsky et al. 2009; Giles et al. 2012). Understanding and modeling the variability of the large-scale gyre circulation and the associated FWC remain a challenging problem.

The availability of FW sources, the strength of Ekman pumping, and interactions with the Atlantic layer are all factors external to the large-scale circulation that have been studied with respect to Beaufort Gyre variability (Proshutinsky et al. 2002; Giles et al. 2012; Martin et al. 2014; Morison et al. 2012; *S*tewart and Haine 2013; Lique and Johnson 2015; Lique et al. 2015). A recent study, however, emphasized internal dynamics of the gyre and demonstrated that the large-scale halocline deepening due to Ekman pumping is counteracted by the cumulative action of mesoscale eddies (Manucharyan and Spall 2016), a dynamical balance similar to that of the Antarctic Circumpolar Current (ACC; Marshall and Radko 2003) or the Weddell Gyre (Su et al. 2014).

Unlike the ACC, which might be in a so-called eddy saturation regime with weak equilibrium sensitivity of its transport and density structure (Böning et al. 2008; Tansley and Marshall 2001; Hallberg and Gnanadesikan 2001; Munday et al. 2013; Abernathey and Cessi 2014), the Beaufort Gyre appears to be highly sensitive to variations in the Ekman pumping (Proshutinsky et al. 2002; Manucharyan and Spall 2016). This implies that a continuing melting of the sea ice and/or changes in the atmospheric winds that alter the Ekman pumping would lead to substantial changes in the halocline depth and hence in the freshwater budget of the entire Arctic Ocean. Since the Arctic Ocean is rapidly changing, it is important to understand mechanisms behind the transient gyre adjustment.

Manucharyan et al. (2016) point out that the eddy field also plays a key role in Ekman-driven gyre variability. Because of the tendency of the eddy streamfunction to oppose the Ekman transport, a large-scale circulation (i.e., averaged over small-scale eddies) has to be a stable dynamical system that equilibrates on a time scale inversely proportional to the eddy diffusivity *K* (Manucharyan et al. 2016). These theoretical conclusions rely upon the use of a conventional Gent–McWilliams (GM) eddy parameterization (Gent and McWilliams 1990; Gent et al. 1995) that takes the isopycnal eddy thickness fluxes to be proportional to large-scale halocline thickness gradients:

where **v**′ and *h*′ are the velocity and halocline thickness perturbations due to time-dependent motions and the overbar denotes an ensemble average.

The GM parameterization is used in nearly all climate models that do not resolve mesoscale eddies and requires specification of an eddy diffusivity. Marshall and Radko (2003) demonstrate that the equilibrium isopycnal slope *s* of an eddying baroclinic current is directly related to *K*. In addition, Manucharyan et al. (2016) show that the temporal variability of *s* also depends on the sensitivity of eddy diffusivity to halocline slope *dK*/*ds* (see section 4). Thus, climate models that use too small (large) values of eddy diffusivity would result in too deep (shallow) halocline. If the models prescribe too weak (strong) sensitivity of eddy diffusivity to isopycnal slope the temporal variability of the currents would be overestimated (underestimated). The GM parameters *K* and *dK*/*ds* can be either diagnosed from eddy-resolving models or with sufficient coverage inferred through observations.

A key assumption of the GM parameterization, which is to be challenged in this manuscript, is that the eddy fluxes at a particular point in time and space depend only on the large-scale state of the ocean at the same point in time and space. The assumption of temporal and spatial locality can be questioned from both observational and numerical modeling evidence. Observations suggest that mesoscale eddies (in a form of coherent vortices) can persist in the open ocean for years, propagating large distances from their formation regions (Chelton et al. 2011). With respect to the Arctic Ocean, using a numerical model Spall et al. (2008) discuss how shelfbreak eddies propagate away from their formation region and can be further transported by the mean current of the Beaufort Gyre. Manucharyan and Timmermans (2013) discuss a self-propagation mechanism of the submixed layer Arctic eddies that are observed to advect the buoyancy and potential vorticity anomalies up to 500 km away from their presumed formation regions (Timmermans et al. 2008). Thus, it is reasonable to assume that a large-scale eddy field does not only depend on a current state of the ocean but carries a finite memory of its past states and a history of dissipative processes.

We illustrate the persistence of eddy properties schematically in Fig. 1. Consider the equilibration of an ocean that is populated with blue eddies that are less energetic than the equilibrium-type red eddies (Fig. 1, left). By different types of eddies we imply the existence of statistical eddy properties, that is, sizes, eddy kinetic energy, or eddy transport, that are generated by the mean flow. While the mean currents are generating the red-type eddies, the number of blue-type eddies will be decreasing as they are dissipated or absorbed by the mean flow (Fig. 1, center). Eventually, the ocean will be populated only with the equilibrium-type red eddies (Fig. 1, right). Since the eddy transport is associated with an eddy field that has a memory, it is important to understand the time scales of this eddy adjustment process and how it feeds back on the evolution of the large-scale current.

In this manuscript, we use an idealized eddy-resolving model of the Beaufort Gyre (section 2) to present evidence that eddy memory significantly affects the gyre dynamics (section 3). We discuss dynamical constraints that a conventional local-in-time GM parameterization imposes on large-scale ocean dynamics (section 4) and suggest an improvement that includes the eddy memory effect (section 5). Using the new eddy parameterization we bring to light a low-frequency mode of variability of large-scale flows associated with the eddy memory (section 6). Characteristics of the mode and its effects on the variability of the ocean circulation are assessed in section 7. We summarize and discuss implications in section 8.

## 2. Model configuration

We use an idealized model of the Beaufort Gyre in a configuration identical to the one used in Manucharyan and Spall (2016, see their appendix A). The gyre is represented by a cylindrical ocean basin (radius *R* = 600 km, depth 900 m) driven by an anticyclonic surface stress *τ*(*r*) corresponding to a uniform Ekman pumping. The fluid dynamical equations are solved using the MITgcm in its three-dimensional hydrostatic configuration (Marshall et al. 1997). A 4-km horizontal resolution along with variable in-depth vertical resolution between 10 and 60 m is sufficient to permit Rossby deformation-scale eddies (the baroclinic deformation radius is about 20 km in these simulations). The salinity is restored at the edges of the gyre to a fixed salinity profile that consists of a 50-m-deep surface layer of relatively fresh waters of salinity 29, followed by a lower layer of salinity 34 [values are chosen to mimic the hydrography of the Beaufort Gyre (Steele et al. 2001)]. Strong salinity restoring leads to a fixed halocline depth at the gyre boundaries. This configuration, as highlighted in earlier studies (Manucharyan and Spall 2016; Manucharyan et al. 2016), provides an infinite reservoir of freshwater and thus the results here should be treated as an upper bound on the possible fluctuations in the halocline volume. Consistent with Manucharyan and Spall (2016), we include a continental slope: lateral boundaries are vertical down to 300 m, below which the depth increases linearly to the bottom of the basin penetrating 100 km toward the center of the gyre (see Fig. 2b). Note that Manucharyan et al. (2016) considered a flat topography case, which is more applicable to the interior of the gyre. However, we demonstrate in section 8 that the continental slope is essential for the gyre dynamics discussed here.

## 3. Signatures of eddy memory

Spinup simulations are initialized with a horizontally uniform stratification (50-m initial halocline depth) and forced with a spatially uniform, temporally invariant Ekman pumping [corresponding to a linear stress profile , where is the magnitude]. Following the development of the mesoscale eddy field, the idealized Beaufort Gyre model achieves a statistical steady state after several decades (Fig. 2a). Note that this time scale is significantly longer than the 5-yr-long *e*-folding equilibration in a flat basin case considered in Manucharyan et al. (2016). The prolonged equilibration is entirely due to the presence of a continental slope in our simulations. The slope suppresses the development of instabilities and locally reduces the mesoscale eddy diffusivity (Isachsen 2011; Stewart and Thompson 2013) that leads to longer eddy diffusion time scale. Note that, because the gyre dynamics are governed by a halocline thickness diffusion equation, any localized reduction in the eddy diffusivity would have a basinwide influence on the halocline depth and its equilibration time scale.

The equilibrium corresponds to a vanishing residual circulation that is supported by the eddy buoyancy transport that counteracts the Ekman transport. Most of the eddy kinetic energy is concentrated in the upper layer corresponding to the first baroclinic mode (see Fig. 2b). Characteristic eddy velocities in the upper 200 m of the gyre are 0.05–0.1 m s^{−1}. These eddies are about 100–200 km in diameter but are relatively weak compared to more intense but smaller-scale submixed layer eddies (about 20 km in diameter) commonly observed via the ice-tethered profilers (Zhao et al. 2014). Nonetheless, these larger-scale eddies are responsible for maintaining key properties of the halocline such as its depth and adjustment time scale.

In our idealized gyre simulations, the FWC [conventionally defined as a domain-integrated measure of salt content relative to a reference salinity *S*_{ref} (e.g., Proshutinsky et al. 2002)] is directly proportional to the volume of water above the halocline *V* and to the salinity difference Δ*S* across the halocline:

and *h* is the halocline depth. While the equilibration of the gyre FWC and eddy kinetic energy (EKE) can be crudely viewed as an exponential adjustment, a closer consideration reveals important deviations (Fig. 2a). In particular, there exists a significant overshoot of the FWC that lasts for several decades (Fig. 2a, red). The maximum depth of the halocline (at year 30) is up to 10–15 m deeper than its equilibrium value (Fig. 2b). This corresponds to a FWC overshoot of about 2000 km^{3} (Fig. 2a), which is comparable to an observed Beaufort Gyre FWC increase of 3000 km^{3} over the past two decades (Haine et al. 2015).

During the early stages of the spinup, the generated eddies are less energetic as the currents are weak. As a result of their persistence, the halocline deepens beyond its equilibrium because the eddy field is insufficiently energetic to counteract the Ekman transport (see Fig. 2a around model year 10). After the FWC reaches its maximum values and the halocline is at its deepest levels (around model year 20), the eddy field becomes overly energetic, generating excessive thickness fluxes that reduce the halocline depth back to its equilibrium. This cycle is manifested as an overshoot in the EKE that is lagged with respect to the FWC by about 5–7 yr (Fig. 2a, blue). We note that this overshoot behavior is not only associated with the “cold” spinup from the state of rest but also exists for perturbation experiments starting from a fully turbulent equilibrated state (additional experiments are not shown here as they are qualitatively similar).

These lagged overshoots in EKE and FWC are signatures of an oscillatory mode that operates in addition to the exponential decay. Since we observe only one full oscillation before the gyre equilibrates, this mode is heavily damped and needs external forcing to be sustained. Nonetheless, the existence of this damped mode does not fit into our traditional understanding of the mesoscale eddy dynamics as viewed through the lens of the local-in-time GM parameterization.

## 4. Dynamical implications of GM parameterization

Here, we briefly discuss the conventional GM parameterization and a key constraint that it imposes on the evolution of large-scale flows. In particular, we focus on the dynamics of a large class of currents that are forced by transient Ekman pumping with nonzero mean. When the diabatic forcing is small compared to the mean Ekman pumping, the time-averaged state corresponds to a vanishing residual-mean circulation. Transient forcing, however, can produce significant deviations from the equilibrated state.

The GM parameterization operates under the assumption of slowly evolving ocean dynamics such that at any given moment in time and space the mesoscale eddy field can be considered in equilibrium with local, large-scale currents. In particular, it assumes that the eddy streamfunction is proportional to the time-dependent halocline slope *s*(*t*) and an eddy diffusivity *K* that can be slope dependent (Visbeck et al. 1997) but is typically not time dependent. Near the equilibrium halocline slope *s*_{0}(*r*), the perturbation eddy streamfunction can be expressed as

where primes denote perturbations from equilibrium and is a constant-in-time effective eddy diffusivity for the perturbations

Note that the eddy diffusivity for a linearized system is always larger than the background diffusivity because *K* is expected to monotonically increase with *s*. For example, Visbeck et al. (1997) suggest a linear dependence of the eddy diffusivity on halocline slope, which is found to be appropriate for the Beaufort Gyre (Manucharyan et al. 2016). In this case, , a factor of 2 increase. For eddy-saturated currents, like the ACC, the term (*dK*/*ds*)*s*_{0} might provide a dominant contribution to the effective eddy diffusivity .

In cylindrical coordinates the halocline thickness evolution under the GM parameterization obeys a forced diffusion equation:

where *w*_{E}(*r*, *t*) is the time-dependent Ekman pumping, and *r* is the radial coordinate [for a derivation, see appendix A of Manucharyan et al. (2016)]. The eddies act as a thickness diffusion because GM parameterization assumes the eddy thickness flux to be proportional to the halocline slope (). Note the use of instead of *K* in Eq. (5) since it is written for perturbation variables.

Equation (5) provides further insight into the halocline volume dynamics. Consider the least-damped eigenfunction of the diffusion operator on the right-hand side of Eq. (5):

where defines the gyre equilibration time scale (inverse of the smallest eigenvalue). Manucharyan et al. (2016) demonstrate that higher-mode eigenfunctions of this eddy diffusion operator are highly damped (damping time scale increases quadratically with the number of zero crossings) and hence only the least-damped eigenfunction can significantly contribute to changes in the halocline volume. Thus, integrating Eq. (5) over the domain, keeping only a contribution from the least-damped eigenmode, and using Eq. (6), we arrive at

where the overdot indicates the time derivative, and *W*_{E} is the Ekman transport (domain-integrated Ekman pumping projection onto the least-damped eigenmode). Equation (7) demonstrates that under the GM parameterization the large-scale currents are constrained to equilibrate exponentially, that is, no internal oscillations are possible. Numerical simulations, however, demonstrate that an eddying current can significantly deviate from exponential equilibration (Fig. 2), exhibiting oscillatory behavior that cannot be captured by the local-in-time GM parameterization.

## 5. Eddy memory parameterization

Here, we introduce an improvement to the GM parameterization by accounting for the eddy memory and validate its relevance in the eddy-resolving model.

### a. Parameterization

We make a key assumption that the eddy streamfunction has a finite memory of past states:

where *γ* is the eddy memory time scale, and this definition is applicable for sufficiently small perturbations away from the mean state such that does not depend on slope perturbations *s*. In a fully nonlinear formulation, a constant-in-time should be replaced by a slope-dependent eddy diffusivity *K* that can vary in time. The integral form of Eq. (8) implies that the present eddy transport (at time *t*) consists of contributions of past transports (at times *t*′ *< t*) with weights that are exponentially decreasing with time toward the past. Contributions from past transports quantify the eddy persistence (see Fig. 1 and discussions thereof). Note that in the limit of no memory *γ* → 0 or in equilibrium *s*(*t*) = constant, this parameterization [Eq. (8)] is identical to the conventional GM parameterization [Eqs. (3)–(4)].

Since is defined by the mean state and hence is independent of time, our parameterization can be written in the following way:

where we have defined an effective slope that contains memory of past transports and governs the eddy transport. Alternatively, one can define an effective eddy diffusivity such that ; however, the effective diffusivity is not a physically relevant quantity, and it can be inconvenient to use in practice. For example, *K*_{eff} is unbounded when eddies generated in the past are still contributing to the transport while a present slope is negligibly small (i.e., *K*_{eff} → ∞ when *s* → 0). The extremely large *K*_{eff} is misleading since a physically relevant eddy transport () is finite. Hence, *K*_{eff} cannot be directly interpreted as a measure of the mesoscale eddy activity. We thus choose to proceed with our analysis using the effective slope [Eq. (10)] as a physically relevant and singularity-free quantity that defines the eddy thickness transport.

Differentiating Eq. (10) with respect to time leads to an alternative definition of the effective slope

It implies that is exponentially pulled back to the actual slope *s* on a time scale *γ*. Integrating Eq. (11) forward in time is an efficient and simple method of calculating in numerical models as it avoids heavy calculations of the integral in Eq. (10) at every time step. Note that the eddy memory can lead to periods of upgradient thickness transport when and *s* have different signs. However, Eq. (11) implies that the eddy memory does not affect steady or slowly evolving mean currents for which and hence the long-term average of the eddy thickness transport is always downgradient.

The eddy memory parameterization [Eq. (10)] involves an effective slope that is an additional prognostic variable defining the eddy transport. This approach is complementary but different from mixing length arguments (Holloway 1986) that use EKE (or characteristic eddy velocity) as a variable to determine the intensity of the eddy transport [see an example for the ACC in Ferrari and Nikurashin (2010) and Sinha and Abernathey (2016)]. While both and EKE can be linked to the eddy transport using scaling laws, these are conceptually different measures of turbulent flows. The eddy kinetic energy () does not contain information about a correlation between velocity and thickness anomalies that constitute the eddy transport (). Hence, additional arguments are needed to link EKE to eddy transport. However, mesoscale eddies, viewed as coherent vortices, propagate in space and persist in time together with their corresponding isopycnal thickness anomalies (e.g., quasigeostrophic eddy currents explicitly depend on layer thickness anomalies), implying not only a persistence of EKE but also a persistence of a thickness flux. Thus, we consider the eddy streamfunction (a measure of eddy tracer fluxes) as a dynamically relevant quantity that has a finite memory of past ocean states.

### b. Diagnosing eddy memory

We now quantify the impact of the eddy memory as it relates to the eddy thickness transport in the eddy-resolving Beaufort Gyre model. We diagnose time series for the eddy streamfunction by calculating the eddy thickness fluxes during the gyre spinup simulation. Figure 3a demonstrates that and *s*(*t*) when plotted against each other have a relationship characteristic of a spiral sink. First, the eddy field and halocline slope do not equilibrate following a conventional straight line path predicted by the GM parameterization; for example, the blue curve has large deviations from the straight line in Fig. 3a. Instead, away from equilibrium the eddy streamfunction evolves more slowly compared to the halocline slope (blue curve is above the dashed black curve in 3a). Second, the diagnosed and *s* loop approaches the equilibrium in a spiraling trajectory with a quickly decaying radius of the spiral (see the black arrows around 0, 0 in Fig. 3a). Physically, this view is consistent with the isopycnal slope overshooting its equilibrium value (Fig. 2a) because the mesoscale eddy field at that time is not fully developed. At the peak of the isopycnal overshoot, the excessively baroclinic currents start to generate overly energetic eddies that weaken the currents below their equilibrium strength. Such a dynamical behavior is consistent with a damped oscillator.

We do not know the value of *γ* a priori, but we can attempt to infer it by assessing the correlation between the eddy streamfunction and the effective slope defined by Eq. (10). Calculations show that there is indeed an optimal value of *γ* ≈ 6 yr that enhances the correlation (Fig. 3b, blue). Throughout the gyre equilibration is better approximated as a linear function of rather than *s* (Fig. 3a), with the best linear fit of = 150 m^{2} s^{−1} (Fig. 3a, black dashed).

The inferred eddy diffusivity is smaller than the value reported in Manucharyan et al. (2016) because of the continental slope that suppresses baroclinic instabilities. The continental slope also plays a major role in enhancing the eddy memory, since in the interior of the gyre *γ* ≈ 2 yr (Fig. 3b, red). This suggests that there might be a relation between the magnitude of eddy memory and the eddy adjustment time scale that is inversely proportional to the eddy diffusivity. Thus, for regions with high eddy diffusion, the eddy memory is small and vice versa.

## 6. Emergence of the eddy memory mode

Now that we have established the physical basis and diagnosed the eddy memory in the gyre, we proceed to reveal the newly emerging dynamics. Although most of the following derivation will focus on the evolution of *V*, the eddy memory affects all related gyre characteristics such as FWC and *s*. Applying our modified GM parameterization, the perturbations in halocline depth *h* with respect to the equilibrium or mean state of the gyre evolve following a thickness diffusion equation:

where we have introduced an effective halocline depth as . The two terms on the right-hand side of Eq. (12) represent the divergence of the eddy thickness flux and the Ekman pumping. Note that we consider axisymmetric solutions in cylindrical coordinates, and all variables in Eqs. (12)–(13) are perturbations from the equilibrium state corresponding to forcing by the mean Ekman pumping.

Combining the equations above to eliminate we obtain

Note that the model-diagnosed eddy memory is spatially inhomogeneous, but for simplicity of the analytical analysis, we are assuming a constant parameter *γ*, which should be interpreted as an effective memory that affects the bulk gyre dynamics. In the absence of forcing (*w*_{E} = 0) this equation describes the equilibration of the gyre by exponentially damped waves. To further illuminate the dynamics, let us consider the evolution of the halocline volume *V*. Domain integrating Eq. (14), keeping only a contribution from the least-damped eigenmode, and using Eq. (6), we arrive at

where the overdot indicates the time derivative, and *W*_{E} is the Ekman transport [as in Eq. (7)].

Equation (15) illuminates the core internal dynamics behind the equilibration of the halocline, an externally forced damped oscillator. The ratio of the eddy memory to the eddy adjustment time scale *γ*/*T*_{E} determines whether solutions are either overdamped (nonoscillatory) or underdamped (oscillatory). Note that in the absence of memory (a limit of *γ* → 0) Eq. (15) becomes identical to an exponential decay equation [Eq. (7)], derived using a conventional GM parameterization. However, if the eddy memory is sufficiently large (*γ* > 0.25*T*_{E}, as shown below) the system oscillates with the frequency *ω*_{0} expressed as

and hence the period would be proportional to the geometric mean between the eddy memory and eddy diffusion time scales. Using our model-based estimates for the Beaufort Gyre (*γ* ≈ 6 yr, *T*_{E} ≈ 10 yr), we obtain a period of the EM mode *T*_{0} ≈ 50 yr. The damping time scale for the oscillations is given by 2*γ* = 12 yr (much shorter than its period), and hence the mode is highly damped, requiring continuous external forcing to be sustained. While the EM mode has a distinct multidecadal period, its amplitude has a significant response to a wide range of forcing frequencies because of its strong damping. We speculate that the transience of the atmospheric Beaufort high pressure system can efficiently energize this mode.

Solving Eq. (15) for the initial conditions from the spinup simulation (shown in Fig. 2a), taking *W*_{E} = 0 since there are no Ekman pumping perturbations during the spinup, and using *γ* = 6 yr (as implied by Fig. 3b, red) significantly improve the theoretical prediction of the numerically diagnosed evolution of FWC (see Fig. 4a). In particular, our new theory captures the amplitude and duration of the overshoot in addition to the overall exponential equilibration. Furthermore, it captures a lag between the peaks in FWC and the eddy transport (observe that is proportional to in Fig. 4a). Since the inclusion of the EM mode dramatically improves the representation of the halocline dynamics, we proceed to explore several of its major implications.

## 7. Role of the EM mode in halocline dynamics

### a. Halocline equilibration

The equilibration of the FWC anomalies is represented by damped oscillations that can be expressed in a form of complex exponentials *V* ~ exp(−*λt*), where the real part of *λ* corresponds to the amplitude decay rate and the imaginary part of *λ* corresponds to the oscillation frequency. Plugging this solution into Eq. (15), we get two possible values:

If *γ* > 0.25*T*_{E}, the solution oscillates as *λ* has an imaginary part. These oscillations decay in amplitude with time and, for arbitrary initial conditions, their decay rate corresponds to the smallest of the real parts of the two characteristic equation roots [Eq. (17)]. The inverse of the decay rate is the time scale of the gyre equilibration *T*_{eq}, which depends on both *γ* and *T*_{E}. In the absence of memory (*γ* → 0), *λ* = 1/*T*_{E} and hence *T*_{eq} = *T*_{E}, consistent with an exponential gyre equilibration that was discussed in section 4. The presence of memory *γ*/*T*_{E} < 0.5 leads to a reduction of the equilibration time making the gyre more stable despite the presence of the oscillations (by increased stability we imply larger decay rates). In fact, having a memory *γ* = 0.25*T*_{0} reduces the equilibration time by a factor of 2. This is a critical damping limit where the gyre approaches equilibrium state in the fastest possible way without oscillating. Common examples of such critically damped systems include door closers seen on many hinged doors or shock absorbers in car suspensions.

Thus, the fundamental gyre dynamics are strongly dependent on the ratio *γ*/*T*. It is an underdamped oscillator if *γ* > 0.25*T*_{E} and is a faster-equilibrating overdamped oscillator for *γ* < 0.25*T*_{E}. For *γ* > 0.5*T*_{E} the equilibration becomes slower as compared to the memoryless limit. In this limit, the prolonged memory results in the eddy field that is in large disequilibrium with the current forcing such that the overshoots and undershoots in the isopycnal slopes are less damped.

Our numerical simulations conducted for a wide range of Ekman pumping forcing suggest that there is a relation between the two time scales. Figure 4b shows that for a wide range of mean Ekman pumping, most of the diagnosed time scales lie close to *γ* = 0.5*T*_{E}, a boundary at which the gyre equilibration time *T*_{eq} = 2*γ* = *T*_{E} [as inferred from Eq. (17)] would be exactly equal to the eddy diffusion time scale. Indeed, stronger Ekman pumping forcing leads to shorter time scale because of the large-eddy diffusivities. However, a strong flow would lead to a faster reduction in the persistence of the eddy field due to enhanced eddy–mean flow interactions. The opposite occurs for weak Ekman pumping. A mechanistic understanding of the parameter regime where the ratio *γ*/*T*_{E} is constant remains an open question.

### b. Enhanced halocline variability

The FWC of the gyre in our surface stress–driven simulations is directly proportional to the halocline volume *V*. Consider now the variability of *V* for the gyre forced by transient Ekman pumping by numerically simulating Eq. (15). We take the parameters *γ* = 6 yr and *T*_{E} = 10 yr as diagnosed from the eddy-resolving model. For simplicity, we represent *W*(*t*) as a white noise process that has equal energy at all frequencies and highlight the impact of the EM mode by comparing a simulation to the case of *γ* = 0.

Figure 5a compares the FWC evolution with and without memory. The amplitude of FWC variations is larger with the EM mode due to the overshoots that are particularly prominent when decadal trends are present. For example, near years 70 and 130 the EM mode gives an additional 2000 km^{3} of FWC anomaly for a gyre that would otherwise have a 4000 km^{3} oscillation in FWC. That is a 50% increase in the amplitude of FWC, comparable to the observed FWC increase of 3000 km^{3} in the Beaufort Gyre (Haine et al. 2015).

It is perhaps more illustrative to assess the effects of EM mode in frequency space. According to Eq. (15), the spectrum of *V* depends on the eddy memory in the following way:

where *σ*^{2} represents the spectral energy of the Ekman transport *W*_{E} (a white noise process with equal energy distribution for all frequencies) and , as defined in Eq. (16). In the memoryless limit (*γ* = 0) we recover an expected red noise spectrum:

For comparison, both spectra (with and without memory) are plotted in Fig. 5b, demonstrating an enhanced energy at all frequencies. The two spectra approach the same values at very low frequencies (for which ) as well as at high frequencies for which Ekman pumping dominates the dynamics and eddies do not play a significant role. Note that the frequency *ω*_{max} of peak EM-mode energy is significantly shifted from *ω*_{0} toward lower values; in particular, Eq. (18) dictates that *ω*_{max} ≈ 0.34*ω*_{0} when *γ* = 0.5*T*_{E}. However, it is the maximum relative increase of spectral energy that occurs at *ω* = *ω*_{0}.

We can assess a total variance of *V* by taking the integral of its power spectral density over all frequencies:

Here, the integral has been calculated exactly via the Cauchy’s residue theorem, making use of the integrand having four simple poles on a complex plane . Equation (20) implies that the EM mode enhances the variance by a fraction *γ*/*T*_{E}. The ratio of the diagnosed values for *γ* and *T*_{E} is *γ*/*T*_{E} = 0.5 ± 0.15 (see Fig. 4b), and hence the expected increase in FWC variance should be 50% ± 15%. Note that the standard deviation is a square root of variance such that the contribution of the EM mode is approximately 0.25*γ*/*T*_{E}. Overall, the standard deviation of FWC time series is about 2000 km^{3} without and 3000 km^{3} with the EM mode. Thus, Fig. 5 together with our analytical calculations demonstrates a clear enhancement of FWC variability due to the EM mode. This extra variance is not accounted for in climate models that implement local-in-time eddy parameterizations.

## 8. Summary and discussion

An Ekman-driven, eddy-resolving model of the Beaufort Gyre was used to assess the large-scale impacts of the eddy memory. The key manifestations of the eddy memory are the overshoots in halocline slope and a lagged development of the eddy kinetic energy (Fig. 4). These features cannot be represented by the conventional GM parameterization that assumes time locality of eddy fluxes.

Overshoots in FWC of the simulated Beaufort Gyre reach 2000 km^{3}, a magnitude comparable to FWC variations observed over the past two decades. Note that because there are no sufficient observations of the eddy field in the Arctic Ocean, previous attempts to explain the gyre variability via Ekman pumping likely carry a significant uncertainty due to the eddy thickness fluxes that are unaccounted for.

Using a transformed Eulerian-mean theory, we diagnose the time-dependent eddy streamfunction and show that it is more closely related to the effective slope that takes into account the history of ocean evolution (Fig. 3a) than to the present value of isopycnal slope *s* (as assumed by the GM parameterization). With Eq. (8), we have introduced an improvement of a GM parameterization by relaxing its key assumption of time locality. The improved parameterization reproduces well the transient behavior of the eddy-resolving gyre model (Fig. 4a).

Our theoretical analysis of the proposed parameterization reveals that the eddy memory leads to an emergence of a decadal variability mode that has a period (approximately 50 yr for the Beaufort Gyre). Despite the EM mode operating on multidecadal time scales, it increases the overall isopycnal slope variance by a fraction of *γ*/*T*_{E} ≈ 0.5 that stays relatively constant for a wide range of mean forcing (Fig. 4b). This suggests that in eddy-dominated flows there might be an inverse relation between eddy memory *γ* and eddy diffusivity (since ).

Note that in this manuscript we have identified the bulk memory of the current as it relates to the cumulative thickness transport of the eddy field. Nonetheless, specific dynamics of individual eddies that can lead to a current having a memory remain unclear. We expect the extent of memory to depend not only on eddy growth rates but also on eddy dissipation rates and on the intensity of the inverse energy cascade. These processes can suppress the eddy transport and affect not only the eddy memory but also the eddy diffusivity.

To emphasize the role of eddy memory, we have made several simplifications. We have used a memory as a parameter characteristic of the entire current. However, the spatial heterogeneity of the eddy diffusivity implies that eddy memory might also be spatially variable. Indeed, Fig. 3b demonstrates that the memory is significantly enhanced over the continental slope, a region with weakened eddy diffusivity. The halocline evolution [Eq. (12)] is valid for a general case of spatially dependent eddy memory, but its analytical treatment is too convoluted to highlight the essential dynamics. Instead, we have simulated the evolution of Eqs. (12) and (13) with an enhanced memory near the coast and confirmed that our key conclusions still hold (not shown). Note that the continental slope occupies only a small portion of the gyre (about 100 km wide); however, the bulk memory that has been diagnosed from FWC evolution is close to the local memory on the continental slope. This implies that the enhanced eddy memory, even in localized regions, impacts the interior gyre dynamics.

Another potentially important factor that was omitted in our theory is vertical diffusion (a diabatic process). Mixing is likely to be important over the continental slope at the boundary where there are freshwater sources. In the Arctic Ocean vertical diffusivity estimates are small *O*(10^{−6} to 10^{−5}) m^{2} s^{−1}, but this weak mixing might still play a role in water mass transformations. In particular, enhanced mixing near gyre boundaries can restrict the availability of freshwater sources and thus limit the temporal variations of the FWC. In our gyre simulations we have used a relatively high estimate of vertical diffusivity (10^{−5} m^{2} s^{−1}), and we have confirmed that the amplitude of the overshoot due to the eddy memory mode increases by a factor of 2 when the vertical diffusivity is reduced to 10^{−6} m^{2} s^{−1} (see Fig. 6). The reasons for such a dramatic suppression of the EM-mode amplitude by vertical mixing remain unclear.

We note that a low-frequency energy enhancement resembling the EM mode has been reported for other turbulent atmospheric and oceanic flows. For example, Thompson and Barnes (2014) highlight a 20- to 30-day periodic variability in the large-scale Southern Hemisphere atmospheric circulation. They suggest it arises because the time rate of change of the eddy heat flux is proportional to the baroclinicity of the flow. Using the similarities between our mathematical formulations we can infer that for the atmospheric flow *γ* ≈ 2 days is of the order of an inverse of the Eady growth rate. In addition, Sinha and Abernathey (2016) demonstrate that the ACC has a maximum response to external forcing at a period of 4 yr. They interpret the ACC behavior from an energetic perspective that also bears mathematical resemblance with our eddy transport description. Assuming the EM mode is pertinent, the 4-yr time scale implies that *γ* ≈ 2 months for ACC [note that from Eq. (18) assuming *γ* = 0.5*T*_{E}]. Consistent with our results, the inferred ACC memory is close to a time lag between EKE and APE evolutions [see Fig. 12 in Sinha and Abernathey (2016)]. Given a zonal current of about 0.1 m s^{−1}, the eddy field with a memory of 3 months can influence ACC dynamics for about 800 km downstream.

The climate modeling community is constantly seeking to improve predictions of the mean climate. However, it is just as important, especially in the context of recent climate change, to simulate and understand low-frequency climate variability, which is largely dictated by ocean dynamics. We have demonstrated here that mesoscale eddies can provide yet another mechanism of low-frequency variability for baroclinic currents. This effect may be amplified by potential feedbacks that involve atmospheric buoyancy fluxes that are in many cases coupled with the ocean dynamics. We thus argue that the implementation of eddy parameterizations that account for eddy memory and an assessment of their implications for coupled climate dynamics is a necessary step forward in climate modeling.

## Acknowledgments

GEM acknowledges the Stanback Postdoctoral Fellowship Fund at Caltech and the Howland Postdoctoral Program Fund at WHOI. MAS was supported by NSF Grants PLR-1415489 and OCE-1232389. AFT acknowledges support from NSF OCE-1235488. GEM thanks Glenn Flierl for discussions during the 2014 Geophysical Fluid Dynamics Summer School at WHOI. The authors acknowledge the high-performance computing support from Yellowstone provided by the NCAR CIS Laboratory, sponsored by the NSF. This work used the Extreme Science and Engineering Discovery Environment (XSEDE; Towns et al. 2014), which is supported by NSF Grant Number ACI-1053575.The manuscript benefited from discussions at the annual Forum for Arctic Modeling and Observing Synthesis (FAMOS) funded by the NSF OPP Awards PLR-1313614 and PLR-1203720.

## REFERENCES

## Footnotes

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