Several methods have been developed to estimate the MLD, including the difference threshold (Kara et al., 2000; Kara et al., 2003; de Boyer Montégut et al., 2004), gradient threshold (Lukas and Lindstrom, 1991), curvature method (Lorbacher et al., 2006), split and merge method (Thomson and Fine, 2003), and hybrid method (Holte and Talley, 2009). After quantitative comparison, in this study the widely adopted threshold method (de Boyer Montégut et al., 2004; Peralta-Ferriz and Woodgate, 2015) was selected. The criterion was critical, as the small vertical fluctuations of the temperature and density within the ML brought confusion when calculating the MLD. A group of criteria was tested, and the results were shown in Supplementary Material VI . The adopted optimal criterion for the MLDd was Δσ = 0.125 kg/m³, and the reference depth was 5 m. The suitable criterion for the temperature difference method was 0.5°C. Visual examination confirmed that the criteria successfully identify the variations of MLD. To distinguish the role of temperature and salinity in the spatial variation of MLD, the MLDd (MLD estimated from the density profiles) and the MLDt (MLD from the temperature profiles) were both calculated, and in the following, MLD specifically referred to MLDd.

For the stratification index (SI) due to temperature (salinity), the density profile p is calculated using the temperature (salinity) profile and a vertically averaged salinity (temperature), as Equation (1) shows.

where V_θ represents the stratification index due to temperature, ρ ( θ , S ¯ ) is the density profile derived from the temperature profile and a vertically averaged salinity, and S ¯ is the vertically averaged salinity. h is the value of MLD here.

The MLD was closely related to the distribution and interaction of water masses. There were several water masses in the Bering Sea and Chukchi Sea, including three in the northern BSs: Anadyr Water (AW), Alaskan Coastal Water (ACW), and Bering Shelf Water. As shown in Figure 2 , the average temperature of the cold and salt AW (S > 32.5, T < 4°C) was approximately 2°C and the average salinity was about 32.6. The average temperature of ACW was about 9°C with the salinity significantly lower than that of AW ( Figure 2 ). The water mass below the surface was termed the cold intermediate layer (CIL) in BSb. The core temperature (the lowest temperature in the vertical section of potential temperature) of CIL was approximately 3°C ( Figure 2 ).

In summer, the water mass in the Chukchi Sea shelf was a product of the mixing of Bering Sea Water, ACW, and Meltwater (Corlett and Pickart, 2017; Pisareva, 2018). Fundamentally, the farther north, the colder and fresher the water mass was (blue circle in Figure 2 ). The cold water mass in the intermediate layer was Winter Water and the deep water was Atlantic Water on the north of the Chukchi Sea slope ( Figure 2 ). Vertical sections of potential temperature and salinity can be found in Supplementary Material VII ( Figure S4 , Figures S6 – S8 ).

In general, the mean MLD in BSb was 17 m (with standard deviation being 1.5 m), greater than 13 m (with standard deviation being 4.2 m) on the BSs. The mean MLD in BSp was 28 m (with standard deviation of 3.3 m), significantly greater than BSb and Bss. The BL section was representative of the main circulation and water masses during the expedition period in both the BSb and BSs. Stations BL01–BL06 were located in BSb, and the MLDs at these stations were all greater than 15 m ( Figure 3A ). The maximum value of MLD at these stations reached 30 m, which was observed at station BL01 near Aleutian Islands. This deep ML was caused by an anticyclonic eddy and will be discussed in Section 3.5.1. In contrast, at stations BL11–BL14 in the west of BSs, the MLD was shallower than 15 m. In addition, the minimum MLD at BL section was 6 m, which was observed at station BL14 on the northwestern BSs. The MLD at station BL07 located in BSp was 25 m. Both of the MLD at stations BL01 and BL07 were markedly larger than other stations in BSb and BSs.

The stations BR00–BR03 were located in BSp. The maximum MLD reached 30 m at station BR03 ( Figure 3B ). The average MLD at stations BR04–BR09 was 14 m in BSs, which was much smaller than the 23 m in BSp. Stations BR10 and BR11 were in the northeastern BSs, and isohalines and isotherms were almost vertical there. Correspondingly, the MLD at stations BR10 and BR11 were greater than those at other BR stations on BSs.

