Submesoscale modulation of deep water formation in the Labrador Sea

Submesoscale structures fill the ocean surface, and recent numerical simulations and indirect observations suggest that they may extend to the ocean interior. It remains unclear, however, how far-reaching their impact may be—in both space and time, from weather to climate scales. Here transport pathways and the ultimate fate of the Irminger Current water from the continental slope to Labrador Sea interior are investigated through regional ocean simulations. Submesoscale processes modulate this transport and in turn the stratification of the Labrador Sea interior, by controlling the characteristics of the coherent vortices formed along West Greenland. Submesoscale circulations modify and control the Labrador Sea contribution to the global meridional overturning, with a linear relationship between time-averaged near surface vorticity and/or frontogenetic tendency along the west coast of Greenland, and volume of convected water. This research puts into contest the lesser role of the Labrador Sea in the overall control of the state of the MOC argued through the analysis of recent OSNAP (Overturning in the Subpolar North Atlantic Program) data with respect to estimates from climate models. It also confirms that submesoscale turbulence scales-up to climate relevance, pointing to the urgency of including its advective contribution in Earth systems models.

www.nature.com/scientificreports/ IC contributes to restratifying the interior of the basin 17 and prevents ice formation in the central portion of the basin in winter. Convective activity to the south of the Greenland tip is controlled by momentum fluxes and wind forcing, while in the central portion of the basin the LSW formation is driven by local surface buoyancy loss and modulated by the atmospheric heat fluxes 18 and by ocean dynamics 19,20 . Long-lived coherent mesoscale eddies populate the basin and modify the heat and salt budgets of the gyre interior. The largest are the Irminger Rings (IR) with a diameter between 30 and 60 km 19 . They form through localized baroclinic instability 21 near a constriction of the isobaths along West Greenland, north of Eirik Ridge 22 , are predominantly anticyclonic and transport off-shore the warm and salty water from the IC. As they approach the convective region, they release the heat in their cores and compensate the heat loss from the surface during wintertime convection 19,23,24 .
Several other factors contribute to the LS hydrography and its variability, including freshwater inputs from Arctic and Greenland Ice Sheet (GrIs) melting, continental runoff and precipitation. Over the recent decades, GrIS mass losses accelerated, especially along west Greenland, resulting in increased freshwater inputs into the adjacent seas. Future projections indicate as highly probable a further exponential acceleration 25,26 . Freshwater fluxes contribute to the upper stratification in the LS 27 , also influencing the marine ecosystem primary productivity [28][29][30][31] . Recently, global climate models suggested that these freshwater anomalies may reduce dramatically the LS convection, weakening the AMOC 32,33 . Regionally focused experiments at 2.5 km horizontal resolution, on the other hand, have shown that most of the surface meltwater runoff from southwest Greenland is not transported offshore the continental shelf 34 , with strength and direction of the winds in August and September determining the offshore transport of GrIs meltwater [34][35][36] . Little is known, however, of the contribution that SMCs may have in offshore advection of heat and freshwater anomalies.
In this work we probe how the LSW formation is represented in a regional model, at mesoscale permitting (15 km, LBR15), mesoscale resolving (5 km, LBR5) and submesoscale permitting (1.7 km and 1 km, LBR1.7 and LBR1, respectively) horizontal resolution, focusing on the contributions of submesoscale circulations. The integrations are performed with and without the GrIs meltwater input.

