Oceanic eddies are an important component in preconditioning the central Labrador Sea (LS) for deep convection and in restratifying the convected water. This study investigates the different sources and impacts of eddy kinetic energy (EKE) and its temporal variability in the LS with the help of a 52-yr-long hindcast simulation of a 1/20° ocean model. Irminger Rings (IR) are generated in the West Greenland Current (WGC) between 60° and 62°N, mainly affect preconditioning, and limit the northward extent of the convection area. The IR exhibit a seasonal cycle and decadal variations linked to the WGC strength, varying with the circulation of the subpolar gyre. The mean and temporal variations of IR generation can be attributed to changes in deep ocean baroclinic and upper-ocean barotropic instabilities at comparable magnitudes. The main source of EKE and restratification in the central LS are convective eddies (CE). They are generated by baroclinic instabilities near the bottom of the mixed layer during and after convection. The CE have a middepth core and reflect the hydrographic properties of the convected water mass with a distinct minimum in potential vorticity. Their seasonal to decadal variability is tightly connected to the local atmospheric forcing and the associated air–sea heat fluxes. A third class of eddies in the LS are the boundary current eddies shed from the Labrador Current (LC). Since they are mostly confined to the vicinity of the LC, these eddies appear to exert only minor influence on preconditioning and restratification.
The Labrador Sea (LS) in the western subpolar North Atlantic is one of the few regions in the world’s oceans where deep convection occurs on a regular basis. Large air–sea heat fluxes, sometimes exceeding 1000 W m−2, during the winter months cause the surface waters to loose buoyancy until the stratification of the water column is overcome and deep convection is triggered (e.g., Marshall and Schott 1999). The strength of the stratification is the limiting factor that usually prevents deep convection and thus a preconditioning of the LS is required. In this context preconditioning refers to a cyclonic large-scale circulation with upward doming isopycnals at its center causing a weakly stratified thermocline and a weak stratification in the near-surface ocean. The preconditioning and the LS’s stratification can be influenced by a variety of processes linked to the atmospheric circulation, such as the strength of the subpolar gyre (Treguier et al. 2005) or Ekman advection (Schulze Chretien and Frajka-Williams 2018). However, over the past two to three decades, numerous studies on the basis of observations and numerical simulations suggested that oceanic mesoscale eddies play a major role in determining both the location and strength of deep winter convection and the subsequent restratification of the convected water mass during spring and summer. The “ocean mesoscale,” according to the quasigeostrophic framework, refers to the horizontal scale that is similar to or larger than the local first baroclinic Rossby radius of deformation (Ferrari and Wunsch 2009; Su and Ingersoll 2016), which is between 5 and 10 km in the LS (Gascard and Clarke 1983; Gelderloos et al. 2011). The Rossby radius has the approximate scale of the most unstable wavelength of baroclinic instability and Rossby numbers of the mesoscale are usually much smaller than one. In the LS, these mesoscale eddies can hinder deep convection by transporting heat and freshwater toward the interior, thereby increasing the stratification in the central LS.
Three types of eddies can be identified in the LS: Irminger Rings, convective eddies, and boundary current eddies.
a. Irminger Rings
Irminger Rings (IR) originate in the West Greenland Current (WGC) at the west coast of Greenland between 61° and 62°N (Heywood et al. 1994; Brandt et al. 2004; Zhang and Yan 2018). Near Cape Desolation (60.44°N, 48.10°W), the WGC encounters steep topography that induces instabilities, leading to the generation of the anticyclonic IR. The nature of these instabilities is not fully understood to date. While earlier studies with ocean general circulation models (OGCM) found barotropic instabilities (BT) induced by the large horizontal shear of velocities to be the main driver (Eden and Böning 2002; Chanut et al. 2008), simulations with idealized numerical models pointed to baroclinic instabilities (BC) near the bottom (Katsman et al. 2004; Bracco et al. 2008) that depend on the horizontal gradient of density and thus, through the thermal wind balance, the vertical gradient of velocity. More recently Luo et al. (2011) also found near-bottom BC to be the major driver in an OGCM, while Zhu et al. (2014) point out the importance of both types of instabilities in their OGCM.
The IR have diameters of 30–60 km, have a relatively warm core (Lilly et al. 2003), and propagate toward the west and southwest with a speed of 3–6 cm s−1 (Katsman et al. 2004; Chanut et al. 2008; Zhang and Yan 2018). Mixing with the surrounding water in the northern LS increases the stratification there, effectively prohibiting deep convection (Chanut et al. 2008). Whether or not IR also aid in rapidly restratifying the central LS after the convection period is a matter of ongoing debate. Some studies find IR to have relatively little effect on the rapid restratification in spring (Chanut et al. 2008; Zhang and Yan 2014) while others point to a shared importance of IR and other sources (Gelderloos et al. 2011) or even attribute a major role in this process to the IR (Katsman et al. 2004). Part of these discrepancies can possibly be explained by the difficulties of many models to simulate the convective region at the right location. Often it is simulated too far north (e.g., Saenko et al. 2014) in the westward path of the IR.
The eddy kinetic energy (EKE) maximum generated near Cape Desolation exhibits a seasonal cycle with a maximum in winter (January–March; Lilly et al. 2003; Brandt et al. 2004; Zhang and Yan 2018) in accordance with the winter maximum of the WGC velocity in that region. On interannual time scales the variations of EKE and IR in the WGC are less well understood. Brandt et al. (2004) speculate that the density gradient between the central LS and the boundary current is responsible for a stronger WGC and thus EKE, which is supported by a model study by de Jong et al. (2016). Zhang and Yan (2018) attribute the changes to the strength of the Subpolar Gyre (SPG) circulation.
b. Convective eddies
Convective eddies (CE) are thought to originate from the rim current along the boundary of the convective patch driven by the large lateral buoyancy gradient between the convected water and its surroundings (Gascard and Clarke 1983; Marshall and Schott 1999). Theoretical and idealized numerical studies show the rim current to be baroclinically unstable and to generate cyclonic and anticyclonic CE (Send and Marshall 1995; Jones and Marshall 1997; Lilly et al. 2003). However, the few CE that have been observed were exclusively anticyclonic (Lilly et al. 2003). Their diameter ranges from 10 to 36 km and maximum velocities of 0.1–0.3 m s−1 are found at middepth around a core containing water from within the convective area, usually cold and fresh compared to the surroundings. In OGCMs with sufficiently high resolution, EKE generation is often found in the central LS (Chanut et al. 2008; Luo et al. 2011; Zhu et al. 2014), however CE have not yet been described in detail in a simulation with realistic geometry and forcing.
Due to their close proximity to the convective region, CE and the associated baroclinic instabilities are supposed to play a role in the quick restratification in spring due to the lateral buoyancy flux they induce (e.g., Jones and Marshall 1997). However, as with IR, their relative importance is debated. Gelderloos et al. (2011) attribute only minor restratifying effects to the CE, while Chanut et al. (2008) find them to be essential for the restratification right after convection. The importance of CE for the central LS is underlined by float observations between 1000- and 1500-m depth showing EKE at these depths to dominate over surface EKE in most of the convection region (Fischer et al. 2018).
The seasonal variability of CE is dictated by the seasonal cycle of convection, making EKE in the central LS strongest in early spring (March, April; Brandt et al. 2004; Luo et al. 2011; Zhang and Yan 2018). However, most studies do not explicitly investigate CE, so that the variability could also arise from mesoscale features propagating into the region (Brandt et al. 2004). The interannual variations in CE generation are thought to be linked to the strength of deep convection (Lilly et al. 2003; Brandt et al. 2004).
c. Boundary current eddies
Boundary current eddies (BCE) are generated in Labrador Sea’s boundary currents, the WGC, and the Labrador Current (LC; Chanut et al. 2008). At the eastern boundary, BCE are generated between the southern tip of Greenland and Cape Desolation. Most of the BCE at the western boundary are generated between 56° and 60°N along the shelf (Eden and Böning 2002; Brandt et al. 2004), with a pronounced maximum in winter when the boundary currents are strongest. BCE originate from baroclinic instabilities of the WGC and LC, are surface intensified, and generally do not propagate far from the boundary (Chanut et al. 2008). However, in the western LS, the convection area is situated close to the LC. An impact of the BCE on the convection could thus be possible.
The mesoscale variability of the LS is far from understood in all its details. While the general properties of the different types of eddies in the LS are well defined, there is still a lack of explanation for some of the generation mechanisms, as well as their relative importance for preconditioning and restratification. We address these questions by combining three important aspects in one simulation: a high-resolution (1/20°) ocean model, in realistic geometry, run over a 52-yr period with forcing from a well-established reanalysis product. Using a realistic instead of an idealized geometry is expected to yield a better representation of the locations and timings of convection and mesoscale activities. The higher resolution compared to earlier studies with OGCMs (e.g., Zhu et al. 2014) improves the representation of the different types of eddies in the LS. The long hindcast simulation enables us to investigate temporal variations beyond interannual time scales. This is important as, especially near strong currents, a large part of the temporal variability of the mesoscale is intrinsic to the ocean, that is, not attributable to a certain atmospheric forcing (Sérazin et al. 2015; Rieck et al. 2018), even at interannual time scales.
After describing the model, data, and methods in section 2, we focus on investigating the properties, generation mechanisms, and potential impact on deep convection of Irminger Rings (section 3a), convective eddies (section 3b), and boundary current eddies (section 3c), summarizing and discussing the results in section 4.
2. Model, data, and methods
This study utilizes a global ocean general circulation model based on the Nucleus for European Modeling of the Ocean (NEMO) code version 3.6 (Madec 2008), developed as part of the DRAKKAR collaboration. The ocean component of NEMO is based on Océan Parallélisé (OPA; Madec et al. 1998), while the sea ice component is the Louvain-la-Neuve Ice Model version 2 (LIM2; Fichefet and Maqueda 1997; Vancoppenolle et al. 2009). The specific configuration used here, VIKING20X, is an extended and updated version of VIKING20. It has been built on the experiences with VIKING20 in simulating various aspects of the subpolar North Atlantic, ranging from the Denmark Strait overflow (Behrens et al. 2017) and the circulation in the Deep Western Boundary Current (Fischer et al. 2015; Handmann et al. 2018) to the flow field along the North Atlantic Current (Mertens et al. 2014; Breckenfelder et al. 2017) and the influence of Greenland meltwater on the Labrador Sea (Böning et al. 2016). VIKING20X consists of a global, orthogonal, curvilinear, tripolar Arakawa-C type grid with a nominal resolution of 1/4° in longitude, which is refined to 1/20° in the Atlantic between 34°S and 70°N using adaptive grid refinement in FORTRAN (AGRIF; Debreu et al. 2008). Both the base model and the refined nest share the same 46 vertical levels with 6-m thickness at the surface, increasing toward a maximum of 250 m in the deep ocean and a partial-cell formulation at the bottom (Barnier et al. 2006). In VIKING20X, the advection of tracers is discretized using a combination of an upstream and a centered scheme, the total variance dissipation (TVD) scheme (Zalesak 1979). The advection of momentum is accomplished by an energy- and enstrophy-conserving second-order centered scheme adapted from Arakawa and Hsu (1990) and modified to suppress symmetric instability of the computational kind (Ducousso et al. 2017).
An atmospheric forcing based on the Coordinated Ocean-Ice Reference Experiments (CORE.v2; Griffies et al. 2009; Large and Yeager 2009) along with the associated bulk formulations is used. The model is spun up from rest with temperature and salinity from climatology (PHC2.1; updated from Steele et al. 2001) and an initial sea ice field from 31 December 1992 from a previous 1/4° simulation. The 30-yr spinup period is forced by the interannually varying atmospheric forcing for the period 1980–2009. A snapshot of the oceanic state at the end of this spinup is then used as initial state for the hindcast simulation from 1958 to 2009. The model simulation applies a weak sea surface salinity restoring (SSSR) to climatology of 33.33 mm day−1, which corresponds to a relaxation time scale of 1500 days over a 50-m surface layer. Convection in the model is parameterized by enhanced vertical diffusion (10 m2 s−1) of tracers and momentum whenever the local Brunt–Väisälä frequency is negative.
Preliminary tests showed the lateral boundary condition to be a critical aspect to the simulation of IR in the WGC. Accordingly, a no-slip condition is implemented in the WGC (59°–62°N, 43°–51°W) to represent the topographic effect on the current necessary to produce instabilities that generate IR, while the model uses a free-slip condition elsewhere. The reliance of the generation of IR on this no-slip condition is exploited to study the impact of IR on deep convection in the central LS. An additional sensitivity simulation differing only in the definition of a free-slip lateral boundary condition in the WGC is forked off the hindcast on 1 January 2007 and run for 3 years with the same atmospheric forcing.
The influence of resolution on CE in the central LS is investigated by comparing the 1/20° behavior to a 1/12° simulation (ORCA12), which is driven by the same forcing as the 1/20° hindcast, for the period 1 January 1980–31 December 1991. Allowing for a short spinup phase, the years 1983–91 are compared to the 1/20° simulation.
EKE is defined to be , where . The terms are the 5-day mean zonal and meridional velocities, respectively. The mean velocities are based on 1-yr averages to avoid interannual current variability from being considered EKE (cf. Penduff et al. 2004; Rieck et al. 2015). The observational estimate for EKE is derived from the gridded, delayed time geostrophic surface velocities made available by EU Copernicus Marine Service (CMEMS) and calculated as described above. Additionally, 1-Hz along-track sea level anomaly (SLA) data developed, validated, and distributed by the Centre for Topographic Studies of the Oceans and Hydrosphere (CTOH)/Laboratoire d’Etudes en Géophysique et Océanographie Spatiales (LEGOS), France (X-TRACK; Birol et al. 2017) are used in the Labrador Sea (TOPEX/Poseidon, Jason-1/2/3, Track 72), to account for the underestimation of EKE derived from gridded SLA products at high latitudes (Zhang and Yan 2018). The SLA is low-pass filtered with a ~40-km cutoff locally estimated scatterplot smoothing (LOESS) filter, and EKE is calculated from the across-track geostrophic velocity, derived from the along-track SLA gradient. To compare the modeled EKE to this observational value, daily mean simulated SLA are subsampled with a nearest-neighbor lookup to the times and locations of the satellite-derived SLA and then treated identically.
An eddy detection and tracking method based on Chelton et al. (2011) and Oliver et al. (2015) is implemented to investigate the number, strength, and other statistics of IR. The method has been adapted to detect and track IR based on the Okubo–Weiss parameter (Okubo 1970; Weiss 1991) with a threshold parameter of 0.2 × 10−9 (e.g., Chelton et al. 2007) derived from simulated horizontal velocities at 191-m depth (at this depth, IR are well isolated and not masked by other, more surface intensified mesoscale features), contrary to the original implementation by Chelton et al. (2011), which is based on sea surface height. Irminger Rings are defined to be anticyclonic eddies originating from 59.5° to 62.5°N and from 46° to 53°W that exist for at least 15 days and are larger than ~300 km2 upon initial detection (depending on their location of generation).
To estimate the respective contribution of baroclinic and barotropic instabilities to the generation of EKE, the energy transfer terms are considered (Wright 1981; Böning and Budich 1992; Zhu et al. 2014). Baroclinic instability can be described in terms of the transfer from mean to eddy available potential energy (MAPE and EAPE, respectively):
where ρ0 = 1025 kg m−3 is the reference density; is the annual mean Brunt–Väisälä frequency , with σ being the horizontally averaged simulated density and z the depth; and and are the horizontal velocity anomalies with respect to low-pass filtered with a boxcar window of length 100 days. Analogously, is the density anomaly. Barotropic transfer from mean kinetic energy (MKE) to EKE is expressed as the work of the Reynolds stresses against the mean shear:
Both terms, and , must be considered approximations in this context, as they only represent energy conversions when integrated over a closed volume. When considering subregions as done in this study, advection of energy by the mean flow or the mean eddy advection and divergence of the pressure work can influence and (Böning and Budich 1992). Nevertheless, the two terms can be used to address the relative importance of baroclinic and barotropic instability, and their temporal and spatial variability is a useful measure to contrast times and regions of higher instability to those of lower instability.
Decadal changes of the circulation of the SPG are indicated by the SPG index calculated following Berx and Payne (2017) and Hátún and Chafik (2018) as the principal component of the first empirical orthogonal function (EOF) of the sea surface height in the subpolar North Atlantic between 40° and 65°N.
The model simulation is well representing the features of the circulation in the subpolar North Atlantic. Figure 1 illustrates the geographical distribution of the near-surface currents in the region of interest. The North Atlantic Current (NAC) along with the Northwest Corner at ~45°W, 50°N stands out as a swath of vigorous eddying motions. The boundary current system in the western subpolar gyre is composed of the East Greenland Current (EGC), the WGC, and the LC. Different types of eddies are expected to be found in the central LS: the Irminger Rings in the northern LS, the convective eddies in the central LS around the region of deep convection, and the boundary current eddies. In the following sections the simulated characteristics of these three types of eddies are studied.
a. Irminger Rings
The 1/20° simulation is well resolving the large anticyclonic (negative relative vorticity) eddies originating in the WGC near Cape Desolation (Fig. 2). Relative vorticity is calculated from 5-day mean velocities as and then normalized by f, the Coriolis frequency, such that the absolute value of is the Rossby number. The WGC can be seen along the west coast of Greenland between Cape Farewell and 61°N with positive vorticity offshore of the current and negative vorticity toward the coast. Additionally, three IR can be clearly observed between 60° and 62°N. Normalized relative vorticity in this snapshot is generally elevated between 59° and 62°N and between 48° and 58°W, with absolute values > 0.4, while it is below 0.1 in the region where convection occurs. The deep convection is indicated by the March mixed layer depth (MLD).
Sections through one of these IR at 61.1°N, 56.5°W in Fig. 2 are shown in Fig. 3a. It compares favorably with observations (Hátún et al. 2007; de Jong et al. 2014), exhibiting a surface intensified structure with velocities well above 0.4 m s−1 in the upper ocean, a warm and saline core with anomalies compared to the surroundings of up to 0.5°C and 0.02, respectively, and a fresh cap in the upper 100–200 m. This vertical structure stems from the WGC, consisting of cold, fresh waters originating in the East Greenland Current on top of warmer, saline waters from the Irminger Current. Relative vorticity normalized by planetary vorticity () exhibits a minimum < −0.3 in the core, surrounded by a maximum of up to 0.2 which shields the ring’s center from the surroundings. Another Irminger Ring (Fig. 3b) has significantly different properties in the upper ~1000 m about 2 months later. While propagating farther south (58.8°N) it has experienced convection and displays nearly homogeneous temperature and salinity in the upper 1000 m, as well as a weaker, more barotropic velocity structure (0.2 m s−1). This compares well with observations of IR in the central LS, which do not show the fresh cap observed farther north (Lilly et al. 2003).
The integrated effect of the IR is shown in Fig. 4a. An EKE (at 94-m depth) maximum of more than 250 cm2 s−2 near the coast of Greenland between 60° and 62°N is extending toward the west and southwest as a result of the IR propagation into the interior. The location of this maximum compares well with EKE calculated from gridded SLA data from satellite altimetry (Fig. 4b). The EKE in the model, however, is about twice as large as in observations. This discrepancy, at least partly, arises due to the inability of gridded SLA data to resolve mesoscale variability at high latitudes (Zhang and Yan 2018). The SLA data used in this study are on a 1/4° grid and thus only resolve eddies with scales larger than ~30 km in the Labrador Sea (Chelton et al. 2011). When comparing EKE based on higher-resolution along-track SLA data to the simulated EKE, the discrepancy is significantly reduced (Fig. S1 in the online supplemental material).
The simulated IR are detected within the region 59.5°–62.5°N, 46°–53°W and are tracked with a detection and tracking algorithm based on Chelton et al. (2011) and Oliver et al. (2015) (see section 2). A total of 1939 (38.8 yr−1) IR are identified between 1960 and 2009. This number is higher than previous estimates based on observations (Lilly et al. 2003) and model simulations (Chanut et al. 2008) that found roughly 30 eddies (cyclonic and anticyclonic) being generated off Cape Desolation. The high number of IR detected in this study can be explained by frequent deformation, splitting and merging of eddies close to the WGC. The tracking algorithm is not able to keep track of the IR under these circumstances and probably detects them more than once. Indeed, the number of IR leaving the region of detection to the west or south is significantly lower at 204 (4.1 yr−1) and in better agreement with Lilly et al. (2003) and Chanut et al. (2008) who find 1–2 IR per year south of 58°N. The simulated IR that propagate offshore have a mean diameter of 36 km (based on their vorticity at 191 m) and a mean lifetime of 147 days while roughly 10% of them exist for longer than 1 year.
The most important effect of the IR is the limitation of deep convection to the southern central LS. While winter heat loss is large over most of the central LS, deep convection, approximated by the mean March MLD, only occurs south of 59°N due to the stratifying effect of the IR farther north, where IR flux buoyancy toward the interior LS (Fig. S2). The location of the deepest MLD around 57°N, 55°W in the model also compares well to MLD from ARGO float observations (Fig. 4) for the period 2000–09, although the extent of the convective area appears to be overestimated, especially northwest of 58°N, 55°W. For the MLD from ARGO, the variable density threshold method depending on the local temperature and salinity is used (Holte et al. 2017). In the cold waters of the LS this variable threshold is close to 0.01 kg m−3. Accordingly, a density threshold of 0.01 kg m−3 is also used in the model simulation.
The hypothesis that IR mainly affect the preconditioning is additionally supported by a sensitivity study that aims to reduce the number of IR shed from the WGC. A free-slip lateral boundary condition is implemented at the west coast of Greenland, whereas in the hindcast simulation a no-slip condition is used. Free-slip reduces the lateral and vertical shear near the boundary, and thus significantly weakens the EKE generation by the local baroclinic and barotropic instabilities in the WGC (Figs. 5a,c,d). As a result, almost no IR are shed in this simulation and in consequence a reduced stratification in the northern central LS can be seen, with the mean MLD of March 2008 up to 1000 m deeper than in the simulation with IR (Fig. 5b).
Arguably, the elimination of all IR is not a realistic scenario. Nevertheless, Fig. 5 illustrates the impact of the IR on the lateral extent of deep convection in the central LS. It also shows that the response of the northern central LS MLD to a reduced number of IR is rather quick. Thus, a connection between the IR shed from the WGC and the MLD in the northern central LS can be expected even though the changes in EKE due to realistic variations in the forcing are much smaller compared to the sensitivity study above.
1) Temporal variability
We first investigate the seasonal cycle of EKE in the WGC, as this is observationally well constrained and can be used to test the model’s ability to simulate the temporal variability of the instability processes occurring near Cape Desolation. A box ranging from 60° to 62°N and from 49° to 56°W (cf. Fig. 4) is chosen to average the EKE maximum off the WGC. The simulated EKE peaks in February at 200 and 140 cm2 s−2 in the upper (112–322-m depth) and the deeper (382–1655-m depth) ocean, respectively (Fig. 6), which compares well with observations (Brandt et al. 2004) and other model studies (Luo et al. 2011). The two energy conversion terms and described in section 2, averaged from 60.0° to 61.5°N and from 48.0° to 52.0°W, are used to approximate the baroclinic (BC) and barotropic (BT) instability processes. The results gleaned from this analysis are qualitatively independent of the exact location and size of the box chosen to average (not shown). Both BC and BT exhibit a similar seasonal variability in the upper and deeper ocean, peaking in winter. The BC in the deeper ocean is the largest of the instability processes and peaks in January at 4.6 × 10−5 kg m−1 s−3, while BC in the upper ocean has its maximum of 2.6 × 10−5 kg m−1 s−3 in November. Barotropic instabilities in the upper ocean evolve similarly to BC in the deeper ocean with a maximum (3.0 × 10−5 kg m−1 s−3) in January. BT in the deeper ocean is low (0.5–1.0 × 10−5 kg m−1 s−3) throughout the year and does not show a pronounced seasonal cycle. Generally, both routes of energy transfer, from MAPE to EAPE as indicated by BC and from MKE to EKE as indicated by BT, are of importance in the generation region of the IR. While the baroclinic transfer dominates in the deeper ocean, the barotropic transfer appears as the major term in the upper ocean. These findings are in general agreement with previous idealized model studies (Katsman et al. 2004; Bracco et al. 2008) and more recent studies using OGCMs (Zhu et al. 2014), and deviate from earlier model studies (Eden and Böning 2002; Chanut et al. 2008). It should be noted that, while these results indicate the importance of BC and BT for the whole generation process of IR, they do not yield information about the exact instability process triggered near Cape Desolation. A detailed investigation of the WGC’s interaction with the topography and the initial instabilities is beyond the scope of this study.
After establishing the models ability to produce instabilities around Cape Desolation with a seasonal variability that supports the simulated and observed seasonal cycle of EKE, we now focus on variations on longer time scales. As the ocean integrates the high-frequency atmospheric forcing into lower-frequency variability of the general circulation, links between year-to-year variations in the atmosphere and large-scale ocean circulation are hard to detect. In fact, oceanic EKE variability is largely intrinsic, that is, internally driven, in most regions on interannual time scales (Wilson et al. 2015; Rieck et al. 2018), which might explain the difficulties previous studies had relating the interannual EKE variability in the WGC to an external forcing (e.g., Luo et al. 2011). In the following, we thus focus on variations on decadal time scales.
A subpolar gyre index is defined as described in section 2 to indicate periods of stronger and weaker circulation in the North Atlantic subpolar gyre (Fig. S3). Based on this index, the two periods 1966–75 and 1988–97 are chosen, representing weak (SPG−) and strong (SPG+) circulation in the SPG, respectively.
During SPG+, the boundary currents of the LS are significantly stronger compared to SPG− (Fig. 7a). The WGC is additionally more confined to the shelf. The current velocity averaged from 382- to 1655-m depth reduces offshore, while the velocity closer to the shelf increases by more than 10 cm s−1. The WGC, although stronger, appears to be more stable and locked to the topography. More precisely, BT in the upper ocean decreases by more than 5 × 10−5 kg m−1 s−3 between 60.0° and 60.5°N, where the velocity increase is strongest, and increases by a similar amount farther downstream, north of 60.5°N (Fig. 7d). Similarly, BC in the deeper ocean decreases by more than 5 × 10−5 kg m−1 s−3 at the shelf between 60.0° and 60.5°N and increases by more than 10 × 10−5 kg m−1 s−3 farther downstream where the velocity increase is strongest in this depth range (Fig. 7c). The behavior of the instability processes during SPG+ influences the generation of IR and thus EKE. In accordance with BT and BC, EKE decreases by up to 40 cm2 s−2 during strong WGC circulation east of 50°W and south of 61°N, while it increases by more than 60 cm2 s−2 in a large area farther downstream and offshore (Fig. 7), which indicates an increased number of stronger IR propagating into the northern central LS. Indeed, during SPG+ a total of 388 IR are generated of which 3.8 yr−1 leave the formation region, compared to 375 IR during SPG− of which 3.4 yr−1 propagate offshore. Additionally to the arguably small difference in IR numbers, IR during SPG are on average larger (734 km2) and stronger (mean relative vorticity minimum of −1.58 × 10−5 s−1) compared to the SPG− period (717 km2 and −1.51 × 10−5 s−1, respectively).
Three boxes have been chosen to further investigate the processes in the WGC. Box 1 (59.5°–60.7°N, 47.0°–49.3°W) encompasses the region where, despite stronger velocities, the instabilities are reduced during SPG+. Box 2 (60.0°–61.5°N, 49.5°–52.0°W) is the area with stronger instabilities during SPG+ and Box 3 (60.0°–61.5°N, 52.5°–57.0°W) is located offshore in the path of the IR to capture the variations of EKE in response to a stronger SPG. Generally, EKE and BT act in concert in one region, while BC has a larger influence on downstream EKE. For long time scales (time series are smoothed with a 5-yr running mean), the EKE in Box 2 is positively correlated to BT in the same box (0.68) and to BC in the upstream Box 1 (0.72), all correlations are significant at the 95% level. However, neither EKE and BT in Box 2, nor BC in Box 1 are correlated to the SPG index, indicating that, even at long time scales, much of the EKE generation is chaotic and only when averaging over decades, an effect can be seen (Fig. 7). Farther downstream, BC in Box 2 is positively correlated to EKE offshore in Box 3 (0.86) and this EKE in Box 3 is correlated to the large-scale circulation, that is, the SPG index (0.90), as already indicated by Fig. 7. Thus, moving farther downstream from the instabilities in the WGC has an integrating effect, organizing the intrinsic oceanic variations into longer time-scale variability. This complex pattern explains contradictory results with regard to the instability processes in earlier studies. Depending on the exact region and depth range chosen, the interannual and decadal variability of BC, BT, and EKE is different, and additionally does not always have a local expression.
As demonstrated above, the generation of IR in the WGC is strongly connected to the large-scale circulation of the SPG. It controls the extent of the convection region to the north through the stratifying effect the IR exert on the northern central LS. Whether the IR variability is responsible for decadal changes of the northward extent of the convection region cannot be determined in a realistic scenario, because the area where deep convection occurs is to first order controlled by the local atmospheric forcing. Only in idealized cases, where EKE in the WGC is altered dramatically, an effect on the extent of the convection region can be observed (Fig. 5). Note, that there is a significant negative correlation of −0.39 of EKE off the WGC (i.e., IR) and the northward extent of the convection region (approximated by the area, where MLD > 1000 m north of 57°N, between 50° and 60°W) for 1958–88 (not shown), however, there is no significant correlation in the last 20 years.
2) Impact on restratification
A first hint at the role of IR on the restratification of the convection region in spring can be gleaned from Fig. 5. When the number and strength of IR is drastically reduced, the effect on the MLD in the southern central LS is minimal. This suggests that, while the effect on the northern central LS is large, only few IR propagate to south of 58°N and immediately affect the deep convection there. Additional support for this hypothesis is given by Fig. 8. The zonal mean EKE is observed to propagate southward from the generation region. This happens at a speed of ~5 cm s−1, in agreement with previous studies ranging from 3–4 (Chanut et al. 2008) to 6 cm s−1 (Katsman et al. 2004). The southward propagating EKE can be interpreted as IR, however they only intermittently reach south of 58°N, where deep convection regularly occurs. The EKE variations in this region emerge instantaneously in the whole southern central LS, suggesting a local generation. While it is clear that IR that do reach the convection region in spring must in some form influence restratification, the model simulation does not support the idea that they represent a continuous source of restratification, acting in every year. Additionally, EKE off the WGC shows a maximum in winter and the IR need about half a year to reach the convection region (assuming a southward propagation of 5 cm s−1). Thus, the largest impact of IR on the southern central LS is to be expected in late summer, and only for years where IR actually reach the convective area. In fact, in SPG+ years, the MLD in the central LS is deeper compared to SPG− years (Fig. 7b). As in the sensitivity study above (Fig. 5), the effect of IR is restricted to north of 58°N.
b. Convective eddies
A second type of eddy in the LS is the convective eddy. CE are generated along the edge of the mixed patch in the central LS (e.g., Marshall and Schott 1999). Figure 9 illustrates this process exemplarily for the winter 2007/08. In November, before the onset of deep convection, at 989-m depth in the central LS is below 0.05, indicating no substantial mesoscale activity. At the end of March, however, small-scale vorticity structures > 0.2 in amplitude can be seen at the rim of the convective area, indicating the generation of instabilities, while inside the mixed patch, vorticity is still low. In May, much of this variability is still visible, albeit organized in slightly larger structures. Furthermore, the mesoscale structures have populated the whole central LS in May, including the formerly convective region that did not exhibit mesoscale variability in March. Over the summer, the small-scale variations get eradicated and only few mesoscale features remain in October.
The distribution of in the vertical yields further insights into the generation of CE (Fig. 10). In November, before the onset of convection, relative vorticity away from the boundaries and the bottom is low. At the end of March, convection is at its height as indicated by the Ertel potential vorticity  close to zero in large areas. Regarding the CE, the edge of the convective areas at depth is of particular interest. There are five locations in this snapshot (indicated by the black arrows), where large horizontal gradients of PV occur below 1000-m depth which can be identified as boundaries to convective patches. At all five locations, shows clear maximums > 0.1 at the bottom of the mixed layer that extend to more than 500 m deeper. Additionally, there is a local minimum in PV at 56°N, 800-m depth, seemingly disconnected from the mixed patches to the south and north, which is associated with a relative vorticity minimum and further develops into the anticyclonic eddy in Fig. 11.
Thus, some of these small-scale vorticity structures organize into coherent vortices, as observed by Lilly et al. (2003), although the horizontal and vertical resolutions of the model are probably not sufficient to capture the full range of instabilities at the rim of the mixed patch. The simulated CE have a diameter of about 30 km, which is on the upper end of the range given by observations (Lilly et al. 2003). Similarly to the observed lens-shaped vortices, the simulated velocity has a maximum of up to 0.2 m s−1 at depth, below 1000 m (Fig. 11a). The origin of the additional surface-intensified currents depicted in Fig. 11a remains unknown. Whether they are part of the original CE or rather stem from another mesoscale structure, merged with the CE, cannot be determined by careful analysis of the velocity field at the preceding weeks (not shown). The exact depth of the deep core depends on the depth of the mixed layer, as can be seen in Fig. 10. The core of the CE has low PV and carries temperature and salinity characteristics of the mixed patch. In most cases, and also in the exemplary CE in Fig. 11, the core is cold and fresh compared to the surroundings, that is, the waters outside the mixed patch.
1) Temporal variability
To further manifest the importance of local instabilities for EKE in the central LS, its temporal variability is investigated in the following. The seasonal cycle of EKE in the central LS (56.0°–57.5°N, 53.0°–56.0°W) exhibits April maximums of 20 cm2 s−2 for the upper (112–322 m) and 14 cm2 s−2 for the lower (382–1655 m) ocean (Fig. 12). During spring, summer, and early autumn, EKE then decreases gradually toward its minimum of 12 (upper) and 7 cm2 s−2 (lower) in October/November, followed by an increase throughout winter. In contrast, the conversion term for barotropic instability does not show a seasonal cycle and remains negative and small at all depth in all months. Also the conversion term for baroclinic instability in the upper ocean does not exhibit a clear seasonal cycle. Only BC in the deeper ocean has a pronounced maximum of 1.0 × 10−6 kg m−1 s−3 in April, in phase with the EKE, indicating that local baroclinic instability at depth is the major source of EKE at all depths in the central LS, that is, the convection region.
On longer time scales, a similar picture emerges (Fig. S4). The annual mean EKE in the central LS is significantly correlated to BC in the deeper ocean (correlation coefficient 0.86) and to the local winter heat fluxes (0.55) responsible for deep convection (averaged over January–March). Stronger buoyancy loss in winter leads to stronger convection, which causes larger lateral density gradients responsible for stronger baroclinic instabilities and thus EKE. As with the seasonal variability, BT at all depths is small and mostly negative. Thus baroclinic instabilities at depth are also responsible for changes in central LS EKE at interannual time scales. In some years (e.g., 1985 and 2003), the baroclinic instabilities in the upper ocean are significantly stronger than in the deeper ocean, and upper-ocean EKE consequently reacts with an increase. These years are years with weak convection and thus, the processes at work are most likely the same, just shifted toward the surface due to the reduced MLD. Indeed, the depth of the local maximum of BC varies depending on the depth of the mixed layer in each year (not shown).
At decadal time scales (the time series are smoothed with a boxcar window of 5 years length) EKE and BC at depth in the central LS are both significantly correlated to the SPG index (correlation coefficients 0.87 and 0.66, respectively) with 1–2 years lag, as the SPG is driven by the same large-scale atmospheric forcing as the local buoyancy fluxes in the LS. The IR and BCE, as discussed in section 3c, are both forced by the same atmospheric pattern as the deep convection thus some influences of IR and BCE on the decadal variations of EKE in the central LS cannot be ruled out completely. However, the correlation of local BC and EKE points to the conclusion, that the variability of EKE in the central LS is mainly controlled by local forcing at all time scales from seasonal to decadal.
2) Impact on restratification
Convective eddies are the most promising mechanism for rapid restratification of the mixed patch in spring. CE are formed during the convection process directly along the convection area and thus inevitably lead to an exchange of the waters outside and inside the mixed patch. This is supported by a comparison of the seasonal evolution of the potential density (with reference to 1000 m) anomaly averaged over 534–2710 m in the central LS between the simulation with 1/20° resolution and a short simulation with 1/12° (Fig. 13). At 1/12° (roughly 6 km in the central LS), CE are not well resolved and while instabilities still occur along the mixed patch (not shown), they are less energetic and the effect on restratification should be reduced. At the same time, IR are still well resolved at 1/12° and a comparison of these two simulations can be used to approximate the effect of CE on restratification. Figure 13 shows that restratification is indeed much slower in spring right after the convection period, when the horizontal resolution is reduced. With 1/20°, the positive winter density anomaly is halved just 1.5 months after the end of deep convection, compared to 4 months at 1/12°. Later in the year (from June/July on), the density decreases at a similar rate, independent of resolution. Toward the end of the year, the decrease in the 1/12° simulation is larger compared to 1/20°.
The analysis of the eddy buoyancy flux convergence confirms the role of CE in restratifying the central LS. Eddy buoyancy flux convergence is defined as , where u and υ are the zonal and meridional velocities, respectively, and b is the buoyancy. The prime denotes the 5-day mean deviation from a 100-day running mean and the denotes a long-term average. In the 1/20° simulation the eddy buoyancy flux convergence exhibits a clear seasonal cycle, with largest values in March and April, when the CE are expected to most actively restratify the convective area and lead to a rapid decrease in density. Later in the year, the eddy buoyancy flux convergence decreases and thus the restratification is slower. The behavior of the 1/12° simulation is substantially different. There is a strong convergence of buoyancy due to eddy fluxes in March, suggesting that there is some mesoscale activity in the central LS at 1/12° when convection is strongest. From April to July however, the eddy buoyancy flux convergence does not show a consistent behavior that could lead to restratification, with divergence occurring in April and June. This period of weak restratification is followed by large convergence of unknown origin in August and September, inducing the observed rapid decrease in density.
c. Boundary current eddies
Boundary current eddies are continuously generated along the boundary currents of the LS. While BCE exist along the WGC south of the region where IR are generated (not shown), this section is focused on the BCE of the LC, as they are more likely to influence convection due to their proximity to the area of deepest winter mixed layers (Fig. 4). As shown by Fig. 14 for a typical BCE at the western boundary, these eddies are strongly surface intensified and do not reach as deep as the IR from the WGC. The velocity maximum of 0.2 m s−1 is found at the surface. The core of the BCE is highly stratified, with a combination of fresh and cold water on top of warm and saline water, similar to the IR shed from the eastern boundary. However, the more saline water below ~200-m depth is still fresh compared to the surroundings.
To illustrate the seasonal cycle of BCE, the climatological seasonal evolution of EKE is shown with a focus on the LC in Fig. 15. Eddy kinetic energy has a maximum of 50–100 cm2 s−2 in the LC in winter, and decreases throughout spring and summer to a minimum of 10–20 cm2 s−2 and starts to increase again in autumn. This seasonal evolution is observed at all latitudes from 55° to 64°N in the LC and is linked to the varying strength of baroclinic instability (Eden and Böning 2002; Thomsen et al. 2014), which is driven by the seasonal cycle of the boundary current, mostly caused by large-scale variations of the wind field (Brandt et al. 2004; Han et al. 2008). Offshore of the LC, there is a distinct minimum of EKE in all seasons, in agreement with observations (Brandt et al. 2004). Farther toward the interior, EKE increases again due to the IR in the north and locally generated EKE farther south. Especially close to the convection region south of 58°N, the seasonal variations of EKE are confined to a narrow band along the shelf and are not observed to reach into the area of deep convection. Thus, BCE are most likely not contributing substantially to the restratification of the mixed patch in spring. Support for this hypothesis is found in Fig. 5. The reduction of IR in the sensitivity study yields an increased strength of the LC (not shown). This in turn causes more baroclinic instabilities and thus generation of EKE (light red area in Fig. 5a). However, an influence on the MLD of this increased EKE is only found in the far western part of the convection region. At this point it is unclear, whether this is due to the work of the BCE against preconditioning of this area for convection or due to an active restratification of the mixed layer during and after convection.
4. Summary and conclusions
This study investigates different sources of eddy kinetic energy in the Labrador Sea with the help of a 52-yr-long hindcast simulation of a 1/20° nested OGCM. Besides a good representation of the mean near-surface EKE and MLD in the LS (Fig. 4), the model also is able to simulate all three major types of eddies in the LS, Irminger Rings, convective eddies, and boundary current eddies, with characteristics close to observational accounts (Figs. 3, 11, and 14).
Irminger Rings are generated between 60° and 62°N in the WGC along the west coast of Greenland. The observed seasonal cycle of EKE is well reproduced by the model and the associated instability processes are resolved. In concert with the strength of the WGC, the seasonal cycle peaks in winter. The instability processes reach their maximum in January, while EKE, as a resulting quantity, is highest in February. Both, the transfer from mean available potential energy to eddy available potential energy associated with baroclinic instabilities at depth (below ~400 m) and the transfer of mean kinetic energy to eddy kinetic energy associated with barotropic instabilities between 112 and 322 m (Fig. 6), are important for the generation process of IR. However, no definite conclusion can be drawn with respect to the initializing mechanism. The vertical structure of the IR and the importance of the lateral boundary condition suggest a barotropic extraction from the WGC while the sensitivity simulations in section 3a indicate that a reduction in IR is to first order associated with changes in BC (Fig. 5) and the analysis of decadal variations of EKE in the WGC shows EKE to be correlated to upstream BC and local BT which indicates that BC is the triggering mechanism. A more detailed analysis of the current–topography interaction at Cape Desolation and the initially triggered instabilities is required to unravel the complex processes at work in this region.
On longer time scales, the generation of IR and the EKE in the northern central LS are connected to the large-scale circulation of the subpolar gyre, which appears to a large degree to be driven by atmospheric variability. A prevailing positive phase of the Arctic Oscillation facilitates a stronger subpolar gyre, and increases velocities in the boundary currents. In the WGC, this leads to a downstream (northwestward) shift of the occurrence of instabilities and also the generated EKE (Fig. 7). A sensitivity study with a switch from “no-slip” to “free-slip” as a lateral boundary condition along the west coast of Greenland results in a strong suppression of IR generation and demonstrates that their main effect is limited to the northern part of the LS (north of 58°N): under the same atmospheric forcing, a reduction in IR abundance thus alters the preconditioning and results in an extension of the area of deep convection to the north. In a realistic scenario where variations of the numbers and strength of IR shed from the WGC are related to temporal changes in the atmospheric forcing, the consequences are more complex. A stronger forcing (positive AO) leads, simultaneously, to stronger convection due to the local heat loss in the central LS, and to a preconditioning less favorable for convection due to more IR. As the convection is to first order related to local heat loss, the effect of the IR is only minor in a realistic setting.
Due to the limitation of the IR’s effect to the northern central LS, no substantial impact on rapid restratification of the mixed patch in spring is expected. While occasionally there are IR arriving in the central LS in spring that are possibly contributing to the restratification of the freshly convected water masses, this does not occur annually, so that IR are just an additional, intermittent factor in the restratification of the deep winter mixed layer in some years.
The main source of EKE in the central LS and thus of rapid restratification in spring are the convective eddies. Convective eddies emerge along the rim of the mixed patch, due to the baroclinic instabilities triggered by the large lateral buoyancy gradient across the edge of the convection area (Figs. 9 and 12). In agreement with the time of deepest mixed layers at the end of winter/beginning of spring, the seasonal cycle of EKE in the central LS peaks in April. As derived from theoretical considerations (Jones and Marshall 1997; Dewar 2002), the strongest baroclinic instabilities occur near the bottom of the mixed layer. Another sensitivity study with a 1/12° configuration of the same model supports the CE’s role in rapid restratification of the mixed patch. A resolution of 1/12° is not resolving CE and consequently, the restratification is significantly slower (Fig. 13).
The interannual to decadal variability of CE is controlled by the atmospheric forcing. Large buoyancy loss at the surface triggers deep convection and thus the baroclinic instabilities at depth. The local air–sea fluxes and the resulting convection and EKE generation are correlated to the large-scale atmospheric state as expressed, for example, by the AO index.
The third class of eddies in the LS are the boundary current eddies. They are not found to propagate far from the boundary currents, where they are generated by baroclinic instabilities. In the vicinity of the Labrador Current the BCE’s impact is restricted to the far western part of the convection region. Furthermore, it could be overestimated by the model because deep convection appears to be simulated unrealistically close to the boundary current. In an earlier configuration of the model used here, Handmann et al. (2018) find the western boundary current in the LS to be too barotropic, compared to observations. We speculate, that a too barotropic boundary current is likely to produce insufficiently strong baroclinic instabilities and hence fewer BCE. A weak eddy field along the western boundary then partially fails to act against the loss of stratification, facilitating strong convection close to the shelf. The complex interplay between the LC, BCE, and deep convection needs further, detailed investigation.
In conclusion, the Irminger Rings in the northern central Labrador Sea and to a lesser extent the boundary current eddies along the Labrador Current appear as vital components acting against the preconditioning for deep convection. In the mean, they limit the area of deep convection to the southern part of the central LS. During periods of favorable atmospheric conditions, a possible increase of deep convection in the northern central LS due to local air–sea buoyancy loss is counteracted by a strengthened limiting, stratifying effect of IR and BCE. The convective eddies in the central LS, tightly linked to the strength of deep convection, are the main driver of the rapid restratification in spring.
As suggested by this study, a critical aspect of ocean model formulations to successfully simulate the Labrador Sea’s mesoscale eddy field and its impact on deep convection is the combination of high resolution and a realistic location of the area of deep convection. The implementation of an adequate lateral boundary condition (no-slip in this study) in the WGC region is additionally necessary for a correct representation of the Irminger Rings. Idealized model configurations set up to study mesoscale variability in the LS generally fulfill the requirement for a high resolution and an adequate lateral boundary condition. However, deep convection is often represented by a circular mixed patch located rather centrally in the LS (e.g., Katsman et al. 2004). With the area of deep convection located too close to the westward path of the IR, their effect on restratification is likely to be overestimated. In OGCMs with realistic geometry and forcing, the horizontal resolution often is the limiting factor. Generally, a low resolution also implies a misrepresentation of the mixed patch’s location and lateral extent. Without explicitly resolving the IR and BCE, the area preconditioned for convection in the central LS is too large. Furthermore, CE are not resolved and thus do not counteract deep convection or restratify the mixed patch in spring. The combined effect is an overestimation of the extent and depth of deep convection (e.g., Rattan et al. 2010) in most OGCMs with coarse resolutions. Simulations at 1/20° still exhibit biases in the hydrographic properties of the Labrador Sea (Handmann et al. 2018) and additionally are not able to resolve submesoscale processes which could impact the stratification (Fox-Kemper and Ferrari 2008; Su et al. 2018) and convection (Callies and Ferrari 2018). However, the presented study shows 1/20° to be sufficient to simulate most of the mesoscale dynamical features.
We are grateful to three anonymous reviewers for their helpful comments on the manuscript. The model system has been developed by the ocean modeling group at GEOMAR in the framework of the DRAKKAR collaboration. The model simulations were supported by the cooperative programme “RACE - Regional Atlantic Circulation and Global Change” (BMBF Grant 03F0651B). All integrations were performed at the North-German Supercomputing Alliance (HLRN). Data shown in this paper are available at http://data.geomar.de.
Denotes content that is immediately available upon publication as open access.