{'fraction_non_alphanumeric': 0.030726658520807876, 'fraction_numerical': 0.018957090490907783, 'mean_word_length': 4.920531914893617, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The warmer the ocean surface, the shallower the mixed layer. How much of this is true?

Abstract Ocean surface warming is commonly associated with a more stratified, less productive, and less oxygenated ocean. Such an assertion is mainly based on consistent projections of increased near‐surface stratification and shallower mixed layers under global warming scenarios. However, while the observed sea surface temperature (SST) is rising at midlatitudes, the concurrent ocean record shows that stratification is not unequivocally increasing nor is MLD shoaling. We find that while SST increases at three study areas at midlatitudes, stratification both increases and decreases, and MLD deepens with enhanced deepening of winter MLDs at rates over 10 m  decade−1. These results rely on the estimation of several MLD and stratification indexes of different complexity on hydrographic profiles from long‐term hydrographic time‐series, ocean reanalysis, and Argo floats. Combining this information with estimated MLDs from buoyancy fluxes and the enhanced deepening/attenuation of the winter MLD trends due to changes in the Ekman pumping, MLD variability involves a subtle interplay between circulation and atmospheric forcing at midlatitudes. Besides, it is highlighted that the density difference between the surface and 200 m, the most widely used stratification index, should not be expected to reliably inform about changes in the vertical extent of mixing.


Introduction
As greenhouse gases increase the Earth radiative imbalance, the oceans sequester up to 93% of the extra heat in the climate system [Alexander et al., 2013], and as a result, the ocean is warming [Balmaseda et al., 2013a[Balmaseda et al., ,2013bRoemmich et al., 2015]. However, some oceanic regions absorb more heat than others, and ocean warming does not occur homogeneously in space ( Figure 1) and gradually from the surface to the bottom as one might expect from a regular removal of atmospheric heat.
Certainly, the marked seasonality in the heat exchange with the atmosphere (summer heat uptake versus winter heat release) drives the broadly noticed coupled changes at seasonal time-scales between the increase/ decrease in upper water temperature and stratification, as well as their control in the extent of vertical mixing, nutrient supply to the upper layers, and primary productivity [Lozier et al., 2011]. However, as also noticed by Lozier et al. [2011], the long-term variability of these processes is not necessarily dominated by the same driving forces. Indeed, the ocean stratification and mixing control from the overlying atmosphere is representative of processes governing the evolution of parameters within and above the seasonal pycnocline, a feature that is created and destroyed annually by definition. It is the permanent pycnocline which acts as a barrier layer against turbulent mixing  and prevents the ocean from being progressively mixed from the top to the bottom. The properties of the permanent pycnocline are controlled by air-sea exchanges at mid and high-latitudes at mode water formation places and their annual subduction rate since it is this flux that ventilates the permanent pycnocline. Thus, temperature and salinity anomalies generated by winter airsea exchanges at these locations are carried at depth by ocean currents to remote locations where from below they exert control over the permanent pycnocline . As a consequence, unlike the atmospherically driven control on stratification of the seasonal pycnocline from above, the permanent stratification of the water column is sustained from below by ocean circulation.
Bearing that in mind, the assumption that ocean surface warming is unequivocally associated with a more stratified and less mixed ocean seems inaccurate. The existence of such a definitive relation is commonly supported by model outputs projecting increased near-surface stratification and shallower mixed layers under global warming scenarios [Behrenfeld et al., 2006;Polovina et al., 2008;Boyce et al., 2010;Tyrrell, 2011;Gruber, 2011;Capotondi et al., 2012;Schmidtko et al., 2017]. However, CMIP5 models have shown difficulties in reproducing the MLD seasonal cycle and the location of the maximum MLDs in the Southern Ocean [Sall ee et al., 2013]. The seasonal cycle of the mixed layer at midlatitude areas of thick mixed layers is intimately related to, and critically controlling, the ventilation of the permanent pycnocline [Marshall and Nurser, 1993]. This causes any inaccuracy in determining the location of maximum MLDs and the MLD seasonal cycle to misrepresent the thermohaline properties of the permanent pycnocline and its ventilation rates [Sall ee et al., 2013]. These biases are not specific only to the Southern Ocean, and CMIP5 simulations have generally been found affected by large drifts in circulation and stratification from mid to high-latitudes [Russell et al., 2015], which call into question the presumed unequivocal relation between warming, increasing stratification, and shallower MLDs.
It seems evident that increased stratification from surface warming and freshening would hinder the ventilation of the ocean. However, the ventilation of the ocean depends on a second competing force, winddriven divergence that can compensate and even reverse the effects of ocean warming from above causing increased stratification [Russell et al., 2006;Lauderdale et al., 2013]. Apart from the certain warming of the upper ocean over the last 30 years [Alexander et al., 2013], the poleward intensification of the westerly winds in the Southern Ocean is one of the most significant trends in the global climate system [Russell et al., 2006;Sall ee et al., 2010]. In the Northern Hemisphere, a poleward migration of the jet stream (maximum westerlies) has also been found but in this case associated with a weaker and wavier (higher wave amplitude) jet stream [Archer and Caldeira, 2008;Capua and Coumou, 2016]. Both the poleward migration and the wave amplitude and asymmetry changes of the jet stream may result in wind-driven divergence changes [Lauderdale et al., 2013;Landsch€ utzer et al., 2015]. On the other hand, the weakening of the jet stream has been related to a higher frequency of atmospheric blocking anomalies in the Northern Hemisphere [Barnes, 2013;Francis and Vavrus, 2012] which can have severe effects on the upper ocean vertical structure and ocean circulation [Somavilla et al., 2016]. All these changes have an effect on MLD and stratification variability that may oppose the expected one-directional changes due to a regular warming from the overlying atmosphere.
Here we derive from hydrographic data the MLD and stratification trends observed in the upper ocean surface layer. In contrast to model simulations, we find that while SST increases in every area studied at midlatitudes, stratification either increases or decreases, and MLD deepens with enhanced deepening of winter MLDs. We present these results in section 3. The reasons behind these results and their consequences for biogeochemistry are discussed in section 4, focusing on: (1) caveats regarding the assumption that a warmer ocean will be more stratified and the suitability of the density difference between the surface and 200 m, the most widely used stratification index, to represent the tendency of a water column to be mixed (section 4.1); (2) the roles of changes in the buoyancy forcing (densification of the ocean surface due to Figure 1. Global SST Trends. 1981-2015 SST trends from the optimum interpolated SST data set (OISST, https://www.ncdc.noaa.gov/oisst). Location of the oceanographic timeseries HOTS, BATS, and SATS in the study areas of the North Pacific subtropical gyre (sg-NPac) and subtropical gyre North Atlantic (sg-NAtl) as defined in McClain et al. [2004] based on dynamic height, and in the midlatitudes of the Eastern North Atlantic (ml-ENA) as in Somavilla et al. [2016], where MLD long-term trends will be shown later in the paper.
Journal of Geophysical Research: Oceans 10.1002/2017JC013125 cooling and/or evaporation) and wind-driven divergence/convergence (Ekman pumping) on MLD deepening trends (section 4.2); and (3) the potential effects of our findings on primary productivity and ocean oxygenation (section 4.3).

Study Areas and Data Sources
The study focuses on the long-term changes and trends in Mixed Layer Depth (MLD) and stratification in the subtropical gyres and midlatitude regions of the North Atlantic and North Pacific in the areas shown in Figure  1. These areas are of particular interest because ocean warming is observed there and has been related to reduced primary productivity and ocean deoxygenation through increased stratification and shallower MLDs [Behrenfeld et al., 2006;Polovina et al., 2008;Boyce et al., 2010;Talley et al., 2016] and/or they comprise or partially overlap with mode water formation regions. Besides, in each of these areas, long-term oceanographic time series are also regularly maintained providing high-quality hydrographic data, which make feasible assessments of actual long-term changes in ocean stratification and MLD. Oceanographic time-series used are the time series at ocean observatories HOTS (http://hahana.soest.hawaii.edu/hot/, 22:8 N, 158 W) in the North Pacific subtropical gyre (sg-NPac); BATS (http://bats.bios.edu/32:2 N, 64:5 W) in the North Atlantic subtropical gyre (sg-NAtl); and SATS (http://www. boya-agl.st.ieo.es/, 43:8 N, 3:8 W) at midlatitudes of the eastern North Atlantic (ml-ENA). A caveat of using standard oceanographic time-series is that they can be affected by temporal gaps and provide local conditions for only a single spot, and so may not be consistently representative of a wider region. To overcome this issue, we combine them with other data sets that provide greater temporal and spatial coverage (Argo floats, global gridded ocean reanalysis, and satellite data).
To assess changes in MLD and stratification, we use the density profiles from the oceanographic time-series, Argo floats, and ocean reanalysis (ECMWF Ocean Reanalysis System 4, ORAS4) [Balmaseda et al., 2013a]. To evaluate the warming of the ocean surface, data of sea surface temperature (SST) based mainly on satellite data (OISST, https://www.ncdc.noaa.gov/oisst) are used. To estimate the effects of buoyancy fluxes and Ekman pumping on MLD winter deepening, we use surface fluxes from NCEP/NCAR Reanalysis data set [Kalnay et al., 1996] (For further information on these data sets, see section 1.1 in supporting information). Our results would not change by using ECMWF Reanalysis data, since the comparison between NCEP/NCAR and ECMWF reanalysis data show very good agreement (the coefficient of variation of the RMSD is less than 3%).

Stratification and MLD Indexes
Stratification and MLD are among the most useful properties in oceanography because of their relevance both in climate processes and biogeochemical cycles. However, there are no standard criteria for their definitions and hence to establish objective metrics [Gonz alez-Pola et al., 2007;Thomson and Emery, 2014]. Figure 2 provides a view of different methods within a typical midlatitude winter profile and the problem of comparing methods for MLD determination when not perfectly homogeneous mixed layers exist (www. mldidentificationexperiment.com)-in this case due to the existence of both a deep MLD and an evanescent shallow MLD. Figure 2 also shows two extreme hypothetical profiles with the same bulk overall stratification.
The temperature (StrI h ) or density (StrI r ) difference between the sea surface and 200 m is the stratification index that is most frequently employed [Behrenfeld et al., 2006;Capotondi et al., 2012]. The MLD is often estimated using the so-called threshold methods, which define the MLD as the depth at which temperature or density differs from the surface value by a given amount (Dh or Dr h ) ( Figure 2). Alternatively, the use of other more advanced metrics for MLD determination is not rare. Here, we include the measures of MLD obtained following the Gonzalez-Pola et al. [2007] algorithm (MLD GP ) and the Holte and Talley [2009] algorithm (MLD HT ). The Gonzalez-Pola et al. [2007] algorithm is based on the best fit of density profiles to a predefined ''ideal'' functional form (i.e., converting the actual profile to an analytical curve). This fitting depends on various parameters with physical meaning among which is the MLD (MLD GP in Figure 2) and from which also various measures of stratification can be obtained. The most representative are the slope of the permanent pycnocline, denoted by b 3 , and the maximum stratification of each profile, G max , i.e., the maximum density increase per meter in the water column in the profile (Figure 2b). Since both temperature and salinity data are available in the data sets used in this study, we focus on the density-based estimates of these methods.

Trends and Seasonality in the Upper Ocean
Long-term changes and trends in signals like SST, stratification, and MLD having clear, well-defined seasonal signals can be analyzed through well-established statistical tools. If we expect a smooth variation, a good approach is the assumption that the data can be modeled by a periodic signal and optionally a linear trend. This will be the case of SST and StrI r . Through a harmonic fit, a periodic signal can be approximated by: where time is in decimal years. This mixed model corresponds to a linear trend plus annual and semiannual harmonics. The parameters are then fitted by generalized least squares method and the inclusion/exclusion of more parameters in the fitting is determined by a t test [Jenkins and Watts, 1998].
For the MLD a fitting using harmonics is unsatisfactory, as the mean evolution shows a discontinuity from the maximum MLD reached at the end of the winter to the surface at the beginning of the spring ( Figure 3b). Moreover, between the end of the winter and beginning of spring, the MLD can alternate from very deep to very shallow values when incipient stratification is established at the surface ( Figure 3b). Actually, at times a deep but vanishing MLD and a shallow incipient MLD coexist. Since this transition period marks the end of the winter and beginning of the spring, we establish this day as the date of minimum climatological SST when on average the maximum MLD is reached (65.33 days). Centered on the SST Date min , we define after visual inspection the period of time when the MLD is not unequivocally defined with a range of 0:1 year (0:05 year to each side) ( Figure 3b). During this time of the year, in order to avoid underestimation of the climatological maximum MLD due to evanescent very shallow MLDs that occur under transient calm conditions, we retain for the fitting only the maximum MLD reached each year during this period (filled circles in Figure 3b). Using these data, the climatological seasonal cycle of the MLD is fitted by using a polynomial of third degree, as done by Gonzalez-Pola et al. [2007] (the orange line in Figure 3c). If we had used all available data instead, the maximum winter MLD would be biased toward the shallower values in the range of deepest MLDs reached at the end of the winter (Figure 3c). In order to explore trends, this exclusion of transient shallow MLDs does not make a difference.
Climatological seasonal cycles are also commonly obtained from weekly, biweekly, or monthly average values. For signals with a smooth variation like SST or StrI r , the differences among these approaches are not particularly crucial when a large amount of data are available (i.e., both curves almost overlap, Figure  3a), but for the MLD the use of such averages may introduce important deviations and artifacts. As shown in Figure 3c, this approach provides a maximum winter MLD biased toward the shallowest value in the range of deepest MLDs reached at the end of the winter and artificially introduces a shallowing of the MLD during the period when the MLD is not unequivocally defined. Such shallowing as a physical process may occasionally occur, but in general the MLD simply alternates between very shallow and deep values (Figure 3b).
MLD long-term trends are calculated as linear regression of MLD anomaly time-series. This anomaly timeseries is the sequence of observed minus fitted seasonal cycle values. Regarding the ocean ventilation and ocean productivity, the winter conditions are especially relevant, so we also calculate the winter and summer MLD long-term trends. These are calculated as linear regression of winter and summer MLD anomalies against time (e.g., only MLD anomalies during the winter (December, January, and February) or the summer (June, July, and August) are retained for the regression).

Winter MLD Trends and Its Forcing
In addition to the analysis of the long-term changes and trends, we investigate the role of the different drivers of winter MLD variability. Destabilizing buoyancy forcing from cooling and/or evaporation making the ocean surface colder and saltier induces convective mixing of surface water with deeper water. Wind stirring always causes vertical turbulence in the upper mixed layer; however, it is less effective than convective mixing to thicken the mixed layer, and the autumn-winter mixed layer deepening is controlled by the buoyancy forcing [Kantha and Clayson, 2000;Alexander et al., 2000;Somavilla et al., 2011]. Because the ocean mixed layer responds so rapidly to surface densification by buoyancy-forced processes, the deepening of the mixed layer can often be estimated successfully using one-dimensional approaches as the nonpenetrative convection equation [Marshall and Schott, 1999]: where the vertical extent of the mixed layer (h, MLD) results from the integration in time of the buoyancy forcing (B 0 ) assuming the water column stratification (N 2 , Brunt-Vaisala frequency) has a constant value with depth. Using this expression, the expected changes in winter MLD in each of the study areas according to changes in (i) winter buoyancy fluxes estimated from NCEP/NCAR reanalysis data and (ii) stratification observed in oceanographic time-series are calculated (for detailed information, see section 1.2.1 in supporting information).
Horizontal convergences/divergences of the water masses due to wind stress curl generate downward/ upward nonturbulent vertical velocities that pump (Ekman pumping) the upper ocean surface layer. This pumping generates a vertical displacement of the interface that defines the base of the mixed layer, displacement that thickens mixed layers in convergent zones and thins in divergent zones [de Szoeke, 1980;Carranza and Gille, 2015]. The effects of vertical advection due to Ekman pumping on winter MLD are considered here by calculating the vertical velocities (w EK ) resulting from horizontal convergences/divergences of the water masses due to wind stress (s) curl using NCEP/NCAR reanalysis data (positive wEK generates downwelling and thickening the MLD, and negative wEK upwelling and thining of the MLD) (for detailed information, see section 1.2.2. in supporting information).

SST, Stratification, and MLD Seasonal Cycles
As introduced above, our paper focuses on MLD and stratification long-term changes at midlatitude ocean basins of the Northern Hemisphere. Before addressing this main topic, we briefly consider the ability of the different data sources and methods for representing the amplitudes and phases of the MLD and stratification seasonal cycles for the different study areas.
The time series of SST, stratification measured as StrI r , and MLD in the subtropical gyres of the North Pacific (sg-NPac) and North Atlantic (sg-NAtl) and in the midlatitudes of the eastern North Atlantic (ml-ENA) show, as expected, marked seasonal cycles ( Figure 4). For SST and StrI r , the cycles start to increase at the beginning of spring after reaching their annual minima. They continue increasing until the end of the summer when they reach their maxima and decrease in autumn and winter until they reach again their annual minima. Conversely to SST and stratification, the MLD evolution does not have a smooth annual cycle. Its annual cycle starts at the surface at the beginning of spring, remains more or less stable during late spring and summer, and deepens progressively during autumn and winter until the deepest MLD is reached before any distinct signature vanishes. The days of maximum and minimum climatological MLD are, therefore, consecutive. A shallowing of the MLD is not observed as SST and StrI r increase, and the MLD actually deepens slowly during late spring and summer. In Figure 4 only MLD GP is shown but climatological seasonal cycles from MLD HT and MLD Dr show similar results (supporting information Figure S2). This agreement among MLD determination methods is also found for long-term changes presented in the next section. In order to avoid redundancy of results and considering that the Gonzalez-Pola et al. provides information of the upper ocean vertical structure and stratification that will be used in the discussion, we will not show the figures for MLD HT and MLD Dr here. These can be found in supporting information.
The minimum and maximum values of the climatological seasonal cycles of SST from OISST, StrI r and the different MLDs indexes estimated from oceanographic time-series and the dates when they are reached (Date min and Date max ) are summarized in Table 1 (similar tables for estimates from ORAS4 and Argo floats time-series are summarized in supporting information Tables S1 and S2). For SST, the highest temperatures are reached in the sg-NPac (25:58 C) while the lowest temperatures and greatest amplitude of the seasonal cycle are found in the ml-ENA (between 12:66 C and 20:14 C). For stratification measured as StrI r , the greatest values and amplitude are reached in the sg-NAtl (between 0:14 and 2:84 kg m 23 ) and the lowest in the ml-ENA (0:12 kg m 23 ). The climatological MLD seasonal cycles measured as MLD GP from oceanographic time-series range between the surface and approximately 100 m in the sg-NPac, 220 m in the sg-NAtl, and 180 m in the ml-ENA ( Figure 4 and Table 1). Thus, climatologically, the deepest winter MLDs are found in the sg-NAtl.
Despite the different spatial and temporal resolution of oceanographic time-series, ocean reanalysis, and Argo floats, the amplitude of the StrI r seasonal cycles from the different data sources differs by less than 6% (see Table 1 and supporting information Tables S1 and S2) showing also similar phases as seen in Figure 4. For MLDs, just considering the seasonal cycles from observations and ORAS4, their amplitudes differ by less than 4% for MLD GP , 10% for MLD HT , and 5% for MLD Dr . It is interesting to highlight that the climatological maximum winter MLDs estimated with Argo floats are on average 30 m deeper than estimated on oceanographic time-series or ocean reanalysis independently of the method used for the MLD determination.
Date min Date max sg-NPac sg-NAtl ml-ENA  (2003-2015) gives a climatological maximum MLD that is 20 m deeper than for the period 1990-2015. Thus, using the same sampling period, the difference between the climatological maximum winter MLDs from observations and Argo is reduced to 10 m, which is within the expected uncertainty due to the lower vertical resolution of Argo floats. From this, we rule out that the different vertical resolution is responsible for the 30 m deeper climatological maximum MLDs in Argo floats and conclude that it is due to their different, later sampling period.

SST, StrI r , and MLD Long-Term Changes and Trends
Long-term changes and trends are analyzed over SST, StrI r , and MLD anomaly time-series ( Figure 5). Anomaly time-series are obtained as the difference between observed and fitted climatological seasonal cycle values. As shown in the previous section, the different data sources reproduce adequately the climatological seasonal cycles, and so we rely on them for the investigation of low-frequency changes.
Low-frequency changes derived from oceanographic and reanalysis time-series seem consistent, and this inference is supported by similar long-term trends obtained from both data sources. Data from Argo floats are not shown in Figure 5 since they do not span the same sampling period, and the trends obtained are not comparable with those from ocean observations and reanalysis data.
SST warming trends are found in all the three oceanic areas. The warming trends in the sg-NPac and sg-NAtl are of 0:09 and 0:15 C decade 21 , respectively. Both trends are significant at 0.05 significance level and are within the average in the global ocean (0:09-0:13 C decade 21 in the upper 75 m) [Alexander et al., 2013]. In the ml-ENA, the surface warming is less pronounced (0:03 C decade 21 ) due to intense cooling events associated with the occurrence of winter blocking events since 2005 in the area [Somavilla et al., 2009[Somavilla et al., , 2011[Somavilla et al., , 2016] that were not compensated by a slight summer temperature increase. Concurrent with surface warming, stratification measured as StrI r is observed showing increasing trends in the sg-NAtl (10:05 kg m 23 decade 21 ) but also decreasing trends in the sg-NPac (20:16 kg m 23 decade 21 ). The moderate surface warming in the ml-ENA coincides with relatively moderate changes in StrI r (10:004 kg m 23 decade 21 )-1 order of magnitude lower than that observed in the sg-NAtl.
Regarding the MLD long-term changes, MLD deepening trends ranging from 1 to 8 m decade 21 are found at the three sites ( Figure 5) and by all three different MLD determination methods MLD GP , MLD HT , and MLD Dr (Figure 5 and supporting information S2). Mostly due to the alternation of very deep to very shallow values of the MLD between the end of the winter and beginning of the spring (see section 2.3), the anomaly time-series are very noisy. Still the trends are significant both in ocean observations and reanalysis data except for the ml-ENA, and this can be more easily recognized on the MLD time-series shown in Figure 6 with the trend shown by the white line. The deepening trends in the sg-NPc and sg-NAtl are consistent in ocean observations and ORAS4 reanalysis data and significant at 0.05 significance level (Table 2). In the ml-ENA, deepening trends are found in ocean observations by the three MLD determination methods MLD Dh , MLD GP , MLD HT (Table 2), but they are not significant and are not found in ORAS4 reanalysis data.  Separately, due to the crucial role of deep winter MLDs in ocean ventilation and global biogeochemical cycles, we also investigate the winter MLD long-term trends. As shown in Figure 6 and summarized in Table  2, winter MLDs are getting deeper at rates between 2.5 and 5 times faster than the annual trends. The faster winter MLD deepening is found for all three study areas by all three MLD determination methods. However,  The blue and red lines represent, respectively, the winter and summer MLD long-term trends calculated as linear regression of MLD anomalies against time using only the winter (blue circles) or summer (red circles) data. With explanatory purposes, the trends (white, blue and red lines) -although calculated as linear regression of MLD anomaly time-series against time-are made to vertically coincide with the actual average annual, winter or summer MLDs observed at each of these sites (white, blue, and red dots), summing up to the long-term trends the corresponding long-term average value (e.g: winter mean value for winter MLD trends or summer mean value for summer MLD trends). The waters above/below the MLD are shown in light/dark blue.

Journal of
winter MLD deepening trends are again significant only in the sg-NPac and sg-NAtl, not in the ml-ENA data.
In view of Figure 6c, in the ml-ENA rather than a trend, a jump toward much deeper winter MLDs is observed from 2005 onward when intense cooling events occurred in the area [Somavilla et al.,2009[Somavilla et al., , 2016. Unlike the winter MLDs, summer MLDs are deepening 2.5-5 times slower than the annual trends; MLDs are even getting shallower in the case of the sg-NAtl (Table 2). Thus, the winter MLD deepening trends at rates between 10 and 17 m decade 21 dominate the MLD deepening trends.
Argo floats trends cannot be used to provide further support to the trends in oceanographic time-series and ocean reanalysis. However, the finding of climatological maximum winter MLDs on average 30 m deeper when estimated from Argo floats than from the oceanographic time-series or the ocean reanalysis is nevertheless indicative of progressively deeper winter MLDs. As explained in the previous section, this finding cannot be explained by the different vertical resolution of Argo data and seems to be due to its shorter and more recent period of sampling (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015).
Overall, the finding of warmer and deeper MLDs relies on the estimation of MLDs by several objective methods from hydrographic profiles both from long-term hydrographic time-series and ocean reanalysis; this finding is consistent with Argo floats within the period of overlap. Accordingly, it should not be attributable to a particular methodology or only noticeable at specific oceanic stations. The maps of MLD long-term trends in the sg-NPac, sg-NAtl, and ml-ENA in Figure 7 confirm that the MLD deepening trends found at the location of oceanographic time-series by the different data sources and determination methods are representative of larger oceanic regions.

StrI r : A Good Measure of Stratification?
Based mainly on model results, ocean surface warming is commonly associated with increased stratification and shallower mixed layers. In contrast to this assumption, we find that stratification is not unequivocally increasing, and MLD is not shoaling, despite evident and substantial ocean warming.
The results from the climatological seasonal cycles calculated here also do not support the widespread idea that the warmer the ocean surface, the shallower the mixed layer. During spring and summer stabilizing buoyancy forcing from heating-causing the ocean surface to become warmer-stratifies the surface and isolates it from the deeper waters . Thus, stabilizing buoyancy forcing hinders upper ocean mixing, and during the summer months wind-induced mixing controls the vertical extension of the mixed layer [Alexander et al., 2000;Kantha and Clayson, 2000;Somavilla et al., 2011]. Since different processes govern the variabilities of SST and MLD during summer months, MLD does not necessarily shoal as SST and StrI r increase; MLD remains more or less stable during late spring and summer. Consequently, the apparent shallowing of the MLD found in MLD climatological seasonal cycles obtained from weekly, biweekly, or monthly averaging (Figure 3c and section 2.3) cannot be explained by physical process of buoyancy and wind forcing. Such shallowing and climatological maximum MLDs shallower than actually observed-artefact of the  It is true that during autumn-winter destabilizing buoyancy forcing from cooling makes the ocean surface colder and the MLD deeper through convective mixing. However, the comparison of the different climatological seasonal cycles shows that the lowest temperature and the minimum stratification values (StrI r ) at the end of the winter are in the ml-ENA but not the climatological deepest winter MLD which is found in the sg-NAtl. SST is often used as index of physical changes within the surface mixed layer including its vertical extension, having been used to identify permanently stratified zones [Behrenfeld et al., 2006[Behrenfeld et al., , 2008O'Malley et al., 2010; Gonz alez Taboada and Anad on, 2014] characterized by low primary productivity. However, this approach has been reported to be subjected to several caveats, since SST changes cannot be solely attributed to shoaling or deepening of the MLD at synoptic scales and is influenced by mesoscale processes and horizontal advection [Carranza and Gille, 2015;Dave and Lozier, 2015].
Similarly, the finding of surface warming concurrent with decreasing stratification measured as StrI r in the sg-NPac is not surprising ( Figure 5). The commonly involved link between increasing stratification and surface warming trends is supported by the seasonal coupling between stabilizing buoyancy forcing from atmospheric heating and the warming of the ocean surface that stratifies the upper layer. However, stabilizing/destabilizing buoyancy forcing not only occurs due to heating/cooling causing the ocean surface to become warmer/cooler but also due to precipitation/evaporation causing the ocean surface to become fresher/saltier . The study areas are subjected to changes in the evaporationprecipitation balance [Alexander et al., 2013] that can compensate for the expected surface density decrease due to surface warming. This is observed in the sg-NPac where surface density increase (0:08 kg m 23 decade 21 ) due to salinity changes (0:13 decade 21 ) outpace the temperature contribution to making the surface waters lighter (for density and salinity the linear trends are also calculated through a harmonic fit as for SST and StrI r ). At depth, however, salinity is decreasing (20:012 decade 21 ) making the waters at depth lighter (20:02 kg m 23 decade 21 ), decreasing the density difference between the surface and 200 m (StrI r ) thought the ocean surface is warming. On the other hand, even with temperature changes dominating the density changes, the expected finding of increasing stratification concurrent with surface warming assumes that warming is faster at the ocean surface than at depth. The properties at depth are controlled from below through advective pathways of remotely formed water masses . Faster warming at water formation regions and changes in circulation can result in warming trends faster at depth than at the surface as observed in the ml-ENA (not significant trend at the surface versus a warming trend at depth of 0:16 C decade 21 ) [Somavilla et al., 2011]. A somewhat similar situation has been found in CMIP5 models under mitigation scenarios, when a transient imbalance between cooling at the ocean surface and continued warming in subsurface waters takes place [John et al., 2015].
Finally, the finding of increasing StrI r trends concurrent with MLD deepening trends is perhaps more troubling. The index of stratification StrI r shown in the above figures is the most widely used. However, beyond its simplicity, the rationale for its use is questionable. A good stratification index must represent a measure of the tendency of a water column to be mixed. The deepening of the mixed layer is inversely proportional to the density difference between the mixed layer and the water below it (Dr bMLD ) according to the equation for the time evolution of the turbulent kinetic energy (TKE) [Barton et al., 2015;Follows and Dutkiewicz, 2002]. With an identical density difference between the surface and 200 m, and so identical StrI r , will be more difficult to mix a water column with this density difference concentrated at a specific depth (step-like profile in Figure 2, Dr bMLD 5 StrI r ) than a water column with this density difference uniformly distributed from the surface to 200 m (continuous type profile in Figure 2, Dr bMLD 5 StrI r =200 m). In other words, the same forcing from wind and cooling will mix deeper the water column in which the density difference between the surface and the 200 m level is evenly distributed than the one in which the same density difference is concentrated at a specific depth. Thus, a good stratification index requires information about the vertical distribution of the density difference between the surface and the deeper layers that StrI r does not include.
The search for a stratification index that provides better representation of the tendency of a water column to being mixed needs to consider all kinds of stratification situations that can be found in the ocean, from strongly to very weakly stratified areas. Such a specific study is clearly beyond the scope of this study, which is focussed on midlatitude ocean regions. As a preliminary analysis (results not shown), we have examined three stratification indexes that quantify in different ways the vertical distribution of the density difference between the surface mixed layer and the waters below. These alternative stratification indexes are: the maximum stratification of the water column, G max ( Figure 2b); the slope of the steepness, b 3 ( Figure 2a); and the potential energy U criteria that integrates the vertical density distribution [Simpson et al., 1977;Huret et al., 2013]. Note that these different stratification indexes do not measure the same feature of the upper ocean vertical structure. Hence, we do not find agreement of results regarding the stratification seasonal cycles and trends from these indexes, unlike MLD (in which the different determination methods all try to estimate the vertical extent of the mixed layer). G max temporal variability is dominated by its seasonal cycle as for SST and StrI r . On the other hand, b 3 and U show distinct, low-frequency variability resembling that depicted from winter MLDs at the three study areas. While G max and StrI r seem indexes representing the stratification of the seasonal pycnocline, b 3 and U seem to be more representative of the stratification through the permanent pycnocline. Thus, our analysis from midlatitudes does not allow us to identify one stratification index as clearly better than the others. For that purpose, a global study that considers all kind of stratification situations that can be found in the ocean is necessary.

Winter MLD Deepening Trends and Their Forcing
In the previous section, we explained why ocean warming does not necessarily imply an increase of stratification and why changes in stratification measured as StrI r should not be expected to reliably inform about changes in the vertical extent of mixing. However, we still have to provide an explanation for the observed MLD deepening trends. Deepening MLD trends seem to be mostly driven by the faster trends of winter MLDs. We discuss here the roles of changes in the buoyancy forcing (densification of the ocean surface due to cooling and/or evaporation) and wind-driven divergence/convergence (Ekman pumping) on winter MLD Journal of Geophysical Research: Oceans 10.1002/2017JC013125 deepening trends. In order to facilitate the understanding of this discussion, we first briefly describe the limitations of these estimates. A schematic of the processes considered is shown in Figure 8.
The expected changes in winter MLD in every study area due to changes in (1) winter buoyancy fluxes (B 0 ) estimated from NCEP/NCAR reanalysis data and (2) stratification (N 2 ) observed in oceanographic time-series are calculated according to the non-penetrative convection equation [Marshall and Schott, 1999] (equation (2)). Although winter mixing is mostly governed by the accumulated effect of buoyancy (heat) loss throughout the winter, extraordinary convection episodes appear more likely to occur when heat loss is concentrated in intense mixing events, rather than distributed evenly but more weakly throughout the winter [Marshall and Schott, 1999;Vage et al., 2008;Somavilla et al., 2011]. The time integration of B 0 used in equation (2) does not allow us to reproduce specific events deviating largely from the long-term mean value of winter MLDs. However, we are not interested in explaining specific years but the long-term changes, and this expression seems to provide plausible estimations of the MLD long-term mean and trends as argued below.
Additionally, we calculate the enhanced deepening/shoaling of winter MLD that can occur in downwelling/ upwelling regions due to changes in the Ekman pumping. Basically, we calculate the vertical velocities (w EK ) from horizontal convergences and divergences of the water masses due to wind stress curl using the NCEP/ NCAR reanalysis data and integrate throughout the winter the resulting vertical displacement (Dh EK ) that this pumping generates (Figures 8b and 8c). Vertical velocities due to Ekman pumping and the resulting vertical displacement of MLD may seem small and without effect on MLD variability (e.g., daily average of w EK 5 3.5 10 26 m s 21 and Dh EK 5 0.3 m). However, at seasonal time-scales the effects of Ekman pumping on MLD variability are noticeable and when integrated (accumulated) throughout the winter can drive the MLD to be tens of meters deeper than expected from the buoyancy forcing in convergent zones and vice versa in divergent zones [de Szoeke, 1980] (Figure 8a). With all that in mind, the lower plots of Figure 9 show that changes in the buoyancy forcing and stratification included in the non-penetrative equation do not completely explain winter MLD trends since the trends estimated from the buoyancy forcing (red lines in Figures 9h, 9i and 9j) do not coincide with the observed ones (blue lines). The average difference between the line observed and the line estimated from the buoyancy forcing seems to coincide in each study area with the expected integrated vertical displacement due to Ekman pumping (Dh Ek ) (Figures 9a and 9c) (in other words, adding the estimated effect of Ekman pumping displaces the red line to coincide with the observed blue line, Figures 8b1 and  c1). Thus, in the sg-NAtl area the average winter Dh Ek is 223 m (Figures 9c and 9f) coinciding approximately with the downwards displacement of 27 m that would make the mean winter MLD estimated from buoyancy forcing (165 m) to reach the average value found in the observations (192 m) (Figure 9i). In the sg-NPac, the average winter Dh Ek within the HOTS station location is positive (124 m) indicating that divergence is pushing the water column up in contrast to the rest of the subtropical gyre ( Figures  9a and 9e). This value of winter Dh Ek is also in the range of the upward displacement of 37 m that would make the mean winter MLD estimated from buoyancy forcing (135 m) reach the average value found in the observations (98 m) (Figure 9h). In the ml-ENA, observed and estimated trends from buoyancy fluxes almost overlap in the vertical, in agreement with no significant upward or downward Dh Ek in the area (Figure 9c). In the sg-NAtl, the downward velocities are increasing, and so the integrated downward displacement through the winter is also increasing, at a rate of 7:6 m decade 21 (Figures 9d and 9f). This compensates and outpaces the expected shallowing of MLD due to changes in the buoyancy forcing and stratification (21:1 m decade 21 , red line in Figure 9i) enabling the MLD to deepen with time as observed (situation scheme Figure 8b2). In the sg-NPac, upward w Ek are decreasing (Figures 9b and 9e), contributing to the faster deepening trends found in observations (16 m decade 21 ) with respect to the deepening MLD trends expected from changes in the buoyancy forcing and stratification (7:8 m decade 21 ) (situation scheme Figure 8c2). In the ml-ENA, as observed for winter MLDs, rather than a trend a regime shift is observed after 2005 (Figure 9g) in correspondence with the severe changes in the vertical extent of mixing in the area.
Evidence of the relationship between winter MLD variability and changes in that season's air-sea heat and salt exchanges (buoyancy forcing) is common in the literature [Carton et al, 2008;Vage et al., 2008;Somavilla et al., 2011] but the relationship of long-term changes in MLD with changes in the Ekman pumping has not been recognized. However, changes in wind patterns including a poleward migration and wave amplitude and asymmetry changes of the maximum westerlies in both Northern and Southern hemispheres that may result in such changes in wind-driven divergence and Ekman pumping have been extensively reported [Russell et al., 2006;Sall ee et al., 2010;Archer and Caldeira, 2008;Capua and Coumou, 2016]. Similarly, circulation changes in the subtropical gyres estimated from Sea Surface Height (SSH) that are a reflection of Ekman pumping changes have also been described [McClain et al., 2004;Roemmich et al., 2007]. The maps of average winter vertical displacement due to Ekman pumping (Dh Ek ) shown in Figure 9a and c are consistent with those found in the literature for both the Atlantic and the Pacific [Marshall and Nurser, 1993;Qiu and Huang, 1995]. Besides, the estimated deepening trends of the MLD due to Ekman pumping (7 m decade 21 ) are similar to those found for the Atlantic and Pacific from SSH changes. The 10:5 cm yr 21 SSH (g) trends found in these areas would be equivalent to a 5 m decade 21 deepening of the pycnocline (D) according to the expression g 0 D 2 =2f 5gg 2 =2fg 0 where f is the Coriolis parameter and g' the reduced gravity (10 22 m s 22 ) [McClain et al., 2004]. These trends due to Ekman pumping are not usually applied to the MLD. However, since they pump the entire upper layer vertical structure, they inevitably also affect the MLD, as shown in this study, being able to enhance, compensate, and even outpace the effects of ocean warming and increasing stratification on MLD deepening trends.
Overall, changes in the buoyancy forcing and stratification considered in the non-penetrative equation together with wind-driven divergence changes provide an explanation for the finding of deepening winter MLD trends concurrent with surface warming in the three study areas. Applying the nonpenetrative equation for N 2 constant over time (the long-term average value of N 2 in each study area), the estimated winter MLD trends would be of 1:1 m decade 21 in the sg-NPac (versus 7:8 m decade 21 considering changes in stratification); 5 m decade 21 in the sg-NAtl (versus 21:1 m decade 21 considering changes in stratification); and 10 m decade 21 in the ml-ENA versus (9 m decade 21 considering changes in stratification). Thus, both changes in the horizontal convergences/divergences of the water masses due to wind stress curl (circulation) and in the buoyancy forcing contribute to deepening the winter MLDs in all the three areas. Changes in stratification seem to modulate the response of the upper ocean to the atmospheric forcing, enhancing deepening trends in the sg-NPac and damping them in the sg-NAtl.

Biogeochemical Implications
According to our findings, ocean observations do not support an unequivocal relation between ocean surface warming, increasing stratification, and shallower MLDs. The existence of such a link in model simulations under global warming scenarios has been extensively used to explain observed or projected declines in global oceanic oxygen content and primary productivity [Behrenfeld et al., 2006;Polovina et al., 2008;Boyce et al., 2010;Tyrrell, 2011;Gruber, 2011;Capotondi et al., 2012;Schmidtko et al., 2017]. How do our findings reconcile deepening MLDs with the observation of less productive and oxygenated oceans?
Global satellite observations have been widely used to document the existence of an inverse relationship between changes in SST and surface phytoplankton chlorophyll concentrations (a proxy for phytoplankton biomass). However, rather than being a universal relationship, it seems to be restricted to satellite observations for specific areas [see e.g., Ram ırez et al., 2017, supporting information Figure S1] and time-scales of variability. As pointed by Lozier et al. [2011], ocean surface warming, increasing stratification, and reduced primary productivity are linked by the ''straightforward'' relationship existing between the physical Journal of Geophysical Research: Oceans 10.1002/2017JC013125 processes and primary productivity at seasonal time-scales. The problem is that no such relationship exists at longer, decadal time-scales [Lozier et al., 2011]. Besides, it has been shown that global-scale correlations between upper ocean stratification and chlorophyll a are driven by strong associations between the two properties in the central and western equatorial Pacific [Dave and Lozier, 2013], and that these correlations are not due to SST changes but result from the advection of heat, salt, and nutrients in the region [Dave and Lozier, 2015].
In situ primary production observations also do not support such a definitive relationship either: increases of primary production have been found concurrent with surface warming trends at BATS and HOTS oceanic stations [Wallhead et al., 2014;Saba et al., 2010;Corno et al., 2007]. Interestingly, biogeochemical ocean general circulation models were tested to reproduce this finding, and only one of twelve models compared for each site was able to reproduce an increase of primary production concurrent with the surface warming trends [Saba et al., 2010]. The deeper winter-time MLDs and more temperate nature of the ecosystems were proposed as possible explanations for the low skill of the models [Saba et al., 2010].
The origin of these findings of increasing primary production together with surface warming trends from specific oceanic stations [Wallhead et al., 2014;Saba et al., 2010;Corno et al., 2007] has been used to claim that the station data cannot necessarily be extrapolated to larger spatial scales. However, we have demonstrated that the finding of warmer and deeper MLDs in long-term oceanographic timeseries is representative of larger oceanic areas, and so it can also be the case for their effects on the availability of nutrients and primary productivity. The net effect of deepening winter MLDs on the nutrient supply to the euphotic zone enhancing primary productivity depends on the drivers of this deepening. Deepening winter MLDs as a result of intensified buoyancy losses and/or weaker stratification would increase convective mixing and so entrainment of nutrient-rich deep waters into the upper layer. This nutrient enhancement is expected to be accompanied by an increase in primary production. Deepening winter MLDs due to changed wind-driven divergence would not necessarily entrain more nutrient rich deep waters into the upper layer, since the Ekman pumping would push the whole water column down, not solely the MLD. Just the maintenance of the nutrient supply may lead to an increase of primary production in a warmer upper ocean under the assumption of a bottom-up control and nutrient-limited production. However, in this case a dilution effect on phytoplankton due to deeper MLDs should be also considered. In agreement with these expectations, in the vicinity of HOTS oceanographic time-series, where winter MLD deepening trend is more strongly dominated by changes in the buoyancy forcing and stratification, the enhancement of nutrient supply, and primary production is larger than at BATS oceanographic time-series [Saba et al., 2010].
Recently, the discrepancy between satellite and in situ phytoplankton biomass observations at global scales has been examined, showing that contemporary relationships between chlorophyll changes and ocean warming are not indicative of proportional changes in productivity, as light-driven decreases in chlorophyll can be associated with constant or even increased photosynthesis [Behrenfeld et al., 2016]. Still this explanation seems to implicitly require the occurrence of shallower MLDs to cause the light-driven decreases in chlorophyll concurrent with surface warming. As explained through this paper, warmer MLDs are not necessarily shallower, and further investigation is needed.
The global oceanic oxygen content decline detected over the past few decades has been attributed in ocean models, and so in observations, to warming-induced declines in oxygen solubility, and reduced ventilation of deeper waters from enhanced upper-ocean stratification [Helm et al., 2011;Schmidtko et al., 2017]. However, our results do not necessarily contradict the finding of oceans losing breath, because these results may indicate that ventilation increases for warmer waters which have lower capacity to contain dissolved gasses, including oxygen. The oxygen decline in this case should be similar to the oxygen decrease due to warming-induced solubility losses since ventilation is not reduced. Interestingly, Schmidtko et al. [2017] found the oxygen decline in the depth range of mode waters (300-700 m)-to which the midlatitudes areas of the ml-ENA and sg-NAtl would contribute to ventilate-coinciding with the expected decrease solely due to solubility changes. The waters above and below this depth range show oxygen declines larger than the expected decrease due to solubility changes that are attributed to reduced ventilation.

Conclusions
We have shown that surface warming is not definitively linked to a more stratified ocean with shallower mixed layers. The previously assumed relationship between increasing stratification and surface warming assumes that warming is faster at the ocean surface than at depth which is not necessarily the case, as found here for the midlatitudes of the Eastern North Atlantic. On the other hand, even when warming is faster at the surface, salinity changes can compensate for the expected surface density decrease due to surface warming. Moreover, even when the density difference between surface waters and those at depth increases, MLD can get deeper with time. The vertical extent of mixing depends not only on stratification but also on wind-driven divergence and buoyancy losses. These can compensate for and outpace the effects on MLD of ocean warming and increasing stratification. Changes in both wind-driven divergence and in the buoyancy forcing have been found to contribute to deepening the winter MLDs at the three study areas at midlatitudes considered in this work.
Overall, we have found that while SST increases at three study areas at midlatitudes, stratification both increases and decreases, and MLD deepens with enhanced deepening of winter MLDs at rates from 10 to 17 m decade 21 . These results rely on the estimation of several MLD and stratification indexes of different complexity from hydrographic profiles collected by long-term oceanographic time-series, by ocean reanalysis and by Argo floats. Because of the multiple data sources, the apparent deepening of winter MLDs cannot be attributed to a particular methodology or a specific oceanic station. Instead of that, the finding of warmer and deeper MLDs is representative of large oceanic regions and the likely mechanism could operate in other ocean regions [Martinez et al., 2016]. These conclusions may be extended to the expected enhancement of the nutrient availability and primary productivity due to deeper MLDs. Besides, warmer and deeper MLDs may indicate that subsurface waters are ventilated more with warmer waters which have lower capacity to contain dissolved gasses providing an alternative explanation for the oceanic oxygen content decline in some regions as midlatitudes areas of mode water formation.
Finally, it must be repeated that caution must be taken when using the density difference between the sea surface and 200 m (StrI r ) as stratification index. StrI r is the most widely used stratification index, but it does not provides a good representation of the tendency of a water column to be mixed since it does not inform about the vertical distribution of the density difference between the surface and the deeper layers.