Methods
We use the Regional Ocean Modeling System (ROMS) 37,38 in its Coastal and Regional Ocean COmmunity model (CROCO) version 39,40 . ROMS is a free-surface, terrain-following, hydrostatic primitive equation model, configured here loosely following 34 .
Convection in the LS takes place within vertical plumes with radius of order O(1) km and vertically extent up to 2 km 5 that can be simulated directly only using non-hydrostatic models 41 . In ROMS the non-local K-Profile Parametrization (KPP) scheme 42 accounts for the unresolved vertical mixing and realizes the vertical turbulent fluxes as a summation of down-gradient fluxes and non-local contributions. Previous studies have shown that ROMS can properly simulate extent, seasonality and variability of convective episodes whenever the mesoscale circulation of the basin is resolved 18,34,43 .
The horizontal resolution in the study increases from 15 km (LBR15) to 5 km (LBR5), 1.7 km (LBR1.7), and finally 1 km (LBR1), with 30 vertical levels in all cases. We followed the offline-nesting procedure described in 44 as we moved towards decreasing the grid-point size, starting from the mesoscale permitting case (Fig. 1). The topography is derived from ETOPO2 45 . To avoid potential errors associated with the pressure gradient and the s-coordinate horizontal layers 46 but capture as much bathymetric detail as possible, the topography is smoothed using a logarithmic interpolation method 47 to a maximum slope parameter r max = 0.35. The parameter r max is defined as the ratio (r max = Δh/h mean ) of the maximum difference between adjacent grid cell depths and the mean depth at that point. The bathymetric details are important given the topographic control of the boundary current meandering and eddy formation along the West coast of Greenland, around Cape Desolation 21,22,48 .
The model domain in the LBR15 case covers approximately 48° N-66.5° N and 34° W-65° W. All boundaries are open and, at the boundaries, the velocity fields, along with temperature and salinity profiles, are nudged to the Simple Ocean Data Assimilation (SODA) reanalysis version 3.4.2 49 . The circulation in the basin is sensitive to the width, strength and variability of incoming currents, and the choice of SODA is supported by a previous comparison with other reanalysis products 43 . In all runs, the model is forced by daily surface winds stresses and heat fluxes from the ERA-interim product 50 and the surface heat fluxes are corrected towards the NOAA Extended Reconstructed Sea Surface Temperature (ERSST), Version 4 51 , available at (https ://psl.noaa.gov/thred ds/catal og/Datas ets/noaa.ersst .v4/catal og.html). As we increase horizontal resolution, the size of the nested domains decrease and boundary conditions are extracted from the run at the immediately lower resolution, so that LBR1 covers only the Labrador Sea proper (52.5° N-65.5° N and 44° W-64° W) (Fig. 1B). We choose to impose high frequency boundary conditions (3-day averages) to avoid deflections along the southern boundary.
Stratification and horizontal density gradients along the west and east Greenland coasts (west of 47.6° W, east of 44.2° W and south of 66° N) are strongly affected by the seasonal meltwater inflow from the Greenland fjords which is greatest during summer, between July and September. We focus in the period September 2007-August 2013, during which very different momentum and meltwater fluxes conditions were observed. In a first set of runs, the meltwater is included and the input dataset used is that in 34 . The runoff at the 108 major fjords along the whole coast of Greenland is estimated by the Modele Atmospheric Regional (MAR) coupled to the 1-D Surface Vegetation Atmosphere Transfer scheme Soil Ice Snow Vegetation Transfer (SISVAT) 52 . The amount of meltwater introduced from GrIs (Fig. 1C) does not vary among the experiments but the number of source points 'fjords' may differ due to grid size limitations. For example, in the lower resolution case the meltwater is introduced at 80 grid points by summing up the discharge amount of two or three fjords within one CROCO grid cell into one source point. A passive tracer is also released at each source point with constant concentration and in each experiment, except for the highest resolution run. At all 'fjords' , the meltwater and the corresponding passive  Figure S1). Differences in the mean circulation among experiments can be summarized as follows: the intensity of the WGC, EGC and LC currents increases for increasing resolution; while intensity increases, WGC narrows and in the 1.7 km and 1 km resolution runs the intensification is confined to the core of the current that sits offshore the shelf and away from the fjords. The broader WGC flows westward at 60° N following a steady path along the 3000 m isobath in the mesoscale permitting resolution, while the westward veering begins at lower latitudes, south of Cape Desolation, in the LBR5, LBR1.7 and LBR1 cases. In Fig. 2A the variance of the sea surface height averaged over the integration period is presented as a proxy of eddy activity. The maximum variability occurs off the west coast of Greenland, around 61° N and 52° W, and a secondary maximum in the central Labrador Sea is found at about 58° N and 52° W, as in the altimetric data 55 . The strength of the signal in both simulated areas of maximum variability increases with resolution. Surface vorticity differs significantly among simulations, especially in winter, when eddies are more abundant 56 (Fig. 2B). CROCO captures the formation of the IRs 23,48,57 only at the three finer resolutions, 5 km, 1.7 km and 1 km; more intense mesoscale eddies, submesoscale coherent vortices (or SCVs) and vorticity filaments 1 fill the basin in the LBR1.7 and LBR1 simulations, as quantified by the time series of the domain-averaged vorticity (Fig. 2C). SCVs form abundantly where the continental slope meets the shelf which corresponds, on the Greenland side, at the meeting point between the IC and WGC. The time-series of absolute relative vorticity display small and comparable interannual variability in all cases and a seasonal cycle that becomes more pronounced the higher is the resolution, with elevated values from mid-November to May and a less active period during summer and early fall. In the submesoscale permitting cases, the winter amplification follows that found in open ocean studies [58][59][60] . Noticeably, the increase in vorticity towards the end of summer happens not only faster but also slightly earlier in LBR1.7 and LRB1 compared to the mesoscale resolving case (LBR5).
Resolution dependency on deep water formation. Next, we explore the horizontal resolution dependency of the representation of deep convection in the LS, spanning different dynamical ranges, from mesoscale to submesoscale permitting. We quantify the impact of horizontal resolution in the representation of LSW formation computing the mixed-layer depth (MLD) over the basin, defined as the depth at which density differs from the surface by 0.008 kg m −3 , as in 18 . This threshold follows from the observation in 61 that density differ-  Fig. 3A for all simulations. The modelled convective patches are centered at about [57° W-57° N], consistent with observations 6,63 , but their extent differs significantly among runs. Convective area and convective volume decrease as resolution increases. As previous regional simulations have shown 18 , and in agreement with observational and modelling studies 18,64-67 , temporal characteristics of convective activity as initiation, intensity, seasonality and duration are modulated foremost by the strength of the atmospheric heat fluxes, and secondly by the characteristics of the IC that is advected to the interior of the LS basin. From a numerical point of view, the representation of the boundary current system and its instabilities is important for reproducing correctly such stratification. The IRs, in particular, contribute through eddy-induced lateral fluxes and lateral mixing 19 ; however their effectiveness depends greatly on model resolution. In the simulation that barely resolves the Rossby deformation radius of the basin (LBR15), the absence of vigorous mesoscale activity causes a weak stratification of the upper water column in the central portion of the basin, as to be expected, resulting in stronger mixing and a deeper mixed layer patch that extends further north compared to the higher resolution cases (Fig. 3). Resolving the mesoscale circulations (LBR5) induces a nearly 50% reduction of the convective volume produced by LBR15. This is shown in Fig. 3B by the time series of the difference in convective volume between the lowest resolution experiment in the absence of meltwater input (LBR15-NoMW) and all other runs. The meltwater from the GrIS further decreases convection www.nature.com/scientificreports/ by an additional 2.7% (LBR5-MW). As Fig. 3C suggests, in LBR5, this additional reduction takes place in the northeast corner of the convective patch case (thick versus thin blue lines). The contribution of submesoscale advection amounts to another 30% decrease. Differences between LBR1.7 and LBR1, with and without GrIS input, remain linear. This information is further quantified in Fig. 3D, where the convective volume averaged over the convective season (Jan-May) is plotted against the mean value of near surface vorticity averaged over the LBR1 domain. A linear fit describes well the dependence of the time averaged quantities. It is enlightening to visualize the link between meso-and submeso-scale circulations and MLD also through time-snapshots of MLD, shown for various days in February 2009 in Figure S2. Previous simulations 34 at 2.5 km horizontal resolution suggested that salinity anomalies are trapped by the energetic current system along the shoreline and on the shelf, winds are key to exporting the meltwater signal off-shore 34,36 , and the meltwater modulation of deep convection is likely small in current meltwater conditions. This is confirmed by the present experiments, independently of resolution. As the resolution increases towards submesoscale permitting, we observe additional export of meltwater towards the central Labrador Sea. The surface meltwater signal delimits in all cases the edges of the convective area ( Figure S3) but has only a small impact on the overall convective activity for a given resolution ( Table 1). The impact is inversely proportional  www.nature.com/scientificreports/ to model resolution. Whenever the resolution is low, some GrIS meltwater can be advected into the convective area by the broader and more diffused boundary currents.
To verify that the representation of convective activity improves in realism whenever submesoscale processes are included, in Fig. 4A,B we compare the simulated temperature profiles averaged in a region that encompasses the common convective area (MW-experiment) over the integration period, and the same quantity reconstructed from ARGO floats. The average number of ARGO profiles in each month is about 20. The best representation of the mean temperature stratification is provided by LBR1 with LBR1.7 (not shown) being a close second, while a cold bias characterizes the other two resolutions, in line with the greater convective volume. We note the poor representation of the details in the vertical structure found in the observed temperature profiles in the upper 400 m. The limited vertical resolution adopted is likely responsible for it. Indeed, 50 vertical levels would be needed to resolve the first baroclinic mode, with an additional 25 levels per subsequent mode 68 .
The very small differences between the integrations with and without meltwater suggest that the salinity anomalies linked to the GrIs do not contribute significantly (yet) to the stratification of the LS interior. Differences in overall heat transport must therefore drive the resolution dependence by controlling the stratification in Table 1. Convective volume as a percentage to the volume produced in the lowest resolution 15 km case in the absence of melt water input from the GrIS (LBR15-NoMW) calculated over the convective period January to May and 2008-2013.   Figure S5. The spatial structure of both net mean and eddy contribution is far more complex in the three finer resolutions (LBR5, LBR1.7 and LBR1) than in LBR15 where a local maximum cannot be distinguished, but the integral of the net mean contribution where convection can take place-the LBR15 convective patch-does not vary significantly across resolutions. For the eddy component (Fig. 5A), on the other hand, very high values of E are located off the Greenland coast, in the region of eddy formation, in the mesoscale resolving and submesoscale permitting simulations (LBR5, LBR1.7 and LBR1), consistent with 57,69 . Furthermore, in LBR1.7 and LBR1 the convergence of heat by eddy advection is large also in the central Labrador Sea following the IR path, with the maximum values reaching E > 500 W m −2 to the north-northeast of the LBR1.7 and LBR1 convective patches. This heat convergence therefore delimits the convective patch in the submesoscale permitting runs by controlling the stratification in the off-shore portion of the basin. The time-mean eddy contribution to heat advection over the LBR5 convective area is 70 W m −2 and over the same area increases to 83 and finally 90 W m −2 in LBR1.7 and LBR1, respectively.

