Hotspots of Dense Water Cascading in the Arctic Ocean: Implications for the Pacific Water Pathways

We explore dense water cascading (DWC), a type of bottom‐trapped gravity current, on multidecadal time scales using a pan‐Arctic regional ocean‐ice model. DWC is particularly important in the Arctic Ocean as the main mechanism of ventilation of interior waters when open ocean convection is blocked by strong density stratification. We identify the locations where the most intense DWC events occur and evaluate the associated cross‐shelf mass, heat, and salt fluxes. We find that the modeled locations of cascading agree well with the sparse historical observations and that cascading is the dominant process responsible for cross‐shelf exchange in the boundary layers. Simulated DWC fluxes of 1.3 Sv (1 Sv = 10 m/s) in the Central Arctic are comparable to Bering Strait inflow, with associated surface and benthic Ekman fluxes of 0.85 and 0.58 Sv. With ice decline, both surface Ekman flux and DWC fluxes are increasing at a rate of 0.023 and 0.0175 Sv/year, respectively. A detailed analysis of specific cascading sites around the Beaufort Gyre and adjacent regions shows that autumn upwelling of warm and saltier Atlantic waters on the shelf and subsequent cooling and mixing of uplifted waters trigger the cascading on the West Chukchi Sea shelf break. Lagrangian particle tracking of low salinity Pacific waters originating at the surface in the Bering Strait shows that these waters are modified by brine rejection and cooling, and through subsequent mixing become dense enough to reach depths of 160–200 m. Plain Language Summary In this study we explore dense water cascading, а specific type of bottom‐trapped gravity current. This current of very dense waters originates on the shelves, due to winter cooling and sea ice freezing, and slowly propagates downslope to deep waters. It is specifically important in the Arctic Ocean as themainmechanism of deep water mass formation and carbon storage. We use numerical model of the Arctic Ocean to predict preferable locations of the cascading, its intensity mass, and heat fluxes. Our results are in agreement with available very sparse observations. Our model predicts that with sea ice melting, cascading formation will accelerate.


