1 INTRODUCTION

Breaking internal waves can alter temperature and salinity gradients, both horizontally and vertically. In doing so, they drive water-mass exchange and contribute significantly to the diapycnal mixing required to maintain the global ocean circulation (Garrett and Munk, 1979; Wolanski et al., 2004; Xu et al., 2013, 2021). Coverage of mixing estimates either from direct microstructure observations or from parameterization is now becoming sufficient to provide global patterns of turbulent mixing with both spatial and temporal variability (parameterization: e.g., Whalen et al., 2012; observations: e.g., Waterman et al., 2013; Rippeth et al., 2015; Kunze, 2017). Nevertheless, there are gaps in the coverage, with one of the biggest gaps in the Arctic Ocean. This is due to the patchy, episodic, and scarce observations of mixing in the Arctic Ocean (e.g., Sundfjord et al., 2007; Kunze, 2017; Whalen et al., 2018). Consequently, the important role that mixing plays in water mass structure and circulation in the Arctic Ocean is not well known. Addressing this shortcoming and understanding mixing properties of the Arctic regions are clearly necessary for advancing both Arctic Ocean and global climate research. Therefore, our work will focus on two specific issues: the spatial pattern and temporal variability of mixing in the Arctic Ocean and the main driving mechanisms on these mixing properties.

The water column of Arctic Ocean is typically divided into three primary layers, including a near-surface layer with temperatures near freezing most of the year, a cold halocline, Atlantic Water, and a deeper layer (Rudels et al., 1996).The Pacific sector could have more layers due to the intrusion of the Pacific Water. The warm and salty Atlantic Water (AW, 200–900 m) is separated from the freezing near-surface layer by halocline (HL, often occurs in 50–200 m), resulting in strong stratification (e.g., Aagaard et al., 1981; Steel et al., 1995; Shimada et al., 2005; McLaughlin et al., 2009). This makes the Arctic Ocean different to the temperate and tropical ones. However, the stratification of the Arctic Ocean is highly variable, both temporally and horizontally (Rudels and Quadfasel, 1991; Pemberton et al., 2015). The relationship between such variable stratifications and the variability of mixing is still unclear.

The Arctic climate and environment have changed significantly in recent decades, especially sea ice, which has declined in extent, thickness (e.g., Rothrock et al., 2008; Stroeve et al., 2012; Lindsay and Schweiger, 2015), and the percentage of multiyear ice (e.g., Comiso, 2002; Lei et al., 2016; Bi et al., 2020). Simultaneously, sea ice drifted speeds have increased (Hakkinen et al., 2008). This points out that the Arctic Ocean entered a new state (Carmack et al., 2015; Polyakov et al., 2017). However, it is unclear how mixing will be influenced and modulated in this lower sea ice state.

The primary mixing mechanisms in the Arctic are wind, internal tides (e.g., Halle and Pinkel, 2003; Guthrie et al., 2013; Wang et al., 2021; Zhao et al., 2021), and double diffusion (e.g., Robertson, 1995; Ruddick and Gargett, 2003; Timmermans et al., 2008; Shibley et al., 2017; Fer et al., 2020). In the lower latitudes, internal tides can propagate a long-distance with low-mode and redistribute the mixing (e.g., Liu et al., 2017; Guo et al., 2018). Much of the Arctic is poleward of the critical latitude for the generation of a freely propagating internal tide, but there are two main kinds of internal waves that dissipate and drive the mixing in the Arctic: tidal lee-waves (e.g., Rippeth et al., 2017; Fer et al., 2020; Schulz et al., 2021a) and bottom-trapped internal waves (e.g., Falahat and Nycander, 2015; Hughes and Klymak, 2019). Although the tidal mixing is found over a much more limited area of the Arctic Ocean, the averaged heat flux from areas with tidal mixing is comparable to that supported by double diffusion (Fer et al., 2009, Fer et al., 2020). Near-inertial internal waves can be generated by the wind forcing at the atmosphere-ocean interface and can transfer wind energy downward until the breaking depth of internal wave (e.g., Hebert and Moum, 1994; Cole et al., 2014; Fer, 2014; Wang et al., 2014). Some works have provided the evidence of the enhancing near-inertial internal waves with decreasing sea ice (e.g., Lenn et al., 2011; Dosser and Rainville, 2016; Polyakov et al., 2020b). On the other hand, double diffusion and intrusions could play a bigger role in some regions with weak tide- and wind-induced mixing. Compared with other oceans, the Arctic Ocean is always viewed as a quiescent, weakly mixed ocean due to low turbulent energy levels (Rainville and Winsor, 2008). Mixing of the Arctic Ocean has been observed in some previous studies, and showed a highly variable pattern, reflecting the dominance of various mechanisms. For example, diapycnal diffusivities were found to be O(10^-5–10^-4) m²/s in Yermak Plateau (Fer et al., 2010), O(10^-6–10^-4) m²/s in Canadian Basin (Guthrie et al., 2013), and O(10^-6–10^-4) m²/s in the Amundsen Gulf (Chanona et al., 2018), indicating relatively weak mixing and high variability.