In comparison to LBR 15 km (NoMW)
Submesoscale circulations, and specifically mixed-layer eddies, can also contribute to restratifying locally the upper portion of the water column 70 , and this contribution can be parameterized in climate models 71 . This local restratification mechanism can be evaluated through the available potential energy (APE) release in the convective season, from December to April (Fig. 5B). Such release is estimated by the vertical eddy flux of buoyancy, w ′ b ′ , where b is buoyancy, here integrated over the upper 400 m of the water column, given the vertical extent of the submesoscale circulations and the core of the IR (see Fig. 7 below). The APE release is small in LBR15 everywhere but in the convective area, as to be expected, and increases by increasing resolution. Such increase is strong along the western portion of the basin between the 2000 and 3000 m isobaths in all other runs, and also along the IR path and the Greenland slope in LBR1.7 and LBR1. The APE release inside each respective convective area does not vary significantly or linearly between mesoscale resolving and submesoscale permitting runs (149, 146 and 158 × 10 -10 m 2 s −3 from LRB5 to LBR1), but increases steadily by 20 and 30 × 10 -10 m 2 s −3 in the region comprised between the western boundary of the LBR1.7 and LBR1 patches and the 2000 m isobath to the west of it.
Overall, the latitudinal shrinking of the convective patch for increasing model resolution is dominated by the eddy-driven heat convergence, while the longitudinal shrinking is modulated by both mechanisms, with the  While the meltwater signal is confined near the surface and does not influence convection strongly, the advection pathways of the meltwater from Greenland into the interior provide another indication of the processes at play. In the submesoscale permitting cases, advection from the boundary currents is accomplished by the longliving population of anticyclonic vortices and cyclonic vorticity filaments characterized by Ro of order 1. The submesoscale process behind the subsurface intensification of these circulations in LBR1 or LBR1.7 compared to the lower resolution cases is, predominantly, strain-induced frontogenesis 72  www.nature.com/scientificreports/ encompassing the IC core located between 100 and 300 m depth and their generation and overall impacts are independent of the meltwater input. They form in correspondence of the steepening of the bathymetry and result from baroclinic instability of the horizontal shear layers induced by the warm IC that meets the steep continental slope 24 . Resolution influences their characteristics in two ways. First, the steepness and veering of the actual bathymetry are better resolved in the LBR1 (and LBR1.7) than in LBR5 (or LBR15), contributing to the confinement, intensification and increased variability of the IC (Fig. 7C,D). Second, the horizontal shear layer coinciding with a non-zero vertical component of the vorticity tensor that extends into the interior, because the interior mean flow has to be identically zero at the sloping bottom, is better resolved and intensified through straininduced frontogenesis (Fig. 6C,D, Figure S4 and Video VS1) in the submesoscale permitting cases. In sum, in www.nature.com/scientificreports/ LBR1 and LBR1.7 this horizontal shear layer is more intense and has an elevated positive frontogenic tendency. There is indeed a nearly perfect linear relationship also between the frontogenetic tendency along the coast of Greenland at 200 m depth and the convective volume, similarly to what seen for surface vorticity (Fig. 6F). The end result is the formation of eddies of size comparable to the mesoscale case, being the size controlled by the bathymetry 21,24 , but of submesoscale strength (Ro ~ 1) being originated from the submesoscale fronts ( Fig. 7E-H).
The eddies are then surface-intensified by the strong and variable winds. The IRs in the LBR1.7 and LBR1 simulation are not only stronger, but they are also longer living that the LBR5 counterpart, are surrounded by submesoscale filaments, and, together with the filaments, contribute more effectively to the advection of warm water into the center of the LS. They form abundantly especially in late fall and winter, when the winds and the boundary current are stronger and more variable 56 , and they stratify the portion of the basin they travel to, effectively delimiting the convective patch and its volume. In this regard, Figure S3, showing the distribution of passive dye released at the fjords at the various resolution is illuminating. In doing so, they also transport freshwater originated from the GrIS melting, but the freshwater anomaly remains within the stratified region, where convection does not occur anyway, because of the IC heat contribution and the submesoscale eddy-induced restratification.

