This study primarily focuses on interannual changes in CLS SLA and SHA during 1992-2025, a period marked by the most considerable recorded changes in heat content, freshwater accumulation, and winter mixing depth, and coinciding with the satellite altimetry era, and in CLS SHA starting from 1948. By comparing the long-term trends, multiyear cycles, and extremes of SLA, SHA (as well as its thermosteric and halosteric components), and associated hydrographic properties, we investigate the underlying origins and drivers of these changes by sequentially addressing the first-to-fifth questions. The SLA and SHA changes concurrently resolved since 1992 are then placed in the context of SHA variability over the period of 1948-2025. The key goal of this retrospective analysis is to determine how recent conditions differ from those of the past, particularly with respect to the third-to-fifth questions. Furthermore, we assess the contribution of certain vertical layers to the total steric height changes, with particular emphasis on the role of the deep ocean, as asked in the sixth question. Finally, we discuss how air-sea heat exchange and Arctic sea ice meltwater contribute to regional heat and freshwater content changes, and how this knowledge can help us to accurately model thermosteric and halosteric variations, and hence sea level trends, thus addressing the last two questions.
Our SHA analysis is exclusively based on the original (i.e., unaveraged and ungridded spatially and temporally) Argo float and ship-based hydrographic measurements. Analogously, the along-track altimetry data is used to construct the CLS SLA time series. As part of cross-validation, the CLS SLA time series derived from along-track and gridded altimetry (Supplementary Figs. 6 and 7) show similar patterns across seasonal, interannual, and multidecadal timescales throughout the record.
The raised research questions are addressed by studying the interannual changes, trends, and extremes of SLA, SHA, and thermosteric and halosteric heights. Prior to analysis, the respective climatological seasonal cycles are subtracted from individual measurements, which are then filtered, smoothe,d and bin-averaged, as described in the "Methods" section.
Since we analyze anomalies, computed as deviations from the respective climatological seasonal cycles, we show these cycles in Figs. 3 and 4, and Supplementary Figs. 6-10 to facilitate comparisons and summations with the anomaly-derived values (Supplementary Figs. 9 and 10 have seasonal cycles added to anomalies).
The distinct seasonal cycles of the detrended SLA and SHA have considerably different magnitudes (Fig. 3). This difference is highly insightful for the regional freshwater and ocean dynamics, as it explicitly points at a summer mass gain, supported by satellite gravimetric data (Fig. 3, "Methods" section), further discussed in the "Remaining challenges and future steps" subsection.
The thermosteric component dominates the seasonality of SHA (Fig. 4 and Supplementary Fig. 8). The specific volume, temperature, and salinity, and steric, thermosteric, and halosteric height seasonal cycles examined across the 10-900 dbar pressure range (Fig. 4) reveal the pressure-dependence of the respective seasonal amplitudes and phases. While the thermosteric amplitude monotonically decreases with pressure, the overall smaller halosteric amplitude reaches its maximum within the 100-200 dbar range. While the temperature and thermosteric height phases monotonically increase with depth, the salinity and halosteric phases reverse between 10 dbar and 150 dbar.
In the CLS, SLA has varied significantly over the satellite altimetry era (1992-2025), reaching a record low in 1994, a notable local minimum in 2017, and a record high in 2025 (Fig. 5). The full-depth SHA, which shows similar interannual variability, is constructed as explained below.
Since the core Argo profiling depth is limited to 2000 dbar, there are no sufficient measurements to resolve sub-annual hydrographic changes in the deeper layer. Fortunately, the CLS water column below 1500 dbar is dominated by interannual variability (and the deeper the stronger, Fig. 2). Therefore, to sufficiently cover the entire water column to resolve the prevailing timescales, we analyze the 1900-3300 dbar layer separately from the overlaying, 10-1900 dbar layer. This approach facilitates a focused assessment of the deep-layer SHA interannual variability and long-term trends, providing a baseline study for measuring and understanding the contribution of the deep layer to the steric height changes. To represent the full-depth, 10-3300 dbar, water column in our research, we combine the upper, 10-1900 dbar, and deep, 1900-3300 dbar, layer-based estimates.
While having measurably different seasonal cycles (Fig. 3), SLA and full-depth SHA show strikingly similar interannual changes, trends, and extremes. This similarity is evidenced by a high, 0.985, correlation of the yearly averaged SLA and SHA, hereafter called YASLA and YASHA, respectively and a low, 0.67 cm, standard deviation of their difference. Detrending of YASLA and YASHA does not affect the strength of their similarity, giving, respectively, 0.960 and 0.67 cm for the noted metrics. Furthermore, YASHA accounts for approximately 97% of the interannual variability of YASLA before detrending, and 92% after detrending, with individual values within 1.4 cm of each other in 97% of all cases, both before and after detrending. In 2025, both YASLA and YASHA reached record highs, signified by the respective squares in Fig. 5 standing 18.4 cm and 15.0 cm higher than 31 years earlier, at the dawn of the satellite altimetry era. The 34-year, 1992-2025, trends of YASLA and YASHA, 0.308 ± 0.007 and 0.300 ± 0.007 cm/year, respectively, which small difference of ~2.6% is attributed to the water-column mass change. Coincidentally, a 30-year period is commonly regarded as a baseline for climatological normal, signifying this study as the first-of-its-kind climatological assessment of sea level changes in the CLS.
The strong alignment of SLA and SHA underscores the fact that sea level variability in the CLS is essentially (i.e., by 97%) steric, and is driven primarily by full-depth variations of temperature and salinity. This answers the first key question.
Figure 2 provides insight into the role of the multiyear convective cycles in defining and shaping the SHA and, consequently, SLA trends and extremes (Fig. 5). Thermosteric and halosteric components of SHA, shown in Fig. 6, quantify the respective roles of the full-depth temperature and salinity changes driven by both heat and freshwater fluxes and convection.
The answer to the first key question in the previous subsection is supported by the correlation between the altimetric sea level (SLA and YASLA) and the total steric height for the entire water column (SHA, YASHA), and associated statistics. Since all analyzed years, except 1992-1995, the direct effect of winter convection is restricted to the top 2000 db (Fig. 2), in this subsection, we provide analogous statistics for the 10-1900 dbar layer steric height instead of the full-depth one. The correlation and standard deviation of the difference of YASLA and 10-1900 dbar YASHA are 0.961 and 1.08 cm before detrending, and 0.932 and 0.95 cm after detrending, respectively. YASHA explains 92% of the interannual variability of YASLA before detrending and 87% after detrending.
The aforementioned convective intensification-to-relaxation phase change, evident in Fig. 2, had a direct impact on YASHA and YASLA (Figs. 5 and 6 for full-depth YASHA, and Supplementary Figs. 11 and 12 for 10-1900 dbar YASHA), manifesting in their respective trend changeovers around 2018. Stronger convective mixing events produce denser thicker pycnostads (Fig. 2, bottom panel), and, in turn, lower SHA and SLA. In contrast, weaker convections produce lighter, thinner pycnostads, stronger vertical stratification, and, in turn, higher SHA and SLA. This explains why the main driver of convection and deep-ocean cooling -- net winter surface heat loss -- is also in control of YASHA and YASLA through ocean cooling dominating mixed layer density, convection depth, and vertical stratification.
Summarizing the material presented in this subsection, we conclude that recurrent winter convection, with its 2010-2011 and 2021, 2023-2025 lows, and 1994 (discussed later) and 2018 (Fig. 2) highs, through the thermosteric height, dominates YASHA (Fig. 6). This answers the second key question.
On the other hand, driven by the massive upper-layer freshening, the effect of the year-to-year halosteric height changes (our next focus) has recently switched from moderating (i.e., reducing) to reinforcing (i.e., amplifying) the SHA and, consequently, SLA upward trends, earlier influencing the timing of the SHA/SLA fall-to-rise reversal, shifting it from 2018 to a slightly earlier date.
The yearly averaged full-depth halosteric height resides within about half the thermosteric height range (Fig. 6, squares). Besides, the interannual halosteric and thermosteric height changes are often opposite in sign, with the former tending to counterbalance the latter. Can it be, therefore, concluded that the sea level changes driven by thermal expansion are generally half-compensated by concurrent haline contraction? Apparently, it can. Indeed, according to Fig. 6, nearly every year from 1990 to 2015, the full-depth halosteric component consistently counteracts the thermosteric component, offsetting year-to-year thermosteric height changes by about a half. Such persistence of the halosteric-thermosteric counterbalancing showing through 2015 might imply that the partial compensation of thermal expansion by haline contraction would recur, if the situation had not radically changed after that year.
In 2016, the halosteric contributions to the year-to-year steric height changes switched from counterbalancing to reinforcing the respective thermosteric contributions. This answers the third key question. This transition arose as the halosteric height maintained a positive trend throughout the 2011-2025 period, while the thermosteric height rebounded from its decline in response to the 2012-2025 convective cycle's phase change (Fig. 2). The physical processes responsible for the recent transition from halosteric-thermosteric counterbalancing (before 2016) to their allying and covarying (2016-2025) become evident from Fig. 2, Supplementary Figs. 1-4, and the published analysis of the relevant signals and their causes recapped here. The CLS water column steadily cooled as convection deepened between 2011 and 2018. Then, these trends reversed as convection entered a relaxation phase. In turn, salinity, unlike temperature (Fig. 2), and thus the halosteric component, unlike the thermosteric component (Fig. 6), exhibited distinct steady multiyear trends throughout the entire 2011-2025 period, with upper-layer freshening being interrupted only briefly in 2016 and 2017 without disturbing the overall trends. This answers the fourth key question.
As shown in Fig. 2, and Supplementary Figs. 1-4, and is further illustrated and explained in the last, "Linking halosteric height changes to the extreme Arctic sea ice losses", subsection of the "Results" section, in 2025, the upper (15-100 m) and upper intermediate (300-500 m) layer salinities reached record lows, directly caused by advection of anomalous quantities of freshwater produced by extreme Arctic sea ice melt a few years earlier, while the full-depth halosteric height reached its maximum somewhat earlier, during 2022-2023 (this shift is explained in the noted subsection). Driven by the unbalanced freshwater inflow into the Labrador Sea, the halosteric height increased by approximately 5 cm between 2011 and 2023. In contrast, the thermosteric height was lower in 2023 than in 2011, as post-2018 warming could not fully compensate (and thus offset) the 2011-2018 cooling. Overall, the contribution of salinity has amplified that of temperature after winter convection entered the relaxation phase, making the sea level to be higher in 2023 than in 2011 (Figs. 5 and 6). The low-salinity state of the upper 800 m layer acquired in 2023 persisted through most of the two following years, 2024 and 2025, keeping the halosteric height relatively high, although it slightly declined toward the record's end (Fig. 6). Concurrently, between 2023 and 2025, the further shoaling of winter convection to record shallow and sea warming sustained the upward trends of SHA and SLA, which reached record highs in 2025. The sea level rise between 2023 and 2025 was mostly due to the thermosteric component.
Overall, in contrast to moderating (i.e., reducing or partially offsetting) the density responses to both positive and negative temperature changes prior to 2016, the salinity changes of the subsequent years amplified these responses, aiding first the convection-driven downward and then upward thermosteric trends. By reinforcing instead of compensating the thermosteric height changes in the latter years, the halosteric component boosted the temperature-driven sea level rise, creating a marked extreme. Indeed, while through 2023-2025 the thermosteric height stayed 3.6-1.0 cm below its record high that was registered in 2011 (Fig. 6, red squares), the accentuated halosteric-thermosteric relationship changeover from counterdirectional to codirectional played a critical role in the CLS sea level rise to a record high in 2025. This answers the fifth key question.
Moreover, record-breaking rates of the YASHA and YASLA changes over a sliding eight-year interval also fall on 2017-2024. These rates surpass even those for the 1994-2001 period brought up by a rapid recovery from record cold, fresh and dense conditions, alongside record low YASHA and YASLA levels (Figs. 5 and 6). A fundamental difference between these two periods lies in the behavior of the halosteric height: in 1994-2001, it counterbalanced about half of the thermosteric height change, whereas in 2017-2023, it reinforced the other, accelerating the total SHA and hence SLA increase.
The role of the deep layer in the recent SLA trend and extreme becomes evident when comparing the full-depth YASHA (Fig. 5) with that of the upper 1900 dbar layer (Supplementary Figs. 11 and 12). While the 2012 and 2025 upper-layer YASHA were not significantly different, the full-depth YASHA was at a record high in 2025. After two decades, 1992-2011, of relative stability with rather small fluctuations (± 0.6 cm) deep-layer YASHA underwent a positive trend of ~0.22 cm/year, implying a cumulative deep SHA increase of ~2.4 cm throughout 2012-2023 (Figs. 5 and 6). This change exceeds the ~0.9 cm increase in full-depth YASHA over the same period, emphasizing the critical role of the deep layer in raising full-depth SHA, and, consequently, SLA, to their unprecedented highs reached in 2025.
The origin of the positive 2012-2023 deep-layer YASHA trend is revealed through a decomposition of the deep-layer YASHA series into the respective thermosteric and halosteric components and their close co-examination (Figs. 6 and 7). Between 1990 and 2011, these components recurrently rebalanced each other, effectively suppressing development of any significant decadal trend in deep-layer YASHA (Figs. 5 and 6). This delicate equilibrium was irreversibly disrupted in 2012, causing a shift in the thermohaline balance that persisted for eleven consecutive years. Throughout this period, the stronger variability of the deep thermosteric height compared to that of the deep halosteric height dominated the deep-layer YASHA; consequently, a positive trend. This answers the sixth key question.
Recognizing that the halosteric reinforcement of the thermosteric height changes occurring over the period of 2016-2023 played a critical role in the recent sea level rise to a record high (Fig. 6), we extend our analysis back to 1948 in order to find out if similar events ever occurred before 1996.
The ship and Argo float-based hydrographic measurements collected in the CLS from 1948 to 2023 have provided insights into multidecadal cycles of winter convection, water-column cooling and warming, freshening, and salinization. To investigate the impacts of these changes on the regional steric height components, we, upon quality-checking and editing all in situ observations, construct for each year of the 1948-2025 period, meeting the data sufficiency requirement, composite vertical profiles. The steps of this data synthesis are outlined in the "Data Sources" and "Methods" sections, and our earlier publications. The resulting full-depth annual profiles, compiled in Supplementary Fig. 1, fully and accurately capture all major hydrographic developments in the Labrador Sea, including the multiyear convective cycles, since the late 1940s. To optimize quality and comparability of the in situ observations entering the underlying annual state-of-the-ocean profile making process (and hence Supplementary Figs. 1 and Fig. 7), here, the Argo and, when available, Deep Argo profiles only fill in for missing or insufficient ship-based observations (all Deep Argo float profiles have been checked and calibrated to achieve 0.002 °C temperature and 0.002 salinity accuracies as explained in the "Supplementary Information" section of the mentioned study), like in 2017 and 2021, while being left out otherwise. The annual profiles are then used to derive the YASHA time series and its components across pressure levels, spaced at 5 dbar intervals from 200 to 3000 dbar and referenced to 3300 dbar, as shown in Fig. 7 (the 3000-3300 dbar layer is not displayed as changes are weak there).
As seen in Figs. 2 and 4, the seasonal cycle dominates temporal variability in the top 200 dbar layer (e.g., the seasonal cycle accounts for >92% of the total temperature variance at 10 dbar). Furthermore, even with the regular seasonal cycle subtracted from infrequent pre-Argo (1948-2002) ship-based measurements, undersampled irregular seasonal variations may still be present there, contributing spurious signals to yearly averaged values. The upper 200 bar (meter) layer data, prone to residual seasonal aliasing, have been excluded from the pre-2003 compilations shown in Fig. 7. This exclusion has an insignificant effect on the full-depth, 10-3300 dbar, steric, thermosteric, and halosteric heights. The corresponding 10-3300 dbar and 200-3300 dbar height differences, particularly for the steric height, are much smaller than the decadal changes (Fig. 7, top), justifying our choice of the 200-3300 dbar heights for analyzing longer-term variations in CLS steric sea level. Furthermore, while the exclusion of the upper 200 dbar layer slightly reduces the range of height changes, the patterns remain unaffected.
The 200-3300 dbar steric, thermosteric, and halosteric height extremes for the 1948-2025 period are highlighted with triangles in Fig. 7 (top). The timing of these events suggests the following:
[1] A thermosteric high and a halosteric low occurred between 1970 and 1971.
[2] Conversely, a thermosteric low and a halosteric high were recorded in 1994.
[3] The total steric height reached its absolute minimum concurrent with the lowest thermosteric and highest halosteric values, whereas its absolute maximum, achieved in 2024, is comparatively less pronounced an extreme in either component.
[4] While all major hydrographic events (e.g., anomalies, trends, mixing events), occurring in the CLS between 1948 and 2015 had their sustained thermosteric signals counterbalanced and significantly compensated by the respective halosteric signals, the halosteric and thermosteric changes during the 2015-2023 period were of the same sign and had comparable magnitudes, reinforcing each other in SHA and SLA. Notably, the positive coupling of the post-2015 thermosteric and halosteric height changes, coinciding with the most substantial freshening of the upper 700 dbar layer (Fig. 2 and Supplementary Figs. 1-4), challenges the conventional vision of CLS steric height changes based on the assumption of counteraction of halosteric and thermosteric changes and trends through density compensation. It also emphasizes the role of changes in horizontal advection of heat, salt, and, especially, freshwater in accelerating or moderating the sea level rise.
Between 1953 and 1956, we see another instance of the two steric components changing in the same direction. However, taking into consideration the short duration of that event, the sparseness of water sample observations, and the varying accuracy and resolution of the hydrographic measurements of that time, we refrain from making any conclusive statement about the possibility of a halosteric effect reversal during that time (the mentioned sparseness, low quality and low quantity of salinity measurements is reflected in Supplementary Fig. 1in the patchy vertical strips between 1948 and 1975. In contrast, the period of 1990-2019 features uniformly high quality and sufficiency of hydrographic observations, cataloged in Supplementary Table).
Figure 7 also provides detailed insight into the intermediate and deep layers' contributions to YASHA variability. YASHA time series constructed for pressure levels spaced at 5 dbar are brought together in the second panel from the top. This compilation clearly shows where in the water column, when, and how fast each trend developed, weakened, and reversed, and which layer shaped it the most. In most cases, the interannual changes and trends tend to reverse anywhere between the pressure levels of 1250 dbar and 2000 dbar. These reversals, regardless of their exact vertical positions, explain why the trends observed at 200 dbar and 2250 dbar are so profoundly different. Over the 1990-2021 period, the 200 dbar and 2250 dbar YASHA series displayed opposing trends, each punctuated by short-term reversals in the opposite direction. This inverse symmetry of YASHA changes between the upper and deep layers suggests a counterbalancing effect, where upper-layer changes are offset by those in the deep layer, or vice versa. What is particularly important to our study is that in 2022-2024, unlike the previous years, YASHA were consistent across all pressure levels, indicating a positive contribution from the deep layers to the 200-3300 dbar and full-depth steric heights.
Further insights into the impact of the deep layer on upper ocean SHA changes are obtained by comparing thermosteric and (reversely color-coded for easier comparison) halosteric component variations with depth (Fig. 7). Below 1500 dbar, the two components largely counterbalance each other, although the resulting compensation is not full at all times. For instance, the halosteric component dominated deep-layer residuals from 1991 to 1999, while the thermosteric component did so from 2005 to 2016. In 2022-2024, both deep halosteric and thermosteric height anomalies were positive, amplifying the upper-layer signal significantly.
Unlike its deeper counterpart, the water column's segment above ~1500 dbar, throughout its extent, is typically strongly affected by the thermosteric component as the upward halosteric gain remains relatively weak. This promotes faster accumulation of thermosteric anomalous signals toward the surface and their prevalence in full-depth YASHA. Strikingly, here again, 2022-2024 were exceptional. In these years, the upward halosteric gain matched the thermosteric one.
The provided explanation of the deep and intermediate layer contributions to the exceptional CLS sea level rise during 2022-2025 underscores the uniqueness of the event and raises its importance for the subpolar and larger North Atlantic domains, as local deep and full-depth regime shifts are likely to affect broader-scale ocean dynamics and exchanges.
The sea level of the central Labrador Sea (CLS) reached a 78-year record high in 2025. The rapid sea level rise that led to this event was caused by the joint action of mild winters, warm summers, sustained shutdown of deep convection, exceptional upper ocean freshening, deep-ocean halosteric-thermosteric balance shift, and full water-column mass gain. While both deep, 1900-3300 dbar, layer steric height and full-column mass changes contribute to the long-term sea level trend, the interannual-to-decadal variability is predominantly shaped by the upper, 10-1900 dbar, layer thermosteric and halosteric components. Our understanding of the interannual-to-decadal sea level changes is further advanced following the approach presented here for the reconstruction and prediction of the two steric height components. We first demonstrate how realistically detailed and accurate our pseudo thermosteric heights are simulated using the atmospheric forcing data (Fig. 8), and then discuss the predictability of the halosteric height.
The thermosteric component dominates both seasonal (Fig. 4 and Supplementary Fig. 8) and interannual (Figs. 6 and 7) variabilities of SHA, and, consequently, those of SLA (Figs. 3 and 5). The seasonal patterns of temperature and thermosteric height change their shapes with pressure (depth) (Fig. 4), revealing two distinct vital signals with unique vertical penetration routines, pressure-dependent time lags, and appearances -- a deeper winter cooling signal and a shallower summer warming signal. The winter cooling is regulated by Winter Surface Heat Loss (WSHL), which was thoroughly analyzed in ref. and is calculated by integrating all components of surface heat budget over an individually-defined cooling period. In contrast, Summer Surface Heat Gain (SSHG) is calculated for the period when the accumulated net heat flux is directed into the sea. Another characteristic of seasonal warming, used in our work, is Summer Heat Peak (SHP). SHP is averaged over a fixed-length time interval (e.g., 15, 21, 31 days) centered on a daily surface heat gain peak. The surface heat exchange characteristics are further detailed in the "Seasonal air-sea heat exchange metrics" subsection and Supplementary Fig. 13 caption.
To empirically model the observed yearly averaged thermosteric height changes we make three key assumptions:
[1] air-sea heat exchange is the leading factor controlling the thermosteric height,
[2] the interannual changes of WSHL and SSHG have different in magnitude and differently lagged effects on the thermosteric height, and need to be assessed separately, and
[3] the residual impacts of previous cooling and warming events (preconditioning) can be approximated by some asymmetric low-pass filtering (e.g., autoregressive; or, here, with left-side half-triangularly weighted window).
We reconstruct the yearly averaged thermosteric height series by optimally low-pass filtering, scaling, and merging WSHL and SSHG or, alternatively, SHP (SSHG | SHP). All sought parameters are found through iterative approximations aiming to minimize either squared or absolute deviations from the thermosteric heights. The atmospheric variables (e.g., WSHL, SSHG, SHP, and NAO) are low-pass filtered using a left-sided triangular window with weights decreasing linearly backward from the central point and equal zero forward. The best thermosteric reconstruction is achieved with different WSHL and SSHG | SHP filter window sizes.
The reconstructed thermosteric heights closely approach their targets, especially after 2000 (Fig. 8). The reconstruction captures 96% of the observed variance. The strong correlation (0.98) underscores the robustness of the model in simulating thermosteric contributions to sea level changes over the past three and a half decades, accurately tracking both the overall trend and individual cycles, with observed changes replicated with a 1.0 cm accuracy in 27 out of 33 years (~ 83%). The overall level of fit that is achieved by using optimally filtered and weighted (scaled) CLS WSHL and SSHG | SHP time series supports our assumptions and demonstrates the effectiveness of the proposed empirical model, making it suitable for further investigation and interpretation of both atmospheric forcing and signal transfer. This answers the seventh key question.
The difference between the optimally-fitted WSHL and SSHG | SHP left-side triangularly weighted window sizes of 7 and 13 years, respectively, emphasizes the different roles of the previous winter conditions retained by the water column, known as convective preconditioning, and the cumulative effect of summer warming. Indeed, while the winter cooling and mixing are uniquely strong and deep in the Labrador Sea, leaving an immediate trace over a thick layer, the direct effect of summer warming is not that deep. Additionally, anomalous winter cooling situations usually have smaller spatial scales than ocean-wide summer heat waves (in preparation). Therefore, it must have taken a longer time and, possibly, a larger region of influence for SSHG to achieve a sizable effect on YASHA and YASLA. By accumulating its signal over a broader domain, SSHG spreads its influence on the thermosteric height over a longer time, hence a longer memory of SSHG changes in YASHA. Yet, while both WSHL and SSHG take turns driving thermosteric height, and the SSHG changes are more influential there on the longer timescales, WSHL dominates in this linkage as a whole. Based on our results, the WSHL-thermosteric interaction is performed through convection, hence the leading role of convective cycles in sea level variability in the CLS domain.
Notably, while the low-pass filtered winter (DJFM) NAO index shows some limited agreement with the yearly averaged thermosteric height, it lacks the precision and detail captured by the presented heat-based model, highlighting the superior predictive accuracy of the heat-based reconstruction approach.
The significant reductions in Arctic sea ice in 2007, 2012, and 2019-2020 (Figs. 8 and 9; and Supplementary Fig. 2) were each followed by pronounced freshening in the upper layer of the Labrador Sea approximately two years later (Fig. 9). This freshening could also be influenced by a recent shift in the Beaufort Gyre's regime from freshwater accumulation to stabilization, with a potential release phase that may contribute to the latest CLS freshening event. In contrast, the Greenland freshwater flux anomaly, which changes more gradually over time, is unlikely to drive rapid freshening events in the CLS. For this reason, we focus our analysis on the impacts of extreme Arctic sea ice losses as the primary driver of intermittent CLS halosteric height changes.
The arrows in Fig. 8 point from the 2007, 2012, and 2020 extreme winter-to-summer Arctic sea ice reductions to the 2009, 2015, and 2022 CLS halosteric height maxima lagged by 2-3 years. The linkages between the extreme points of these two key ocean state variables are inferred from the timing of the respective local salinity minima in the upper (15-100 m) and upper intermediate (300-500 m) layers of the CLS (Fig. 9). The upper layer salinity minima of 2008, 2013, 2022 and 2025, and, lagged by about a year, corresponding upper intermediate layer salinity minima of 2009, 2015, 2023 and 2025, can be traced back to the extreme Arctic sea ice losses of 2007, 2012, 2020 and 2024. The latter are inferred from increased year-to-year Arctic sea ice extent and volume reductions, similarly showing in the time series of August-October means and winter-to-summer reductions of both sea ice metric (Fig. 9). Whenever this happens, an increased quantity of meltwater is passed to the Arctic outflow to be delivered to the Labrador Sea through the Davis Strait (and, possibly later, through the Denmark Strait) to reduce salinity of its upper layer.
In contrast to the 2007, 2012, and 2020 extreme releases of Arctic sea ice meltwater being closely followed by the respective CLS halosteric height maxima, the meltwater release of 2024 was trailed by a halosteric height decline. Even though the halosteric height may respond with a rise later on, 2025 stands out as the year when the top 500 m layer of the Labrador Sea experienced the largest (most massive) freshening in the instrumental oceanographic record (Figs. 2 and 9; and Supplementary Figs. 1-4). According to the salinity sections (Supplementary Figs. 3 and 4), this, as well as earlier freshening events, entered the CLS from the Labrador slope side.
The seeming contradiction between the anomalous sea freshening and halosteric height decline between 2014 and 2015 is resolved by considering the broader impacts of the increased Arctic meltwater influx on the full-depth water column. Indeed, through inhibiting winter convection and consequently reducing its depth to 500 m, the sustained freshening of the upper layer led to salinification of the deep intermediate (800-1500 m) layer between 2022 and 2025 (Fig. 9). As winter convection shoaled, the amount of freshwater that reached the deep intermediate layer had reduced, shifting the freshwater balance toward the saltier deep intermediate waters entering the CLS from the neighboring SPNA basins. This resulted in a sustained salinification of the deep intermediate layer, which during 2023-2025 competed with the freshening of the upper and upper intermediate layer over the halosteric height dominance. Obviously, the halosteric height decline seen during 2014-2015 explicitly points to the last-year winner of this everlasting competition.
However, one might question the lack of substantial annual reductions in Arctic sea ice preceding the 1996-2005 period, during which CLS halosteric height was decreasing alongside rising salinity. This period coincided with exceptionally strong and sustained deep convection from the late 1980s to the mid-1990s, which infused the CLS water column with an estimated 7-meter freshwater equivalent. As the sea entered a convective relaxation phase in 1996, freshwater began to discharge from its intermediate layer into the broader North Atlantic, resulting in a gradual reduction in halosteric height during this relaxation phase.