In the Arctic Ocean, near-inertial waves can also be generated by the motion of drifting sea ice or icebergs. The absence or presence of sea ice can profoundly affect the capability of the kinetic energy of nearinertial waves in driving vertical mixing (Lenn et al., 2011). As the sea ice extent declines, direct interaction between sea surface and wind is increasing (Manucharyan et al., 2016; Armitage et al., 2020). This is exacerbated by the recent enhancement of windstorms and sea ice drift speeds. Evidence underscores the growth of near-inertial motions in early winter after freezing (Merrifield and Pinkel, 1996; Halle and Pinkel, 2003). Dosser and Rainville (2016) suggested the increase in near-inertial waves through estimates of the near-inertial wave amplitudes in the Canadian Basin, the most significant decline region of sea ice extent. Thus, it is reasonable to anticipate a future enhancement of turbulent mixing in the upper ocean. However, some prior, direct microstructure observations portray continued quiescence in the Arctic Ocean interior (Lincoln et al., 2016). Whether the increasing internal wave generate stronger turbulence or not, and what their impacts on the changes of turbulent mixing and turbulent heat flux in the Arctic remain to be answered.

Stratification is also a critical factor in turbulent mixing. The sea ice will be melted when the heat is transferred from the warm AW to the typically cooler surface mixed layer (SML) by mixing. However, the mixing is usually suppressed by strong stratification. The impediment to this mixing is the cold halocline layer (CHL, Lenn et al., 2009). The recent reductions in sea ice coverage and a thinner CHL with weaker stratification will enhance the transfer of AW heat into the upper layer and significantly decrease the rate of sea ice formation in winter. This was found in the Eurasian Basin (Polyakov et al., 2017; 2020a). On one hand, a viewpoint suggested that this represents an essential step towards a new Arctic climate state (Carmack et al., 2015; Polyakov et al., 2017). On the other hand, McLaughlin et al. (2004) suggested that a strengthened stratification resulting from relatively fresh water from the Pacific Ocean at mid-depth (60–220 m) would limit the exchange between SML and AW. This view mainly focused on the Pacific sector. Chanona et al. (2018) suggested that turbulent mixing might still be restrained by strong stratification. They found that the strong stratification over the slope of Canadian Basin played an important role in mixing. It is still unknown that whether the mixing properties in the Atlantic sector and the Pacific sector are similar or not.

In recent decades, some studies focused on the mixing trend in the Arctic Ocean interior due to climatic variability of internal waves and stratification and the implications of these trends on the heat budget and sea ice predictions (e.g., Lique et al., 2014; Chanona et al., 2018). However, it was not clear whether internal waves significantly contribute to the mixing in the Arctic Ocean. A possibility to address this issue and discuss the mixing and its modulations in the Arctic Ocean are provided by Ice-Tethered Profiler (ITP)'s data set.

The spatial pattern and temporal variability of diapycnal diffusivities and dissipation rate in the Arctic were examined in this paper, and we reviewed some possible modulating factors for turbulent mixing, including wind, sea ice, stratification, and near-inertial waves. This paper is organized as follows. In Section 2, we described the ITP data set and our method of inferring the associated wave-driven turbulent diffusivity and kinetic energy dissipation rate. Our results, including horizontal maps, vertical structure, and seasonal and long-term variability of turbulent mixing of the Arctic, are shown in Section 3. In Section 4, we discussed the apparent roles of wind, sea ice, stratification, and near-inertial waves, which underpin the observed spatial and temporal patterns. In Section 5, the spatial pattern and temporal variability of heat flux are discussed. In addition, in Section 6, we summarized our works and conclusions.

2 DATA AND METHOD 2.1 Ice-tethered profiler (ITP)

The primary dataset in this study was ice-tethered profiler (ITPs). The ITP data set was used to examine diapycnal diffusivities and dissipation rates over several main regions of the Arctic, including the Nansen, Amundsen, and Canadian Basin along with other seas. The data set covered 15 years from winter of 2004 to fall of 2019 (Fig. 1). ITPs are anchored in perennial ice floes by a surface buoy, and report potential temperature and conductivity at 1 Hz, as they profile from 7–750-m depth. The dataset we used in this paper is ITP level 3 product, which is provided by the ITP group and the Woods Hole Oceanographic Institution (www.whoi.edu/itp; Krishfield et al., 2008). Processing steps include identification and removal of unphysical data, corrections for pressure and sensor response behavior. The vertical sampling interval of the ITPs after averaging is 1 m. Typically ITPs measure the water column twice a day, with profiles beginning at 00:00 and 06:00 UTC.

We examined 79 798 ITP potential temperature and salinity profiles, discarding 18 078 invalid profiles (data missing or incomplete) and 34 276 profiles due to the existence of double-diffusive staircases (Fig. 1). Double-diffusive staircases are a widespread phenomenon, and generally common in the central Canadian Basin, Amundsen Basin, and Fram Strait (Timmermans et al., 2008; Bebieva and Timmermans, 2017; Shibley et al., 2017; Polyakov et al., 2019). In this paper, we focus on the region where the fine scale parameterization method is robust—most of the Eurasian basin and the low latitude region of the Canadian Basin. These available profiles were used to estimate diaypcnal diffusivities and dissipation rates in the Arctic Ocean. The temporal coverage of these data spans 15 years, with over 1 000 usable number of profiles in each year, with the exceptions of 2005 (105 profiles), 2006 (182 profiles), and 2018 (609 profiles).