Introduction
Gradual reduction in Arctic sea ice cover in the 1990-2000s in response to global warming probably acted as a trigger for the accelerated transition toward a seasonally ice-free Arctic Ocean (Kwok et al., 2009). Enhanced ocean-air energy exchange, driving vertical mixing in the ocean, is a fundamental consequence of increasing exposure of the sea surface to atmospheric forcing. Multiple examples indicating that this effect is already at work at the margins of the Arctic Ocean deep interior have recently been reported (Ivanov et al., 2016;Ivanov & Repina, 2019;Lind et al., 2018;Polyakov et al., 2017;Timmermans et al., 2018).
One of the crucial processes contributing to the upper ocean to abyssal exchange in the Arctic Ocean is shelf convection, also known as dense water cascading (DWC), a specific type of bottom-trapped gravity current (Shapiro et al., 2003). At high latitudes, the dense waters feeding DWC are formed on shallow shelves in the cold season due to heat loss to the atmosphere, resulting in freezing and subsequent brine injection into the water column. Preconditioning of DWC happens if the convection reaches the bottom, forming a well-mixed water mass denser than the waters downslope. In this case, dense water propagates mainly along-slope (due to Coriolis forcing) and downslope (due to bottom friction) as a bottom boundary current, mixing with the ambient waters en route. Descending DWC waters interact with bottom topographic features such as sills and canyons (Chapman & Gawarkiewicz, 1995;Kämpf, 2005), with large-scale currents (Aksenov et al., 2011;Ivanov et al., 2015) and form meanders and eddies due to baroclinic instabilities (e.g., Tanaka & Akitomo, 2001). After reaching the deep basin floor or the neutral density level (where the density in the plume equals the density of the ambient waters), the DWC plume spreads as either a bottom layer (Ivanov et al., 2004) or intrudes in the adjoining water column as a lens shape with specific thermohaline characteristics (Koenig et al., 2018).
In the Arctic Ocean, DWC is particularly important because of the strong stratification in the upper water column, making the ocean essentially two layered (e.g., Proshutinsky & Johnson, 1997). Extremely strong density gradients between the fresh upper layer and the underlying waters inhibit open ocean convection from reaching deep water. This makes DWC the only process capable of deep ventilation of the interior Arctic basins. Specifically, in Arctic and sub-Arctic regions, cascades may replenish water in the halocline (Aagaard et al., 1981), in the intermediate warm and salty Atlantic layer (Ivanov & Golovin, 2007), and in the deep layer (Falina et al., 2012;Quadfasel et al., 1988). Indirect evidence exists that DWC may contribute to the replenishment of the bottom water mass in the Canada Basin, thus making it saltier and warmer than the bottom water in the Eurasian Basin (Rudels, 1986). However, no direct measurements supporting this concept have been reported. The geography of cascading around the Arctic Ocean, as revealed by historical hydrographic data in the twentieth century, is nonuniform. The highest density of DWC events ever captured by direct measurements was observed in the East-Atlantic sector of the Arctic Ocean (see Figure 1 and Ivanov et al., 2004). Several DWC sites were also observed at the North American continental slope (see Figure 1 and Ivanov et al., 2004). However, the broad shallow shelves between the Laptev Sea and the Bering Strait are not suitable for dense water formation due to the extreme freshening of surface waters by river runoff. Recent high-resolution observations by submarine gliders captured submesoscale lenses filled with relatively cold and fresh waters in the core of the West Spitsbergen Current, west of Svalbard (Koenig et al., 2018). These lenses are highly likely to contain shelf waters that had cascaded downslope and intruded into the West Spitsbergen current at their neutral density levels.
With the transition toward seasonal Arctic sea ice cover, the DWC spatial distribution is expected to change. While the Eurasian Arctic seas progressively evolve to a permanently ice-free and well-mixed Atlantic-dominated hydrographic regime (Lind et al., 2018;Polyakov et al., 2020), the North American Arctic shelf seas acquire more pronounced seasonal features, with extended ice-free areas in summer (Timmermans et al., 2018) while remaining fully ice covered in winter. This change theoretically creates favorable prerequisite conditions for stronger salinization of water on shelves after the onset of freezing and for enhanced polynya activity in winter in the regions historically predisposed to polynyas (Cavalieri & Martin, 1994;Winsor & Björk, 2000). Ivanov and Watanabe (2013) hypothesized an increased rate of cascading waters in the future and found support for this hypothesis in an idealized modeling study and a global model, comparing 1986-1987 (more ice) and 2006-2007 (less ice) conditions in the Laptev Sea (Ivanov et al., 2015).
Being a localized and sporadic process and occurring mostly under sea ice in winter, DWC is hard to observe. So far, only 27 locations around the Arctic Ocean have been identified where DWC events were confirmed at some times within the entire period of instrumental observations (see Figure 1). To the best of our knowledge, there are no well-justified estimates of the net DWC fluxes on a pan-Arctic scale, except a very approximate one by Shapiro et al. (2003). In this study we make a first attempt to fill this gap by using model-based estimates of DWC and comparing them with other cross-shelf mass and heat flux estimates and their trends.
This study is performed within the framework of the Forum for Arctic Ocean Modeling and Observational Synthesis (https://famosarctic.com/teams/Cascading) and aims to evaluate the ensemble mean cross-shelf volume, heat, and salt cascading fluxes, along with their uncertainties. An ensemble of global and Arctic models with different resolutions (from 2 to 13 km in the Arctic) and with different types of vertical coordinates presents contributions from two Nucleus for European Modeling of the Ocean framework for ocean climate, research and operational oceanography (NEMO) model configurations. In this study, we introduce and test methods of flux evaluations at observed locations of DWC and estimate the possible strength of the cascading compared with other components of cross-shelf flux. We perform the first flux evaluation using the NEMO-shelf model with 1/4°spatial resolution, explicitly simulated tides, and well-resolved benthic layers (Luneva et al., 2015). We also use output from a global NEMO model with 1/12°spatial resolution but without tides and coarser vertical resolution to track the spread of the Pacific waters in the Arctic Ocean (Kelly et al., 2018).
The paper is organized as follows. In section 3, we describe existing analytical estimates of DWC fluxes and methods to evaluate DWC and cross-shelf fluxes and discuss the importance of vertical and horizontal resolutions in simulations of DWC. Section 4 gives details on model setup and numerical experiments. In section 5, we present results for the entire Arctic Ocean, while section 6 focuses on case studies in the Beaufort and Chukchi Seas. Section 0 summarizes the findings and discusses the results in the context of the recent changes in Arctic sea ice.

Scaling and Theoretical Studies
The depth range to which dense water descends downslope depends on many factors. Horizontal density contrast across the shelf edge and the steepness of the slope are the major parameters that determine the speed of dense water propagation along the slope, the so-called Nof (1983) speed:  Table 1. Contours of 300-m isobaths, used in analysis, are marked by color. The general bathymetry is defined by the associated color scale.
where g is the gravity constant, δρ/ρ o is the relative density anomaly in the dense plume compared with ambient waters, f is the Coriolis parameter, and s = ∇ H is slope of bathymetry, H. The pace of cross-slope "spillage," which is proportional to the Nof speed, together with density stratification in the ambient water, controls the ultimate depth of the dense plume spreading (Shapiro et al., 2003). In a laterally homogeneous 2-D case (with no variations along the shelf break), the dynamics of the downward-propagating dense front, in a steady state, are evaluated from a "1½-layer" model (Shapiro & Hill, 1997, their equation 22). That study incorporates bottom topography, Earth rotation, internal and bottom friction, and entrainment, as well as externally imposed pressure gradients, and defines the downslope-propagating velocity of a front as a fraction of Nof velocity which, for the case of a fully developed plume with maximal upstream thickness, is Following Ivanov et al. (2004), the observationally estimated Nof velocity is in the range of 2 × 10 −3 to 0.35 m/s, which gives the downward (cross-slope) velocity 4 × 10 −4 to 0.07 m/s. Theoretical estimate of the length scale of the thickness of a cascading plume is about 1.8h eb , where h eb is bottom Ekman depth, which is defined for a flow with constant turbulent viscosity (A v ) as The resulting cascading flux is Equation 4 is consistent with the laboratory experiments by Condie (1995) and Shapiro and Zatsepin (1997). Wobus et al. (2011), in a model with fine vertical resolution, reproduced the laboratory results. This suggests hydrostatic models are capable of reproducing the physics of cascading even on a steep slope, which is essential in the context of the current study.
On sufficiently steep slopes with large density gradients, the steady-state solution does not exist. Huthnance (2009) examined an initial value problem analytically and found that steady state exists when the dimensionless parameter If this condition is violated, the cascading flux accelerates with time and depth, and the thickness of the cascading tongue becomes much larger than the Ekman depth.

Processes Affecting Cascading
In the real ocean, cascading fluxes are affected by ambient stratification, ice conditions, along-slope geostrophic currents, and upwelling and downwelling caused by local and nonlocal winds and by eddies (Harden et al., 2014;Huthnance, 1995). If the cascading velocity is strong enough, Kelvin-Helmholtz instability will result in the entrainment of the ambient fluid, changing the density and thickness of the plume. At the edge of a plume, entrainment is likely in the benthic boundary layer. Mesoscale eddies may appear as a result of DWC flow instability and potential energy release, thus leading to increased cross-isobath water exchange (Platov & Golubeva, 2019). Ivanov and Watanabe (2013) and Ivanov et al. (2015) suggested that a large seasonal variability of sea ice conditions on the shelves (ice-covered shelves in winter and ice-free in summer) would result in a stronger heat loss during ice-free period and brine rejection during new ice formation, thus leading to a stronger buoyancy loss and stronger cascading ( Figure 2a). More extensive ice cover prevents heat loss and new ice formation, resulting in less dense water forming on the shelves. With the summer ice edge retreating beyond the outer continental shelf and the pack ice transforming in the Marginal Ice Zone, due to break up by waves and wind over the shelf slopes (Aksenov et al., 2017), dense water forms closer to the shelf break. Offshore winds expose shelf waters to densification due to cooling and new ice formation, providing further preconditioning for cascading ( Figure 2a).

Figure 2.
Schematics of the Arctic Ocean stratification, cascading process, and processes affecting cascading. Ice is indicated by cyan color, mixed upper waters remain uncolored, cold fresh halocline is shown by navy blue, pink is warm and salty Atlantic origin water (AW), magenta shows cold, salty Arctic abyssal waters, and green represents the cascading plume. (a) Ice concentration and offshore winds. Wind blowing offshore results in larger heat uptake from the ocean on the shelf, stronger ice formation, and brine rejection. A cascading plume (shown by green) forms on the shelf break, descends along the shelf break to the depth of ambient neutral density, and propagates along it. (b) Cyclonic winds move surface waters onshore, resulting in compensating downwelling, offshore fluxes in the benthic layer, and support descending cascading flux. (c) Anticyclonic winds produce offshore fluxes in the surface layers and compensating upwelling-onshore fluxes. This might result in uplift of warm and saltier AW, cooling, and ventilation of which would support the formation of the relatively dense waters in the plume and thus deeper cascading.
In the Arctic Ocean, geostrophic currents closely follow the contours of f/H in a cyclonic manner, with the shelf to the right of the current (e.g., Aksenov et al., 2011). This causes a benthic Ekman drain of the shelf waters directed offshore, supporting downward water fluxes, while anticyclonic along-slope currents generate upslope Ekman flux (e.g., Bacon et al., 2014). In contrast, the effects of the wind and surface Ekman currents on DWC are not well understood. Basin-scale cyclonic winds drive water onshore in the surface Ekman layer, resulting in compensating downwelling in the benthic layer ( Figure 2b). Therefore, cyclonic winds can enhance DWC. Basin-scale anticyclonic winds drive surface waters offshore, creating compensating upwelling in the benthic layer, thus impeding DWC (Figure 2b). Following this mechanism, we might expect a negative correlation between offshore surface Ekman flux (SEF) and DWC fluxes, although the details depend on the shelf topography that controls the pathways of compensating upwelling waters. On the other hand, anticyclonic winds in winter drive compensating flux of denser waters on the shelf (Proshutinsky & Johnson, 1997). This increases the ventilation and mixing of warm Atlantic waters with newly formed DWC water masses ( Figure 2c). Cabbeling can also increase the density of the mixed waters. This effect can be compared with the preconditioning of deep convection in the Labrador Sea, where doming of isopycnals toward the surface supports ventilation and enhances deep convection (Marshall & Schott, 1999).
Wind-driven circulation can also change the freshwater pathways from the rivers along the Siberian and Canadian coasts (Timmermans et al., 2011), affecting stratification on the shelves. Local effects of winds may also force cascading waters across the shelf and off the shelf through the generation of Ekman downwelling cells (Harden et al., 2014). Examining DWC in the East Greenland shelf, the so-called Spill Jet Current, Pickart, Torres, and Fratantoni (2005) found that the combined effect of cyclonic barrier winds and eddies results in the strongest spilling events. The observed near-bottom intensified offshore currents exceeded 0.25 m/s, which is substantially more than the cross-stream velocities estimated from the steady-state theory of Shapiro and Hill (1997) and by Ivanov et al. (2004).

Flux Evaluation
To examine cascading, its preconditioning, and other components of cross-shelf exchange, we use model output (5-day means) for variables θ, S, u, v, w, τ ! s , and τ ! b representing, respectively, potential temperature, salinity, velocity, and ocean shear stresses on the surface (ocean-atmosphere or ice-ocean). Bottom shear stress τ ! b , if not present in the model outputs, is calculated as where u ! b is the velocity at the closest liquid grid point near the bottom and C D is the quadratic bottom drag coefficient. Here we calculate the latter from the logarithmic law: where κ = 0.4 is the Von Karman constant, Δz is the near-bottom-layer thickness, and z 0 = 0.001 m is the bottom roughness. The tidal component of bottom stress is not considered in this formula, which could result in an underestimation of bottom shear stresses by 50% (Luneva et al., 2015).
All fluxes are calculated in the normal and along-slope directions and integrated over the particular depths relevant to each process considered. The modulus of the unit normal vector is Then the normal and tangential components of velocity flux are Here, u n is taken as positive offshore, and u l is positive while following contours of f/H, with shallow water on the right-hand side. We define ds l as an element of the isobath and ds n as an element of the contour directed downslope:

Preconditioning to Cascading
Here we describe conditions appropriate for the initiation of a cascading event. For water to be able to cascade downslope in the benthic layer, the mixed layer depth (MLD) on the shelf-slope must reach the bottom: where MLD is defined by the density criteria, with threshold value of Δρ ρ 0 ¼ 10 −4 , where Δρ is the density difference between the ocean surface and the bottom. If Condition 9 is satisfied, we define the mean salinity over the depth of mixed waters from sea surface ζ, to the bottom as and density as Here, T is the period of temporal averaging (5 days in our case) over local variable t′ ∈ t − T 2 ; t þ T 2 ! , and S is salinity. To separate well-mixed zones from the stratified areas, we set S ini (x, y, t) = 0, ρ ini (x, y, t) = 0 if 9 is not satisfied. The temperature of the mixed waters is typically at the local freezing point. If ρ ini is larger than that of water downslope and the MLD reaches steep enough slopes, water will start to propagate along and downslope with the speed in the plume defined by the Nof velocity (if there are no other effects). We take a threshold value of the Nof velocity, defined by Equation 1, as 2 × 10 −3 m s −1 , consistent with observations (Ivanov et al., 2004). The condition for cascading then becomes u n > 0; δρ · s=ρ 0 < −2 × 10 −8 : This is a weak condition because it corresponds to small density differences, only twice the threshold of mixed waters (|δρ/ρ 0 | > 2 × 10 −4 ) and to weak slopes, s = ∇ H > 10 −4 , corresponding to 1 m per typical dx~10 km in our simulations.
We apply this condition to the adjacent model cells, where the density difference may not be large. The Nof velocity scale is relevant to the density difference across the cascading front, that is, between the maximum density inside the cascading plume and the density at the same depth in the undisturbed by cascading ambient water. The spatial scale of this front is expected to be much larger than the separation between adjacent model cells.

Offshore Cascading Fluxes
In practice, cascading fluxes are hard to separate from dense water flows of other origin. We identify a cascading flux as that where the denser water lies upslope and propagates downslope, subjected to condition 11. The cascading flux is then a mass flux, F n m , integrated over the layer z ∈ (−H, −Z*), that is, from the bottom to the upper plume boundary Z*, where condition 11 is just satisfied.
Corresponding salinity and heat fluxes are Here, θ ref and S ref are reference potential temperature and salinity, and c p is specific heat. F n m is always positive.
The choice of the reference values for the heat and freshwater fluxes calculations is not unique and depends on the focus of the analysis (e.g., see Bacon et al., 2015). For the purpose of this study, we use 5 years of mean benthic values to evaluate cascading fluxes with reference to the ambient waters. The relative salt flux anomaly (in freshwater volume flux units) is evaluated with reference to the 5-year mean salinity at the bottom, Finally, the height of the cascading plume is defined as We also calculate theoretical estimate of the cascading flux, for conditions when all other forcing (along-slope currents, Ekman drains, and stratification of surrounding waters) are absent, given by Equation 3 (Shapiro & Hill, 1997).
As a theoretical estimate of DWC flux, we use Equation 3, V down , defined by 2, and for a density contrast, δρ/ρ 0 , we take the maximal value in the dense water plume.

Ekman Surface and Benthic Fluxes
The theoretical estimates of the Ekman wind driven flux in steady state follow Gill (1968) and are independent of the choice of vertical eddy viscosity: where τ sl = τ sx n y − τ sy n x and τ bl = τ bx n y − τ by n x are along-slope components of the surface and bottom shear stresses, applied to the ocean. Bottom shear stresses act in the opposite direction to the current, and steady-state Ekman flux is perpendicular to the shear stresses. In the Eurasian Arctic, where the near-bottom circulation is cyclonic (follows f/H contours; Aksenov et al., 2011;Gavilan Pascual-Ahuir et al., 2020), the bottom Ekman flux (BEF) is positive and offshore. If the wind circulation is anticyclonic, shear stresses applied to the ocean are negative, and the SEF is positive and offshore.
In a layer of constant viscosity A v , a theoretical estimate of the surface Ekman depth, h es , depends on the turbulent viscosity and the Coriolis parameter: and for a flow with constant shear, as often used in the atmospheric boundary layer (Wimbush & Munk, 1970): Note that, unlike the Ekman fluxes ( (14a) and (14b)), estimates of Ekman depths ( (15a) and (15b)) are approximate and valid (and correspond to MLD) only in the unstratified case.
The benthic Ekman depth is an important parameter in this study. Indeed, the theoretically predicted thickness of a cascading plume is about twice the Ekman depth. Thus, the capability of the vertical discretization of the model to resolve the benthic Ekman depth at least in the initial stages of cascading is crucial in this study. Theoretical estimates of Ekman depth depend on the type of flow, presence/absence of tides (Prandle, 1982;Soulsby, 1983), and stratification. The commonly used Ekman depth boundary layer formulation, similar to the surface one 15b, was found to be typically smaller by a half than that observed in the ocean and shelf seas (Thorpe, 1988). Soulsby (1983) proposed (assuming C D = 1.9 × 10 −3 ) In stratified conditions, Weatherley and Martin (1978) suggest In our study we will use (15b) and (15c) for our estimates of surface and benthic boundary layers. Recalculating the expression with C D = 2.5 × 10 −3 gives is the absolute value of velocity at the grid cell adjacent to the bottom, that is, twice as large as In addition to the theoretical expressions ((14a) and (14b)), we estimate surface and bottom cross-shelf fluxes in the Ekman layers (thereafter SEFm and BEFm): Unlike the steady-state, balanced theoretical estimates of cross-shelf Ekman fluxes ( (15a) and (15b)), SEFm and BEFm ((16a) and (16b)) also include cross-slope component of the geostrophic currents, for example, induced by mesoscale eddies, and contributions of cascading flux. Note that cascading depth and surface Ekman depths could be similar, especially at the preconditioning phase of cascading when all waters are well mixed from the surface to bottom: In this case we expect strong interactions of benthic-, surface-, and density-driven fluxes.

Effects of Vertical and Horizontal Resolutions
In the previous section we emphasized the importance of resolving the benthic boundary layers in the simulation of cascading. Horizontal resolution is another important issue, specifically related to the resolution of topographical features such as canyons and troughs, as these are preferred pathways for cascading plumes (Ivanov & Golovin, 2007;Kämpf, 2005). If the propagation speed of dense water is sufficiently high, eddies associated with baroclinic instability of the current and the interaction with eddies external to the plume may also play a role (Harden et al., 2014;Magaldi & Haine, 2015).
There are few studies examining the effects of resolution on the simulation of cascading. In idealized simulations of frictional downslope flow, Winton et al. (1998) found that a z-coordinate model gave results consistent with theory when the (1) horizontal resolution was high enough to resolve the slopes and (2) the vertical grid spacing resolved several points within the bottom frictional layer. Shapiro et al. (2013), on examining different vertical discretizations, found that their new hybrid scheme gave the best results both in reproducing cascading compared with z coordinates and reducing pressure gradient errors compared with terrain-following s coordinates (Song & Haidvogel, 1994). This hybrid scheme, combined with the NEMO s-coordinates scheme, masking the very steep slopes (Madec & NEMO Team, 2016, Figure 4.5d), was applied only on the shelf, with z-partial steps elsewhere. A recent study by Bruciaferri et al. (2018) examines the effects of vertical resolution on an idealized numerical experiment of cascading in the Black Sea, imitating conditions of the Shapiro and Hill (1997) theory. Their model has a very fine horizontal resolution (~3.5 km, which is 7 points over a baroclinic Rossby radius of 20 km). However, they were unable to reproduce the analytical solution with z-partial steps vertical coordinates while successfully simulating it with hybrid coordinates on the shelf (the same type of coordinates as used by Shapiro et al., 2013, Luneva et al., 2015, and in current study). They found that using z coordinates resulted in spurious entrainment and a slowdown of the cascading plume by two orders of magnitude.
Several studies have examined similar issues of vertical coordinate type and resolution in the simulation of overflows (Colombo, 2018(Colombo, , 2020Ezer & Mellor, 2004;Legg et al., 2009). Unlike cascading (see Table 1 from Legg et al., 2009), overflows are powerful flows with velocities~0.5-1.5 m/s, large plume thickness (~100-300 m), and propagate long distances (~1,000 km) across the ocean. Hence, unlike in cascades, nonlinearity (large Froude numbers), entrainment due to Kelvin-Helmholtz instability, and formation of eddies are crucial in overflow dynamics. In idealized overflow simulations, on constant slope, Ezer and Mellor (2004) examined the impact of both vertical and horizontal resolutions. They found that using terrain-following coordinates, locally or globally, significantly improves the representation of overflows: Eddies were reproduced in the model with a sigma-coordinate vertical grid and 10-km horizontal resolution (and low horizontal diffusivity), while in the z-level calculations, eddies develop only at a 2.5-km horizontal resolution. Increasing horizontal and vertical resolutions in the z-level grid converges the results toward those obtained by a much coarser horizontal resolution sigma-coordinate model. Similarly, a recent study by Colombo et al. (2020) found that using local terrain-following coordinates in the simulation of the Denmark Strait overflow with a 1/12°model (~6-to 8-km resolution) and 75 hybrid z-s-levels is equivalent to their best z-coordinate simulations with 1/60°resolution (1-2 km) and 150 vertical levels. All the studies discussed demonstrate that refining the vertical resolution at some point compensates for coarseness of horizontal resolution.

Pan-Arctic Model
In this study we use a medium-resolution (6-18 km) pan-Arctic regional NEMO model with explicitly resolved shelf processes (http://www.nemo-ocean.eu/). The horizontal resolution of the baseline model corresponds to a 1/4°Ocean General Circulation Model (ORCA) global ocean model, but the vertical resolution is much finer across the shelf-slope, with terrain-following coordinates on the shelf, refining at the bottom to resolve the benthic layer. The model setup is similar to the regional pan-Arctic set up of Luneva et al. (2015, using NEMO v3.2) but featuring a number of improvements described below. A more-recent version of the NEMO model, v3.6 (Madec, 2008;O'Dea et al., 2017), is used in this study.
NEMO v3.6 has been modified and tested to include shelf sea and shelf-break processes (Luneva et al., 2015;O'Dea et al., 2017), explicitly simulating tides and resolving benthic layers with hybrid (terrain followinggeopotential) vertical coordinates, and hereafter referred to as NEMO-shelf. To resolve tides, a nonlinear free-surface formulation with variable volume and mode-split explicit time stepping has been employed. However, it is noted that NEMO-shelf does not yet have the option for the low-diffusive, vertical piecewise parabolic advection scheme, which was available in NEMO v3.2 (used in Luneva et al., 2015).
The regional, pan-Arctic configuration considered here (see domain in Figure 3a) employs the mesh and topography from the global 0.25°NEMO tripolar grid (Madec, 2008) but reindexed to create a seamless rectangular model grid. Two open liquid boundaries are located in the Pacific Ocean (at 55°N) and in the Atlantic Ocean (near 48°N). The horizontal resolution corresponds to 18 km at the southern periphery of the domain, 12-13 km in the deep Arctic Ocean, 10 km on the Siberian shelf, and reduces to 6 km in the Canadian Arctic Archipelago. So, in these regions the horizontal resolution is comparable to models frequently used in shelf-ocean exchange studies (e.g., Holt et al., 2009). The model is eddy-and internal-tide permitting in the deep part of the Arctic Ocean but too coarse to resolve eddies and internal tides on the 10.1029/2020JC016044

Journal of Geophysical Research: Oceans
shelf, where the first baroclinic radius Rossby is below 5 km. Vertical coordinates are prescribed by a generalized hybrid vertical coordinate , with finer resolution than that used by Luneva et al. (2015). The model has 74 vertical levels: 35 terrain-following s levels in the upper 300 m, condensing toward surface and bottom, and 40 partial-step z levels below. The s levels allow the resolution of the benthic boundary layer on the shelf with layer thicknesses near the bottom varying from 0.3 m in shallow water to 3 m at 300 m and naturally resolve the topographic slope (Winton et al., 1998). Vertical mixing/viscosity is defined by the k-ε closure with the Canuto et al. (2001)

Journal of Geophysical Research: Oceans
The initial and open boundary conditions for temperature, salinity, currents, and sea surface elevation and initial conditions for sea ice fraction and thicknesses are derived from 5-day mean fields from the global ORCA 1/12°NEMO integrations; see Aksenov et al. (2016) and section 3.2 for details. Our regional model is driven by geopotential tidal forcing with 15 constituents and lateral boundary conditions with nine tidal harmonics from TPX07.2, a 1/4°resolution inverse barotropic tidal model (Egbert & Erofeeva, 2002). Surface forcing is taken from DRAKKAR (2012) forcing set DFS5.1.1 reanalysis 3-hourly mean fields at 0.7°horizontal resolution. The forcing fields are interpolated on the model grid and used to calculate surface fluxes to the ocean and sea ice at each model time step following the CORE bulk formulae (Large & Yeager, 2004). The river runoff forcing data set is a combination of Greenland glacial freshwater runoff (Bamber et al., 2012) and Dai et al. (2009). The latter is interpolated on the model grids by the Total Runoff Integrating Pathways procedure (Oki & Sud, 1998).
The numerical simulations begin on 15 September 1979 with the model being spun up for the first 6 years. The analysis is performed for the following 25 years of the simulations, 1986-2010. No assimilation or restoring is used in the model.

Particle Tracking With Global High-Resolution Model
The global high-resolution eddy-permitting/resolving 1/12°ocean model (horizontal resolution of~3 km in the Arctic Ocean and~8 km globally) is used to drive Lagrangian particles, tracking the propagation of the Bering Strait waters and cascading in the Chukchi Sea (Kelly et al., 2018). The model is a z-level global ocean NEMO model (Madec & NEMO Team, 2016) coupled with the same version of the LIM2 ice model as the regional model described above. It simulates vertical mixing using the turbulent kinetic energy mixing scheme (e.g., Blanke & Delecluse, 1993). This global model does not include tides and, with full-depth z-partial steps, does not resolve the benthic boundary layers and shelf slopes as accurately as the regional model: At 100-m depth, the global model vertical grid is about 10 m and, at 300-m depth, is about 30 m.
For this Lagrangian analysis we use the Ariane numerical particle tracking tool developed at the Laboratoire de Physique des Océans in Brest, France (Blanke & Raynaud, 1997). The method calculates particle positions on the faces of the grid cells and solves the volume continuity equation on a C-grid, preserving local threedimensional nondivergent flow: It is, therefore, capable of accurate volume tracing. The Ariane tool has been extensively used and validated for the Arctic (e.g., Popova et al., 2013) and globally (e.g., Van Sebille et al., 2018).
Model particles are released during the summer of 2001 (Kelly et al., 2018), with 100 particles initially seeded at the ocean surface over a regular 10 × 10 km grid in the Bering Strait. Releases take place every 10 days between the start of June and end of October (a total of 15 releases) to give a total of 1,500 trajectories. All particles are tracked for a full year and allowed to be advected horizontally and vertically, with particle positions recorded daily. The experiment design is the same as that used by Kelly et al. (2018). While having the advantage of being computationally fast with a good level of numerical accuracy in representing pathways of the water masses, the off-line particle tracking methods nevertheless do not include vertical mixing (see, e.g., Kelly et al., 2020, for further discussion and references). However, using the full 3-D velocity field allows us to track details of the mean vertical motion, which are key for the present study. Here we use the eddy-permitting/admitting model with~3-to 4-km resolution in the area of the analysis, which accounts for the horizontal mixing primarily due to eddies. The Rossby radius in the area is~5-10 km (Nurser & Bacon, 2014).
Both the global and regional model configurations are extensively validated in the Arctic against observations of hydrography, satellite dynamical topography and ocean circulation, and sea ice; see Luneva et al. (2015), Aksenov et al. (2016), Armitage et al. (2016), and Kelly et al. (2018Kelly et al. ( , 2020 and Supporting Information S1.
In this paper we focus on results from the regional model with a comparatively coarse horizontal resolution but using high-resolution hybrid terrain following vertical coordinates on the shelf and explicitly resolved tides. Neither of these features are available in the 1/12°global model also used in this study. We expect that, particularly in the Central Arctic, cascading currents will be weak and plumes will be relatively thin in the vertical, and so the resolution of the benthic boundary layers is more important than capturing the eddying of cascading plumes.

Journal of Geophysical Research: Oceans
Moreover, the horizontal resolution of the regional model corresponds to the 1/4°ORCA global ocean model, which is currently used as a state-of-art ocean component for forecasting and climate coupled models (Storkey et al., 2018). Hence, analysis of this regional configuration is of particular interest, to examine the model skill in simulating DWC in a configuration relevant to future climate simulations in the Arctic and globally.

Results
We leave the detailed analysis of cascading fluxes in each of the Arctic domains to the intercomparison "cascading team" project and restrict this study to an overview. Figure 3 shows the maximum of "preconditioning" over 25 years, that is, the maximum salinity of well-mixed waters, defined by Equations 9, 10a, and 10b. These are the regions where cascading waters originate. Typically, these source areas are shallower than 200-m depth. Twenty-five-year mean dense water descending fluxes (only fluxes exceeding 1 m 2 /s are shown) are adjacent to these areas but do not necessarily directly fit the edge of the "preconditioning zone." Note that the extensive spatial coverage of the dense water flux in Figure 3 does not mean that these events are repeated at the same place and the same year. Some of these expansive cascading locations could be the signature of the path of single or multiple cascading events. Animations (see supporting information for preconditioning and cascading flux) show intermittent sporadic events propagating along the slope and then descending when they reach steep slopes. These animations are consistent with the unsteady theory of Huthnance (2009) Table 1).
A preconditioning zone and intensive DWC locations are simulated on the shelf break of the Franz Josef Land Archipelago, corresponding to observed DWC in Location 16 ( Figure 1 and Table 1).
The model predicts the intense dense water descending fluxes originating in St. Anna Trough, located in the Kara Sea (Figure 3b). These dense waters propagate cyclonically along the Severnaya Zemlya shelf break to the Laptev Sea, passing the East Severnaya Zemlya coastal polynya (Aksenov et al., 2010(Aksenov et al., , 2011Bareiss & Görgen, 2005;Dmitrenko et al., 2009) Thus, the model qualitatively predicts DWCs in the areas where cascading has been observed, shown in Figure 1. For quantitative estimates, we chose the 300-m isobath, since our model resolves the benthic boundary layer at depths shallower than this. The coarser z-partial step grid cells deeper than 300 m are unlikely to resolve the veering in the benthic Ekman layer and so not capture the correct balance of different physical factors in a cascading plume.
The along-slope distribution of cascading flux around the 300-m isobath (Figure 4a) is consistent with the spatial map in Figure 3. The strongest fluxes, reaching 3-4 m 2 /s, occur on the Baffin Bay eastern slope, the Irminger Sea, and the Kara and the Barents Seas. Spatial and temporal variability is very strong, and the standard deviations exceed mean values in most places (Figure 4d). The SEFm is the dominant flux in the open-water regions, which are free of ice for most of the time (Figure 4b; note the difference in scale compared with Figure 4a). The DWC is the strongest component of cross-shelf exchange in the Arctic Ocean interior (sectors from 6 to12 in Figure 1, thereafter the "Deep Arctic"). The standard deviation of SEFm is Table 1 Coordinates and Locations of Observed DWCs (Figure 1) and  In the Barents Sea inflow, there are small positive heat flux anomalies and small negative freshwater flux anomalies (i.e., implying a buoyancy flux of the wrong sign, assuming that cascading waters are denser). We attribute this to uncertainties in the definition of the reference temperature/salinity over the contour. Specifically, salinity and temperature in the DWC plumes vary with distance from the bottom, and the thickness of the plume varies in time and space. We take here the 5-year average bottom salinity/temperature as reference values. Mean salinity over the thickness of the DWC plume could be smaller, resulting in errors in the definitions of DWC heat and salt fluxes in regions with high temporal variability of water mass properties, such as the Barents Sea inflow.
Integrated over each sector (14 sectors are shown in Figure 1, and the 15th sector, the Deep Arctic, is the sum of 6-12), time-dependent Ekman fluxes correlated with theoretical values of SEF (BEF), given by ( (16a) and (16b)) and model estimates of the Ekman depths ((15a) and (15b)) SEFm (BEFm) with correlations above 0.95 for SEFs and 0.74-0.97 for BEFs, except for Baffin East slope.
The 25-year mean cascading flux, following our estimates, is the dominant flux of cross-shelf exchange in the Arctic Ocean, reaching 1.3 Sv in the Deep Arctic ( Figure 5a) and a maximal annual mean of 1.6 Sv. For SEF, the 25-year mean flux in the Deep Arctic is 0.85 Sv (65% of DWC; Figure 5b) versus 0.58 Sv of benthic drain (58% of DWC; Figure 5c).
Theoretical estimates of cascading flux, based on Equation 3 and constrained by conditions in Equation 11, exceed the direct model estimates by a factor of 10-30. This is not surprising considering the idealized conditions used to derive the theory (Shapiro & Hill, 1997) based on laboratory experiments (Shapiro & Zatsepin, 1997) and idealized model experiments (Wobus et al., 2011) that consider the constant inflow of dense water on a cone with constant slope. Moreover, the estimate of downslope velocity in Equation 2 is valid for a steady propagating front. In the real ocean, the sporadic and unsteady blobs of dense water, formed during cold outbreaks over relatively short periods of time (2-3 months), propagate over irregular slopes under the influence of Ekman fluxes and geostrophic currents and eddies. The rate of dense water production depends on the conditions on the broad shallow slopes, where the Nof velocity is much smaller.
In the Deep Arctic, SEFs are reduced compared with less ice-covered areas, due to reduced atmospheric exposure. Positive offshore flux is a signature of anticyclonic winds (Figure 5b). The good agreement between magnitudes of simulated flux in the surface Ekman layer and theoretical estimates (SEF and SEFm) and very high correlations show that the expression in Equation 15c for Ekman layer depth in the upper ocean is reasonably applicable and the geostrophic component of currents mostly follows bathymetry contours.
The most frequent month of cascading events over the 300-m isobath (Figure 5d) is calculated as follows: Each winter, from September to May, we define the month when cascading reaches a maximum for each sector and then average this maximum cascading flux over all years. The preferred month of cascading varies with region from late December to early February (Figure 5d). In general, more frequent cascading occurs from middle to late January. The year of maximal cascading shifts toward 2000-2007 in the Deep Arctic ( Figure 5e). This period is characterized by a long-lasting period when anticyclonic atmospheric circulation dominates (Proshutinsky et al., 2015).
The flux integrated over the entire deep Arctic Ocean has very strong synoptic, seasonal, and interannual variability (Figure 6), reaching maximum values of 7 Sv for SEF, 3.5 Sv for cascading, and 1.5 Sv for the BEF. After filtering high-frequency synoptic variability (using 85-day running averages, approximating seasonal means), both SEF ( Figure 6a) and BEF (Figure 6b) theoretical and simulated fluxes are very close. The seasonal peaks in BEF precede the peaks in DWC flux (Figure 6b). However, this may be a signature of enhanced benthic currents in summer and DWC in winter. After filtering the seasonal component (365-day running average; Figure 6c), there is no correlation evident between BEF and DWC. However, it is still not clear whether dense water flows contribute to cyclonic benthic currents or whether a benthic drain enhances cascading. All the yearly fluxes are positive, which is not surprising considering that the sectors are not closed, both laterally and vertically (i.e., only the boundary layers are considered). Moreover, these vertical boundary layers intersect in wintertime as the water column becomes well mixed during preconditioning. Thus, the cascading layer can include both surface and bottom Ekman layers (Figure 4). SEF and BEF layers can also overlap under the conditions of strong winds.
DWC fluxes are relatively small up to 1990: With an 85-day running average, they are less than 1 Sv. After 1990, with the acceleration of summer ice melt, DWC fluxes grow to a local maximum of 2.5 Sv during 1990-1995. This is a period when the cyclonic phase of atmospheric circulation dominates (Proshutinsky et al., 2015). Positive correlations between the seasonal cycles in SEF and DWC are evident for the period 1990-1995. This is followed by a 5-year period of negative correlations. The interannual variability (Figure 6c) does not show any obvious correlations. Different factors, discussed in section 3 (such as ice conditions, river runoff pathways, deflection by wind-driven circulation, and upwelling/downwelling regimes) are, presumably, location dependent and play conflicting roles on the feedbacks between net SEF and cascading.
Finally, both SEF and DWC have positive trends with time of 0.023 and 0.0175 Sv/year, respectively. Both trends may be explained by the sea ice decline, which enhances the transfer of momentum to the water compared with ice covered conditions, and also by the increase of new ice production on shelves, thus enhancing shelf convection.

Cascading in the Canadian Basin
Calculated flux values give a general estimation of the scale of mass and heat exchange by different processes between the shallow shelves and deep ocean. Although such estimation is very useful for understanding the significance of DWC in shelf-basin exchange on the pan-Arctic scale, it does not say much about the specific cascading events contributing to this integral flux. In this section we scrutinize time records of model output 10.1029/2020JC016044

Journal of Geophysical Research: Oceans
at several cross-slope sections around the Canada Basin shown in Figure 7, where cascading has been reported in historical data (Ivanov et al., 2004). A characteristic feature of the upper ocean layer at the margins of the Beaufort Gyre is the inflow of the Pacific-origin water, which enters through the Bering Strait. Thermohaline properties of the Pacific Water change seasonally between being warm and fresh in summer Pacific Summer Water (PSW) and cold and salty in winter Pacific Winter Water (PWW). At the entrance to the Canada Basin deep interior, the Pacific Water inflow splits into three branches: Barrow Canyon in the west (Pickart, Weingartner, et al., 2005), the Central Channel , and Herald Canyon in the east . As will be shown later, the distance of the area of interest from the Bering Strait plays an important role in the DWC features and frequency. These case studies are now considered in some detail.

Beaufort Sea (Mackenzie River) Shelf
The shelf area off the Mackenzie River in the southeastern Beaufort Sea (Sector 7 in Figure 1, Locations 23 and 24) was reported by Ivanov et al. (2004) as the area where the fresh water over the shelf is occasionally substituted by relatively salty upwelling-origin waters. Freezing of this water leads to excessive salinization in the water column, causing surface-to-bottom convection and subsequent motion of the dense plume toward the shelf break. The broad and flat shelf (about 150 km to the steep slope) requires a substantial buoyancy deficit in the dense water to reach the shelf break without notable dispersion, in order to cascade downslope. Direct measurements in winters 1981 and 1988 revealed large areas of surface water to the east of the Mackenzie mouth with salinity exceeding 33 Practical Salinity Units (PSU)-the value typical for 150-m depth in the adjacent deep basin. According to the assessment by Melling (1993), the mean volume transport of this cascade to the halocline layer in the long term is about 0.04 Sv.
The time-depth-distance plots (Figures 8a-8c) show the variations of the modeled near-bottom thermohaline properties for the entire model run. The characteristic feature of this section is two different regimes of the bottom water state. The threshold between these two regimes is in the middle of the record. Before 1998, the water on shelf is fresh and light almost all the time (except for 1988). Starting from 1998 until the end of the record, the frequency of winter events, when salty dense water occupies the entire shelf down to the shelf break between 66-and 164-m depth, substantially increases. The highest salinity/density on shelf, reaching those values in the adjacent basin at a depth between 164 and 362 m, is seen in winter 2008. Note that the winter 2008 followed the second successive summer minimum ice cover in the Arctic Ocean and particularly in the Canada Basin (Perovich et al., 2008). In winter 1988, the model data show a local extremum (on an interannual scale) of bottom salinity and density. This extremum coincides with

10.1029/2020JC016044
Journal of Geophysical Research: Oceans the observed data, cited above. However, the modeled bottom salinity (about 31.5 PSU) is about 2 units less than that measured (over 33 PSU). On the other hand, the fact that the model correctly represents the case of favorable prerequisite conditions for cascading in winter 1988 implies that the tendency of generally increasing DWC potential in the 2000s is correctly captured by the model. The latter corroborates the general projection of intensification of cascading under conditions of increased seasonality of the Arctic sea ice cover (Ivanov & Watanabe, 2013).

Barrow Canyon
There is observational evidence of dense water downslope spilling at the southern flank of the Barrow Canyon (Cases 57 and 58 from Ivanov et al., 2004, here shown at Locations 21 and 22 in Figure 1 and Table 1). The relatively high density of this water mass permits the occasional penetration below the upper halocline in the Canadian Basin (Weingartner et al., 1998). Examples confirming this were reported in observations by Pickart, Weingartner, et al. (2005) and Hirano et al. (2016). DWC in the Barrow Canyon may be considered a rather rare event, at least in the past. For instance, even a dedicated survey in 1986-1987, specifically aimed at measuring the DWC down the canyon, failed to capture any signs of the process (Aagaard & Roach, 1990).
The time-depth plot of calculated near-to-bottom thermohaline properties for the entire model run is presented in Figures 9a-9c. Regular seasonal cycling of temperature and salinity in the water column is clearly visible down to the 80-m depth. Such regularity of a strong seasonal signal is not observed on the previously discussed Mackenzie River shelf (see Figure 8). The probable cause of this is the stronger impact of the alternating PWW and PSW, due to the proximity of the Barrow Canyon to the Bering Strait.
The shelf area between Cape Lisburne and Point Barrow along the Alaska coast, including Barrow Canyon, is the domain of the Barrow Coastal Polynya (Winsor & Björk, 2000). In winter, the ocean surface is rapidly cooled to the freezing point. Subsequent brine rejection and associated salinification of the water column inside the polynya provide the prerequisite conditions for DWC. However, as previously mentioned, the rare occurrence of cascading events in the past points to some concurrent processes blocking DWC development despite generally favorable background conditions. According to Hirano et al. (2016), a potential process is the upwelling of warm water forced by the persistent northeasterly winds. Hirano et al. (2016) argue that a rapid salinity drop after the termination of upwelling typically prevents preconditioning of shelf convection by the process suggested by Aagaard et al. (1981), thus explaining the rarity of cascading events in this region.
The time record of model salinity (Figures 9b and 9e) shows salinity exceeding 32.8 PSU over the shelf between 16 and 70 m during the winters of 2000. In 2000, the time intervals of extremal salinity did not exceed 7-10 days. In 2007, the residence time of extremely salty water on shelf was between mid-March and the end of June, thus pointing to a high probability of DWC. The vertical distribution of model potential density and salinity at the transect across the Barrow Canon in February and April 2007 is shown in Figures 10a and 10b. The formation of a compact dense "pool" of water near the bottom and to the north of the canyon in February is probably the result of polynya formation. The model does not resolve the details of polynyas well; however, it shows repeatable appearance of ice concentrations below 0.9 (Figures 10c and 10d) at the edge of the coastline and near Barrow canyon in the period from January to April 2007. Polynya formation is further corroborated by the freezing point temperature and local salinity maximum from surface to bottom. By mid-April, the most salty and dense water accumulates along the center of the canyon, acquiring potential for further motion toward the mouth of the canyon as described by Pickart, Weingartner, et al. (2005).

Western Chukchi Sea
The vast and shallow shelf adjacent to the Bering Strait on the Arctic side is the pathway of the central branch of spreading Pacific-origin waters . Observational evidence of a cold dense water "pool" sitting on the outer shelf in this region is provided by Ivanov et al. (2004); see Case 20 in Figure 1 and in Table 1. The main difference in terms of dense water formation at this site, when compared to the previous two, is the very flat, extended, and shallow seabed. This fact explains why the preconditioning of DWC in this region substantially differs from those described above: They had initial replacement of  Figures 8 and 9), is the frequent occurrence of separate salty/dense periods in the upper shelf region (excluding the shelf break). This fact may be interpreted as the signature of localized surface-to-bottom Journal of Geophysical Research: Oceans thermohaline convection, forming solitary dense water pockets that decay before reaching the shelf break. The reason behind this decay is the shallow slope and associated small Nof speed, implying that moving dense water patches have time to decay by horizontal and vertical mixing before they descend into the deep basin. However, there are three instances (1996, 1998, and 2008), when DWC into the deep basin down to a depth of 150-180 m can be anticipated on the basis of the model. The vertical transect of salinity and potential density across the shelf and the upper slope of the central Chukchi Sea is shown in Figure 11d. The pattern of density contours with maximum values on shelf near the bottom shows cascading of cold salty water from 60-80 to 100-120 m along the slope. The local maximum of salinity within the entire depth of the water column at 40-60 m indicates the source area of dense water formation, possibly linked with polynya formation.
Other evidence of DWC in this area was obtained by Conductivity-Temperature-Depth (CTD) measurements in July 2002 and in September 2004 (Pickart et al., 2009;Pickart, Weingartner, et al., 2005). Vertical transects of potential temperature and salinity, the latter roughly coinciding with the pattern of potential density for water near the freezing point, based on July 2002 data and on the model results for the same date, are presented in Figure 12. Both the observation and model results show a plume of cold and relatively fresh water that descends from the shallow shelf (~20-50 m) down to the deeper water (~100-120 m). Model simulations emphasize increased vertical mixing in the descending plume (Figures 12c and 12d), with diffusivity exceeding the background value by an order of magnitude. The plume meets the salty and warm Atlantic origin water at the intermediate depth. This warm and salty water sets the vertical limit for the descending plume, which then detaches from the seabed and spreads to the east along isobaths within the 80-to 120-m layer (Pickart, Weingartner, et al., 2005; Figure 12a). The modeled vertical bound of DWC spreading in this area is rather stable. The lower limit, as was captured by the model in March 2008 (see Figure 11d), is the same as was measured in summer 2002.

Journal of Geophysical Research: Oceans
Animations of water properties-temperature, salinity, density, and vertical diffusivity along 165.5°W, located just one grid point apart from (Pickart et al., 2009) transect-are shown in Movie S3 and illustrate the subsequent stages of cascading process across the shelf: (i) uplifting of the warm salinity tongue (December 2001), (ii) formation of the blob of cold and saltier waters below the fresh and cold subsurface layer at~73°N, and (iii) the ventilation of the isopycnal of σ = 26 kg/m 3 above the outer shelf. In July, a large tongue of cold water propagates slantwise below the newly formed pycnocline between the isopycnal layers of σ = 26 − 27 kg/m 3 , while mixing and cascading of denser waters is still evident along the shelf break. In December 2002, water structure nearly repeats the transect in December 2001.
Examples of the reproduction of observed cascading events, at the same location and time (see also Figure S3), gives us confidence that, despite the modest horizontal resolution, this model performs satisfactorily in the direct simulation of cascading. Thus, estimates of cascading fluxes, discussed previously (see section 5), can be seen as reliable first-order estimates.

Cascading and the Pathways of Pacific Waters
To obtain more spatial details on the DWC event in July 2002 in the Chukchi Sea and to examine the spread of Pacific Water in the Canadian Basin, we analyze the pathways of trajectories originated in the Bering Strait with the use of the output of the 1/12°global model (see section 3.2). Figure 13 shows the 1/12°model temperature (a) and salinity (b) fields in the Canadian Basin at~100 m, as well as the depth of trajectories of the Pacific Water tracers. The temperature and salinity fields are evidence of the cold (at about freezing point) and saline (up to 33.5) waters cascading from the shelf following the canyons (Herald, the Barrow Canyons, and the smaller canyons, e.g., a deep extension of Central Channel) and then spreading eastward (clockwise, anticyclonically) and offshore in the Canada Basin. This pattern has been suggested by the observational transects completed on the western Chukchi shelf (e.g., Pickart et al., 2009Pickart et al., , 2010. The Barrow Canyon also has a signature of westward propagating cold waters at 100 m. 10.1029/2020JC016044

Journal of Geophysical Research: Oceans
The particle tracking experiments show Pacific water crossing the shelf in several branches during 3-4 months. The water meanders on the shelf and then descends to a depth of 160-200 m in the deep Canada basin (Figures 13c and 13d), moving clockwise. Descending water was formed during winter, thus the deepening stage would correspond to middle and late summer (July-September). Note that the deepest trajectories propagating across the Chukchi shelf, and nearly corresponding to the depth of the shelf at these locations (Figure 13c), coincide with the maximum of MLD salinity seen in Figure 7 (see also Movie S1). The locations of cold and saline water at the depth of 100 m in July 2000 (Figures 13a and 13b) correspond well with the preconditioning shown in Figure 7. Thus, both regional and global models, with different horizontal and vertical resolutions, demonstrate similar results. The question that arises is does the upwelling observed at the 166°W transects and at the other locations support cascading and mixing of waters of Pacific and Atlantic origin? If so, is it a repeating feature? Autumn upwelling, associated with storms crossing the western Chukchi and Beaufort Seas, is a common feature of the Canadian Arctic dynamics (Pickart et al., 2009). Indeed, all the Lagrangian trajectories originated on the surface on the Bering Strait, tracing the waters with very low salinity. However, after water modification due to brine rejection and cooling, or mixing, the trajectories become dense enough to reach depths below 160-200 m. To support this conclusion, we present an animation (Movie S4) of potential temperature, density, and vertical diffusivity, simulated with the regional pan-Arctic model, at a transect along the 100-m isobaths (contour and key locations are shown in Figure 7). This animation demonstrates yearly repeatable upwelling of warm dense waters onto the shelf during October-November. At some locations, strong diapycnal mixing, exceeding the background by 3-4 orders of magnitude, results in the cooling of the mixed waters and the ventilation of isopycnals.
The Barrow and Herald canyons are suggested to be the main gateways for the passage of Pacific waters (Pickart et al., 2009;Spall et al., 2008; Figure 1), which is supported by modeled trajectories (Figure 13).
To examine the interaction of upwelling and cascading in more detail, we consider time slices at transects along Barrow and Herald Canyons, shown in Figure 7.
Strong upwelling of warmer and saltier waters is evident in the Barrow Canyon in February 2002 (Figures 14a and 14b), penetrating up to a depth of 50 m and distance of 200 km from the base of Barrow Canyon. These waters are separated from the cold (T~−1.9°C, −1.5°C) and fresh surface (S < 29.5 − 30 PSU) waters by a strong halocline. Just a week later, cooling and ventilation start at the shallowest part of the canyon, accompanied by strong vertical mixing in the intermediate layers (Figures 14c and  14d). In early March (Figures 14d and 14e), mixed waters on the shelf become colder and fresher but less dense (σ ≈ 26 − 26.5 kg/m 3 ) compared with the stratified waters in the deeper part of the basin (at 100 m; σ ≈ 26.5 − 27 kg/m 3 ), while cooling by 0.3-0.5°C of the layer (σ ≈ 26.5 − 27.5 kg/m 3 ) at the shelf break  (Figures 14g and 14h), with a cascading plume penetrating to a depth of 150 m on the steep slope and a vertical diffusivity in the plume exceeding 10 −3 m 2 /s (Figure 14h). Cooling due to entrainment and vertical mixing is extended up to a depth of 180 m, corresponding to the σ ≈ 27.5 kg/m 3 isopycnal and S ≈ 33.5 PSU isohaline, which is 0.5 − 0.8 PSU more saline than the mixed waters on the shelf. Note that, according to Movie S3 and Figure 13, these waters propagate westward being replaced by the cascading shelf waters. Two weeks later (2 June; Figures 14i and 14j), local upwelling uplifts warmer saltier waters back to a depth of 80-100 m, accompanied by strong vertical mixing (diffusivity 10 −2 m 2 /s). At this time a strong summer seasonal halocline starts to form, preventing the mixing between subsurface and benthic layers. Ten days later, in the middle of June (Figures 14k and 14l (Figure 15a; at the distancẽ 120-150 km) forms a blob of water denser than σ ≈ 27 kg/m 3 (Figure 15c). A cold front propagates downslope and reaches a depth of 100 m by the middle of June. The propagation of well-mixed waters and shape of the cascading tongue is highlighted by the high vertical diffusivity contour (blue solid lines in Figures 15b,  15d, 15f, and 15h). At this time, the benthic and surface mixed layer are separated by a newly formed halocline due to melting of ice at the surface, which is in accordance with Movie S3.
Thus, our simulations support the hypothesis that upwelling and vertical mixing with uplifted Atlantic Waters are likely mechanisms amplifying the formation of heavy dense waters and presumably increasing ventilation of intermediate and deep waters. Deeper (500-1,000 m) and stronger (>1 m 2 /s) DWC fluxes are simulated in the Chukchi Plateau (Figures 3a and 3b).

Summary
The process of DWC, a specific type of topographically trapped flow originating on the shelf, plays a particularly important role in the Arctic Ocean. It is the principal mechanism of ventilation whereby intermediate and deep waters are replenished. DWC events can also contribute to the export and storage of carbon from the shelf as a constituent of the "shelf carbon pump" (Holt et al., 2009) and have been shown to contribute to the modification of Atlantic water (Ivanov & Golovin, 2007). The export of water from the shelf can induce a compensating flow of the deep water onto the shelf (Ivanov & Golovin, 2007;Kämpf, 2005), with the potential to provide a supply of nutrients to the upper ocean and to affect primary production. DWC events are hard to observe due to their sporadic nature and the technical difficulties of obtaining wintertime observations under sea ice. Until now, to our knowledge, only 27 locations of cascading in the Arctic Ocean have been reported in the literature. There have also been few regional modeling studies examining cascading events in the Arctic Ocean: Previous examples include Ivanov and Golovin (2007), Ivanov et al. (2015), Marson et al. (2017), , and Magaldi and Haine (2015).
In this work a pan-Arctic ocean-ice coupled model, with moderate horizontal grid spacing (10-15 km) and fine vertical discretization (0.3-3 m) resolving shelf slopes, is used to examine the ability of 3-D general circulation models to reproduce these observed cascading events and to diagnose their importance on a pan-Arctic scale. Given the choice between fine vertical and fine horizontal resolution for this simulation, we here opted to focus on fine vertical resolution to capture the dominant benthic layer processes at the expense of topographic detail and their interaction with the mesoscale eddies. Using a 30-year simulation, we analyzed model output to detect the most favorable locations of DWC (see Figure 3). We found that most of the locations of the observed DWC events (cf. Figure 1) are reproduced by the model.
The shelf-slope areas around the Canadian Basin provided the focus for a detailed comparison of model results with observations. The present-day sea ice seasonal cycle in the North American Arctic shelf seas, with extended summertime ice-free areas and complete ice coverage in winter, favors pre-conditioning of 10.1029/2020JC016044 Journal of Geophysical Research: Oceans DWC events. The results reveal that the model, despite its relatively coarse horizontal resolution, has the capability to simulate DWC events, or at least intense shelf convection, in most of these cases. Therefore, the trend of increasing cascade intensity under conditions of enhanced seasonality of the Arctic sea ice cover, revealed in the model simulation, may be taken as a plausible projection for the near future.
Examining the conditions and preconditions of DWC in the western Chukchi Sea (Figures 12-14), we found that active DWCs are preceded by autumn upwelling and uplifting of the Atlantic warm and saltier waters on the shelf. Upwelling associated with storm tracks is a frequently observed feature in this region (Pickart et al., 2009). The simulation shows the presence of active diapycnal mixing and cooling of uplifted waters, resulted in the formation of denser waters, able to descend to depths as large as 160-200 m. Therefore, these simulations support the hypothesis by Aagaard et al. (1981). The Pacific water Lagrangian experiments (see sections 3.2 and 5.4 and Figure 13) suggest a counterclockwise initial spread of the cascading waters next to the shelf slope, with a change of pathways to clockwise in the deep Canada Basin, at about 160-to 200-m depth. The time scale for the winter Pacific water to cross the shelf and the deep Arctic Ocean is about 3-4 months.
We provide the first model-based estimates of the various components of shelf-ocean transport fluxes in the deep Arctic Ocean and seas of the Northern Atlantic Ocean connected to the Arctic ( Figure 6): DWC and SEF and BEF. Specifically, the 25-year mean exchange in the "Deep Arctic" (Regions 6 to 12 in Figure 1) across the 300-m isobath is dominated by a DWC flux of 1.3 Sv compared with 0.85 Sv for SEF and 0.58 Sv for BEF. Thus, the simulated DWC flux is comparable to the total Bering Strait inflow (0.8-1.3 Sv). Finally, in most sections along the 300-m isobath, the most intense cascading happens during mid-January to February ( Figure 6). This would imply that the cascading events detected in observations, mostly in ice-free conditions in September-October, underestimate the actual intensity of such events.
Simulated cascading fluxes bring colder waters to the ocean interior, with 25-year mean maximal negative vertically integrated heat flux anomalies (per unit length along the 300-m isobath) exceeding 4-10 × 10 6 W/m in the Irminger Sea, Barents Seas, and near the North-East Greenland shelf break (Figure 4g).
Interannual variability of DWC flux in the "Deep Arctic" is strong, with relative anomalies reaching 50% (~1.6 Sv in 1993vs. 0.8 Sv before 1990and in winter 2004 (Figure 7c), and does not correlate directly with the phase of Arctic Ocean Oscillation index (Proshutinsky et al., 2015). However, the SEF and DWC have positive trends with time, 0.023 and 0.0175 Sv/year, presumably caused by summer ice decline (Figures 7a and 7b).

Data Availability Statement
Baseline model is a regional configuration of NEMO (Nucleus for European Models of the Ocean), at Version 3.6 stable (Madec & NEMO team, 2016). Model code (Version 3.6 stable) is available from the NEMO website (www.nemo-ocean.eu). Model output, used in this study, is now under archiving process at British Oceanographic Data Centre (BODC, https://www.bodc.ac.uk/).