The western BS section was under the influence of AW in the northwestern Bering Sea shelf, and the eastern BS section was under the influence of ACW in the northeastern Bering Sea shelf. Therefore, the BS section represented the MLD under the influence of the advection of these two water masses. As the water column at BS01–BS03 was well mixed, the MLDs were larger than 35 m ( Figure 3C ). The MLD in the northeast was much smaller and was in the range of 6–13 m. The MLD at stations BS04 and BS05 was small, because the interaction of the light ACW and dense AW stratified the upper ocean there.

In general, the MLD in CSs had an overall trend of deepening from south to north ( Figure 3D, F ). The MLD increased from 5 m in the southern CSs to 10 m in the northern CSs, which covered a distance of about 880 km ( Figure 3D ). The MLD at Section R was in the range of 4–12 m. The MLD at R06 was shallower than the other stations along Section R ( Figure 3D ). A northward increase was also observed in the BT section. The maximum MLD was 20 m at station BT15. The MLD at stations BT13–BT16 was greater than 15 m and was also greater than the MLD in CSs ( Figure 3F ).

Though Sections M and BT were at the same latitude and close to each other, the MLD at Section M was in the range of 5–10 m, much shallower than that at Section BT ( Figure 3E ). The water mass along Section M was colder and fresher than that along Section BT. It seemed that Section BT was more easily influenced by the warm and salt water through the Bering Sea shelf, while Section M was more influenced by the sea ice melt water.

The MLD and MLDt had little difference in the southern Bering Sea, while the MLDt was generally larger than the MLD in the northern BSs and the Chukchi Sea ( Figure 3 ). The MLD and MLDt had little difference at stations BL02–BL05, but the difference approached 4 m at station BL06 ( Figure 3A ). It should be noted that the MLDt (14 m) was much larger than the MLD (6 m) at station BL14. The average difference d_MLD (MLD minus MLDt) was 0.51 m in the southern Bering Sea (including BL01–BL06, BR00), while the d_MLD was −3.25 m in the northern BSs and the Chukchi Sea. The largest d_MLD was less than 4 m in the southern Bering Sea, while the largest difference between MLDt and MLD was greater than 11 m in the Chukchi Sea.

The spatial variation of MLD was closely related to the water masses. The larger vertical density gradient and shallower MLD at stations BS04 and BS05 were related to the westward advection of the warm and fresh ACW. The northward increase in the MLD in the Chukchi Sea was accompanied by a high meridional gradient of salinity and temperature. This was the result of the advection of low-salinity water generated from the melting of sea ice in summer in the Chukchi Sea. Around station R06, the warm and fresh ACW, with shallower MLD, was diverted onto the western CSs due to Ekman transport under the influence of anomalously strong northerly winds (Pisareva et al., 2015). The same situation has been observed in September 2004 (Pisareva, 2018), 2008 (Linders et al., 2017), 2009 (Pisareva et al., 2015), and 2012 (Pisareva, 2018). The shallower MLD at station M was also related to the deviation of ACW.

Although the MLD increased in CSp as that in BSp, there was a difference between them. Remarkably, from the ocean basin towards the continental shelf, the isotherm and isohaline tended to be parallel to the continental slope in the Chukchi Sea, while they tended to be perpendicular to the continental slope in the Bering Sea. The MLD variation in CSp was related to the low-salinity water generated from the melting of sea ice in summer and topographical constraints. However, in BSp the isotherm, isohaline, and MLD were mainly affected by BSC.

The stratification index (SI) was calculated over the whole depth of the ML. The value of SI was O (1,000 J/m²) in BSb and the southern BSs, including BL01–BL13, as shown in Figures 4A, B . In BSb, eastern BSs, and western BSs, the SI was dominant by salinity. This was different from the research of Ladd and Stabeno (2012), which concluded that temperature and salinity both influence the stratification of the Eastern Bering Sea shelf with their relative importance varying spatially and temporally. This was because their measurements were influenced by ACW with high temperature and low salinity. In the northeastern BSs, due to the high salinity of AW, the SI was significantly larger ( Figure 4C ). In the northwestern BSs and the Chukchi Sea, the SI was smaller ( Figures 4C–F ). The SI showed a trend of northward decreasing and was dominated by salinity. The contribution of the temperature to stratification was too weak and could be ignored for our observations. Therefore, it was reasonable to assume that the characteristics of the ML were dominated by the low-salinity water generated from the melting of sea ice in the Chukchi Sea and the northern BSs in the summer of 2019. In addition, the farther north, the greater the difference between the MLDt and MLD was in the Chukchi Sea, and salinity played a more important role in determining the MLD. This was consistent with the research of Johnson et al. (2012), which showed that the seasonal variation of the ML in the Arctic was dominated by salinity.

