# A Three-way Balance in The Beaufort Gyre: The Ice-Ocean Governor, Wind Stress, and Eddy Diffusivity

Edward W. Doddridge, Gianluca Meneghello, John Marshall, Jon Scott,
Camille Lique

► To cite this version:

Edward W. Doddridge, Gianluca Meneghello, John Marshall, Jon Scott, Camille Lique. A Three-way Balance in The Beaufort Gyre: The Ice-Ocean Governor, Wind Stress, and Eddy Diffusivity. Journal of Geophysical Research. Oceans, 2019, 124 (5), pp.3107-3124. 10.1029/2018JC014897. hal-04202278

HAL Id: hal-04202278

https://hal.science/hal-04202278v1

Submitted on 6 Oct 2023

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

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

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

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

JGR Oceans

RESEARCH ARTICLE

10.1029/2018JC014897

**Special Section:**

Forum for Arctic Modeling and
Observational Synthesis (FAMOS)

2: Beaufort Gyre phenomenon

**Key Points:**

* We combine eddies and the Ice-Ocean Governor into a single theory of Beaufort Gyre equilibration

* Our analytical theory can predict equilibria and adjustment timescales

* The equilibrium state of the Beaufort Gyre is set by a three-way balance between the Ice-Ocean Governor, wind stress, and eddy diffusivity

**Correspondence to:**

E. W. Doddridge,
ewd@mit.edu

**Citation:**

Doddridge E. W., Meneghello, G., Marshall, J., Scott, J., & Lique, C. (2019). A Three-way balance in the Beaufort Gyre: The Ice-Ocean Governor, wind stress, and eddy diffusivity. *Journal of Geophysical Research: Oceans, 124*, 3107–3124.
https://doi.org/10.1029/2018JC014897

Received 3 JAN 2019
Accepted 27 MAR 2019
Accepted article online 1 APR 2019
Published online 16 MAY 2019

# A Three-Way Balance in the Beaufort Gyre: The Ice-Ocean Governor, Wind Stress, and Eddy Diffusivity

Edward W. Doddridge¹ ID, Gianluca Meneghello¹ ID, John Marshall¹, Jeffery Scott¹, and Camille Lique² ID

¹Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA, ²Laboratoire d'Océanographie Physique et Spatiale, Université de Bretagne Occidentale, CNRS, IRD, Ifremer, IUEM, Brest, France

**Abstract** The Beaufort Gyre (BG) is a large anticyclonic circulation in the Arctic Ocean. Its strength is directly related to the halocline depth, and therefore also to the storage of freshwater. It has recently been proposed that the equilibrium state of the BG is set by the Ice-Ocean Governor, a negative feedback between surface currents and ice-ocean stress, rather than a balance between lateral mesoscale eddy fluxes and surface Ekman pumping. However, mesoscale eddies are present in the Arctic Ocean; it is therefore important to extend the Ice-Ocean Governor theory to include lateral fluxes due to mesoscale eddies. Here, a nonlinear ordinary differential equation is derived that represents the effects of wind stress, the Ice-Ocean Governor, and eddy fluxes. Equilibrium and time-varying solutions to this three-way balance equation are obtained and shown to closely match the output from a hierarchy of numerical simulations, indicating that the analytical model represents the processes controlling BG equilibration. The equilibration timescale derived from this three-way balance is faster than the eddy equilibration timescale and slower than the Ice-Ocean Governor equilibration timescales for most values of eddy diffusivity. The sensitivity of the BG equilibrium depth to changes in eddy diffusivity and the presence of the Ice-Ocean Governor is also explored. These results show that predicting the response of the BG to changing surface forcing and sea ice conditions requires faithfully capturing the three-way balance between the Ice-Ocean Governor, wind stress, and eddy fluxes.

## 1. Introduction

The Beaufort Gyre (BG) is a large anticyclonic surface circulation in the Arctic Ocean. In the Canada Basin, surface stress drives Ekman convergence, which leads to the accumulation of low-salinity surface waters in the BG (Proshutinsky et al., 2009, and references therein; Figures 1a and 1b). This low-salinity water sits above a layer of warmer saltier water from the Atlantic and isolates the sea ice from the warmth below (Davis et al., 2016). When this low-salinity surface water is exported to the subpolar North Atlantic, it can potentially affect regional and global climate by influencing the rate of deepwater formation and, therefore, the strength of the Atlantic Meridional Overturning Circulation (Jahn & Holland, 2013; Proshutinsky et al., 2009; Stouffer et al., 2006).

Previous research hypothesized that the accumulation of freshwater driven by Ekman pumping is balanced by the rectified effect of mesoscale eddies (see, e.g., Davis et al., 2014; Manucharyan & Spall, 2016; Manucharyan et al., 2016). This can be diagnosed as a balance between two gyre-scale overturning circulations: the wind-induced circulation, which gathers surface waters toward the center of the gyre steepening isopycnals, and an eddy-induced circulation, which acts to flatten these isopycnals (Doddridge et al., 2016; Marshall & Radko, 2003). The eddy-induced overturning circulation opposes the wind-driven overturning circulation, thereby halting the accumulation of freshwater and equilibrating the BG. However, an alternative mechanism for equilibrating the BG has recently been proposed: the Ice-Ocean Governor (Meneghello, Marshall, Timmermans, et al., 2018; Meneghello, Marshall, Campin, et al., 2018). The Ice-Ocean Governor relies on the fact that the relative velocity between the sea ice and the ocean controls the surface stress at the ice-ocean interface; if the ocean and the ice are moving at the same speed, there is no transfer of momentum, and the gyre is equilibrated. This mechanism is dubbed the Ice-Ocean Governor by analogy with mechanical governors that regulate the speed of engines and other devices through dynamical feedbacks (see, e.g., Maxwell, 1867). The ocean velocity has long been included when calculating the ice-ocean stress in numerical models (see, e.g., Hibler, 1979), but appreciation of the importance of this effect has been very recent
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) The Beaufort Gyre can be clearly seen in the accumulation of fresh water content (colors) and the doming of dynamic ocean topography (contours). Colors show mean fresh water content for the years 2003 to 2013 (inclusive), relative to a reference salinity of 34.8, from an extension of the data set described in Proshutinsky et al. (2009), while contours show mean dynamic ocean topography between 2003 and 2014 from Armitage et al. (2016). (b) Salinity at 75°N from the World Ocean Atlas climatology (Zweng et al., 2018). (c) Schematic of the three-way balance: wind stress (blue arrow) and the Ice-Ocean Governor (orange double-headed arrow) contribute to Ekman pumping, and the residual between these two is balanced by eddy fluxes and diapycnal mixing (horizontal and vertical squiggly gray arrows, respectively). The red and purple segments indicate the radius of the gyre, *R*, and the thickness of the stratified halocline, *δ*, respectively.</figcaption></figure>

(Dewey et al., 2018; Kwok & Morison, 2017; Meneghello et al., 2017; Meneghello, Marshall, Timmermans, et al., 2018; Meneghello, Marshall, Campin, et al., 2018; Zhong et al., 2018) following the publication of a data set that estimates sea surface height and surface geostrophic velocity in sea ice-covered areas (Armitage et al., 2016).

As shown by Meneghello, Marshall, Campin, et al. (2018), the Ice-Ocean Governor predicts that the BG adjusts to step changes in surface forcing with an equilibration timescale of 1–2 years under current conditions. This is noticeably faster than the decadal adjustment timescales estimated when eddy diffusivity controls the equilibration in both idealized models (Davis et al., 2014; Lique et al., 2015; Manucharyan & Spall, 2016; Manucharyan et al., 2016) and general circulation models (Condron et al., 2009; Stewart & Haine, 2013). The equilibration timescale is a crucial parameter for understanding the dynamics of the BG since it determines the length of time over which the gyre retains the imprint of past surface forcing (see, e.g., Johnson et al., 2018; Marshall et al., 2017). This memory of past forcing has important implications for the storage and release of the low-salinity surface waters held in the BG.

The Ice-Ocean Governor does not require mesoscale eddies to equilibrate the BG. However, mesoscale eddies are commonly found in the Arctic Ocean (Dmitrenko et al., 2008; Hunkins, 1974; Manley & Hunkins, 1985;
Meneghello et al., 2017; Timmermans et al., 2008; Zhao et al., 2014, 2016). It is therefore important to extend the theory of the Ice-Ocean Governor presented by Meneghello, Marshall, Campin, et al. (2018) to include the effect of eddy diffusivity. A schematic of these three processes, wind stress, the Ice-Ocean Governor, and eddy fluxes, is shown in Figure 1c.

In the following section we derive an equation for the halocline depth anomaly across the BG that includes both mesoscale eddies and the Ice-Ocean Governor. We then present analytical solutions for the equilibria and the temporal evolution of this equation. In section 3 we present solutions from highly idealized numerical simulations that support the predictions of our analytical theory; increasing the eddy diffusivity shortens the adjustment timescale of our idealized BG. The inclusion of eddy diffusivity also reduces the halocline depth anomaly across the gyre. In section 4 we present an analysis of a coupled ocean-sea ice general circulation model of the Arctic and show that, consistent with our theory, the Ice-Ocean Governor reduces the sensitivity of the equilibrium gyre depth to the parameterized eddy diffusivity. We then present our conclusions in section 5.