Discussion
Our findings support the existence of a direct link between submesoscale instabilities and the Labrador Sea contribution to the global meridional overturning circulation. Model results reveal surprisingly large differences in the representation of Labrador Sea Water formation for varying resolution, with a 50% reduction in the modeled volume of convected waters (calculated over the area where convection reaches depths greater than 1000 m) when mesoscale advection is included, and over 80% reduction if submesoscale processes are accounted for. Current state-of-the-art climate models do not fully resolve the mesoscale dynamics at high latitudes, and indeed tend to greatly overestimate the formation of Labrador Sea water 16 . Interestingly, the model predicts a linear relationship between the amount of LSW formed, and the mean near-surface vorticity and/or the frontogenetic tendency along the west coast of Greenland, opening the possibility to a physically based parameterization of this contribution.
In the submesoscale permitting simulations, mesoscale anticyclonic eddies with submesoscale characteristics (in primis a Rossby number of order 1) and submesoscale cyclonic vorticity filaments form in correspondence of the intense horizontal shear layers induced by the Irminger Current interacting with the continental slope. These long-lived eddies (longer-lived and more intense than their counterpart in the mesoscale-resolving simulation) carry heat from the Irminger Current to the interior, effectively delimiting where deep convection takes place. They form independently of the near surface stratification along the Greenland coast, their generation being controlled by baroclinic instability associated with the bathymetry and not by the presence/absence of freshwater inputs from the Arctic and Greenland Ice Sheet. At the same time, the large number of submesoscale eddies formed locally in the region bounded by the 2000 m and 3000 m isobaths in the western portion of the basin controls locally the stratification, limiting the lateral extension of the convective patch towards the Labrador shelf.
This work provides physical context to the recent attribution of the differences in sea-level rise between the penultimate and last deglaciation to differences in subsurface warming in the North Atlantic subpolar gyre 75 . It also helps understanding the lesser role of the Labrador Sea in the overall control of the state of the MOC argued through the analysis of recent OSNAP (Overturning in the Subpolar North Atlantic Program) data with respect to estimates from climate models 19,76 , while indicating that the representation of the interannual variability of the LSW formation can be captured independently of resolution.
Ultimately, these findings call for observational efforts aimed at carefully monitoring heat content and trends of the Irminger Current and point to the need for better investigating the rich complexity of interactions and feedbacks between processes localized at the ocean boundaries (in this case along the continental slope of West Greenland) and the global ocean, as they can affect the climate trajectory of our planet. They also call for parameterizations in Earth system models that account for the advective role of coherent vortices generated by submesoscale instabilities in regions where submesoscale processes cascade to the larger (from basin to global) scales. For the LSW this can be easily achieved by accounting for the linear relationship between convective volume and surface mean vorticity or frontogenetic tendency. Predicting the future evolution of the meridional overturning circulation and its persistence as global climate change progresses will hinge upon our ability to model correctly the response of a multiscale environment.