2.2 Fine-scale parameterization

To estimate diapycnal diffusivity k_z and turbulent kinetic energy dissipation rate ε, a fine-scale parameterization was implemented. This method has been widely applied to measure turbulent mixing in global ocean, and applied to the Arctic in recent years (e.g., Fer et al., 2010; Guthrie et al., 2013; Lique et al., 2014; Kawaguchi et al., 2016; Chanona et al., 2018). Chanona et al. (2018) have verified the applicability of the strain-based fine-scale parameterization method based on the CTD dataset in the sub-region of Arctic Ocean, which motivated us to reveal a more comprehensive view of mixing properties in the whole Arctic Ocean, using the ITP data with a broader spatial coverage. This method is based on the hypothesis that energy can be transported from large scales to small scales, at which waves break due to shear or convective instabilities through weakly nonlinear interactions between internal waves (Kunze et al., 2006).

The dissipation rate ε can be expressed as

where N₀=5.24×10^-3/s and ε₀=6.73×10^-10 W/kg, according to the Garrett-Munk (GM) spectrum (Garrett and Munk, 1979). is an averaged stratification of each vertical profile segment of interest, which can be calculated by N²=-gdρ/ρdz. R_ω represents the shear/strain variance ratio. The global averaged R_ω is 7 (Kunze et al., 2006), which is generally used in this method (Whalen et al., 2018, You et al., 2021), and some works have used 7 to estimate the dissipation rate in the Arctic (Dosser et al., 2021). Considering R_ω in the Arctic is larger than that in lower latitudes (for example, R_ω are 4–12 in Guthrie et al. (2013); R_ω are 3–50 in Fine et al. (2021)), in this work R_ω is fixed at 11 (Fer et al., 2010). The estimated dissipation rates are compared with the observations in Appendix. The function h(R_ω)accounts for the frequency content of the internal wave field and is calculated according to

The function L(f, N) corrects for a latitudinal dependence, where f is the local Coriolis frequency, f₃₀ is the Coriolis frequency at 30°, and N is the vertically averaged buoyancy frequency for the segment.

where 〈ξ_zGM²〉 and 〈ξ_z²〉 are the strain variance from the Garrett-Munk (GM) spectrum (Gregg and Kunze, 1991) and observed strain variance, respectively. The 〈〉 indicated integration over specified vertical wavenumbers of internal wave fields on each profile. Strain ξ_z was calculated from each segment as

where N² was derived from 1 dbar observed potential temperature, salinity, and pressure, N_ref² was fitted to the observed N² profile using a smooth piece-wise quadratic fitting with a window of 20 m. By applying a fast Fourier transform (FFT) on 75% overlapping 128-m segments along each vertical ξ_z profile, we computed the spectra S_str(k_z) and integrated to determine the strain variance. According to the strain spectrum of the Arctic, k_zmax and k_zmin corresponded to 0.1 and 0.02 cpm, respectively (Chanona et al., 2018). Substituting 〈ξ₂^z〉 into Eq.1 ultimately yields 24-m resolved vertical profiles of ε. Here the integrated constant 0.2 was used to estimate the mixing based on the Argo data (Whalen et al., 2015). Some research found that changing this constant only has a minor effect on mixing amplitudes and does not affect our conclusions (Wu et al., 2011).

The turbulent dissipation rate ε has been compared with direct observations, and the error of estimated ε has been discussed in Appendix. The diapycnal diffusivity k_z was estimated from the dissipation rate ε using the Osborn relation (Osborn, 1980) in the regimes of turbulent mixing (Dosser et al., 2021). The mixing efficiency Γ is set as 0.2 in this study (Gregg et al., 2018), which is generally used in this method

