How subsurface and double-core anticyclones intensify the winter mixed layer deepening in the Mediterranean sea

. The mixed layer is the uppermost layer of the ocean, connecting the atmosphere to the subsurface ocean through atmospheric ﬂuxes. It is subject to pronounced seasonal variations: it deepens in winter due to buoyancy loss and shallows in spring while heat ﬂux increase and restratify the water column. A mixed layer depth (MLD) modulation over this seasonal cycle has been observed within mesoscale eddies. Taking advantage of the numerous Argo ﬂoats deployed and trapped within large Mediterranean anticyclones over the last decades, we reveal for the ﬁrst time this modulation at a 10-day temporal scale 5 and free of the smoothing effect of composite approaches. The analysis of 16 continuous MLD time series inside 13 long-lived anticyclones at a ﬁne temporal scale brings to light the importance of the eddy preexisting vertical structure in setting the MLD modulation by mesoscale eddies. Extreme MLD anomalies of up to 330m are observed when the winter mixed layer connects with a preexisting subsurface anticyclonic core, greatly accelerating mixed layer deepening. The winter MLD sometimes does not achieve such connection but homogenizes another subsurface layer, then forming a multi-core anticyclone with spring 10 restratiﬁcation. A MLD restratiﬁcation delay is always observed, reaching more than 2 months in 3 out the 16 MLD timeseries. The water column starts to restratify outside anticyclones while mixed layer keeps deepening and cooling at the eddy core for a longer time. These new elements provide new keys for understanding anticyclone vertical structure formation and evolution

pointed out by Villas Bôas et al. (2015) and Hausmann et al. (2017) -is their focus on surface-intensified eddies with the most coherent surface signature. Indeed the relation between eddy SSTA and SSH amplitude strongly relies on the hypothesis of a surface-intensified structure and Assassi et al. (2016) showed that it should not be the case for subsurface anticyclones.
Subsurface eddies are mesoscale structure were the density anomaly (compared to the outside-eddy density profile) is overlaid 60 by an anomaly of opposite sign. For instance a subsurface anticyclone has a lighter core at depth overlaid by a negative density anomaly near the surface. Following the thermal wind equation, depth of maximal geostrophic speed is below the surface. The isopycnals and isotherms doming above a subsurface anticyclone core could greatly impact the upper layer stratification, and subsequently the inside-eddy mixed layer dynamics. In the South China sea and still using the composite method, He et al. (2018) found the anticyclones to be predominantly subsurface eddies.They also observed a linear trend between the subsurface 65 temperature anomaly and SSH on an annual average, and an eddy-induced MLD anomaly but based on monthly climatology and did not find a relationship with MLD. A more detailed comparison with more observational data and especially better temporal resolution is then lacking.
The Mediterranean sea is an interesting region to study eddy influence on MLD. At first, due to repeated oceanographic 70 campaigns the density of in-situ measurements is relatively high and in particular several campaigns specifically targeted longlived mesoscale anticyclones in both Western and Eastern basins. Without aiming for exhaustiveness, one can in particular list: EGYPT (Poulain et al., 2006), BOUM (Moutin and Prieur, 2012), 'Eye of the Levantine' (Hayes et al., 2011), PROTEVS (Garreau et al., 2018) and PERLE (Ioannou et al., 2019). Additionally several Argo profiling floats were launched inside eddies and remained trapped for a long time (Ioannou et al., 2020), altogether allowing in the Mediterranean sea to accurately follow 75 particular eddies and go beyond the averaged composite vision. Moreover data from these programs were often only analyzed in the scope on the campaigns, and an eddy study with a larger statistical focus is still lacking. A second relevance for Mediterranean eddies is the variety of mesoscale structures in terms of dynamics, from intense Algerian and Ierapetra anticyclones needing cyclogeostrophic corrections (Ioannou et al., 2019) to subsurface eddy with strong density anomalies but weak SSH signature (Hayes et al., 2011). Moutin and Prieur (2012) also showed the vertical structure, in temperature and salinity, of 80 mesoscale eddies to be very different from one basin to another. Barboni et al. (2021) showed the marked subsurface difference between a new anticyclone detached from the coast compared to an offshore structure having been tracked more than a year.