Peralta-Ferriz and Woodgate (2015) found that the overall ML gradually became shallower from 1978 to 2012 by studying the seasonal and inter-annual changes of MLD in the Arctic region, including the northern Chukchi Sea in winter. However, MLD in the Bering Sea and summer Chukchi sea was not presented in their study. According to our observations, the MLD in 2019 was shallower, and the temperature of the ML was higher than the mean temperature of the other five (1999, 2003, 2010, 2012, and 2014) years along the BL section in the Bering Sea ( Figures 5A, B ). The mean MLD was 17 m (standard deviation being 5 m) in 2019, while the mean MLD was 23 m (with standard deviation being 5 m) for previous years mentioned above. The shoaling of the ML was more evident in the Bering Sea basin than that in the Bering Sea shelf. It should be noted that, although warmer ML was found along the whole BL section ( Figure 5B ), the shallower MLD does not occur in BSp and the continental slope north of Aleutian Island (BL01) ( Figure 5A ). This shallower MLD was accompanied by the warming of 7the surface layer and CIL. The minimum temperature of the CIL water mass in BSb showed an increasing trend: 0.54°C, 0.94°C, 0.82°C, 0.69°C, 1.99°C, and 2.50°C for 1999, 2003, 2010, 2012, 2014, and 2019, respectively. The ML was saltier in 2019, which led to denser ML in the Bering Sea shelf ( Figures 5C, D ). On the contrary, the ML was fresher in the Bering Sea basin in 2019, which led to lighter ML ( Figures 5C, D ).

The MLD in 2019 was shallower than the mean MLD of 1999, 2003, 2010, 2012, 2014, and 2017 along the R section in the Chukchi Sea ( Figure 5E ). The mean MLD was 8 m (standard deviation being 2 m) in 2019, while the mean MLD was 19 m (with standard deviation being 5 m) for previous years mentioned above. This was accompanied by the warming of the ML ( Figure 5F ). The ML showed a trend of cold and freshening northward along the R section in all these years, including 2019. The warming and freshening around 69°–70°N may be caused by the westward advection of ACW, and the ML showed a deepening there. The warming, freshening, and shoaling of the ML around“69°–70°N”indicated by our observations were not identified by WOA2018. Salinity dominated the spatial fluctuation of the density for most of the year ( Figures 5G, H ).

As mentioned in Section 3.2, the large MLD at BL01 in the northern continental slope of the Aleutian Islands was related to the anticyclonic eddy along the Aleutian Islands, as shown in Figure 6A . The maximum speed radius of the eddy next to BL01 was derived from SSH datasets according to the method developed by Dong et al. (2009) and Nencioli et al. (2010). The amplitude, maximum speed radius, and effective radius of the eddy were 0.30 m, 49 km, and 83 km, respectively. The distance between the eddy center and BL01 was 40 km, which indicated that the station was influenced by the eddy, especially by the dynamics adjacent to the eddy periphery ( Figure 6B ). Pumping exerts its impact on the upper layer inside the eddy, while multiple processes, including advection and submesoscales, are involved outside the eddy periphery (Wang et al., 2018). Hence, the large MLD at BL01 was caused by this anticyclonic eddy through the downwelling effect. This coincides with the conclusion that anticyclones deepen the MLD in the research of Gaube et al. (2019). The MLD at BL01 was 30 m, significantly larger than that at BL02, which was 19 m ( Figure 3A ). The upper-ocean current velocity at BL01 was approximately 0.2 m/s, while it was less than 0.1 m/s in BSb according to the ADCP observations. There was enhanced shear at the base of the ML at BL01 ( Figure 6C ). Under the interaction of ocean currents and the Aleutian Islands, the eddies near the Aleutian Islands were relatively active, which result in deep MLD.