## 2. Governing Equation

If we consider an idealization of the BG on an f plane in which the effect of bathymetry is ignored, then we can construct an analytical equation for the depth of an isohaline just below the fresh surface layer, $h_{halo}$. The temporal evolution of $h_{halo}$ is governed by

$$ \frac{\partial h_{halo}}{\partial t} = \nabla \times \frac{\tau}{\rho f} - \kappa \nabla^2 h_{halo}, \quad (1) $$

where $\tau$ represents the surface stress, $\rho$ is the background oceanic density, $f$ is the Coriolis parameter, and $\kappa$ is the eddy diffusivity. The second term on the right-hand side represents mesoscale eddies and vertical diffusivity acting to flatten density surfaces, such that

$$ \kappa = \kappa_h + \frac{R^2 \kappa_v}{\delta^2}, \quad (2) $$

where $\kappa_h$ is the isopycnal thickness diffusivity, $R$ is a representative horizontal length scale, $\kappa_v$ is the diapycnal diffusivity, and $\delta$ is a representative vertical length scale for the halocline. Multiplying $\kappa_v$ by $R^2/\delta^2$ converts the vertical diffusivity into an equivalent horizontal diffusivity. Equation (1) is equivalent to the nonlinear version of equation (22) from Manucharyan et al. (2016). The surface stress term can be decomposed into ice-driven and wind-driven components (Meneghello, Marshall, Timmermans, et al., 2018), written as

$$ \tau = \alpha \rho C_{Di} |\mathbf{u}_i - \mathbf{u}_s|(\mathbf{u}_i - \mathbf{u}_s) + (1 - \alpha) \rho_a C_{Da} |\mathbf{u}_a|(\mathbf{u}_a), \quad (3) $$

in which $\alpha$ is the sea ice concentration; $C_{Di}$ is the ice-ocean drag coefficient; $\mathbf{u}_i$ is the ice velocity; $\mathbf{u}_s$ is the surface ocean velocity, which means that $\mathbf{u}_i - \mathbf{u}_s$ is the ice-ocean relative velocity, $\rho_a$ is the atmospheric density, $C_{Da}$ is the atmosphere-ocean drag coefficient, and $\mathbf{u}_a$ is the surface wind velocity. Because the wind velocity is much larger than the ocean velocity, equation (3) neglects the ocean velocity when computing the air-ocean stress.

We will henceforth treat $\mathbf{u}_i$ and $\mathbf{u}_a$ as if they are independent parameters. In reality the sea ice velocity, $\mathbf{u}_i$, and the wind velocity, $\mathbf{u}_a$, are linked since the wind imparts momentum to the sea ice. However, the exact relationship between the wind velocity and the ice velocity is complex: the relationship depends on the internal stresses in the ice field, the ocean velocity, the geometry of the surrounding coastlines, and the timescale over which the velocities are compared (Cole et al., 2017; Hibler, 1979; Martin et al., 2014). It is therefore important that the chosen parameter values be physically consistent.

To ensure that our equation remains tractable when examining the balance between ice-mediated Ekman pumping and mesoscale eddies, we consider a highly idealized configuration when constructing our analytical model, that is nevertheless similar to configurations employed in the literature (see, e.g., Davis et al., 2014; Yang et al., 2016). We will assume an f plane, a perfectly circular gyre, and define the depth anomaly, $h$, as the difference between the isohaline depth, $h_{halo}$, at the edge of the gyre and at the center of the gyre. While our configuration is similar to others that have been used in the literature, it nevertheless ignores processes that may affect the BG. These include interactions between the surface gyre and the flow in the
Atlantic Water layer (cf. Lique et al., 2015), interactions between the surface flow and the bathymetry (Nøst et al., 2008; Regan et al., 2019; Yang et al., 2016), and the possibility for surface buoyancy forcing to modify the circulation (Howard et al., 2015). It should also be remembered that the surface forcing experienced by the BG varies on a wide range of timescales and length scales; as well as synoptic variability, there is a large seasonal cycle in sea ice concentration (Simmonds, 2015), and interannual variability in both the winds and the ice cover (Simmonds, 2015; Thompson & Wallace, 1998). We do not represent this variability in our analytical solution. Instead, we use parameter values that represent long-term means and leave the impact of short term variability and the seasonal cycle for a future contribution.

The absolute value in equation (3) makes finding a general analytical solution difficult. However, we can begin by considering two cases: if the circulation of the BG is slower than the ice speed, then $u_i - u_s > 0$ and the ice is a source of momentum for the ocean, while $u_i - u_s < 0$ means that the ice is removing momentum from the oceanic circulation. We therefore analyze solutions for both $u_i \ge u_s$ and $u_i \le u_s$. The scenario in which $u_i = u_s$ can be derived from either case. These two cases are also discussed by Dewey et al. (2018), albeit with a focus on synoptic variability in the atmosphere.

Following Meneghello, Marshall, Campin, et al. (2018), we will henceforth approximate the surface ocean speed, $u_s$, with the geostrophic speed of the gyre, $u_g = g'h/fR$ in which $g'$ is the reduced gravity across the halocline and $R$ is the radius of the gyre. Using the above idealizations means that equation (1) can be written as

$$ \frac{1}{\xi} \frac{dh}{dt} = \frac{\alpha C_{Di} |u_i - \frac{g'h}{fR}| (u_i - \frac{g'h}{fR})}{fR} + \frac{(1-\alpha) C_{Da} u_a^2 \rho_a}{fR\rho} - \frac{\kappa}{R^2} h, \quad (4) $$

where $\xi \ge 1$ is a scaling factor that converts between a volume flux and the isohaline depth anomaly and depends on the initial shape of the isohaline (see Appendix A), $u_i$ is the ice speed, and $u_a$ is the wind speed. We assume that both $u_i$ and $u_a$ represent the rotational component of the ice and wind velocities averaged over the gyre region. The derivatives have been simplified by evaluating radial derivatives as the difference between the center and the edge of the gyre divided by the radius of the gyre, while azimuthal derivatives have been assumed to be zero. Equation (4) highlights the titular three-way balance; to reach equilibrium the three terms on the right hand side must sum to zero. At that point there is a three-way balance between the sea ice-ocean stress modulated by the Ice-Ocean Governor, the wind stress on the ice-free portions of the gyre, and eddy diffusivity.

## 2.1. Equilibrium Solutions

We now begin our analysis of equation (4) by seeking equilibrium solutions, $h_{eq}$, and determining their stability. To obtain equilibrium solutions we set $dh/dt = 0$ and solve for $h_{eq}$.

### 2.1.1. No Sea Ice Cover: $\alpha = 0$

In the absence of sea ice the surface stress does not depend on the state of the underlying ocean; because the ocean velocity is negligible in comparison to the wind velocity, the momentum forcing is completely determined by the imposed forcing fields. It is worth noting that this scenario is equivalent to neglecting the ocean geostrophic velocity in the calculation of the ice-ocean stress since it removes the feedback between the state of the gyre and the strength of the surface forcing. Both of these scenarios, removing the ice or neglecting the ocean velocity when calculating the ice-ocean stress, explicitly exclude any influence of the Ice-Ocean Governor on the solution. We therefore expect the gyre to equilibrate via eddy diffusivity. With $\alpha = 0$, equation (4) becomes a linear ordinary differential equation. As such, the solutions are straightforward. The equilibrium solution is

$$ h_{eq} = \frac{RC_{Da}u_a^2\rho_a}{\kappa f \rho} \equiv \frac{R\tau_a}{\kappa f \rho}, \quad (5) $$

which corresponds to a balance between wind-driven Ekman pumping and eddy fluxes (cf. Davis et al., 2014; Manucharyan & Spall, 2016; Marshall & Radko, 2003). The predicted isopycnal depth anomaly, $h_{eq}$, tends to infinity as the eddy diffusivity, $\kappa$, tends to 0; without eddy diffusivity or the Ice-Ocean Governor there is no mechanism through which the gyre can equilibrate.

If we follow Manucharyan and Spall (2016) and assume that $\kappa$ is proportional to the square of the isopycnal slope, $\kappa = k(h_{eq}/R)^2$, where $k$ is a constant eddy efficiency coefficient, then we find that the isopycnal depth
anomaly is given by

$$h_{eq} = R \left( \frac{\tau_a}{k f \rho} \right)^{1/3}. \qquad (6)$$

Since the freshwater content (FWC) is proportional to the isopycnal depth anomaly in idealized models in which the salinity of the surface layer does not vary (Meneghello, Marshall, Campin, et al., 2018), equation (6) is consistent with Manucharyan and Spall (2016), who show that, without the Ice-Ocean Governor, the FWC of their idealized BG is proportional to the cube root of the surface stress.

### 2.1.2. With Sea Ice: $\alpha \neq 0$