All these structures should provide different examples of eddy-MLD interactions.
In the Mediterranean sea there is in addition a strong asymmetry between cyclones and anticyclones, remarkable in lifetime 85 difference (Mkhinini et al., 2014). The deformation radius in the Mediterranean sea is indeed about 8 to 12 km (Kurkin et al., 2020), and cyclones are less stable when greater than the deformation radius and more subject to external shear (Arai and Yamagata, 1994;Graves et al., 2006). This leads to cyclones being predominantly below the effective resolution of SSH products about 20km (Stegner et al., 2021). As a consequence anticyclones are coherent larger vortices, while cyclones in the Mediterranean sea as detected by altimetry are instead cyclonic gyres bounded by topography or hydrographic fronts such as the 90 Ligurian, South-Western Crete or Rhodes gyre (Stegner et al., 2021). MLD evolution inside these cyclonic gyres was already surveyed because of their importance for biological production, in particular with the development of BGC-Argo (d 'Ortenzio et al., 2021;Taillandier et al., 2022). Apart for specific campaigns, Mediterranean anticyclones remain poorly analyzed despite being more coherent, and statistical comparison based on vertical profiles is lacking, with the noticeable exception of the BOUM campaign surveying 3 anticyclones across the Mediterranean in 2008 (Moutin and Prieur, 2012). 95 This paper aims to study the temporal evolution of the mixed layer inside a wide diversity of long-lived anticyclones in the Mediterranean sea compared to the evolution of the background MLD. The goal is to quantify more precisely the local impacts of individual eddies on the winter mixed layer deepening. The paper is organized as follows. Section 2 describes the eddy detection and tracking algorithm and the in-situ profiles database. Section 3 details the methodology used to compute the 100 MLD, colocalize profiles and eddies in order to quantify accurately the MLD anomalies induced by individual eddies. In the section 4 we analyze the MLD evolution at anticyclone cores, provide statistical analysis over the variety of structures surveyed and discuss the impact of complex vertical eddy structure on winter mixed layer deepening. Finally in the last section 5 we discuss the possible physical drivers and implication for these MLD anomalies. is a mixed velocity-altimetry approach, its relies on using primarily streamlines from a velocity field and identifying possible eddy centers computed as maxima of local normalized angular momentum (Le . It was successfully used in several regions of the world ocean in altimetric data (Aroucha et al., 2020;Ayouche et al., 2021;Barboni et al., 2021), high 110 frequency radar data (Liu et al., 2020) or numerical simulations (de Marez et al., 2021). From 1 January 2000 to 31 December 2019, AMEDA is applied on AVISO sea surface height (SSH) delayed-time product at a resolution of 1/8°with daily output.
From 1 January 2020 to 31 December 2021, AMEDA is applied on AVISO SSH near-real-time day+6 product (Pujol, 2021), at the same spatial and temporal resolution. In each eddy single observation (one eddy observed one day), AMEDA gives a center and two contours. The 'maximal speed' contour is the enclosed streamline with maximal speed (i.e. in the geostrophic 115 approximation, with maximal SSH gradient) ; it is assumed to be the limit of the eddy core region where water parcels are trapped. The 'end' contour is the outermost closed SSH contour surrounding the eddy center and the maximal speed contour ; it is assumed to be the area of the eddy footprint, larger than just its core but still influenced by the eddy shear (Le . AMEDA gathers eddy observations in eddy tracks, allowing to follow the same structure in time and space, sometimes over several months. The eddy tracks collection in the whole Mediterranean sea constitutes the DYNED Atlas database (Stegner 120 et al., 2019), and is available online (for the years 2000 to 2019) at : https://www1.lmd.polytechnique.fr/dyned/. From 2000 to 2021, a total of 7038 (respectively 8890 ) anticyclonic (cyclonic) eddy tracks were retrieved. The asymmetry in eddy numbers is driven by a lifetime difference, anticyclones living noticeably longer, an asymmetry even more marked in the Levantine basin (Barboni et al., 2021).
A climatological database is created collecting in-situ profiles from the Coriolis Ocean Dataset for Reanalysis (CORA).
Delayed-time (CORA-DT, Szekely et al. (2019b)) profiles are recovered from 2000 to 2019 (113486 profiles), and nearreal-time (Copernicus-NRT, Copernicus (2021)) profiles are recovered from 2020 to 2021 included using the history release (22821 profiles). In addition, some profiles prior to 2020 are not yet released in CORA-DT but available in Copernicus-NRT.
This happens in particular when the salinity sensor of an Argo float has abnormal values but the temperature is still correct 130 (by visual inspection and correct quality flag). As the MLD computation can be performed on a temperature profile alone, profiles were also retrieved in NRT mode after careful check, described in Appendix A. This provides an extra 20746 profiles from 2000 to 2019. Spotted duplicates between CORA-DT and Copernicus-NRT are retrieved only from CORA-DT. Complete database then accounts for 157053 profiles in total. for density and 0.2°C for temperature, based on a reference depth at 10 m to avoid diurnal heating at the surface. In the Mediterranean sea d 'Ortenzio et al. (2005) used this methodology for a 0.5°resolution MLD climatology. Houpert et al. (2015) updated it with 8 supplementary years of data, but opted for a 0.1°C temperature threshold. This more restrictive criterion enables to 140 reduce the difference between the MLD computed on temperature profiles with the one on density profiles. Gradient methods are looking in a similar way for critical gradients as an indicator of the mixed layer base. Typical gradient threshold values in use are 2.5 × 10 −2 • Cm −1 for temperature profiles and range from 5 × 10 −4 to 5 × 10 −2 kg.m −4 for potential density profiles (Dong et al., 2008). Mixed gradient and threshold methods were also developed (Holte and Talley, 2009). Here we aim to capture as accurately as possible the MLD evolution, which can vary on timescales shorter than a month. More specifically we 145 observed in several cases that the threshold method (with criteria ∆σ = 0.03 kg.m −3 and ∆T = 0.1°C) can miss the mixed layer and return the main thermocline instead (see Fig.1a). The main thermocline is indeed characterized by a small jump in potential density (or in temperature) but a significant peak in the gradient profile, and happens mostly in the spring, probably due to a start of restratification quickly mixed. To capture such small-scale restratification events, we built the following methodology combining both threshold and gradient approaches. Using the following thresholds : ∆σ = 0.03 kg.m −3 and ∆T = 0.1 150°C , we derive a first estimate of the MLD. If it is shallower than 20 m, we take it as our estimate of the MLD. Otherwise, we apply a three-point running average to remove small spikes and compute the gradient using a second-order centered difference.
From the subsurface (20 m) up to the first MLD estimate, we apply a gradient method with the given gradient thresholds : If the gradient fails to exceed the threshold within the given depth range, then the first MLD estimate is kept.

155
Threshold and gradient methods are limited by their dependence on the criterion values which can have strong influence on the MLD estimate. The relatively low gradient thresholds chosen here appeared to be necessary to catch the MLD in some of our profiles as higher thresholds would return the main thermocline (see Fig.1a). A sample of 400 randomly-picked profiles colocated inside eddies was used for validation. We chose this validation dataset with profiles inside eddies as it is our main 160 focus. Wrong detection on double gradient profiles inside eddies were found to be quite large, exceeding 100m sometimes. On these 400 profiles, 22 (5.5 %) of them were identified as double gradient profiles, resulting in an overestimated MLD when derived with the threshold method. Moving to our methodology, this issue is now only encountered for 2 profiles (0.5 %). However, with the gradient method comes some issues on profiles with small residual spikes despite the applied smoothing.
For 2 profiles, the gradient method returned wrong MLD detection where the threshold method was correct. However, the 165 gradient method was found overall to be more accurate on estimating the MLD.
Moreover the chosen thresholds should return similar estimates between a MLD obtained from temperature profile (MLD T ) and one from a potential density profile (MLD σ ). Potential density is a better estimate of the stability of a layer, and thus MLD σ should give a more reliable value. However, salinity (and hence density) suffers from data holes, representing about 170 15 % in our dataset. Temperature profiles then offer a good alternative in evaluating the MLD providing MLD T gives similar a estimate than MLD σ . MLD T and MLD σ difference histogram is shown in Fig.1b: the gradient method appears to slightly reduce this difference with 64 % of the profiles leading to the same estimate and 94 % to less than a 30 m difference compare to 62 and 93 % respectively for the threshold method. MLD is then computed on the density profile, and on the temperature if no density is available. Figure 2. Profiles colocalisation with eddy contours for an 'inside-anticyclone' profile (panels (a) to (e)), an 'ambiguous' profile (panels (f) to (j)) and an 'outside-eddy' profile (panels (k) to (o)). Profile cast position is assumed fixed and compared to eddy contours at D-2, D-1, D, D+1, D+2, D being the profile cast date.

Eddy colocalisation and background estimate
In order to characterize the impact of anticyclonic eddies on the MLD seasonal evolution and spatial gradient, we need to accurately colocalize in situ profiles with eddy observations. However due to altimetric product interpolation and disparate satellite tracks, SSH-based contours can vary a lot in size and position, making a single eddy observation less reliable in the Mediterranean sea (Amores et al., 2019;Stegner et al., 2021). Therefore we colocalize eddy observations and in situ profiles 180 at ±2 days. Assuming a profile position fixed at cast date D, it is then labeled as 'inside-eddy' if it remains inside the maximal speed contours of the same eddy at D-2, D-1, D, D+1 and D+2 (at least 4 contours out of 5). This 4-out-of-5 threshold avoids to neglect a colocated profile when the eddy contour is not available for just one day (see Fig.2b). For the same purpose, hereafter the eddy center and the distance of a profile to the eddy center are averaged at ±2 days.

185
AMEDA also gives for each observation the last closed SSH contour (see section 2.1), inside which there is still an impact of the eddy shear, but outside of the maximal speed contour the water particles are not assumed to be trapped. The area between the maximal speed and last closed SSH contours are then considered as an intermediate zone to be discarded. Consistently with the 'inside-eddy' definition, we label as 'outside-eddy' only profiles staying outside any eddy contours at ±2 days of its cast date. Any profile being neither 'inside-' or 'outside-eddy' is considered as ambiguous and discarded. From 2000 to 2021, out 190 of 157053 profiles retrieved in the Mediterranean sea, 104787 are labelled 'outside-eddy', 7939 are 'inside-anticyclone', 14919 'inside-cyclone', and the remaining 29410 'ambiguous' profiles are removed from this analysis. This asymmetry between anticyclones and cyclones sampling is also due to heterogeneous oceanographic surveys (Houpert et al., 2015), in particular the numerous glider missions in the Gulf of Lion, a cyclonic gyre with no large anticyclones (Millot and Taupier-Letage, 2005). 'ambiguous' one  and an 'outside-eddy' one . For this particular 'inside-anticyclone' profile, the maximal speed contour was missing at day D-1, but available the other days, and the profile was indeed cast close to the eddy center.
To follow the accurate evolution of the MLD inside anticyclones, we need a reference for comparison: an unperturbed, local and time-coincident ocean state without eddies, hereafter called 'background'. This outside-eddy background differs from a We therefore chose ∆day = 10 days and ∆y = 1 year as day and year intervals in order to capture MLD variation as short as 215 possible, which is crucial for parameters varying quickly such as the MLD. For the two earliest recorded events (Mersa-Matruh 1 and 2 in 2006 and in 2008, see Table1), ∆y is set to 2 years because no background MLD was available otherwise. Choosing ∆y = 1 allows to have accurate eddy induced anomalies without being corrupted by interannual variability of temperature and salinity fields, which can be marked in the Mediterranean, in particular in the Eastern basin (Ozer et al., 2017) and with a significant warming trend (Parras-Berrocal et al., 2020).

MLD evolution function fit
To describe more objectively the MLD seasonal evolution in the background, we performed a function fit using the Python optimization routine scipy.optimize.curve_fit. MLD data points are selected per 5 days time steps. Background MLD is fitted by a skewed Gaussian, t back max being the time when deepest MLD (M LD back max ) is reached, σ and τ respectively the deepening and restratification timescales : This fit captures the background MLD evolution, somehow smooth, typically with a sharper restratification than deepening (τ < σ). However this is not sufficient for the anticyclonic core MLD evolution that can have more abrupt variations, then calling for a more complex fit with two deepening timescales σ 1 and σ 2 : To fit accurately the MLD evolution, and in particular the maximal depth reached, data are fitted with weights proportional to their depth. Because it is difficult to have long and continuous time series, data are often missing on the previous or next summers. To ensure physical behavior, fit is forced back to 10m on the edges, miming summer stratification. The MLD anomaly (M LD anom ) is defined as the difference between the fitted background and anticyclonic core MLD. M LD anom is a function of time but reaches its maximum (∆M LD) at almost the same time than the absolute anticyclonic core MLD, as the latter 235 has more amplitude than the background one. At last an advantage of the scipy.optimize.curve_fit routine is to provide the parameter covariance matrix, and hence an error estimate taking the square root of the covariance matrix diagonal (Bevington et al., 1993). It can happen for the covariance matrix to have very large values, in this case we used an upper uncertainty of ±30m for ∆M LD and ±20 days for t AE max . A fit illustration is provided in Fig.3 for the Pelops 2 event in 2016 (see Table1), with dots as the real MLD data and fits in continuous lines, background in black and anticyclonic core in red.

240
Using the fit routine, maximal MLD anomaly is then estimated for this event to ∆M LD = 127 ± 13 m. On can also notice an absence of coincidence between deepest inside-eddy and background MLDs. Following previous notation, we can then define a restratification delay of the anticyclonic core MLD, used throughout this study : ∆τ = t AE max − t back max . In the example shown in Fig.3, ∆τ = 26 ± 11 days.

245
Several long-lived anticyclones are tracked for several months, recording up to 16 winter mixed layer deepening events at their core. In order to investigate the relation between the MLD evolution and the vertical eddy structure, we plot together the timeseries of MLD and vertical temperature gradients inside the eddy core. Two different MLD temporal patterns are observed, depending on whether or not the current winter mixed layer reaches the subsurface anticyclone core. This core is constituted by a preexisting homogeneous layer, and in the following, we define as 'homogenized' a layer with temperature 250 gradient constantly below 2.5 × 10 −3°C .m −1 in absolute value.

Winter deepening connecting preexisting subsurface core
Very deep mixed layer can be observed in several anticyclones when the MLD erodes the inside-eddy stratification and abruptly connects with a subsurface homogenized core, an event hereafter called a 'connecting' MLD. An example of this evolution is described below with a long-lived Eratosthenes anticyclone during winter 2008-2009. Its temporal evolution from August 2007   255 to April 2009 is shown in Fig.4, and listed hereafter and in Table1 with name 'Eratosthenes (ERA) 1'. This kind of anticyclone, also called 'Cyprus eddy' or even 'Shikmonah gyre', are large mesoscale structures with almost stationary position south of Cyprus island in the Levantine basin, extensively studied with several CTDs (Brenner, 1993;Krom et al., 1992;Moutin and Prieur, 2012), gliders and Argo floats deployments (Hayes et al., 2011). The anticyclonic density anomaly is characterized on average by a deep (about 400m) and extremely warm temperature anomaly (up to +2°C at 400m) (Moutin and Prieur, 2012; 260 Barboni et al., 2021), sometimes with a strong salt anomaly (Hayes et al., 2011). Thus temperature profiles are considered as a (respectively (f)) shows profiles corresponding position on a map with same color code, together with the eddy maximal speed contour (dark green green shape) and eddy footprint (outermost closed SSH contour, light green shape). Bathymetric data are from ETOPO1 (Smith and Sandwell, 1997). good estimate for relative density and temperature gradient for stratification.  (Fig.4a). This restratification occurring at a different time outside-and inside-eddy -with a delay of about two months -leads to the noticeable situation measured on 1 March 2009 : the insideanticyclone profile is warmer than its environment at depth (100 to 350m deep) but homogenized in its upper part, whereas background profiles are stratified with positive temperature gradient. Such geometrical configuration leads to an anticyclone negative temperature (and hence positive density) anomaly from 50m to the surface, compared to the stratified outside-eddy 285 profile. Such positive density anomaly above the eddy core is then a clear signature of subsurface anticyclone (Assassi et al., 2016).
Over the whole 2008-2009 winter, background MLD barely reached 60m whereas the anticyclonic core MLD went down to 350m. This intense deepening at the anticyclone core is due to the preexisting subsurface eddy, made of a well mixed layer of 290 few hundred meter depth below the summer stratification. When the winter mixed layer deepens, it reconnects to this deep subsurface core and leads to a rapid and strong MLD increase in comparison to the eddy background. This MLD temporal pattern is characteristic of a 'connecting' event, observed 10 times in our analysis and throughout the Mediterranean sea (see section 4.3).

Winter deepening not connecting preexisting subsurface core 295
Conversely for some other anticyclones, it can clearly be seen that the subsurface temperature anomaly does not connect with the winter mixed layer and remains unperturbed and homogeneized at depth. Such event is hereafter called a 'non-connecting' MLD. Figure 5 shows the evolution of another Eratosthenes anticyclone living from 2009 to early 2012, with two recorded anticyclonic core MLD deepening in 2010 and 2012 (respectively listed in Table1 as 'ERA2' and 'ERA3'), with same color codes than Fig.4, with profiles at 20 March 2010 and 15 June 2010. As several profiles were located at the same time inside 300 the anticyclone, they are shown in orange line on Fig.5c and 5e (respectively orange dot for MLD on Fig.5a). The red line highlights only the profile with closest distance to the eddy center, assumed to be more representative of the eddy core.
Similarly to the 'ERA1' event in 2009 described above, a thick and deep subsurface anomaly forms primitive eddy core in late 2009, as an homogeneous layer from 250 to 400m deep ( green bar on Fig. 5c), reaching an anomaly about +2.5°C at 305 400m. However the anticyclonic core MLD did not deepened in the winter 2009-2010 below 150m, only forming a second homogeneous layer above. This constitutes a second surface core, still separated from the primitive core by a temperature stratification, revealed by a temperature gradient continuous in time (Fig.5b). On the vertical profile on 20 March 2010 (Fig.5c) about 1°C temperature jump remains between the two core, forming a double-core anticyclone. In June 2010 (Fig.5e), this second homogeneous layer is itself covered by the spring restratification, then forming what is also referred to as 'thermostad ', 310 or 'mode-water eddy' in the literature (Dugan et al., 1982). Thanks to the trapped Argo floats remaining near the eddy core for months, both cores can be tracked until August 2010 as separated in subsurface.
Such 'non-connecting' winter MLD inside anticyclone reveals the possibility of a persisting separation between a primitive subsurface anticyclone core and the new homogeneous layer formed by the current winter mixing, then constituting a double-315 core anticyclone. The example showed in Fig.5 occurred on a Eratosthenes anticyclone, in the same region than anticyclone in which a 'connecting' MLD example was previously shwon in Fig.4. This MLD pattern is not limited to this example and was also observed 5 times in other regions, as detailed in the next section. Consequences of the formation of a double-core anticyclone are discussed in Sect.5.2 hereafter, with another remarkable example in Fig.10.

Inside-anticyclone MLD statistics
From 2000 to 2021, thanks to extensive Argo deployments sampling eddies, 16 winter MLD deepening events were accurately recorded with vertical profiles in 13 mesoscale eddies, 10 being 'connecting' events, 6 'non-connecting' ones. Several structures were indeed surveyed over 2 winters (see Fig.5 and 10). For each event, the fitting method detailed in Sect.3.3 is applied and parameters are reported in Table1 together with eddy characteristics : eddy SSH amplitude, maximal speed V max and 325 maximal speed radius R max . Eddy measurements are estimated by the mean from November to March of the corresponding winter. Figure 6 shows the location of each structure, which actually corresponds to a type of long-lived structures already identified in the literature (Millot and Taupier-Letage, 2005;Hamad et al., 2006;Budillon et al., 2009;Barboni et al., 2021),  10c) or a new one (see Fig.5e and 10d). In all winter deepening events listed on Table 1, such homogeneous layers of least 50m were indeed visible on vertical profiles. For Tyrrhenian sea anticyclones (TYR1 & 2) with stronger salinity influence (Budillon et al., 2009), homogeneous layers with density gradient below 5.0 × 10 −4 kg.m −4 are also visible. One can also notice that 340 'non-connecting' events are quite common, but double-core structures should be even more frequent. Indeed a 'connecting' event can occur inside a double-core structure and reconnect only the homogeneous core formed in the previous winter but not the deepest anomaly, as shown later in Fig.10b-e. In other worlds the proportion of 'non-connecting' events in Table 1 and Fig.6 should be considered as a lower bound for double-core structures, revealing their high occurrence. and its SSH amplitude, using regional average and monthly climatology. We previously showed that MLD anomaly varying over very short timescales can produce sharp MLD gradient and anomalies reaching several hundreds meters, not captured by smoothed composites. The relation between MLD anomaly and eddy amplitude is tested on Fig.7a  anomaly for 1cm eddy amplitude). This proposed relation is obviously not verified, the deep MLD observed in Mediterranean anticyclones exceeding by far the relation. On the opposite, deepest MLD anomalies seem to be observed in the eddies with weakest SSH signature. Although surprising at first sight, this trend might be explained by the fact that deepest MLD can be observed when mixed layer abruptly connects to an anticyclone deep homogeneous core, in subsurface and hidden by a strong seasonal thermocline. This is in particular the case for 'ERA1' shown in Fig.4, having an extreme ∆M LD deeper than 300m 355 but the lowest SSH signature in Table 1. The relation between the MLD anomaly and the eddy Rossby and Burger numbers are also tested in Fig.7b  MLD anomalies appear for various eddy intensity and size and for both 'connecting' and 'non-connecting' events. One can only notice that 'connecting' events pull MLD deeper in general, and that these events are slightly more observed in large structures (small Bu). Remote-sensing measurements are then hard to link with observed eddy-induced MLD anomalies. On the opposite the diversity of vertical structure shown in this study (Fig.4, 5 and 10) suggests that eddy vertical structure might 365 have more influence, and previously proposed linear relation seem to apply mostly for surface-intensified structures. A new and important observation is that MLD inside anticyclones tends on average to clearly restratify later than the neighbouring background. It was shown for two individual events in Fig.4 (ERA1, 'connecting') and Fig.5 (ERA2,'non-connecting'), 370 but it is statistically robust in Table1 : average t back max is 22 days and average t AE max is 49 days, meaning restratification usually begins in the second half of February in anticyclones, on average one month later than outside-eddy. Restratification delay ∆τ can reach two months in some case : 67 days for ERA2 or 74 days for MM6 (see Fig.9b). Figure 8a shows the relation between ∆τ and ∆M LD and reveals that no clear trend can be identified alone : deep MLD anomalies are observed when the anticyclone MLD restratify early (low ∆τ ) or later (large ∆τ ). However when distinguishing 'connecting' and 'non-connecting' 375 events, a linear trend appears separately : MLD anomalies go deeper as ∆τ increases and for similar ∆τ value 'connecting' events go deeper. Linear fit is performed separately for both types, shown by dashed line on Fig.8 : for each day of continued MLD deepening inside anticyclones, a 'connecting' (respectively 'non-connecting') MLD gets about 3m deeper with correlation coefficient of 0.766 (1m deeper with correlation 0.604). This trend is logical, as a later restratification (large ∆τ ) lets the MLD deepens longer and hence leads to larger ∆M LD.