The relationship between MLD and buoyancy flux as well as momentum flux was explored through correlation and linear regression analysis. The momentum flux and average buoyancy flux caused by sea surface net heat flux and freshwater flux from July 1 to September 8 are shown in Figure 7 . The MLD had a positive correlation with the momentum flux, and the correlation coefficient was 0.60 at 99% significant test level. The kinetic energy input into the ocean due to the large wind speed enhanced the sea surface turbulent mixing. Although this pattern can be found at most stations, it was not always valid, e.g., the larger MLD at BS01, BS02, and BS03 was related to the vertical homogeneous Anadyr Water there (as indicated in Figure S4C in Supplementary Material VII ). Anadyr Water was cold and vertically homogeneous. As mentioned in Section 4.1, the deepening of the ML at station BL01 resulted from the anticyclonic eddy. The correlation coefficient between MLD and buoyancy flux was −0.30 at 98% significant test level.

The MLDs at BR02 and BR03 (in the central Bering Sea) were significantly large under the combined effect of the buoyancy flux and the momentum flux. Though not as significant as BR02 and BR03, the MLDs at BT13, BT14, BT15, and BT26 were still larger than other stations nearby in the Chukchi Sea. As momentum flux played a more important role than buoyancy flux, the small momentum flux led to small MLD, such as R02 and R03 (southern Chukchi Sea). This convinced the modeling work by Hordoir et al. (2022) that wind led to the spatial variation of MLD in ice-free conditions in the Arctic.

However, multiple linear regression showed that the combined effect of momentum flux and buoyancy flux only explained 12% of the variance of MLD, and 11% of the variance of MLD was attributed to momentum flux (p = 0.01). So, there were more important factors to be explored, such as stratification to be discussed in the next section.

Controlling effects of stratification on MLD were estimated through the approach proposed by Peralta-Ferriz and Woodgate (2015) in the study on pan-Arctic surface mixed layer properties. In their study, to investigate reasonable relationships between MLD and stratification, an empirical expression was proposed as follows:

with the exponents α β, and constant C (in S.I. units) to be estimated from multiple regression of observational data. W was wind magnitude (m/s). Δρ (kg/m³) was the stratification for each cast and it was quantified as the difference between the average density of a 20-m-thick layer below the ML and the average ML density.

Based on observational data in 2019, the detailed fit results are shown in Table 2 . Besides the combined effect of wind and stratification, the wind-only and stratification-only exercises were also conducted by assuming β or α to be zero.

In both combined and wind-only exercises, α was ~0.5. In both combined and stratification-only exercise, β was ~0.5, larger than the results from Pollard et al. (1973) (β ~0.25) and Manucharyan (2010) (β ~ 0.33). This suggested that MLD in this region was more dependent on stratification. The estimated value of α and β here were similar to those estimated by Peralta-Ferriz and Woodgate (2015) in the pan-Arctic. This indicated the relationship between MLD and stratification in the Bering Sea and Chukchi Sea had similar patterns with that in the pan-Arctic. The combined effect of wind and stratification explained 57% of the variance of MLD, much larger than the combined effect of momentum flux and buoyancy flux in Section 4.2. This was smaller than the mean value (67%) from the pan-Arctic sample in the study of Peralta-Ferriz and Woodgate (2015). This indicated that there were other factors needed to be considered in this marginal sea. The stratification explained 40% of the variance of MLD, much larger than wind (14%), and this suggested that in order to understand the variation of MLD in this region, it was necessary to focus on investigating processes affecting upper ocean stratification, i.e., sea ice melt water in the next section.

The less sea ice volume resulted in less melt water in summer ( Figure 8 ), which led to the saltier and shallower ML in 2019 than previous years. As the mean wind speed did not change significantly in 2019 compared to previous years ( Figure S9 in Supplementary Material VIII ), it was not related to the shallower of the ML. It should be noted that the salinity of the water in BSs was larger than the climatology ( Figure 5C ), while it was not in CSs ( Figure 5G ). This may be linked to the inflow of freshwater due to the increasing net glacial ablation in the Gulf of Alaska watershed (Danielson et al., 2020). The surface warming in the Chukchi Sea was related to the regional air-sea heat flux and Arctic amplification (Danielson et al., 2020). Recent trends to earlier ice melt and increasing heat flux through the Bering Strait have also impacted the Chukchi Sea (Woodgate, 2018). The warming of CIL may be related to the air temperature warming in the previous winter and the processes in the Bering Sea (Overland et al., 2012).