The sea ice concentration that we use represents a long-term average across a substantial seasonal cycle. Despite declining sea ice concentrations and thicknesses in the Arctic, the wintertime sea ice coverage remains substantial (Simmonds, 2015). Sea ice is generally considered to be in free drift at concentrations below 0.8 (Martin et al., 2014). Despite this, we expect the Ice-Ocean Governor to be important even when the annual mean $\alpha$ is less than 0.8, provided that the wintertime ice is not in free drift.

Retaining the sea ice concentration, $\alpha$, in our analytical equation for isopycnal depth anomaly means that the equation is a nonlinear ordinary differential equation. The equilibrium solution to equation (4) with $u_i \ge u_g$ is

$$h_{eq} = \frac{u_i f R}{g'} + \frac{f^3 R}{2 \alpha g'^2 C_{Di}} \left( \kappa \pm \sqrt{\kappa \left( \kappa + \frac{4 \alpha C_{Di} u_i g'}{f^2} \right) - \frac{4 \alpha (1 - \alpha) C_{Da} u_a^2 \rho_a C_{Di} g'^2}{f^4 \rho}} \right), \quad (7)$$

which collapses to the solution predicted by Meneghello, Marshall, Campin, et al. (2018) when $\alpha = 1$ and $\kappa = 0$.

In this limiting case, with $\alpha = 1$ and $\kappa = 0$, we can greatly simplify equation (7) to become

$$h_{eq} = \frac{u_i f R}{g'}. \qquad (8)$$

This can be rearranged to give $u_i = h_{eq}g'/fR$, which shows that the ice velocity is equal to the geostrophic velocity, $u_g$. Under the assumption of total ice cover, $\alpha = 1$, and no eddy diffusivity, $\kappa = 0$, the equilibrium is reached when the geostrophic velocity and ice velocity are equal. At this point the Ice-Ocean Governor has reduced the surface stress to zero and completely switched off the Ekman pumping.

Using realistic parameter values (Table 1) and $\alpha = 1$, we obtain the solution for equilibrium isopycnal depth anomaly shown in Figure 2, which plots the equilibrium depth anomaly from equation (7) as a function of $\kappa$ and highlights the presence of two solution branches for $\kappa > 0$, as well as a saddle-node bifurcation at $\kappa = 0$. The unstable upper branch is not physical. The solution shown in Figure 2 is the special case in which $\alpha = 1$, and hence, the saddle-node bifurcation occurs at $\kappa = 0$ m²/s; in general, the bifurcation will occur at a positive, nonzero, value of $\kappa$.

If we set $u_i \le u_g$, then the equilibrium solution is given by

$$h_{eq} = \frac{u_i f R}{g'} - \frac{f^3 R}{2 \alpha g'^2 C_{Di}} \left( \kappa \mp \sqrt{\kappa \left( \kappa - \frac{4 \alpha C_{Di} u_i g'}{f^2} \right) + \frac{4 \alpha (1 - \alpha) C_{Da} u_a^2 \rho_a C_{Di} g'^2}{f^4 \rho}} \right), \quad (9)$$

which has the same form as the solution in equation (7), but the signs after the first term are swapped. This solution is only possible when $\alpha < 1$, since it requires the wind stress over the ice-free areas to accelerate the gyre beyond the ice speed.

The presence of a square root in equations (7) and (9) leads to a saddle-node bifurcation in both solutions. One consequence of these saddle-node bifurcations is that it is possible to find solutions for both $u_i > u_g$ and $u_i < u_g$ over a range of parameter values. However, in this overlap region one of the solutions violates the assumed sign of $u_i - u_g$ (see Appendix B for further discussion of the saddle-node bifurcations). Therefore, although there are two mathematically stable solutions, only one of them is physically valid.

Our discussion so far highlights the existence of two possible regimes. In the first regime the ice speed is greater than the current speed, which means that the ice imparts momentum into the ocean. In this
<table><caption>**Table 1**
Default Parameter Values Used to Represent the Beaufort Gyre</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>f</td><td>Coriolis parameter</td><td>1.4 × 10<sup>-4</sup></td><td>s<sup>-1</sup></td></tr><tr><td>α</td><td>Ice fraction</td><td>0.87</td><td></td></tr><tr><td>C<sub>Di</sub></td><td>Ice-ocean drag</td><td>5.5 × 10<sup>-3</sup></td><td></td></tr><tr><td>u<sub>i</sub></td><td>Ice speed</td><td>8</td><td>cm/s</td></tr><tr><td>g'</td><td>Reduced gravity</td><td>6.22 × 10<sup>-2</sup></td><td>m/s<sup>2</sup></td></tr><tr><td>ρ</td><td>Ocean density</td><td>1,028</td><td>kg/m<sup>3</sup></td></tr><tr><td>C<sub>Da</sub></td><td>Atmosphere-ocean drag</td><td>1.25 × 10<sup>-3</sup></td><td></td></tr><tr><td>u<sub>a</sub></td><td>Wind speed</td><td>4</td><td>m/s</td></tr><tr><td>ρ<sub>a</sub></td><td>Atmospheric density</td><td>1.25</td><td>kg/m<sup>3</sup></td></tr><tr><td>R</td><td>Radius of the gyre</td><td>340</td><td>km</td></tr><tr><td>κ</td><td>Eddy diffusivity</td><td>300</td><td>m<sup>2</sup>/s</td></tr></tbody></table>

scenario to Ice-Ocean Governor acts to reduce the surface stress received from the ice and the eddy diffusivity balances the sum of the ice-ocean stress and the air-ocean stress. In the second regime the air-ocean stress over the ice-free sections of the gyre accelerates the geostrophic speed beyond the ice speed such that both the ice-ocean stress and the eddy diffusivity work together to balance the air-ocean stress. Between the two regimes there is a point at which the ice speed and the current speed are equal. In this scenario the eddy diffusivity balances the air-ocean stress, and the Ice-Ocean Governor has equilibrated the gyre such that there is no net transfer of momentum from the ice to the ocean. While the second regime, in which $u_i < u_g$, is likely important over the seasonal cycle and during the transient response to changing surface forcings, it is an implausible equilibrium state since the large fraction of open water required for the wind stress to directly accelerate the gyre means that the sea ice is likely to be in free drift (Martin et al., 2014), and the Ice-Ocean Governor will therefore be ineffective. Dewey et al. (2018) present a discussion of these regimes in the BG. Here we extend the concept to include a rigorous dynamical framework and explore the implications of this framework for the BG.

For the parameter values that represent a reasonable approximation of the BG (Table 1), we find that the regime change occurs at $κ ≈ 125 m^2/s$, which is at the low end of recent eddy diffusivity estimates of 300 ± 200 m^2/s from idealized modeling (Manucharyan et al., 2016) and 100–600 m^2/s from observations (Meneghello et al., 2017). This suggests that the ice cover is currently a net source of momentum for the BG.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Black lines show the equilibrium solutions to equation (4), for the case in which $u_i > u_g$, given by equation (7). The lower branch (solid) is stable, while the upper branch (dashed) is unstable and nonphysical. The background colors represent $dh/dt$. In the red region $dh/dt$ is positive, and so a trajectory started below the solid line will move upward with time, toward larger values of $h$. In the blue region $dh/dt$ is negative, so a trajectory started above the solid line will move downward toward smaller values of $h$. Therefore, the lower branch (solid black line) represents stable equilibrium solutions, and the upper branch (dashed black line) represents unstable equilibrium solutions. The parameters are set to the values shown in Table 1, except for $α$ which is set to 1.</figcaption></figure>
## 2.2. Time-varying Nonlinear Solutions

We now seek time-varying solutions to equation (4). Since equation (4) is a nonlinear ordinary differential equation, the solution is not straightforward. With assistance from symbolic algebra software (Wolfram-Alpha, 2018) we arrive at the following solution for $h(t)$ when $u_i \ge u_g$

$$h(t) = h_{eq-} + \frac{f^3 R^3 k_1}{t_{3wb} \alpha g'^2 C_{Di} (e^{\xi t/t_{3wb}} + k_1)}, \quad (10)$$

where $h_{eq-}$ is the negative branch of the equilibrium solutions shown in equation (7),

$$t_{3wb} = \frac{R^2}{\sqrt{\kappa \left( \kappa + \frac{4\alpha C_{Di} u_i g'}{f^2} \right) - \frac{4\alpha(1-\alpha)C_{Di}C_{Da}u_g^2\rho_a g'^2}{f^4\rho}}}, \quad (11)$$

is the equilibration timescale from our three-way balance theory, $k_1$ is an arbitrary constant of integration whose value is defined by the initial condition, and the rest of the symbols are as previously defined. If we set $h(0) = 0$, meaning that we begin from a state of rest with a flat halocline, then $k_1$ is given by

$$k_1 = \frac{-h_{eq-} \alpha g'^2 C_{Di}}{h_{eq-} \alpha g'^2 C_{Di} + f^3 R \sqrt{\kappa \left( \kappa + \frac{4\alpha C_{Di} u_i g'}{f^2} \right) - \frac{4\alpha(1-\alpha)C_{Di}C_{Da}u_g^2\rho_a g'^2}{f^4\rho}}}. \quad (12)$$

In this scenario the value of $k_1$ is constrained to fall within the open interval $(-1, 0)$. Since $\alpha, \kappa = 0$ fundamentally change the structure of equation (4), this time-varying solution is only valid for $0 < \alpha \le 1$ and $\kappa > 0$. For solutions in the absence of eddy diffusivity we refer readers to Meneghello, Marshall, Campin, et al. (2018), and to Manucharyan and Spall (2016) for eddying solutions in the absence of the Ice-Ocean Governor. From equation (10) it can be seen that as $t \to \infty$ the time-varying solution asymptotes to $h_{eq-}$ regardless of the value of $k_1$. Therefore, the asymptotic value is independent of the initial conditions, as expected. Furthermore, the scaling factor $\xi$ appears only in the exponent of the solution, validating our ansatz that $\xi$ is related to the initial shape of the isopycnal layer and that it does not affect the equilibrium solution.

One of the consequences of the nonlinear nature of equation (4) is that the timescale in the solution depends on the forcing, via $u_i$ and $u_g$, and the flow. We are therefore unable to define a single characteristic timescale for the equilibration of the BG. This stands in stark contrast to the eddy equilibrated gyre proposed by Davis et al. (2014), which equilibrates with a characteristic timescale $t_{eddy} = R^2/\kappa$ regardless of the strength of the forcing. It is possible to obtain an adjustment timescale for an eddy equilibrated gyre that depends on the strength of the forcing if the eddy diffusivity is assumed to depend on the isopycnal slope. For example, Manucharyan and Spall (2016) find that their eddy equilibrated gyre adjusts faster when the surface forcing is stronger.

We obtain two equilibration timescales for our analytical solution: one for each sign of $u_i - u_g$. The equilibration timescale for our analytical model is never longer than $t_{eddy}$, provided that we restrict our attention to the valid solutions (Figure 3a, solid blue and orange lines). We find that our equilibration timescale increases rapidly in the vicinity of the regime change. At the exact diffusivity associated with the regime change, the equilibration timescale reaches a maximum and is the same as the eddy equilibration timescale. The equilibration timescale for the pure Ice-Ocean Governor solution (Meneghello, Marshall, Campin, et al., 2018) does not depend on, or indeed consider, the eddy diffusivity. As such, it is a horizontal line in Figure 3a. The equilibration timescale for our three-way balance theory is longer than the pure Ice-Ocean Governor timescale for eddy diffusivities below approximately 500 m²/s, which indicates that the value computed by Meneghello, Marshall, Campin, et al. (2018) is likely a lower limit for the equilibration timescale of the BG. It is also possible to consider the equilibration timescale as a function of the ice speed. With an eddy diffusivity of $\kappa = 300$, the equilibration timescale shows a distinct peak at $u_i \approx 3$ cm s⁻¹ (Figure 3b), which is where the regime shift from $u_i < u_g$ to $u_i > u_g$ occurs. Our three-way balance theory predicts that the equilibration timescale is a nontrivial function of the forcing and the flow state. This implies that the extent to which the BG remembers past forcing depends on the current and historical state of the gyre as well as the forcing history (cf. Johnson et al., 2018; Marshall et al., 2017).

In the following sections we use a range of numerical simulations to explore the effect of eddy diffusivity and ice velocity on the equilibrium thickness anomaly and equilibration timescale.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Equilibration timescales using parameter values from Table 1 as a function of (a) $\kappa$, and (b) $u_i$ for each sign of $u_i - u_g$ ($> 0$ in blue and $< 0$ in orange), a purely eddy equilibrated gyre (green line), and a gyre equilibrated purely by the Ice-Ocean Governor (cyan line). At the exact diffusivity or ice speed associated with the regime shift the equilibration time of our three-way balance solution is equal to the eddy equilibration timescale. The solid blue and orange lines show the valid solutions, while the dashed lines show the solutions which violate the assumed sign of $u_i - u_g$.</figcaption></figure>

### 3. Idealized Simulations

To explore the impact of the three-way balance between the Ice-Ocean Governor, wind stress, and eddy diffusivity on the evolution and equilibrium of the BG, we require a model in which eddy diffusivity is a parameter, not an emergent property of the flow. Because the equations of motion in 1.5 layer models explicitly preclude baroclinic instability, they are ideally suited to this task. We parametrize the effects of mesoscale eddies as a down-gradient thickness flux using a Laplacian operator and a fixed thickness diffusivity. For these idealized simulations we assume complete ice cover, setting $\alpha = 1$. The surface stress is then modeled as a quadratic drag from an imposed ice-drift field identical to the one used by Meneghello, Marshall, Campin, et al. (2018) in which the ice velocity increases with distance from the center of the gyre up to some maximum value and then decreases to 0 near the edge of the domain. The simulations are labeled by the maximum ice speed in the surface forcing; the gyre-averaged ice speed is approximately half the maximum value. The simulations use a Cartesian grid with 4-km horizontal resolution. The vertical resolution of our model is highly idealized; it consists of a single active layer above an infinitely deep quiescent abyss. Numerical models that make this simplification are known as “reduced gravity” models and have been used to explore a wide range of dynamics (see, e.g., Davis et al., 2014; Février et al., 2007; Johnson & Marshall, 2002; Stern, 1998). We use Aronnax v0.2.0 (Doddridge & Radul, 2018) to perform the idealized simulations discussed here and the configuration scripts can be found at https://doddridge.me/publications/DMMSL2019/. The equilibrium state of these simulations is a circular anticyclonic gyre, with layer thickness increasing toward the center of the gyre (Figure 4a).

To test our hypotheses regarding the parameter sensitivity of the BG, we run an ensemble of simulations in which we vary both ice speed and eddy diffusivity systematically. We begin by attempting to replicate the eddy-resolving MITgcm based results of Meneghello, Marshall, Campin, et al. (2018) with our reduced gravity model. We run an ensemble of simulations in which the ice speed is varied while the eddy diffusivity is set to 0. This provides a means to determine whether our highly idealized numerical model is appropriate.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Eulerian-mean meridional velocity (colors) and layer thickness (contours, with a contour interval of 10 m) from the final 10 years of the $\kappa = 0 \text{ m}^2 \text{s}^{-1}$ $u_i = 8 \text{ cm/s}$ simulation. (b) Temporal evolution of the isopycnal depth anomaly for a range of imposed ice velocities with the thickness diffusivity, $\kappa$, set to 0. (c) The normalized isopycnal depth anomaly shows that the adjustment timescale depends on $u_i$. (d) Equilibrium isopycnal depth anomaly plotted against $u_i$ (dots), with the prediction from equation (7) using the values from Table 1 and $\alpha = 1$ shown by the black line.</figcaption></figure>
We find that the imposed ice velocity has a strong effect on the equilibrium value and the timescale of the initial adjustment (Figure 4). Meneghello, Marshall, Campin, et al. (2018) suggest that the gyre should equilibrate faster as the ice velocity increases; a prediction that is supported by our simulations (Figure 4c). Equation (7) predicts that in the absence of eddies the equilibrium isopycnal depth anomaly should increase linearly as the ice velocity is increased. Our reduced gravity simulations support this prediction (Figure 4d), in accord with the results of Meneghello, Marshall, Campin, et al. (2018).

To explore the effect of eddy diffusivity on the equilibrium value, we now hold the ice velocity fixed and vary the eddy diffusivity; Figure 5a shows the evolution of idealized reduced gravity simulations with an imposed maximum ice velocity of 8 cm/s and a range of thickness diffusivities. As predicted in section 2 both the equilibrium depth anomaly (Figures 5a and 5c) and the equilibration timescale (Figures 5b and 5d) depend on the eddy diffusivity; a larger eddy diffusivity leads to a smaller halocline depth anomaly and a faster equilibration timescale. The theoretical predictions for equilibrium depth anomaly, equation (7), and equilibration timescale, equation (10), closely match the simulations. Despite the substantial simplifications used to derive our analytical expression, it retains the essential physics required to explain the evolution and equilibration of an idealized gyre in a nonlinear ocean model.

## 4. Analysis of a General Circulation Model

The numerical results presented so far rely on highly idealized simulations. Both the domain and the dynamics were simplified in order to provide a clean test of the Ice-Ocean Governor in the presence of (parameterized) eddy diffusivity; the analytical and reduced gravity models rely on many of the same assumptions, including a prescribed ice velocity. In this section we analyze the output of a suite of simulations using a coarse resolution ($\Delta x, \Delta y \approx 36$ km) regional coupled ocean-sea ice general circulation model of the Arctic. These simulations are substantially more complex than the reduced gravity simulations presented in section 3 and include a sea ice model, historical forcing fields, and realistic bathymetry. Because of the added complexity and fidelity, these simulations provide a much more rigorous test of the theoretical predictions presented in section 2. Our ensemble of Arctic simulations do not resolve mesoscale eddies and therefore use the Gent-McWilliams mesoscale eddy parameterization (Gent & McWilliams, 1990; Gent et al., 1995). This allows us to systematically vary the eddy diffusivity while holding all other parameters fixed. The Gent-McWilliams eddy parameterization has been shown to produce more realistic solutions when used with a variable eddy diffusivity coefficient (see, e.g., Danabasoglu & Marshall, 2007; Ferreira et al., 2005; Gent, 2016). However, our analytical solution requires a single value for the eddy diffusivity, so we restrict ourselves to prescribing a single constant value of $\kappa_{GM}$ between 50 and 2000 m$^2$/s for each ensemble member. A more complete description of the model setup can be found in Appendix C.

In keeping with the literature, we define our BG region as a latitude-longitude box from 70.5°N to 80.5°N and 130°W to 170°W excluding grid points where the ocean is shallower than 300 m (Proshutinsky et al., 2009; Figure 6a, red contour). The isohaline depth anomaly across the BG is defined as the difference between the average depth of the 34.5-psu isohaline surface at the edge of the gyre region, shown by the 70-km-wide gap between the white and red contours, and its maximum depth within the gyre region.

The climatology and variability of our MITgcm simulations closely resemble those of the Arctic Ocean. Focusing on the ensemble member with $\kappa_{GM} = 500$ m$^2$/s, we find that the spatial distribution of FWC (Figure 6a), the climatological salinity along 75°N (Figure 6b), and the FWC of the BG region (Figure 6c) all compare favorably with the observations. The model captures not only the mean FWC, but also the inter-annual variability. In particular the large increase in FWC circa 2007 is present in the simulation results. These comparisons indicate that our model produces a very good approximation of the BG. We now seek to assess the predictive power of our analytical theory using these simulations.

In contrast to the pure Ice-Ocean Governor solution, both the three-way balance solution and the eddy-equilibrated solution predict that the isopycnal depth anomaly will decrease as $\kappa$ increases. However, these solutions predict very different relationships between $h_{eq}$ and $\kappa$. Figure 7 shows the diagnosed isohaline depth anomalies from the simulations plotted against $\kappa$. Our analytical prediction requires a single value for each of ice speed, wind speed, the reduced gravity, and ice fraction which we diagnose from ensemble means and imposed forcing fields as 8.7 cm/s, 2.5 m/s, 0.0287 m/s$^2$, and 0.87, respectively. The rest of the parameters remain unchanged from those shown in Table 1. With these parameter values, the regime change is expected to occur at $\kappa \approx 20$ m$^2$/s, suggesting that the ice speed will be faster than the geostrophic
<figure><img src="image_7.png" /></figure>

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

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

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

**Figure 5.** Time evolution of the halocline depth anomaly within the idealized Beaufort Gyre for a range of eddy diffusivities with $u_i$ = 8 cm/s. (a) The isopycnal depth anomaly against time. (b) The normalized isopycnal depth anomaly shows that as eddy diffusivity increases the equilibration timescale shortens. (c) Equilibrium depth anomaly from the simulations plotted against analytical prediction. (d) Equilibration timescale from the theory with $\xi$ = 2.5 plotted against timescales extracted from the simulations.
<figure><img src="image_11.png" /><figcaption>**Figure 6.** (a) Arctic Ocean freshwater content (FWC) for our $\kappa_{GM} = 500 \text{ m}^2/\text{s}$ ensemble member in August averaged over the years 2003 to 2013 (inclusive), relative to a reference salinity of 34.8 psu (colors). The FWC is averaged over the same years as the observational estimate shown in Figure 1a. The red contour shows the Beaufort Gyre region as defined by Proshutinsky et al. (2009). The area between the white and red contours is used to calculate the isohaline depth at the edge of the gyre. (b) Salinity at 75°N from the $\kappa_{GM} = 500 \text{ m}^2/\text{s}$ ensemble member, averaged over the final 25 years of the simulation (1989–2013, inclusive), for comparison with the long-term observational climatology presented in Figure 1b. (c) Monthly mean FWC in the Beaufort Gyre region defined by Proshutinsky et al. (2009; red contour in a) from the $\kappa_{GM} = 500 \text{ m}^2/\text{s}$ ensemble member (blue line) and summertime FWC from observations (an extension of the data set described in Proshutinsky et al., 2009; orange line, error bars represent root-mean-square error estimate.).</figcaption></figure>

current in all of the ensemble members we analyze. Diagnosing the ice and ocean speeds from the simulations confirms this prediction (not shown). Figure 7 shows the simulation data and the three analytical predictions for isohaline depth anomaly in the BG: one for a gyre equilibrated solely by the Ice-Ocean Governor, one equilibrated only by eddy diffusion, and finally, a gyre equilibrated by our three-way balance theory. The three-way balance prediction most closely matches the simulations, while there is a substantial mismatch between the simulations and the pure Ice-Ocean Governor and eddy-equilibrated predictions. The pure Ice-Ocean Governor represents an upper bound: the isohaline depth anomaly across the BG cannot be larger than this value for a given ice speed, ice concentration, and wind speed. Our three-way balance theory also represents an upper bound, albeit a far more stringent upper bound, since the processes that we neglected in its derivation will generally act to remove momentum from the gyre and restrict the accumulation of freshwater, thereby reducing the isohaline depth anomaly. These processes likely include the influence of bathymetry through bottom drag and topographic steering (Nøst & Isachsen, 2003).
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Isopycnal depth anomaly from a series of Arctic simulations. The green circles represent the results from an ensemble of MITgcm simulations in which $\kappa_{GM}$ has been varied, while the lines represent the theoretical predictions. The solid red line represents the three-way balance prediction, the solid orange line is the pure Ice-Ocean Governor prediction ($\kappa = 0 \text{ m}^2/\text{s}$), and the dashed blue line represents the predicted solution for an eddy equilibrated gyre in the absence of the Ice-Ocean Governor. The surface stress for the dashed blue line has been calculated using the ensemble average ice-ocean relative velocity in equation (3). The gray shaded region on the x axis highlights the estimated range of 100–600 $\text{m}^2/\text{s}$ for eddy diffusivity in the BG (Meneghello et al., 2017).</figcaption></figure>

While it is possible to tune our analytical prediction to match the MITgcm simulations, we believe that the difference between the analytical prediction and the simulations contains information about the importance of the processes we have neglected. Our results suggest that these processes reduce the isohaline depth anomaly across the BG by approximately 20 m. It also worth noting that the eddy-only prediction, equation (5), cannot be tuned to match the simulation results, which strongly suggests that the BG in these simulations is not equilibrated solely by the eddy parameterization.

This ensemble of simulations provides strong evidence that the Ice-Ocean Governor is active in this model, even in the presence of substantial (parameterized) eddy diffusivity. The simulations also show that including the Ice-Ocean Governor is crucial for predicting the relationship between eddy diffusivity and the equilibrium isopycnal depth anomaly of the BG.

## 5. Conclusions

We have extended the theory presented by Meneghello, Marshall, Campin, et al. (2018) to describe the evolution of the BG in the presence of both the Ice-Ocean Governor and mesoscale eddies. Our analytical model provides predictions for both the equilibrium state and the adjustment timescale as a function of the wind speed, sea ice speed, sea ice fraction, and eddy diffusivity. In the limiting case with no sea ice (or equivalently, neglecting the ocean velocity when calculating the ice-ocean stress), and therefore no Ice-Ocean Governor, we are able to recover previously published relationships between surface forcing and the state of the BG (Davis et al., 2014; Manucharyan & Spall, 2016). Retaining the effect of the Ice-Ocean Governor in our analytical model substantially alters the predicted equilibrium state and equilibration timescale of the BG. The combination of both eddy diffusion and the Ice-Ocean Governor leads to a shorter equilibration timescale than eddy diffusion acting alone. The equilibration timescale from our three-way balance theory is usually longer than the pure Ice-Ocean Governor timescale predicted by Meneghello, Marshall, Campin, et al. (2018) that neglected the effect of mesoscale eddies. Given the ubiquitous nature of mesoscale eddies in the Arctic Ocean, this suggests that the timescale calculated by Meneghello, Marshall, Campin, et al. (2018) is likely a lower limit for the equilibration timescale of the BG.

We test our analytical theory against output from a hierarchy of models consisting of an ensemble of simulations using a highly idealized 1.5 layer reduced gravity model and an ensemble of simulations using a general circulation model with a coupled sea ice model, and realistic bathymetry and forcing. We find startlingly good agreement between the reduced gravity simulations and the analytical predictions, providing very strong evidence that our three-way balance theory accurately represents the physics controlling the evolution of a circular gyre on an *f* plane. While our theory provides a better estimate for the MITgcm simulations than either the eddy equilibration theory or the pure Ice-Ocean Governor, a nonnegligible mismatch remains. The theory overestimates the isohaline depth anomaly. This is likely due to the processes that were neglected in the derivation of our analytical model, for example, the influence of bathymetry, and spatial and temporal variability of wind and sea ice. Our analytical theory also neglects other circulations in the region of the BG, for example, the lateral flow of Pacific winter water (Zhong et al., 2019).

Our analytical model equilibrates faster when the ice-ocean relative velocity is larger, and when the eddy diffusivity is large. Since the sea ice cover undergoes a substantial seasonal cycle, it is likely that the dependence of the adjustment timescale on the ice-ocean relative velocity influences the seasonal adjustment of the BG. Such an analysis is beyond the scope of this contribution, but represents an exciting avenue for future work.
A concept similar to our three-way balance between the Ice-Ocean Governor, wind stress, and eddy diffusivity in the BG was proposed by Dewey et al. (2018). However, they focus on the effect of short term wind variability and argue that if the sign of the curl of the surface stress is negative approximately half of the time, then the gyre is stable. This assertion would hold if the distribution of the curl of the surface stress were symmetric about zero, but there is no a priori reason to expect such a distribution. Ultimately, it is the integrated momentum transfer that is important, not just the sign of the curl of the surface stress.

Our analytical theory describes the evolution of the isohaline depth anomaly across the BG, which is the dynamically relevant portion of FWC and represents the wind-driven variability in FWC. If the isohaline depth at the edge of the gyre is constant, then changes in the dynamically relevant FWC are equivalent to changes in the total FWC. Most idealized models of the BG force the isohaline depth at the edge of the gyre to remain constant by relaxing the stratification to a set profile (see, e.g., Davis et al., 2014; Manucharyan & Spall, 2016; Yang et al., 2016). This means that changes in total FWC and dynamically relevant FWC are equivalent in these models. The extent to which this equivalence holds in the BG depends on the processes that set the stratification at the edge of the BG and the availability of liquid freshwater.

The concentration and thickness of sea ice affect its dynamics; when the sea ice concentration drops below 80%, the ice is expected to be in free drift (Martin et al., 2014). Ice thickness also affects whether the ice is in free drift, since thin ice is less able to sustain the internal stresses that prevent free drift. The sea ice concentration we consider represents a long-term mean over a substantial seasonal cycle. Therefore, we expect there to be periods during which the sea ice is not in free drift, even for average concentrations well below 80%.

We have presented a rigorous dynamical framework that may be used to predict the response of the BG to changing ice and atmospheric conditions. Our results suggest that both eddy diffusivity and the Ice-Ocean Governor are important in setting the structure of the BG. In particular, our results show that the Ice-Ocean Governor reduces the sensitivity of the equilibrium depth of the BG to eddy diffusivity, especially for the small eddy diffusivities estimated by Manucharyan et al. (2016) and Meneghello et al. (2017). Determining the relative importance of each process in the BG represents a crucial next step in improving our understanding of the circulation in the Arctic Ocean and is the subject of a forthcoming observationally focused companion paper. As Arctic sea ice coverage continues to decline in a warming world, we expect wind stress over the ice-free portions of the gyre and eddy diffusivity to become increasingly important in the three-way balance. However, provided that winter sea ice coverage remains extensive, the Ice-Ocean Governor will continue to be central to the dynamics of the BG for the foreseeable future.

## Appendix A: Temporal Scaling Factor

At equilibrium we expect the average isohaline slope across the gyre to be determined by a three-way balance between the Ice-Ocean Governor, wind stress, and eddy diffusivity. However, the initial shape that the surface layer takes as it is inflated is determined by the distribution of Ekman pumping, not by ocean dynamics. If the layer inflates as an inverted cone, then

$$ \begin{gathered} \frac{d}{dt} \frac{\beta R^2}{3} h = \iint w_{Ek} dA = \overline{w_{Ek}} \beta R^2 \\ \Rightarrow \frac{dh}{dt} = \frac{3}{2} \overline{w_{Ek}}, \end{gathered} \quad (A1) $$

where $\overline{w_{Ek}}$ is the area-averaged Ekman pumping, but, if the isopycnal layer inflates as half an ellipsoid, then

$$ \begin{gathered} \frac{d}{dt} \frac{2\beta R^2}{3} h = \iint w_{Ek} dA = \overline{w_{Ek}} \beta R^2 \\ \Rightarrow \frac{dh}{dt} = \frac{3}{2} \overline{w_{Ek}}. \end{gathered} \quad (A2) $$

By comparing equations (A1) and (A2), we can see that $\xi$ varies with the initial shape of the layer; the sharper the cusp at the center of nascent gyre, the greater the value of $\xi$.
<figure><img src="image_13.png" /><figcaption>**Figure B1.** The location of the saddle-node bifurcation depends on sea ice concentration for both solutions: $u_i > u_g$ in blue and $u_i < u_g$ in orange. The space between the two curves is a region in which both solutions exist. Despite the existence of two solutions in this region of parameter space, only one solution is valid.</figcaption></figure>

## Appendix B: Location of Saddle-Node Bifurcations in the Solutions

The presence of a square root in equations (7) and (9) leads to a saddle-node bifurcation in both solutions. These bifurcations are located at

$$ \kappa_{crit+} = +\frac{2}{f^2}\sqrt{\frac{\alpha C_{Di} g'^2}{\rho}(\alpha C_{Di} u_i^2 \rho + (1-\alpha)C_{Da} u_a^2 \rho_a)} - \frac{2\alpha C_{Di} u_i g'}{f^2} \quad (B1) $$

$$ \kappa_{crit-} = -\frac{2}{f^2}\sqrt{\frac{\alpha C_{Di} g'^2}{\rho}(\alpha C_{Di} u_i^2 \rho - (1-\alpha)C_{Da} u_a^2 \rho_a)} + \frac{2\alpha C_{Di} u_i g'}{f^2}, \quad (B2) $$

where $\kappa_{crit+}$ is the location of the saddle-node bifurcation in solutions with $u_i > u_g$, and $\kappa_{crit-}$ is the location for solutions with $u_i < u_g$. Figure B1 plots these two values as a function of the ice fraction, $\alpha$, using values from Table 1 for the remaining parameters. Between the two curves solutions to both equations (7) and (9) exist. Therefore, Figure B1 shows that for all $\alpha < 1$ there exists a region in which both solutions exist. However, in this overlap region one of the solutions violates the assumed sign of $u_i - u_g$ (Figure B2). Therefore, although there are two mathematically stable solutions in this region of parameter space, only one of them is physically valid.

<figure><img src="image_14.png" /><figcaption>**Figure B2.** Diagnosed geostrophic velocity from the two equilibrium solutions: $u_i - u_g > 0$ in blue and $u_i - u_g < 0$ in orange. The horizontal black line shows the prescribed ice velocity; any portion of the blue curve above this line or the orange curve below it violates the assumed sign of the relative velocity and is nonphysical. The nonphysical portion of each curve is shown by the dashed lines. The parameters are set to the values shown in Table 1, apart from $\kappa$, which forms the x axis.</figcaption></figure>
## Appendix C: Numerical Configuration for MITgcm Simulations

We construct a model of the Arctic using MITgcm (Marshall, Hill, et al., 1997, Marshall, Adcroft, et al., 1997) to solve the hydrostatic primitive equations of motion in a domain with realistic bathymetry and forcing. The simulation is integrated on the Arctic cap of a cubed-sphere grid, permitting relatively even grid spacing throughout the domain and avoiding polar singularities (Adcroft et al., 2004). The Arctic face comprises 210 × 192 grid cells with a mean horizontal grid spacing of 36 km. A conservative, fully nonlinear free surface is employed in conjunction with a real freshwater flux (Adcroft & Campin, 2004; Campin et al., 2004). There are 50 vertical levels ranging in thickness from 10 m near the surface to approximately 450 m at a maximum model depth of 6,150 m. We use the same blend of Smith and Sandwell (1997) and General Bathymetric Charts of the Oceans (GEBCO) 1-arc min bathymetry as Marshall et al. (2017). The nonlinear equation of state of Jackett and Mcdougall (1995) is used. Vertical mixing follows Large et al. (1994) with a background diffusivity of $6.0 \times 10^{-6}$ m$^2$/s. A seventh-order monotonicity-preserving advection scheme (Daru & Tenaud, 2004) is employed. A mesoscale eddy parameterization in the spirit of Gent and McWilliams (1990) and Gent et al. (1995) is used with the eddy diffusivity set to a constant value across the entire domain. The ocean model is coupled to a sea ice model described in Losch et al. (2010) and Heimbach et al. (2010). The model was forced by the JRA-25 (6 hr, 1°) reanalysis for the period 1979–2013, using bulk formulae following Large and Yeager (2004). Initial conditions for the ocean are from the World Ocean Circulation Experiment Global Hydrographic Climatology (annual mean, 1990–1998 from Gouretski & Koltermann, 2004). Open boundary conditions on *S*, *T*, *u*, and *v* were employed using “normal-year” conditions averaged from 1992 to 2002 derived from an ECCO climatology (Nguyen et al., 2011); flow boundary conditions were modified in the North Atlantic-Arctic passage to balance freshwater inputs from the JRA climatology and Dai and Trenberth (2002) river data set into the model domain.

To compare our analytical theory and the numerical simulations, we require long-term means from the simulations. These means are calculated by averaging the final 25 years of the simulations and therefore represent the climatology of our ensemble members between 1990 and 2014. To convert from two-dimensional fields to the scalar values required by our theory, we average the fields over the BG region defined by Proshutinsky et al. (2009).

## Acknowledgments

E. W. D., G. M., J. M., and J. S. gratefully acknowledge support from NSF Polar Programs, both Arctic and Antarctic. J. M. acknowledges support from the MIT-GISS collaborative agreement. C. L. received funding from the French LEFE/INSU program (project FREDY). Readers wishing to reproduce our numerical experiments can obtain the configuration and simulation scripts for our reduced gravity simulations (https://doddridge.me/publications/DMMSL2019/). The configuration of our MITgcm simulations is described in detail in Appendix C with citations to published descriptions of each model component and input field.

## References

* Adcroft, A., & Campin, J. M. (2004). Rescaled height coordinates for accurate representation of free-surface flows in ocean circulation models. *Ocean Modelling*, 7(3-4), 269–284. https://doi.org/10.1016/j.ocemod.2003.09.003

* Adcroft, A., Campin, J.-M., Hill, C., & Marshall, J. C. (2004). Implementation of an atmosphere ocean general circulation model on the expanded spherical cube. *Monthly Weather Review*, 132(12), 2845–2863. https://doi.org/10.1175/MWR2823.1

* Armitage, T. W. K., Bacon, S., Ridout, A. L., Thomas, S. F., Aksenov, Y., & Wingham, D. J. (2016). Arctic sea surface height variability and change from satellite radar altimetry and GRACE, 2003–2014. *Journal of Geophysical Research: Oceans*, 121, 4303–4322. https://doi.org/10.1002/2015JC011579

* Campin, J. M., Adcroft, A., Hill, C., & Marshall, J. C. (2004). Conservation of properties in a free-surface model. *Ocean Modelling*, 6(3-4), 221–244. https://doi.org/10.1016/S1463-5003(03)00009-X

* Cole, S. T., Toole, J. M., Lele, R., Timmermans, M.-L., Gallaher, S. G., Stanton, T. P., et al. (2017). Ice and ocean velocity in the Arctic marginal ice zone: Ice roughness and momentum transfer. *Elementa Science of the Anthropocene*, 5, 55. https://doi.org/10.1525/elementa.241

* Condron, A., Winsor, P., Hill, C., & Menemenlis, J. D. (2009). Simulated response of the Arctic freshwater budget to extreme NAO wind forcing. *Journal of Climate*, 22(9), 2422–2437. https://doi.org/10.1175/2008JCLI2626.1

* Dai, A., & Trenberth, K. E. (2002). Estimates of freshwater discharge from continents: Latitudinal and seasonal variations. *Journal of Hydrometeorology*, 3(6), 660–687. https://doi.org/10.1175/1525-7541(2002)003h0660:EOFDFCi2.0:CO2

* Danabasoglu, G., & Marshall, J. C. (2007). Effects of vertical variations of thickness diffusivity in an ocean general circulation model. *Ocean Modelling*, 18(2), 122–141. https://doi.org/10.1016/j.ocemod.2007.03.006

* Daru, V., & Tenaud, C. (2004). High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations. *Journal of Computational Physics*, 193(2), 563–594. https://doi.org/10.1016/j.jcp.2003.08.023

* Davis, P. E. D., Lique, C., & Johnson, H. L. (2014). On the link between arctic sea ice decline and the freshwater content of the Beaufort Gyre: Insights from a simple process model. *Journal of Climate*, 27(21), 8170–8184. https://doi.org/10.1175/JCLI-D-14-00090.1

* Davis, P. E. D., Lique, C., Johnson, H. L., & Guthrie, J. D. (2016). Competing effects of elevated vertical mixing and increased freshwater input on the stratification and sea ice cover in a changing Arctic Ocean. *Journal of Physical Oceanography*, 46(5), 1531–1553. https://doi.org/10.1175/JPO-D-15-0174.1

* Dewey, S., Morison, J., Kwok, R., Dickinson, S., Morison, D., & Andersen, R. (2018). Arctic ice-ocean coupling and gyre equilibration observed with remote sensing. *Geophysical Research Letters*, 45, 1499–1508. https://doi.org/10.1002/2017GL076229

* Dmitrenko, I. A., Kirillov, S. A., Ivanov, V. V., & Woodgate, R. A. (2008). Mesoscale Atlantic water eddy off the Laptev Sea continental slope carries the signature of upstream interaction. *Journal of Geophysical Research*, 113, C07005. https://doi.org/10.1029/2007JC004491

* Doddridge, E. W., Marshall, D. P., & Hogg, A. M. (2016). Eddy cancellation of the Ekman Cell in subtropical gyres. *Journal of Physical Oceanography*, 46(10), 2995–3010. https://doi.org/10.1175/JPO-D-16-0097.1
Doddridge, E. W., & Radul, A. (2018). Aronnax: An idealised isopycnal ocean model. *Journal of Open Source Software*, 3(26), 592. https://doi.org/10.21105/joss.00592

Ferreira, D., Marshall, J. C., & Heimbach, P. (2005). Estimating eddy stresses by fitting dynamics to observations using a residual-mean ocean circulation model and its adjoint. *Journal of Physical Oceanography*, 35(10), 1891–1910. https://doi.org/10.1175/JPO2785.1

Février, S., Sirven, J., & Herbaut, C. (2007). Interaction of a coastal Kelvin wave with the mean state in the gulf stream separation area. *Journal of Physical Oceanography*, 37(6), 1429–1444. https://doi.org/10.1175/JPO3062.1

Gent, P. R. (2016). Effects of Southern Hemisphere wind changes on the meridional overturning circulation in ocean models. *Annual Review of Marine Science*, 8(1), 79–94. https://doi.org/10.1146/annurev-marine-122414-033929

Gent, P. R., & McWilliams, J. C. (1990). Isopycnal mixing in ocean circulation models. *Journal of Physical Oceanography*, 20(1), 150–155. https://doi.org/10.1175/1520-0485(1990)020<0150:IMIOCM>2.0.CO;2

Gent, P. R., Willebrand, J., McDougall, T. J., & McWilliams, J. C. (1995). Parameterizing eddy-induced tracer transports in ocean circulation models. *Journal of Physical Oceanography*, 25(4), 463–474. https://doi.org/10.1175/1520-0485(1995)025h0463:PEITTli2.0.CO;2

Gouretski, V. V., & Koltermann, K. P. (2004). The World Ocean Circulation Experiment (WOCE) global hydrographic climatology (Tech. Rep) Hamburg, Germany: Bundesamtes fur Seeschifffahrt und Hydrographie. Retrieved from http://rda.ucar.edu/datasets/ds285.4/

Heimbach, P., Menemenlis, D., Losch, M., Campin, J. M., & Hill, C. (2010). On the formulation of sea-ice models. Part 2: Lessons from multi-year adjoint sea-ice export sensitivities through the Canadian Arctic Archipelago. *Ocean Modelling*, 33(1-2), 145–158. https://doi.org/10.1016/j.ocemod.2010.02.002

Hibler, W. D. (1979). A dynamic thermodynamic sea ice model. *Journal of Physical Oceanography*, 9(4), 815–846. https://doi.org/10.1175/1520-0485(1979)009<0815:ADTSIM>2.0.CO;2

Howard, E., Hogg, A. M., Waterman, S., & Marshall, D. P. (2015). The injection of zonal momentum by buoyancy forcing in a Southern Ocean model. *Journal of Physical Oceanography*, 45, 259–271. https://doi.org/10.1175/JPO-D-14-0098.1

Hunkins, K. L. (1974). Subsurface eddies in the Arctic ocean. *Deep Sea Research and Oceanographic Abstracts*, 21(12), 1017–1033. https://doi.org/10.1016/0011-7471(74)90064-3

Jackett, D. R., & Mcdougall, T. J. (1995). Minimal adjustment of hydrographic profiles to achieve static stability. *Journal of Atmospheric and Oceanic Technology*, 12(2), 381. https://doi.org/10.1175/1520-0426(1995)012<0381:MAOHPT>2.0.CO;2

Jahn, A., & Holland, M. M. (2013). Implications of Arctic sea ice changes for North Atlantic deep convection and the meridional overturning circulation in CCSM4-CMIP5 simulations. *Geophysical Research Letters*, 40, 1206–1211. https://doi.org/10.1002/grl.50183

Johnson, H. L., Cornish, S. B., Kostov, Y., Beer, E., & Lique, C. (2018). Arctic Ocean freshwater content and its decadal memory of sea-level pressure. *Geophysical Research Letters*, 45, 4991–5001. https://doi.org/10.1029/2017GL067860

Johnson, H. L., & Marshall, D. P. (2002). A theory for the surface Atlantic response to thermohaline variability. *Journal of Physical Oceanography*, 32, 1121–1132. https://doi.org/10.1175/1520-0485(2002)032<1121:ATFΤSA>2.0.CO;2

Kwok, R., & Morison, J. (2017). Recent changes in Arctic sea ice and ocean circulation. *US CLIVAR Var*, 15(3), 1–6. https://doi.org/10.5065/D6833QQP

Large, W. G., McWilliams, J. C., & Doney, S. C. (1994). Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. *Reviews of Geophysics*, 32(4), 363–403. https://doi.org/10.1029/94RG01872

Large, W. G., & Yeager, S. G. (2004). Diurnal to decadal global forcing for ocean and sea-ice models: The data sets and flux climatologies (NCAR Tech. Note, TN-460+ST, 105 pp.). Boulder, CO: National Center for Atmospheric Research. https://doi.org/10.5065/D6KK98Q6

Lique, C., Johnson, H. L., & Davis, P. E. D. (2015). On the interplay between the circulation in the surface and the intermediate layers of the Arctic Ocean. *Journal of Physical Oceanography*, 45(5), 1393–1409. https://doi.org/10.1175/JPO-D-14-0183.1

Losch, M., Menemenlis, D., Campin, J. M., Heimbach, P., & Hill, C. (2010). On the formulation of sea-ice models. Part 1: Effects of different solver implementations and parameterizations. *Ocean Modelling*, 33(1-2), 129–144. https://doi.org/10.1016/j.ocemod.2009.12.008

Manley, T. O., & Hunkins, K. (1985). Mesoscale eddies of the Arctic Ocean. *Journal of Physical Research*, 90(C3), 4911–4930. https://doi.org/10.1029/JC090iC03p04911

Manucharyan, G. E., & Spall, M. A. (2016). Wind-driven freshwater buildup and release in the Beaufort Gyre constrained by mesoscale eddies. *Geophysical Research Letters*, 43, 273–282. https://doi.org/10.1002/2015GL065957

Manucharyan, G. E., Spall, M. A., & Thompson, A. F. (2016). A theory of the wind-driven Beaufort Gyre variability. *Journal of Physical Oceanography*, 46(11), 3263–3278. https://doi.org/10.1175/JPO-D-16-0091.1

Marshall, J. C., Adcroft, A., Hill, C., Perelman, L., & Heisey, C. (1997). A finite-volume, incompressible Navier Stokes model for studies of the ocean on parallel computers. *Journal of Geophysical Research*, 102(C3), 5753–5766. https://doi.org/10.1029/96JC02775

Marshall, J. C., Hill, C., Perelman, L., & Adcroft, A. (1997). Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. *Journal of Geophysical Research*, 102(C3), 5733–5752. https://doi.org/10.1029/96JC02776

Marshall, J. C., & Radko, T. (2003). Residual-mean solutions for the Antarctic Circumpolar Current and its associated overturning circulation. *Journal of Physical Oceanography*, 33(11), 2341–2354. https://doi.org/10.1175/1520-0485(2003)033h2341:RSFTACi2.0.CO;2

Marshall, J., Scott, J., & Proshutinsky, A. (2017). Climate response functions for the Arctic Ocean: A proposed coordinated modelling experiment. *Geoscientific Model Development*, 10(7), 2833–2848. https://doi.org/10.5194/gmd-10-2833-2017

Martin, T., Steele, M., & Zhang, J. (2014). Seasonality and long-term trend of Arctic Ocean surface stress in a model. *Journal of Geophysical Research: Oceans*, 119, 1723–1738. https://doi.org/10.1002/2013JC009425

Maxwell, J. C. (1867). On Governors. *Proceedings of the Royal Society of London*, 16, 270–283. https://doi.org/10.1098/rspl.1867.0055

Meneghello, G., Marshall, J., Campin, J.-M., Doddridge, E. W., & Timmermans, M.-L. (2018). The Ice-Ocean Governor: Ice-ocean stress feedback limits Beaufort Gyre spin-up. *Geophysical Research Letters*, 45, 11,293–11,299. https://doi.org/10.1029/2018GL080171

Meneghello, G., Marshall, J. C., Cole, S. T., & Timmermans, M. L. (2017). Observational inferences of lateral eddy diffusivity in the halocline of the Beaufort Gyre. *Geophysical Research Letters*, 44, 12,331–12,338. https://doi.org/10.1002/2017GL075126

Meneghello, G., Marshall, J. C., Timmermans, M.-L., & Scott, J. (2018). Observations of seasonal upwelling and downwelling in the Beaufort Sea mediated by sea ice. *Journal of Physical Oceanography*, 48(4), 795–805. https://doi.org/10.1175/JPO-D-17-0188.1

Nguyen, A. T., Menemenlis, D., & Kwok, R. (2011). Arctic ice-ocean simulation with optimized model parameters: Approach and assessment. *Journal of Geophysical Research*, 116, C04025. https://doi.org/10.1029/2010JC006573

Nost, O. A., & Isachsen, P. E. (2003). The large-scale time-mean ocean circulation in the Nordic Seas and Arctic Ocean estimated from simplified dynamics. *Journal of Marine Research*, 61(2), 175–210. https://doi.org/10.1357/002224003322005069

Nost, O. A., Nilsson, J., & Nycander, J. (2008). On the asymmetry between cyclonic and anticyclonic flow in basins with sloping boundaries. *Journal of Physical Oceanography*, 38(4), 771–787. https://doi.org/10.1175/2007JPO3714.1

Proshutinsky, A., Krishfield, R., Timmermans, M.-L., Toole, J., Carmack, E., McLaughlin, F., et al. (2009). Beaufort Gyre freshwater reservoir: State and variability from observations. *Journal of Geophysical Research*, 114, C00A10. https://doi.org/10.1029/2008JC005104
Regan, H. C., Lique, C., & Armitage, T. W. K. (2019). The Beaufort Gyre extent, shape, and location between 2003 and 2014 from satellite observations. *Journal of Geophysical Research: Oceans*, 124, 844–862. https://doi.org/10.1029/2018JC014379

Simmonds, I. (2015). Comparing and contrasting the behaviour of Arctic and Antarctic sea ice over the 35 year period 1979-2013. *Annals of Glaciology*, 56(69), 18–28. https://doi.org/10.3189/2015AoG69A909

Smith, W. H. F., & Sandwell, D. T. (1997). Global sea floor topography from satellite altimetry and ship depth soundings. *Science*, 277(5334), 1956–1962. https://doi.org/10.1126/science.277.5334.1956

Stern, M. E. (1998). Separation of a density current from the bottom of a continental slope. *Journal of Physical Oceanography*, 28(10), 2040–2049. https://doi.org/10.1175/1520-0485(1998)028<2040:SOADCF>2.0.CO;2

Stewart, K. D., & Haine, T. W. N. (2013). Wind-driven Arctic freshwater anomalies. *Geophysical Research Letters*, 40, 6196–6201. https://doi.org/10.1002/2013GL058247

Stouffer, R. J., Yin, J., Gregory, J. M., Dixon, K. W., Spelman, M. J., Hurlin, W., et al. (2006). Investigating the causes of the response of the thermohaline circulation to past and future climate changes. *Journal of Climate*, 19(8), 1365–1387. https://doi.org/10.1175/JCLI3689.1

Thompson, D. W. J., & Wallace, J. M. (1998). The Arctic oscillation signature in the wintertime geopotential height and temperature fields. *Geophysical Research Letters*, 25(9), 1297–1300. https://doi.org/10.1029/98GL00950

Timmermans, M.-L., Toole, J., Proshutinsky, A., Krishfield, R., & Plueddemann, A. (2008). Eddies in the Canada Basin, Arctic Ocean, observed from ice tethered profilers. *Journal of Physical Oceanography*, 38(1), 133–145. https://doi.org/10.1175/2007JPO3782.1

Wolfram-Alpha (2018). Wolfram-Alpha. Wolfram-Alpha LLC. Retrieved 2018-11-15, from https://www.wolframalpha.com/input/?i=dh%2Fdt+%3D+a%h%5E2+%2B+b*h+%2B+c

Yang, J., Proshutinsky, A., & Lin, X. (2016). Dynamics of an idealized Beaufort Gyre: 1. The effect of a small beta and lack of western boundaries. *Journal of Geophysical Research: Oceans*, 121, 1249–1261. https://doi.org/10.1002/2015JC011296

Zhao, M., Timmermans, M.-L., Cole, S., Krishfield, R., Proshutinsky, A., & Toole, J. (2014). Characterizing the eddy field in the Arctic Ocean halocline. *Journal of Geophysical Research: Oceans*, 119, 8800–8817. https://doi.org/10.1002/2014JC010488

Zhao, M., Timmermans, M.-L., Cole, S., Krishfield, R., & Toole, J. (2016). Evolution of the eddy field in the Arctic Ocean's Canada Basin, 2005–2015. *Geophysical Research Letters*, 43, 8106–8114. https://doi.org/10.1002/2016GL069671

Zhong, W., Steele, M., Zhang, J., & Cole, S. T. (2019). Circulation of Pacific winter water in the Western Arctic Ocean. *Journal of Geophysical Research: Oceans*, 124, 863–881. https://doi.org/10.1029/2018JC014604

Zhong, W., Steele, M., Zhang, J., & Zhao, J. (2018). Greater role of geostrophic currents in Ekman dynamics in the Western Arctic Ocean as a mechanism for Beaufort Gyre stabilization. *Journal of Geophysical Research: Oceans*, 123, 149–165. https://doi.org/10.1002/2017JC013282

Zweng, M. M., Reagan, J. R., Seidov, D., Boyer, T. P., Locarnini, R. A., Garcia, H. E., & Smolyar, I. (2018). World Ocean Atlas 2018, Volume 2: Salinity In A. Mishonov (Eds.), *NOAA Atlas NESDIS 82* (40 pp.). Retrieved from https://data.nodc.noaa.gov/woa/WOA18/DOC/woa18_vol2.pdf