380
In order to analyze the MLD evolution together with the mixed layer cooling, Fig.9 shows in the upper panel the MLD fit (see  (2021), the anticyclone core is indeed about +1°C warmer than its background, and continues to cool for a while. The inside-eddy maximal MLD is reached around 1 March 2016 (20 March 2021) or +60 (+80) days on Fig.9c (d). Few weeks after, both background and anticyclonic core mixed layers start to warm again around 1 April (in both PEL2 and MM6), but background MLD started to restratify 1 month earlier in PEL2 (2 months earlier in MM6).

390
Although it is hard to infer a mechanism from few observations, it seems that the beginning of restratification outside-eddy does not mean the mixed layer is warmed up again and outside-eddy mixed layer remains cold. The restratification delay seems to be the consequence of a maintained cooling of the initially warmer anticyclone core. Summer heating seems on the other hand to begin at the same time inside-and outside-eddy. Possible mechanism driving this sustained mixing at the anticyclone core are discussed later in section 5.3. An important observation is also that temperature difference between anticyclone core 395 and the background is on the order of +1°C while MLD deepens, but almost vanishes (or even get slightly negative) when the mixed layer warms again. Although sparse, these in situ observation are in total agreement with observed eddy SSTA switch by Moschos et al. (2022) from winter warm-core anticyclones to predominant cold-core anticyclones with spring restratification in the Mediterranean sea.
5 Discussion on physical drivers and perspectives 400

MLD anomaly scaling
We clearly identified the distinction between 'connecting' and 'non-connecting' events as a more important driver than other eddy parameters such as eddy amplitude, surface intensity or size (see Fig.7-8), and this might explain the difficulty to find a general law for any eddy-induced MLD anomaly. Indeed 'connecting' deepening mixed layers seem limited by the bottom of the preexisting subsurface homogeneous core to which they connect (example in Fig.4e), whereas 'non-connecting' ones 405 by definition do not go deep enough and are then expected to be limited by the heat loss, likely also influenced by the eddy.
The other important parameter is the restratification delay (∆τ ), measuring how long the anticyclone continues to deepen MLD or not, and which eventually scale with the maximal MLD anomaly. In a remote sensing perspective, both parameters seem very hard to assess without in situ profiles inside the eddies. Examples shown in this study ( Fig.4, 5 and 10 below) showed the complexity induced by possible connection with previous subsurface anomaly and more generally the key role 410 of the anticyclone vertical structure, that was totally smoothed in previous composite studies. Relationship between eddyinduced MLD anomaly and satellite measurements are definitely more complex. However as theorized by Assassi et al. (2016), detecting remotely information on the eddy vertical structure could be possible, in particular distinguishing the subsurface or surface-intensified nature by comparing eddy signature in SSH and SST. For instance ERA1 event in Fig.4 is an almost textbook case of a subsurface anticyclone with isopycnals doming leading to a cold eddy-SSTA.

Double-core eddy formation
High occurrence of 'non-connecting' events (crosses on Fig.6) are very interesting as they show the formation of double-core anticyclones through winter deepening of the surface layer above a preexisting density anomaly. Double-core eddies were often surveyed in the world ocean (Lilly et al., 2003;Belkin et al., 2020), including in the Western Mediterranean sea (Garreau et al., 2018). Despite various propositions (see e.g. Belkin et al. (2020) for a list ), no clear formation mechanisms emerged. Several 420 studies focused on the so-called 'vertical alignment' of two eddies with different densities in experimental works (Nof and Dewar, 1994), observations (Lilly et al., 2003) or modelling (Trodahl et al., 2020). Interestingly Lilly et al. (2003) observed well in the Labrador sea that double-core anticyclones mostly consist of convective lenses formed in different winters, the heat flux interannual variability leading to different density anomalies, but they explained it with eddies formed separately which later aligned. There were nonetheless previous observations of a second lighter core generated above a preexisting anticyclone.

425
Thanks to repeated XBT transects, Nilsson and Cresswell (1980) surveyed such phenomenon in an anticyclone detached from the East Australian Current, caused by winter heat loss. Bosse et al. (2019) surveyed this in the Lofoten eddy with winter convection, but through glider sections spaced in time, then with a temporal resolution on the order of the month. More recently Meunier et al. (2018) explained the formation of a double-core Loop Current Eddy by winter diabatic processes. However this case is different from the Mediterranean anticyclones ; as the Loop Current Eddy consists of an advection of a large structure 430 of Caribbean waters into the Gulf of Mexico experiencing different surface fluxes with more heat loss and precipitation than the area where they originate. These diabatic processes by surface winter mixing result in a fresher shallower core above a saline core of subtropical under waters. Moreover Meunier et al. (2018) explained quantitatively the observed anomaly against regional average of atmospheric fluxes, whereas in our study the differential MLD evolution between the eddy core and the background ( Fig.4a and Fig.5a) suggested fluxes variations at the scale of the eddy.

435
What drives the formation of double-core structures should be further investigated, but one could expect the interannual variability of heat fluxes to be the main driver. This was already suspected by Lilly et al. (2003) (although for them it was for separate eddies) and Moutin and Prieur (2012). A winter with strong heat loss is expected to deepen MLD a lot, including inside-eddy, and a subsequently warmer winter could then not achieve to deepen the MLD as much. This mechanism drives 440 mode-water formation, and it was already shown in other regions, mostly the Atlantic Ocean, that eddies could modulate mode-water formation (Dugan et al., 1982;Chen et al., 2022). Such hypothesis could also explain the high occurrence of 'non-connecting' events in the Mediterranean sea, this region being known for a high interannual variability of winter heat loss. Pettenuzzo et al. (2010) indeed found maximal winter heat loss to vary by 20% to 30% (in regional monthly average), and a plausible connection with the Northern Atlantic Oscillation (NAO). This interannual variability of the heat flux was 445 already shown to influence deep convection in the North-Western Mediterranean sea (L'Hévéder et al., 2013), and it could then be expected a higher occurrence of double-core anticyclones due to a stronger Mediterranean sea stratification in a warming climate (Somot et al., 2006).
An important consequence in the formation of this lighter core in a 'non-connecting' winter deepening is that the second 450 core is separated from the surface by a thinner seasonal stratification. The next winter is then likely to connect again the new mixed layer with the upper core, while possibly keeping the primitive deeper core untouched. Such kind of interaction from one winter to another was observed in the Ierapetra eddy and is presented in Fig.10 (with same color code as in Fig.4a-f and 5a-f). The Ierapetra eddy is a recurrent long-lived and intense anticyclone formed South-East of Crete (Theocharis et al., 1993;Lascaratos and Tsantilas;Ioannou et al., 2017) and recently surveyed by the PERLE 1 and 2 campaigns (Ioannou et al.,455 2019; Durrieu de Madron and Conan, 2019). Similarly to Eratosthenes anticyclones previously shown, the density anomaly is mostly driven by a warm core, allowing to use temperature profile as a proxy for stratification (Ioannou et al., 2017). Figure 10 shows the Ierapetra anticyclone formed in autumn 2016. The first winter 2016-2017 turned out to be a 'non-connecting' event ('IER1' in Table1). Indeed in March 2017 a preexisting subsurface homogenized layer remains between 350 and 450m, below the maximal anticyclonic core MLD of 220m (Fig.10c) and with about +2°C temperature anomaly. From April to Decem-460 ber 2017, summer heating restratified the upper layer and let below a second homogeneous layer between ∼ 100 and 200m deep. The primitive core remained homogeneized at depth and separated by a temperature gradient throughout the summer (Fig.10b). In January 2018 (Fig.10d), inside-eddy vertical profile shows the mixed layer deepening at 120m -already deeper than the background MLD -and the double-core structure is still retrieved.At last at the end of February 2018 (Fig.10e), the MLD completely eroded the seasonal stratification and connects the current MLD with previous winter subsurface core, then 465 reaching about 280m. The winter 2017-2018 is then a 'connecting' event ('IER2' in Table1). The timeseries is interrupted inside the anticyclone, but Argo floats are again colocalized in May 2018, and despite some variability a temperature gradient continuously separates the two cores between 200 and 250m (see Fig.10b). IER2 was then a 'connecting' event on a doublecore structure. The primitive anticyclonic core was not mixed but remained homogeneized at depth. These data bring to light a possible formation process of a double-core anticyclone through winter convection, and also documents for the first time the 470 fate of the formed subsurface anomaly, which can be tracked up to the next winter when it was mixed again.

Physical drivers
The observed importance of restratification delay ∆τ should also have underlying physical mechanisms. Prolonged MLD deepening and cooling inside-eddy (see examples in Fig.9) leads to the extreme MLD anomalies sometimes larger than 300m, 475 and hence marked MLD gradient occurring at the scale of the eddy radius or shorter. Indeed Gaube et al. (2019) found anomalies on the order the mesoscale but in a composite vision and for large eddies compared to the deformation radius, MLD gradients in eddies should occur on shorter scales (Meunier et al., 2018). Such marked MLD gradients should trigger mixed layer instabilities leading to restratification (Boccaletti et al., 2007;Fox-Kemper et al., 2008), which calls for mechanisms sustaining the mixing inside-eddy during the restratification delay. It should be noted that a homogenized layer itself does not 480 proof an active mixing, but still reveal the absence of restratification. We also interestingly noticed that in several cases Argo floats remained well in the anticyclone core during the MLD deepening phase, but often leave the eddy soon after, maybe a signature for mixed layer instabilities impacting the eddy. The first mechanism explaining longer mixing in anticyclones would be an eddy modulation of air-sea fluxes by eddy-induced SSTA. Villas Bôas et al. (2015) observed such eddy-modulation on air-sea sensible and latent heat fluxes, but in regions of energetic surface-intensified eddies, with very warm anticyclones (in particular the Algulhas current retroflexion). For subsurface anticyclone, the eddy-induced SSTA is on the opposite expected to be weakened ( see the example of the cold-core anticyclone shown in Fig.4e and the study of Assassi et al. (2016)), and this mechanism might then not be the most important. MLD deepening enhanced in anticyclones could be explained by other eddy retroactions than on the heat fluxes, a possible mechanism being the eddy-induced Ekman pumping (Stern, 1965;Gaube et al., 2015) or enhanced mixing in anticyclones due to near-inertial waves trapping (Kunze, 1985).

Impact on eddy dynamics
Connecting events also raise interesting questions on the consequence of such mixing of deeper subsurface anticyclone core, in particular the role of inside-eddy convection on the eddy dynamics itself. Studies in the literature mostly focused on winter convection inside cyclone because of the preconditioning with isopycnals doming at their center (Legg et al., 1998;Legg and 495 McWilliams, 2001). Such phenomenon should also applied for subsurface anticyclone due to the surface isopycnals doming and subsequent stratification weakening (Assassi et al., 2016). The coincidence of observed multiple 'connecting' winters in long-lived anticyclones like the Mersa-Matruh and Eratosthenes structures suggests a possible mechanism regenerating these structure, and maybe explaining the extremely marked cyclone-anticyclone lifetime asymmetry in the Levantine basin (Mkhinini et al., 2014;Barboni et al., 2021). Interestingly Brenner (1993) already proposed winter cooling as possible mecha-500 nism explaining the sustained lifetime of the anticyclone surveyed south of Cyprus. The other structure calling for comparison is the Lofoten eddy in the Sea of Norway, and the Rockwall Trough eddy offshore Ireland, two long-lived deep anticyclonic structures. Winter convection was observed inside the core of the Lofoten eddy, and once thought to help regenerating the structure (Ivanov and Korablev, 1995;Köhl, 2007;Bosse et al., 2019). Double core formation was also observed in the Lofoten eddy (Bosse et al., 2019). Recent numerical studies showed this regeneration was primarily driven by merging of smaller structures 505 (Köhl, 2007;Trodahl et al., 2020), however de Marez et al. (2021 showed that winter time convection eased this merging process by deepening of the eddy core. Merging of eddies detached from the coast towards an offshore anticyclonic attractor being also observed in the Levantine basin (Barboni et al., 2021), this could provide another explanation to long-lived Mediterranean anticyclones. Cylone-anticyclone asymmetry might not have just one mechanism, as other arguments were already proposed.
Anticyclones have indeed a larger radius and are more coherent.

Biological impacts inside anticyclones
'Connecting' winter mixed layer in the Eratosthenes anticyclone was already observed by Krom et al. (1992) with a biogeochemical focus in 1989 (there called 'Cyprus eddy'). They measured a February inside-anticyclone MLD of 450m (compared to 200m) at the eddy boundary, with later spring restratification in May. Their temperature profile clearly corresponds to a 'connecting' event, with even deeper MLD than in our study. Comparing also nitrates, phosphates and chlorophyll, they showed that 515 chlorophyll production was about 30% more abundant at the eddy core while relatively similar to Levantine basin average at its edge. The main nitracline was consistently measured below the winter mixed layer also at the eddy core, then around 450m.
While spring phytoplankton bloom occurred at the surface, they observed a mixed homogeneous layer remaining aphotic, between the euphotic zone (the upper 120m) and this deep main nitracline (called 'decomposition zone' in Krom et al. (1992)).
They observed instead here from February 1989 to January 1990 an increase of both nitrates and phosphates. Consequently in 520 the eddy a second nitracline formed at the bottom of the euphotic layer, at approximately the same level than a summer deep chlorophyll maximum. Similarly in another Eratosthenes structure in July 2008, Moutin and Prieur (2012) also estimated the maximal mixed layer in the previous winter to have reached 396m and let only a stratified thermocline below, another clear description of a 'connecting' winter mixing. Moutin and Prieur (2012) also observed in the eddy a main nitracline at roughly 400m depth, and a second one around 100m together with a deep chlorophyll maximum. The 100-400m depth zone remained 525 with low mineralization and high value of dissolved organic matter (DOC). From these observations, it seems that a 'connecting' deep MLD induces a strong nutrient input to the euphotic layer, but letting in summer an homogenized aphotic subsurface layer with high DOC export at depth.
None of the two studies mentioned above observed the case of a 'non-connecting' MLD frequently observed in our study 530 and also in an Eratosthenes anticyclone (see Fig.5). However Moutin and Prieur (2012) discussed this possibility : the winter MLD does not reach the main nitracline/phosphacline, it would keeps the upper layer away from the deep nutriment source.
The whole system would then evolve towards an ultra-oligotrophic system because of nutrients being very weakly injected to the euphotic layer. This is expected to be in particular the case when a primitive subsurface core does not connect to the surface for several winter, such as the example of the Ierapetra anticyclone in Fig.10. The high frequency of 'non-connecting' 535 anticyclone MLD observed in our study then suggest that anticyclones in the Eastern Mediterranean sea are ultra-oligotrophic systems more frequently than previously thought. Temporal evolution of such 'non-connecting' events with biogeochemical instruments such as BGC-Argo would be interesting to follow this analysis.

Conclusions
In this paper we were able to analyze, thanks to a combination of satellite observations and numerous in-situ data, several time 540 series that finely describe the evolution of the winter mixed layer in the core of Mediterranean anticyclones. We even succeeded in following, for the same long-lived anticyclone, the evolution of its MLD over two consecutive years. This allowed us to quantify extreme anomalies induced by mesoscale eddies in the mixed layer, which would have been smoothed in a standard composite analysis. Indeed, we observed that the winter mixed layer can go down to 380 m in the core of Levantine Basin anticyclones, while the surrounding background MLD does not go deeper than 80 m or 100m. 545 We also observed a time lag of several weeks, and sometimes up to two months, in the spring restratification between the core of these deep anticyclones and the background sea, revealing that MLD temporal evolution is not uniform. Indeed, when the later restratifies due to rising temperature of the atmosphere, the core of these mesoscale anticyclones which are warmer continues to deepen and to cool. This time lag induces very strong spatial heterogeneities of the MLD in the eastern Mediter-550 ranean Sea during the early spring, with observed maximal MLD ranging from 50 to 330m.
We showed that this localized deepening of the MLD is controlled by the vertical structure of these eddies. When the surface mixing layer connects with the subsurface core of preexisting anticyclones, a rapid deepening of the surface mixed layer is observed. Conversely, when the surface mixed layer does not connect with the subsurface core, a double-core eddy is formed.

555
Connection or not with preexisting subsurface core prove to be more relevant to describe MLD deepening than other eddy parameters such as SSH amplitude or size. MLD anomalies was observed to linearly increase with restratification delay, but increasing roughly 2 to 3 times faster for 'connecting' MLD than 'non-connecting' one.
These extreme MLD deepenings in anticyclone cores reveal complex and rich interaction between surface and subsurface 560 of the eddies. Connection between the mixed layer and subsurface anomalies provide a way to propagate heat at depth while mixing in winter, which consequences that remain to be investigated. These winter deepening inside anticyclones could also play a role in sustaining the extremely long-lived anticyclone in the Eastern Mediterranean. MLD anomalies in cyclonic eddies remain to be investigated, and an open question would be to know if a restratification delay could also be observed in cyclones.  Appendix A: In situ profile checking methodology In both CORA-DT (Szekely et al., 2019b) and Copernicus-NRT (Copernicus, 2021) datasets vertical profiles data coming from XBT, CTD, glider and profiling floats are collected by selecting files with respective data type codes XB, CT, GL, PF. When a profile from 2000 to 2019 was available in both DT and NRT mode, it was retrieved from the CORA-DT dataset which performs more quality checks (Szekely et al., 2019a). Selection was done with the following steps, separately for temperature 765 and salinity, and when available, 'ADJUSTED' properties are collected : -Position and date quality control (QC) flags equal to 1,2 or 5 and position not on land.
-Temperature data below 12°or above 35°are discarded, salinity data below 30 PSU or above 42 PSU are discarded. These 770 parameters are specific to the Mediterranean sea.
-When both temperature and salinity are available, density is computed using the TEOS-10 equation (McDougall et al., 2009) from the Python package gsw (https://teos-10.github.io/GSW-Python/) -Profiles are linearly interpolated on the same vertical grid, with 5m grid step from 5m to 300m depth and 10m grid step from 300m to 2000m. Maximal gap allowed is 20m, and profiles with gaps are discarded.

775
-Profiles with temperature jumps higher than +6°C or -2°C (positive upwards) between two grid points are discarded, as assumed to be unrealistic. This is required in particular to filter out noisy XBT profiles.
-After these steps, only profiles with more than 40 data on the interpolated grid between 50 and 400m are kept. and (b) ∆day (day interval). (c) Sensitivity to the MLD computation method. The background MLD method used throughout this study is ∆day = 10 days, ∆y = 1 year and the gradient method (common navy blue line on panels a-c).