Double-diffusive staircases with O(1–10)-m vertical length scales are widespread in deep basins of the Arctic Ocean. The presence of double-diffusive staircases in the profiles will pollute the strain rate estimates from this fine-scale parameterization method. However, the double-diffusive intrusion scale is generally smaller than 5-m scale of the strain calculation (Guthrie et al., 2013). Double-diffusive staircases will interfere with some results of estimated mixing and vertical heat fluxes. Thus, we checked the temperature profiles and discarded profiles with obvious double-diffusive staircases (e.g., Fig. 2a & d) to minimize this effect. Here, we defined doublediffusive staircases as the dθ/dz much less than (Timmermans et al., 2008). The double-diffusive staircases were identified when the difference of dθ/dz between two adjacent data points in each profile is less than a threshold value (such as dθ/dz < 0.08 . This threshold value was determined to be most appropriate for each profile (Shibley et al., 2017; Polyakov et al., 2019).

3 MIXING IN THE ARCTIC OCEAN 3.1 Arctic ocean stratification

Two representative potential temperature and salinity profiles taken from the Canadian Basin and Eurasian Basin, representing the Pacific (Fig. 3a) and Atlantic (Fig. 3b) sectors, respectively, show clear differences for the two sectors. In the Pacific sector, fresher Pacific water intruded at 50–100 m, strengthening the stratification of the HL and deepening the maximum heat core. Moreover, in the Atlantic sector, water column was mainly composed by a cold and fresh near-freezing surface water, a HL at ~20–150 m and a warm and salty Atlantic Water (AW) at ~300 m with a potential temperature maximum. In the Pacific sector, the averaged HL depth was 20–200 m and averaged AW depth was 200–650 m, whereas, in the Atlantic sector, the averaged CHL depth was thinner 20–150 m and averaged AW shallower 150–550 m.

3.2 Domain-wide statistical distribution

The mixing varies spatially and temporally all over the world. First, the basin-wide statistics of two mixing parameters, diapycnal diffusivity (k_z) and dissipation (ε), and their vertical profiles are presented. Subsequently, the spatial and temporal distributions of these parameters are investigated, both for the entire Arctic Ocean and for the Atlantic and Pacific sectors. Finally, the two parameters are compared between the two sectors.

Statistical distributions of the diapycnal diffusivities and dissipation rates in the entire domain were considered as the primary important feature of mixing for the entire Arctic Ocean, including both the Atlantic and Pacific sectors (Fig. 4). Since the Fram Strait differs from most of the rest of the Arctic Ocean, with extremely strong turbulent mixing by internal tides, it was not included here. The averaged diffusivity and dissipation rate for the entire Arctic were approximately O(10^-5) m²/s and O(10^-9) W/kg (geometric mean), corresponding to standard deviations of O(10^-4) m²/s and O(10^-7) W/kg, respectively. Here, the geometric mean was chosen to avoid the error caused by several extremely large or small values. These diffusivities and dissipation rates indicated a significant high inhomogeneity, both spatially and temporally. Moreover, mixing in the Arctic Ocean, except for the Fram Strait, is weaker than in mid- and lower latitudes. As the depth-averaged results at the whole estimated depth (~100–650 m), diapycnal diffusivities in the Atlantic sector and Pacific sector were 4×10^-5m²/s and 8×10^-6 m²/s, respectively. These differ by more than 0.5 order of magnitude, showing an obvious disparity in mixing properties. The differences of dissipation rate ε were relatively smaller, 1×10^-9 W/kg in the Atlantic sector and 7×10^-10 W/kg in the Pacific sector. Compared with middle or lower latitude diapycnal diffusivities and dissipation rates, the diapycnal diffusivities and dissipation rates of Atlantic sector were on par with them and only the Pacific sector was significantly lower (e.g., Kunze, 2017; Whalen et al., 2018), some details of the spatial pattern and mechanisms will be discussed later. Differences in diffusivities and dissipation rates were hypothesized to result from the different stratification and water masses conditions in the two sectors. In the mid-2010s, the stratification of the upper layer in the Atlantic sector, especially in Eurasian Basin, has been found to keep the weakening trend (Polyakov et al., 2017; Polyakov et al., 2020; Schulz et al., 2021b). In the Pacific sector, the stratification was strong with thick CHL (Toole et al., 2010; Jackson et al., 2012) in early 2010s, and then it was decreasing slightly. The declining sea ice trend was significant in the summer for the Pacific sector, but in winter for the Atlantic sector. The depth of AW in the Pacific sector was deeper than that in the Atlantic sector (Fig. 3). These differences indicated some potential factors leading to the mixing discrepancy between the two Arctic Ocean sectors.

3.3 Horizontal spatial patterns of Arctic Ocean mixing

The horizontal maps of diapycnal diffusivities and dissipation rates in the Arctic Ocean (Fig. 5) were bin averaged into 1° latitude × 3° longitude grid cells of water columns in an upper layer and a deeper layer. Here the upper layer was from 100 (minimum depth estimated) to 200 m, and the deeper layer is 200–600 m, which was fixed according to the averaged profile in entire region (Section 3.1).

A significant difference between the Pacific and the Atlantic sectors for k_z and ε was found at the depths of both two layers. k_z and ε in the depth range 100–200 m (Fig. 5a & c) in the Atlantic sector were higher than these in the Pacific sector by roughly an order of magnitude. k_z and ε at these depths were highest in Fram Strait, second highest in the Amundsen Basin and lowest in the Canadian Basin. Moreover, k_z in some areas always covered by ice were smaller than summer ice free areas, which was consistent with previous research (e.g., Rippeth et al., 2015; Zhong et al., 2018). In the deeper layer, the horizontal maps of k_z and ε were similar to that for the depth of 100–200 m, but with smaller differences between both sectors. k_z and ε strengthened significantly in a portion of the Canadian Basin that was always covered by sea ice. Meanwhile, the topography of this region is rough due to the presence of Canadian Basin slope, which might be able to generate the bottom-trapped waves and trigger mixing in a deeper layer. The magnitude of k_z and ε even climbed to the averaged level of the Atlantic sector (Fig. 5b & d), indicating that the existence of multiyear sea ice is a possible factor modulating mixing in the deeper ocean.

Further discrepancies between the two main layers in same regions were seen in Fig. 5. In the Atlantic sector of the Arctic, the maximum of k_z and ε were found in the upper layer. While in the Pacific sector, the maxima occur in AW and were more than 1 order of magnitude larger than in the upper layer for several regions. This shows different vertical structures of k_z and ε between the two sectors of the Arctic Ocean. The specifics of the vertical structure are discussed in the next section.

3.4 Vertical structure of k_z and ε in different basins

The averaged vertical structures and corresponding standard deviations of k_z, ε, and stratification N also differed between the two sectors (Fig. 6). In the Pacific sector (mainly the Canadian Basin), it was found that k_z at the depth of HL was very small, increasing with depth and peaked at 300–400 m, which coincided with the bottom of thermocline (Fig. 6a). Below the thermocline, it decreased with depth (Fig. 6a). The vertical structure of ε was similar to k_z, but its decline was more rapid, particularly in deeper depths (Fig. 6b). Its maximum of k_z occurred at the bottom of the HL, which was attributed to the strong HL stratification (Fig. 6c) suppressing turbulent mixing in the upper layer.

The difference between the vertical structures of k_z and ε on the Atlantic sector (mainly Eurasian Basin) compared to the Pacific sector were that the maxima of k_z and ε in the Eurasian Basin were both found in upper layer. The upper layer stratification in the Atlantic sector is about a half of that in the Pacific sector. A weaker stratification enables the turbulent mixing to overcome the stratification and drive ocean mixing. Below that, they decreased with depth, with the fastest decline between 300–500 m. Such features were probably due to the weaker stratification in HL in the Atlantic sector. In the Pacific sector, the N was stronger and enhanced at ~250 m due to the Pacific Water Intrusion and thicker HL. While in the Atlantic sector, N was obviously weaker than that in the Pacific sector. In addition, N decreased rapidly with depth in the upper ocean, enabling the internal waves to easily drive the mixing. This is one possible reason for different mixing mechanisms between the two sectors.

3.5 Long-term variability and seasonal cycle

The thickness and extent of Arctic sea ice declined rapidly from 2005 to 2019 (e.g., Rothrock et al., 2008; Lindsay and Schweiger, 2015). This can affect the ocean circulation, internal tides, mixing, and water properties. We focus on the temporal variability of k_z and ε in two sectors (Atlantic and Pacific) and two layers. The long-term change of k_z and ε are provided in Figs. 7–8.

In the Pacific sector, there is no significant trend of mixing in the entire period (Dosser et al., 2021), whereas the k_z (Fig. 7a–b) decreased from 2005 to 2011, with rates for the upper and deeper layers of 6.4×10^-7 m²/(s·a) and 6.3×10^-6 m²/(s·a), respectively. The downtrend of the deeper layer was stronger and more obvious. The annual decreasing rates of ε (Fig. 8a–b) for the two layers were 4.4×10^-10 W/(kg·a) and 3.8×10^-10 W/(kg·a), with the upper layer downtrend more obvious. After 2011, the yearly mean values of k_z and ε in the upper layer were relatively stable until 2016, and slightly increased during 2017–2018. Meanwhile, the yearly mean values of k_z and ε in AW were not stable, but gradually increased after 2012, with rates of 3.0×10^-6 m²/(s·a) and 3.2× 10^-10 W/(kg·a), respectively. These recent magnitudes nearly reached the pre-2005 levels of k_z and ε. In addition, more significant monthly fluctuations of k_z and ε were found in monthly averages (Figs. 7–8), indicating a more obvious seasonal signal.

The annual averaged k_z and ε in the Atlantic sector increased significantly during all the periods we focused on (Figs. 6c–d, 7c–d), which is consistent with Dosser et al. (2021). In spite of missing data for some months and some years due to relatively scarce observations in the Atlantic sector, the available data was sufficient to support our analysis of annual and seasonal variations. The annual increasing rates of k_z in the upper and deeper layers were 2.1×10^-6 m²/(s·a) and 1.2×10^-5 m²/(s·a), respectively. The rates of ε were 5.8×10^-9 W/(kg·a) and 2.7×10^-10 W/(kg·a), respectively.

To analyze the change of seasonal feature of wave-driven mixing, we partitioned the time series into before and after 2011 (Fig. 9). The monthly averaged k_z in the Pacific sector before 2011 were larger than after 2011. The seasonal feature of k_z in the Pacific sector was insignificant before 2011 in the deeper layer, but slightly increased in summer in the upper layer. After 2011, the k_z in the Pacific sector was larger in winter in both two layers, and also slightly enhanced in summer in AW (Fig. 9a). The seasonal characteristic of k_z in the Atlantic sector (Fig. 9b) was enhanced significantly in summer after 2011. Compared to that before 2011, the enhanced k_z in summer is obvious in the upper layer. Note that shifting the time series partition one year each way will not significantly alter the results (not shown).

4 THE CONTRIBUTIONS OF VARIOUS FACTORS ON MIXING

Mixing is generated by a wide variety of processes, each at its own scale. Here we investigate the roles of two major factors: near-inertial internal waves and stratification.

4.1 Geographical inhomogeneity

The near surface 10-m wind speeds (ERA-5 10-m winds from European Centre for Medium-Range Weather Forecasts reanalysis) were used to determine the 15-year averaged wind pattern in the Arctic (Fig. 10a). Histogram of estimated k_z corresponding to the regions of weak and strong winds is shown in Fig. 10b, which shows that the k_z over the region of strong winds are ~1 order of magnitude greater than it over the region of weak winds. This indicates that the wind drives the mixing significantly in the Arctic. According to the wind speed, the region corresponding to the lowest 40% (< 5.9 m/s) wind speed is defined as weak winds, and it corresponding to the top 40% (> 6.1 m/s) wind speed is defined as strong winds.

The minimum sea ice boundaries for summer and winter are marked in Fig. 10c according to the lowest sea ice concentration (https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html) in each season. It can be noted that the extent of sea ice declined most significantly in the west of the Pacific sector, while the concentration declined most in the Atlantic sector. The statistics of estimated k_z in both regions suggested that k_z in the region always covered by ice were higher 0.5-1 order of magnitude than regions of summer melting ice (Fig. 10d). This clearly indicates that sea ice impacts ocean mixing.

Breaking of Internal tides and near-inertial waves are considered as two dominant factors of mixing in mid-latitude (e.g., Xu et al., 2014, 2016, 2021; Wang et al., 2018; Chang et al., 2019). Due to the weak internal tides in the Arctic deep basins, here we consider only the contributions at near-inertial frequencies. The amplitudes of near-inertial waves were estimated from isopycnal displacements taken from the ITP dataset (Dosser et al., 2014; Method). This resultant estimation from each profile (Fig. 11a) illustrated that stronger near-inertial waves were found in the Atlantic sector, with amplitudes about twice as large as those in the Pacific sector. Stronger near-inertial waves were mainly found in ice-free areas during summer, indicating ice-free surfaces were more productive for the generation of winddriven near-inertial waves. The Arctic Ocean was divided into a weak near-inertial wave field and a strong near-inertial wave field according to amplitudes in the lowest 40% (< 1 m) and top 40% (> 1.5 m, black line, Fig. 11a). These results suggest that k_z in the area with strong near-inertial waves was stronger by 0.5–1 orders of magnitude than with the weak near-inertial waves, indicating an important impact of the nearinertial waves on turbulent mixing in the Arctic Ocean.

However, the different results in the sectors were interesting. The Atlantic sector was typical for the entire region. k_z in the region of strong near-inertial waves was larger. However, in Pacific sector, k_z in the region of strong near-inertial waves was smaller, which contradicts the idea that internal waves drive turbulent mixing. One possible explanation is that other factors drive mixing in this region, such as lee waves and bottom-trapped waves. It is because the region with larger diapycnal mixing corresponds to the slope of Canadian Basin.

Stratification also influences turbulent mixing. The buoyancy frequency N is a measure of the stratification. The maximum N for each profile N_max was generally found in the near-freezing surface layer (Fig. 12). This layer was influenced by sea ice formation and decay. The Pacific sector (mainly Canadian Basin) had relatively stronger stratification compared to the center and Atlantic sector. This stratification difference affected turbulence. The inhibition of stratification on the Arctic Ocean mixing was proved in Fig. 12b. The mean value of k_z under the weaker stratification is larger by 1 order of magnitude than that under the stronger stratification.

To determine where the turbulence dominates ocean mixing, we computed the dimensionless buoyancy Reynolds number Re_b. When the Re_b is too small, the Osborn relation might need to be modified (Bouffard and Boegman, 2013). The maximum of each profile is shown in Fig. 12c. Re_b was used to judge whether the turbulence was sufficiently energetic to overcome the strong stratification and generate turbulent mixing (Ivey et al., 2008). Turbulent mixing is generally considered to be suppressed when Re_b < 20. Profile locations with Re_b > 20 are marked in red in Fig. 12c, and profile locations with Re_b < 20 are marked in blue, indicating the turbulence was suppressed. It can be noted that most of the Atlantic sector was dominated by turbulent mixing with significant red dots, while only a few profiles with Re_b > 20 were found in the Pacific sector, showing weaker turbulent mixing in that region. The depth of Re_b maximum for Re_b > 20 are provided in Fig. 12d. This shows the depth dominated by turbulent mixing and is consistent with maps of stratification. The turbulent mixing generally occurred above 200 m in the Atlantic sector, but deeper in the Pacific sector, suggesting that most turbulent mixing in the Pacific sector was suppressed due to strong stratification.

4.2 Contribution of time-varying k_z

The effects of the aforementioned influence factors on the temporal variability of turbulent mixing during nearly 15 years are discussed here (Fig. 13). The variation of diapycnal diffusivities shows a difference between both two sectors of the Arctic (see Section 2.5). The corresponding sea ice concentration in the Pacific sector and the Atlantic sector (Fig. 13c–d) both decreased markedly. The decline rates of the Pacific sector were 0.82%/a in summer and 0.09%/a in winter. Compared to the Pacific sector, the Atlantic sector decreased more rapidly, especially in winter, with decline rates of 0.60%/a in summer and 0.39%/a in winter. This winter decline rate in the Atlantic sector is much larger than that in the Pacific sector.

The amplitude of near-inertial waves increased in both the two main Basins during the last 15 years (Fig. 13e & f). The averaged amplitude in the Atlantic sector was larger than in the Pacific sector. There are two obvious growth spikes in 2010 and 2014 on the Atlantic sector, which by 2017 increased by 1 m compared with 2006. Moreover, a significant seasonal cycle of near-inertial waves with stronger in summer and weaker in winter was found in the Atlantic sector. The amplitude of near-inertial waves of the Pacific sector remained at ~1 m prior to 2014, with the exception of a growth spurt in 2013. After 2014, it increased significantly, reaching about 2 m in 2019. This is consistent with the increase of k_z in the upper layer, which indicates that the increasing amplitude of near-inertial waves will enhance turbulent mixing in the upper layer, but it still controlled by stratification. While during this period, the trend of k_z in the Atlantic sector increased with the significant enhancement of near-inertial wave amplitudes and constantly weakening stratification. The increases of k_z in 200–600 m is also consistent with the changing nearinertial waves, despite the fact that the near-inertial waves mainly trigger mixing in the upper layer. One possible reason is that the weakened stratification might provide conditions that the near-inertial internal waves are able to reach deeper layers.

It can be reasonably inferred that the control factors of two basins in both sides of the Arctic Ocean differed. The k_z of the Atlantic sector (mainly Eurasian Basin) is controlled by internal-waves, because of the weak stratification (Figs. 12a & 13h). The significant increase of k_z was consistent with the increasing internal waves (Fig. 13f), even though we only considered near-inertial waves. The k_z of the Pacific sector (mainly Canadian Basin) is controlled by both strong stratification and increasing near inertial internal waves. Strong stratification in this area suppressed wave-driven mixing in previous period, but such limitation of stratification on mixing began to change due to the weakening stratification (after 2011) and growing near-inertial waves. In later period, the k_z is gradually increasing after 2012 as shown in Fig. 13a.

Moreover, it can be noted that the decline of sea ice extent and sea ice free range in summer (Fig. 10c) were significant in the Pacific sector. But the decline of sea ice concentration is more significant in the Atlantic sector (Fig. 13c–d). Such faster decline of sea ice provides conditions for the increase of nearinertial waves, and the "Atlantification" weakens the stratification (Polyakov et al., 2017). This may be a reasonable explanation for why mixing in the Atlantic sector tends to be dominated by turbulent mixing.

The scatter plots and the linear fitting results show the distinct modulation effects on mixing between the two sectors (Fig. 14). The linear fitting slopes of the diffusivities and near-inertial waves amplitude are 0.036 in the Atlantic sector and 0.096 in the Pacific sector. It can be seen that the response of mixing to the near inertial waves is more significant in the Atlantic sector than in the Pacific sector. The linear fitting slopes of diffusivities and stratification are -0.379 (the Atlantic sector) and -0.392 (the Pacific sector). The mixing responses to the stratification in the two sectors are about the same magnitude, despite the much weaker stratification in the Atlantic sector.

From our results, it can be found that the mixing in the Atlantic sector is increasing in the last decade, which is corresponding to the increasing near-inertial waves and decreasing stratification. And the most part of the Atlantic sector is Re_b > 20, indicating a wave driven mixing regime (Dosser et al., 2021). Whereas in the Pacific sector, k_z first decreased and then increase, with a large part of the Pacific sector is Re_b < 20, indicating the strong stratification suppressed the turbulent mixing. However, in recent years, the increased k_z and near internal waves might weaken the inhibition effect of stratification.

5 TURBULENT HEAT FLUX

Utilizing the estimated k_z, we further calculate the vertical turbulent heat flux, which will influence the upper ocean heat content and impact sea ice and climate.

where ρ, C_p, k_z, and dT/dz are density, specific heat, diffusivity, and the potential temperature gradient, respectively. The heat exchange between the warmer AW and the colder, fresher, upper ocean is of major concern. The averaged thermocline depth of 1° latitude × 3° longitude grid cells were examined (Fig. 15b).

The vertical heat fluxes at the depth corresponding to the thermocline (Fig. 15b) are given in Fig. 15a. The boundary between AW and the surface water occur at the thermocline, so it was the depth of major interest for this process. The maxima of vertical heat flux are generally found in this depth. The largest vertical heat fluxes were found in the Atlantic sector. The AW heat was constantly released upward after it entered the central Arctic Ocean through the Fram Strait. The AW heat core constantly deepened from the Atlantic sector to the Pacific sector, which decreased the available heat for further upward transfer. This result indicated that in the Atlantic sector the heat transfer from AW to the surface occurred more easily and resulted in more sea ice melting. However, other factors might be more important for the heat to melt sea ice in the Pacific sector, such as a stronger stratification, limiting the heat transport to the surface, and double-diffusion, contributing a considerable part of heat flux in the Canada Basin (e.g., Smith et al., 2018; Timmermans et al., 2018), and the heat reservoir in the Pacific Summer Water and Near Surface Temperature Maximum layer. Some other potential factors are the intrusion of the Pacific Water or the increasing of sea ice albedo. In the Pacific sector, it appears to be more difficult for the AW to transfer heat to the surface and directly affect the upper ocean.

Temporal and spatial averaged vertical structures of F_h and potential temperature gradients in the two main basins are provided in Fig. 15c–f (here results are shown for 2006, 2008, 2010, 2012, 2014, 2016, and 2018). In the Pacific sector, both the depth of maximum potential temperature gradient and the depth of maximum of corresponding heat flux deepened with time. A sub-peak in the upward heat flux was found at ~150 m after 2010. We cannot determine whether such peaks will overcome the strong stratification in the upper ocean and transfer heat to sea surface. In addition, the heat flux in the Pacific sector decreased in post 2010, which might due to the decreasing of temperature gradient. In the Atlantic sector, our results showed that the maxima potential temperature gradients were found in the upper ocean (the minimum depth of estimated k_z restrains our focus depth). The AW heat probably affected the sea surface and sea ice, because the maxima always occurred in the upper ocean. The decreasing sea ice in the Atlantic sector will become more significant, and the strong stratification at the upper layer will gradually weaken due to the increasing diapycnal heat (Polyakov et al., 2020a). Such amplification of upward heat flux will enhance the rate of warming in the Arctic region and contribute to polar amplification (Beer et al., 2020).

6 CONCLUSION

In this study, we estimate diffusivity k_z, and dissipation rate ε, in the Pan-Arctic Ocean for water depths between ~100–700 m from high resolution observations of ITP using a strain-based finescale parameterization method. A 3-dimensional wavedriven mixing structure and temporal variability of diffusivity and dissipation rate are given in this paper (Figs. 5–9). It can be noted that the mixing properties for both sides of the Arctic Ocean are apparently different. The diapycnal diffusivity k_z and dissipation rate ε in Atlantic sector were 3.5×10^-5 m²/s and 1.2× 10^-9 W/kg, respectively; while those in the Pacific sector were lower, 9.2×10^-6 m²/s and 7.4×10^-10 W/kg, respectively. The interannual variability of k_z and ε is different in the two sectors. A constantly increasing trend of k_z was found in the Atlantic sector, whereas k_z in the Pacific sector weakened first and then increased (see Section 3.5).

Potential influence factors for mixing in the Arctic Ocean include wind, sea ice, near-inertial waves and stratification. Their roles were discussed (Figs. 10–12, Section 4.1). The estimated amplitude of near-inertial waves and stratification in the Arctic Ocean show two modulations of turbulent mixing in both sides of the Arctic Ocean. The Atlantic sector shows a mixing mechanism dominated by internal waves. However, the mixing in the Pacific sector is limited by stratification and not controlled by internal waves. This difference resulted from the different stratification in the sectors. This viewpoint was supported by the pattern of the buoyancy Reynolds number, Re_b (Fig. 12a).

Based on the above-mentioned viewpoint, we analyzed the temporal variability of k_z, F_h, and relative influence factors (Figs. 13 & 14, Section 5). The mixing in the Atlantic sector was always controlled by internal waves, and k_z and F_h constantly increased with the increasing near-inertial waves. The mixing in the Pacific sector may be in a period of change from being limited by stratification to being dominated by internal waves. However, we cannot predict whether turbulent mixing will dominate in this sector in the near future. The temporal variability of k_z increased after 2011, which may because near-inertial waves gradually overcame strong stratification and drove turbulence after 2011. The heat flux in Pacific sector is decreasing. In addition, the depths of heat flux maxima gradually descending. Thus, it is more difficult for the AW heat to be rapidly released to the upper ocean. Moreover, a sub peak of heat flux was found in the HL during the last 10 years. How such peak will affect the upper ocean is still uncertain.

Due to the limitations of parameterization method, these results and conclusions mainly focused on wave-driven mixing in the Arctic Ocean. The diffusivities and heat flux caused by double diffusion were out of the scope of present work. Our result shows that double-diffusive staircases are generally found in the central Canada Basin (Fig. 1), which probably contributes to mixing here, although recent research have shown that the contribution of doublediffusive mixing is small (Boog et al., 2021) on the global ocean. In addition, other types of internal waves can also cause mixing in some rough topography and ice environment (e.g., Schulz et al., 2021b; Zhang et al., 2021). Further question remains as to the mixing properties of the Pan Arctic. Future work will take double diffusion and other internal waves into account and focus on the transition of mixing properties.

7 DATA AVAILABILITY STATEMENT

The Ice-Tethered Profiler data sets were made available by the Ice-Tethered Profiler Program (Krishfield et al., 2008; Toole et al., 2011) at the Woods Hole Oceanographic Institution (http://www.whoi.edu/itp). ERA-5 10-m winds speed data are from European Centre for Medium-Range Weather Forecasts reanalysis (https://www.ecmwf.int/en/forecasts/datasets). The sea ice concentration are available at https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html.