Weddell Sea open-ocean polynyas have been observed to occasionally release heat from the deep ocean to the atmosphere, indicating that their sporadic appearances may be an important feature of high-latitude atmosphere–ocean variability. Yet, observations of the phenomenon are sparse and many standard-resolution models represent these features poorly, if at all. We use a fully coupled, synoptic-scale preindustrial control simulation of the Energy Exascale Earth System Model (E3SMv0-HR) to effectively simulate open-ocean polynyas and investigate their role in the climate system. Our approach employs statistical tests of Granger causality to diagnose local and remote drivers of, and responses to, polynya heat loss on interannual to decadal time scales. First, we find that polynya heat loss Granger causes a persistent increase in surface air temperature over the Weddell Sea, strengthening the local cyclonic wind circulation. Along with responding to polynyas, atmospheric conditions also facilitate their development. When the Southern Ocean experiences a rapid poleward shift in the circumpolar westerlies following a prolonged negative phase of the southern annular mode (SAM), Weddell Sea salinity increases, promoting density destratification and convection in the water column. Finally, we find that the reduction of surface heat fluxes during periods of full ice cover is not fully compensated by ocean heat transport into the high latitudes. This imbalance leads to a buildup of ocean heat content that supplies polynya heat loss. These results disentangle the complex, coupled climate processes that both enable the polynya’s existence and respond to it, providing insights to improve the representation of these highly episodic sea ice features in climate models.
Shortly after passive microwave satellite observations began in the 1970s, a large region of open ocean appeared for three consecutive winters in a typically ice-covered region of the Weddell Sea (Carsey 1980). These anomalous ice-free regions enclosed by the winter sea ice pack have since come to be known as open-ocean polynyas. The initial observation in the Weddell Sea ranged in size from 2 to 3 × 105 km2 during 1974–76 and was maintained by a heat supply from sustained ocean convection reaching depths of up to 3000 m (Gordon 1982). Reanalysis-based reconstructions of the Weddell polynyas have shown that the tapping of the deep-ocean heat reservoir, combined with intense air–sea interaction over the anomalously ice-free ocean, delivered large quantities of heat to the atmosphere (Moore et al. 2002). More recent wintertime observations of smaller polynyas over the Maud Rise seamount during 2016 and 2017 illustrate that these sea ice features can be influenced by high-latitude atmospheric circulation anomalies and strong polar cyclones over the Weddell Sea (Francis et al. 2018; Jena et al. 2019; Cheon and Gordon 2019; Campbell et al. 2019). While both the 1974–76 and 2016–17 events demonstrate connections between deep convection and interannual high-latitude climate variability, observations of open-ocean polynyas in the Southern Hemisphere high latitudes are still sparse. Important questions remain regarding how polynya heat loss may have influenced the atmosphere before satellite observations began, as well as how their seemingly rare appearances may relate to current and future high-latitude climate.
Recent multidecadal trends in the Southern Ocean are in stark opposition to the precipitous sea ice decline occurring simultaneously in the Arctic. From 1979 to 2011, Antarctic sea ice extent exhibited a positive trend and sea surface temperatures (SSTs) around the continent decreased (Fan et al. 2014). The surface cooling trend reversed in 2016, associated with anomalous atmosphere–ocean variability (Meehl et al. 2019; Wang et al. 2019). This trend reversal seems to coincide with the reappearance of an open-ocean polynya in winters of 2016 and 2017, although a clear, physical link between these phenomena has not yet been established. The warming response to anthropogenic greenhouse gas forcing is expected to be delayed in the near-surface Southern Ocean south of the Antarctic Circumpolar Current (Armour et al. 2016). However, the mechanisms underlying the observed SST cooling trend remain unclear; in the region, many historical hindcast simulations feature a warming bias relative to observations during this recent time period (Kostov et al. 2018). This hemispheric asymmetry necessitates a better understanding of Southern Ocean variability, including variability associated with polynya heat loss, that could be operating alongside global anthropogenic forcing.
Both natural climate variability and anthropogenic trends have been suggested to influence ocean, ice, and atmospheric conditions that, in turn, control polynyas and deep convection in the Southern Ocean. For instance, under preindustrial conditions, several climate model simulations feature natural oscillations in deep convection on centennial time scales (Latif et al. 2013; Martin et al. 2013; Cabré et al. 2017). When convection is absent, annual-mean surface air temperature and SST decrease across the Southern Ocean. The absence of polynyas in recent decades may also be related to externally forced variability. For instance, the Southern Ocean surface has undergone a recent freshening trend, which inhibits convective mixing and polynyas by stratifying the water column. The freshening has been related to increased melting from Antarctic glaciers (Bintanja et al. 2015), an intensification of the global hydrological cycle (Durack et al. 2012), and increased precipitation from a positive trend in the southern annular mode (SAM) since 1979 (Gordon et al. 2007). However, the connection between SAM variability and high-latitude precipitation has since been questioned (Spensberger et al. 2020). Furthermore, recent observations have suggested that positive SAM anomalies can instead promote deep convection through enhanced wind-driven upwelling of high-salinity water in the Weddell Sea (Campbell et al. 2019). The connection between observed, externally forced surface freshening and the reduction of deep convection is supported by trends in phase 5 of the Coupled Model Intercomparison Project (CMIP5). Of the 36 CMIP5 models examined in De Lavergne et al. (2014), the 25 models that feature Southern Ocean deep convection become increasingly stratified as the sea surface freshens under anthropogenic forcing. Convection ceases fully in all of the initially convective models by the year 2300.
Despite these insights into high-latitude ocean dynamics, until very recently fully coupled models have lacked the resolution to explicitly resolve open-ocean polynyas that are spatially or temporally consistent with observations. This shortcoming is an important obstacle to understanding how the ocean, atmosphere, and sea ice interact in the Southern Hemisphere. Models vary widely in their representation of Southern Ocean deep convection, both spatially and temporally (De Lavergne et al. 2014), which may influence the degree to which sea ice and SST respond to anthropogenic freshwater forcing (Armour and Bitz 2015). Inconsistent representation of deep convection may also affect atmospheric variability, since the amount of anomalously open water during austral winter will impact the magnitude of air–sea heat flux (Moore et al. 2002). Our study seeks to improve understanding of coupled high-latitude climate processes by placing emphasis on polynyas and their interaction with the atmosphere. We analyze these interactions in the context of the long-term buildup of the deep ocean heat reservoir, which preconditions the Southern Ocean for deep convection.
Advances in state-of-the-art modeling capabilities have provided new opportunities to study polynya heat loss. Recently, a fully coupled synoptic-scale preindustrial control simulation of the Energy Exascale Earth System Model version v0.1 (E3SMv0) has effectively resolved Weddell Sea open-ocean polynyas that are consistent with observations in their spatial extent and location. Prior research with this model has addressed the preconditioning requirements for Maud Rise polynyas (Kurtakoti et al. 2018), while results from a similar Community Earth System Model (CESM) configuration (Small et al. 2014) have addressed the local atmospheric response (Weijer et al. 2017). Our research builds on this previous work through a larger-scale investigation: we use polynya heat loss as a metric to investigate high-latitude climate variability through coupled ocean–ice–atmosphere interactions. We draw connections between polynya heat loss variability and three important climate processes in the Southern Hemisphere middle and high latitudes: 1) surface temperature and heat flux anomalies, 2) the zonal wind pattern over the Weddell Sea, and 3) the poleward transport of heat and its partitioning between the atmosphere and the ocean. We show that all three processes are significantly impacted by polynya heat flux variability, highlighting the importance of correctly resolving these highly episodic features as models gain higher resolution.
We address the challenge of determining causality in coupled systems, where feedback processes are often present. Along with impacting the atmosphere, polynyas can themselves be influenced by atmospheric processes. For instance, idealized simulations have induced Weddell Sea polynya formation by applying additional westerly wind stress forcing, resulting in a spinup of the Weddell Gyre (Cheon et al. 2014, 2018). Such modeling experiments preclude a wind response to polynyas. We instead assess coupled processes through the statistical framework of Granger causality. This statistically robust tool can identify causal relationships in any time series, observational or simulated. It also accounts for memory in the climate system, an important consideration when performing lagged regressions with geophysical time series (McGraw and Barnes 2018). In this framework, the word “causality” strictly refers to quantifying statistical predictability, rather than true causality. The method is nonetheless useful for disentangling drivers from responses in ocean–ice–atmosphere interactions. Following prior work in the Arctic (Matthewman and Magnusdottir 2011; Kretschmer et al. 2016; Samarasinghe et al. 2018), we analyze the degree to which Weddell Sea polynyas drive and/or respond to climate variability in a comprehensive setting, enabling new insights into the elusive role that Weddell Sea convection can have in the high-latitude climate system.
A high-resolution configuration of the E3SMv0 model (E3SMv0-HR) was run for 131 years under fixed radiative forcing conditions representative of the year 1850. The ocean component is the Parallel Ocean Program, version 2 (POP2), configured at a nominal 0.1° resolution. It is coupled to a sea ice model component at the same resolution, the Los Alamos Sea ice Model, version 4 (CICE4). The atmosphere component is the Community Atmosphere Model, version 5 (CAM5), configured at a nominal 0.25° resolution. The land component is the Community Land Model, version 4.5 (CLM4.5), configured using the same grid as the atmosphere model. All fields considered in this study were saved as monthly averages. Notably, we find that a companion run at 1° ocean and atmosphere resolutions does not simulate strong open-ocean convection.
In E3SMv0-HR, we define an open-ocean polynya as a fully enclosed region of ice-free ocean (sea ice fraction less than 15%) that persists through all 3 months of maximum sea ice extent in the Southern Hemisphere: August, September, and October (ASO). By this definition, 52 out of 127 years (41%) in the E3SMv0-HR simulation feature open-ocean polynyas, a more common occurrence than what is seen in the observational record. That said, episodes of major Weddell Sea polynyas are interspersed by periods of stronger stratification that persist for many decades, matching the post-1976 drought of Weddell Sea polynyas. Our simulation therefore differs from the convecting CMIP5 models that exhibit chronically weak stratification in the Weddell Sea, where mean September mixed-layer depths often exceed 2000 m (De Lavergne et al. 2014). Both the location and spatial extent of the simulated polynyas agree well with observations to date. The polynyas form exclusively in the Weddell Sea, centered around the location of the Maud Rise seamount (65°S, 2.5°E), with a median size of 1.71 × 105 km2, approximately 68% of the mean polynya area observed during the 1974–76 event. Smaller open-ocean polynyas with an area on the order of 104 km2, comparable in size to the open-ocean polynyas observed in 2016 and 2017, appear sporadically in the first four decades of the simulation (e.g., years 24–27, 31–35). The simulation’s larger polynyas propagate westward through the Weddell Sea, occurring consecutively in three distinct episodes associated with deep ocean convection, each lasting 10–15 years. These episodes also feature patches of open water that extend to the sea ice edge, forming large embayments. These features do not meet the conventional definition of an open-ocean polynya; embayments are larger, not fully enclosed by sea ice, and have not been seen observationally in austral winter. However, in E3SMv0-HR, they appear to be driven by the same convective process and release heat to the atmosphere in a similar manner as polynyas. Thus, given the necessity of continuous time series and the utility of a binary classification, we classify the 52 years containing polynyas and/or embayments as “polynya years” through visual inspection, and classify the remaining 75 years as “non-polynya years.”
After identifying open-ocean polynyas, we relate their occurrence to several metrics of the high-latitude energy budget. Upward surface heat flux F for a given year is defined as the ASO averaged sum of the turbulent and radiative heat fluxes: net shortwave radiative flux (FSNS), net longwave radiative flux (FLNS), latent heat flux (LHFLX), and sensible heat flux (SHFLX) in units of watts per square meter:
where SHFLX, LHFLX, and FLNS are positive upward, FSNS is positive downward, and the LHFLX term is modified to account for the latent heat of fusion from frozen precipitation. Accordingly, the sum of these terms is positive when heat is lost from the ocean and taken up by the atmosphere. The value of F is calculated for each grid cell. We isolate the seasonal quantity since open-ocean polynyas are only active during the months of maximum sea ice extent in the Southern Hemisphere.
The mean state of poleward heat transport in the atmosphere (AHT) and ocean (OHT) is quantified as a function of latitude. We follow the convention of northward transport defined as positive, and southward transport defined as negative (e.g., Trenberth and Caron 2001). AHT is calculated as the meridional cumulative integral of the top-of-atmosphere radiative heat flux FT (positive downward) and surface heat flux F (positive upward):
where both flux terms are calculated as zonal means (denoted by square brackets), ϕ is latitude, and Re is Earth’s radius. On the time scales relevant to our study, we assume atmospheric heat storage is negligible.
While the effective ocean heat transport is commonly calculated as a residual between the total and atmospheric transport, it does not allow for a proper diagnosis of heat transport and storage. In particular, the time scale of variability in OHT changes when considering heat storage. This consideration is important for time series analysis of the heat transport components, as well as diagnosing the relative response times of the atmosphere and ocean to a perturbation, such as the occurrence of an open-ocean polynya. Therefore, rather than use the cumulative sum formulation of Eq. (2), we explicitly calculate OHT at each latitude as
where the integration is carried out over the total ocean depth (zbot < z < 0) and longitude λ, ρ is seawater density (103 kg m−3), cp is specific heat of seawater at constant pressure (3850 J kg−1 K−1), θ is the ocean potential temperature, and V is the ocean meridional velocity. The product θV is a diagnostic variable that is calculated at every model time step and saved as monthly averages.
To supplement our analysis of ocean heat transport, we additionally infer ocean heat content south of a latitude band as the difference between ocean heat transport (OHT, positive northward) and the meridional cumulative integral of surface heat flux F (positive upward), where both terms are integrated over time:
Finally, we diagnose two-way predictive causality between polynya heat loss and the climatic processes of interest, combining a vector autoregressive model framework with Granger causality. The general equation for a bivariate (n = 2), lag-1 vector autoregressive (VAR-1) model can be written as
where yn,t are stationary time series of endogenous predictor/response variables, cn are the model intercepts, an,n are regression coefficients, and en,t are independent, identically distributed error terms, calculated as the residual difference between observed and expected values of yn,t. Although VAR models can generally accommodate more than two variables, we solely conduct bivariate analyses, as in Eq. (5). Other causal inference techniques are better adapted to connecting larger numbers of variables in a single model (Runge et al. 2019) and our bivariate approach offers the advantage of simple, straightforward interpretability.
VAR model analysis requires stationary time series. To fulfill this requirement in each variable, we subtract the climatological mean for each month in the seasonal cycle and divide by the entire time series’ standard deviation. Any long-term linear trend is removed. After applying these changes, each time series represents normalized climatological anomalies, and parameters an,n and cn can be solved using ordinary least squares estimation. As a simple example, Eq. (5) has a maximum lag, or lag order, of one time step. However, in our analysis, the optimal maximum lag order is determined statistically by the Akaike information criterion (AIC). The AIC is a measure of relative log-likelihood, and is calculated for each model over 1, 2, …, 12(tmax/100)1/4 possible maximum lag orders, where tmax is the length of the time series. The lag order that minimizes the AIC is selected as the best fit and used for subsequent analysis. Our formulation in Eq. (5) focuses on lagged relationships and does not allow for instantaneous connections.
Following the framework of Samarasinghe et al. (2018), we then use our VAR model to test for Granger causality between each variable, in both directions. To accomplish this task, we designate a restricted model and unrestricted model for each response variable in yn,t, utilizing different subsets of the parameters in the prefitted VAR model [Eq. (5)]. The restricted model is an autoregression of the response variable against itself [Eq. (6)], and the unrestricted model is a multiple linear lagged regression of the response variable against both itself and the predictor variable [Eq. (7)]:
where i denotes individual lags, and the single and double primes on the lhs distinguish the respective restricted and unrestricted predicted values from observed values. In the case of Eq. (7), y2(t) can be said to Granger cause y1(t) if two conditions are fulfilled:
At least one regression coefficient of the predictor y2(t) must be significant according to a two-sided t test.
The unrestricted model must be able to explain more variance in y1(t) than the restricted model, according to an F test on the residuals.
If the null hypothesis of no Granger causality is rejected at the 5% level, the sign and magnitude of the predicted anomalies in y1(t) can be found using the sum of y2(t)’s regression coefficients. As an example, let Eq. (7) be produced from a VAR-12 model with i = 1 month and a1,2,1 + a1,2,2 + ⋅⋅⋅ +a1,2,12 = 0.5. If Granger causality exists, then a one standard deviation anomaly in y2(t) predicts a 0.5 standard deviation anomaly in y1(t). To test Granger causality in the other direction, one simply exchanges the place of y1(t), y2(t), and their coefficients in Eqs. (6) and (7). We present results from bivariate VAR models regressing polynya heat loss with surface temperature, zonal near-surface wind, mixed layer salinity, atmospheric heat transport, and ocean heat transport.
a. Changes in atmospheric surface properties
We analyze changes in the ocean-to-atmosphere surface heat flux F by calculating the ASO composite difference between polynya years and non-polynya years (Fig. 1a). Poleward of the ASO climatological sea ice edge, the Weddell Sea features a significant increase in upward surface heat flux during polynya years, with an average 12 W m−2 increase in the region where open-ocean polynyas occur. This region is denoted by the polynya mask in Fig. 1 and subsequent visualizations, calculated as the area where an open-ocean polynya occurs at least once in the simulation (embayments not included). The mask’s surface area (3.9 × 106 km2) is larger than the simulation’s largest polynya (6.6 × 105 km2). Accordingly, the mask’s shape reflects the year-to-year spatial variability in large polynyas as they propagate westward from the Maud Rise region.
The increase in heat loss over the Weddell Sea is also accompanied by a decrease of similar magnitude, north-northeast of the climatological sea ice edge. These two opposing changes have statistical significance primarily in the Atlantic and Indian Ocean sectors. The mechanism underlying this poleward shift in the surface heat flux, as well as its relationship to meridional heat transport, are analyzed in sections 3b and 3c.
To compare how the turbulent and radiative components comprising F [Eq. (1)] contribute to these changes, we next calculate the composite difference for the turbulent sensible and latent heat fluxes only (Fig. 1b). Both the magnitude and spatial pattern of the composite difference are similar to what is seen in Fig. 1a, indicating that the turbulent heat fluxes are the primary contributors to polynya heat loss. Sensible and latent turbulent fluxes each contribute roughly 10 W m−2 to the average anomaly over the polynya mask in Fig. 1b. The small differences that exist between Figs. 1a and 1b are explained by enhanced shortwave absorption during polynya years, which damps ocean heat loss south of the sea ice edge. As turbulent fluxes transfer heat from the ocean to the atmosphere, sea ice fraction and surface albedo decrease in and around the polynya, thus allowing more shortwave radiation to be absorbed at the surface. At all locations, the difference in net longwave radiative flux is negligible. This somewhat counterintuitive result is consistent with prior work using a similar model. Weijer et al. (2017) found that an increase in upward longwave radiation over a polynya is balanced by an equivalent increase in downwelling longwave radiation, which is associated with a warmer, moister atmosphere over polynyas and an optically thicker cloud deck.
Figure 1c shows the composite difference in ASO-average surface temperature between polynya years and non-polynya years. In the months of maximum sea ice extent, there are large, significant temperature increases of 10–12 K over the entirety of the Weddell Sea region where polynyas occur. In general, the largest differences occur poleward of the sea ice edge, but smaller increases appear throughout the Southern Ocean, extending northward up to 40°S. These changes are also statistically significant in the smaller, annual-mean composite difference (not shown), where the maximum temperature change is 3 K. In the Weddell Sea, this temperature change is consistent with prior analysis of open-ocean convection in a lower-resolution simulation (Cabré et al. 2017). In both the ASO average and the annual average, significant differences in surface temperature are more spatially extensive than the significant differences in ocean heat loss. Furthermore, polynya years consistently show an increase in surface temperature in all regions where a significant difference is found.
The polynya mask region straddles the border between the prevailing midlatitude westerlies (positive U) and the much weaker polar easterlies (negative U). Figure 1d shows that, in the Weddell Sea sector, both midlatitude westerlies and polar easterlies undergo a significant increase in speed during polynya years. It is worth noting that large-scale, albeit not statistically significant, composite differences in wind speed also exist outside the Weddell Sea. In polynya years, midlatitude westerlies exhibit a poleward shift in the east Pacific and Atlantic sectors and an equatorward shift in the Indian and west Pacific sectors.
b. Granger causality between polynya heat loss, surface temperature, wind, and mixed layer salinity
In the previous section, we determine concurrent relationships between polynya heat loss, surface temperature, and near-surface zonal wind. However, the system is complex and strongly coupled, obscuring the direction of causality. In this section, we use a combination of bivariate vector autoregressive (VAR) models and Granger causality to explore the time-lagged causal relationships that can exist between these processes. The spatial maximum of ocean heat loss F in the Weddell Sea motivates the additional definition of polynya heat loss Fp as the spatial average of F within the polynya mask. Since the maximum increase is centered over the polynya mask, Fp can be used as a time series that characterizes the heat loss attributable to open-ocean polynyas. Autocorrelation and coupling between surface temperature and polynya heat loss are explored with a Granger–VAR approach. In each grid cell, we design a bivariate VAR model between surface temperature TS and our index time series Fp. Here, both time series have monthly resolution. After building the statistical models, a Granger causality test is conducted in both directions for each model and each grid cell, providing a spatial perspective on the coupling between the two quantities.
To aid our evaluation of the Granger–VAR analysis, we first examine the temporal relationship between the polynya heat loss and surface temperature time series (Fig. 2a). Along with polynya heat loss, we define an accompanying time series of high-latitude surface temperature Ts,hl as the spatial average of all grid cells poleward of 55°S, including Antarctica. Since most of the temperature differences between polynya years and non-polynya years occur poleward of 55°S (Fig. 1c), it is a useful boundary for assessing changes in the Southern Hemisphere equator-to-pole temperature gradient. The terms Ts,hl and Fp are correlated throughout the entirety of the simulation (r = 0.6), and positive anomalies in both time series peak during polynya years. During the longest consecutive appearances of open-ocean polynyas, Ts,hl and Fp will increase for several years at a time, with the temperature time series sustaining a positive anomaly until the polynyas cease appearing. When these positive anomalies are sustained, peaks in Ts,hl typically lag peaks in Fp by 2–3 years.
The Granger causality tests in each grid cell are largely consistent with this initial examination: increases in Fp predict a positive anomaly in surface temperature ([Fp → Ts], Fig. 2b). For each grid cell where Granger causality exists, the AIC selects maximum lag orders ranging from 12–24 months, with a mean of 19 months. The lags explaining the most variance are at 1–3 and 12–14 months. These peaks reflect the immediate response to polynya heat loss, as well as the tendency for open-ocean convection to occur for multiple consecutive winters. In other words, if polynya heat loss causes positive surface temperature anomalies in one year, it is more likely to do so again the next year. The lagged regressions between the response term (Ts) and predictor term (Fp) are relatively constant across the 12–24-month maximum lag order. By contrast, the autoregressions of surface temperature against itself are large in the first 1–2 months, and then decay to zero. Therefore, the effect of polynya heat loss on surface temperature has a longer time scale than the memory of surface temperature. The total predicted anomalies in surface temperature are highest over the western Weddell Sea and uniformly positive throughout the rest of the Weddell Sea sector. Most of the changes are confined poleward of the ASO sea ice edge. The magnitude of the predicted anomalies is consistent with the annual mean composite difference between polynya years and non-polynya years, but smaller than the ASO-composite difference (Fig. 1c). The smaller predicted anomalies reflect that the Granger–VAR analysis utilizes all months in the year, not just the season of maximum polynyna heat loss. Figure 2c shows the results of the Granger causality test in the other direction: [Ts → Fp]. Temperature has less power in predicting polynya heat loss than vice versa, with the exception of the polynya mask region itself. This difference confirms that the anomalous open water during polynya years is heating the atmosphere through a deep ocean heat source.
Next, coupling between polynyas and near-surface zonal wind U is quantified by constructing additional bivariate VAR models in each grid cell: ([Fp, U]; Figs. 3a,b). As with [Fp, Ts], these models also use monthly time stepping intervals, with a similar 12–24-month maximum lag order selected by the AIC. The spatial pattern of Granger causality in the Weddell Sea shows that polynya heat loss and Weddell Sea zonal wind are coupled over this approximately 2-yr window of predictability. In this bivariate VAR model, enhanced subpolar westerly winds are a response, rather than a driver, of polynya heat loss (Fig. 3a). The [Fp → U] causality test predicts enhanced subpolar westerlies from 60° to 50°S, with similar meridional boundaries as the polynya mask region. This regionally confined intensification of the westerly wind is accompanied by an intensification of the easterlies directly south of the polynya mask. Taken together, this meridionally confined region of predictability is consistent with the zonally asymmetric wind anomaly shown in Fig. 1d: an enhanced cyclonic wind circulation in the Weddell Sea and a locally enhanced temperature gradient between the polynya mask and the surrounding seasonal sea ice region. The westerly wind anomaly is also spatially correlated with the region of enhanced ocean heat uptake in Fig. 1a. We consider the possibility that westerly wind anomalies promote northward Ekman advection and anomalous cooling of the surface ocean north of the polynya on yearly time scales. However, this mechanism should be accompanied by negative subpolar surface temperature anomalies and northward sea ice expansion, which is not seen in the E3SMv0-HR simulation. Instead, we hypothesize that the anomalous ocean heat gain north of the ice edge is promoted by advection of anomalously warm and moist air that has passed over a polynya, and by enhanced heat transfer rates due to larger wind speeds.
Granger causality in the [U → Fp] direction is found through a large extent of the Southern Ocean. Polynya heat loss is predicted by circumpolar easterly wind anomalies near the Antarctic continent, while also being predicted by westerly wind anomalies occurring in more northern latitudes of the Southern Ocean’s Atlantic and Indian sectors (Fig. 3b). Although the westerly wind anomalies are not circumpolar in nature, the spatial pattern of Granger causality is nonetheless suggestive of the SAM’s negative phase: a northward shift in the westerly wind belt. The average VAR model lag order for this test is 18 months, leading us to hypothesize that the equatorward shifted westerlies (prolonged negative SAM; Fig. 3b) precondition the Weddell Sea for open-ocean convection during non-polynya years. We seek to identify the long-term preconditioning mechanism by calculating the contemporaneous relationship between Southern Hemisphere surface conditions and monthly Weddell Sea wind stress curl anomalies, but during non-polynya years only (Fig. 4). The mean Weddell Sea wind stress curl is negative, so a positive anomaly indicates a weaker, less cyclonic wind stress curl, and vice versa. Our use of this metric is motivated by a prior examination of E3SMv0-HR; Kurtakoti et al. (2018) found that non-polynya years generally feature positive wind stress curl anomalies in the Weddell Sea, which then become negative immediately prior to open-ocean convection. We find that these positive wind stress curl anomalies are associated with increased westerlies equatorward of 60°S, consistent with a negative phase of the SAM (Fig. 4a). Therefore, we verify that the simulated Weddell Sea wind stress curl anomalies are negatively correlated with the SAM index, as has been shown in observations (Campbell et al. 2019). During the months of maximum sea ice extent (August, September, and October), the wind stress curl anomalies are positively correlated with sea ice fraction in the Weddell Sea, while being negatively correlated with upward surface heat flux F over the seasonal sea ice region (Figs. 4b,c). For a 1σ positive wind stress curl anomaly, upward surface heat flux reductions in the Weddell Sea are on the order of 10 W m−2, while in the same areas sea ice fraction increases in the range of 5%–10% per grid cell. The prolonged negative SAM identified in Fig. 3b is therefore associated with reduced surface heat loss in non-polynya years. As will be shown in section 3c, this reduced surface heat loss leads to a buildup of high-latitude ocean heat content that Granger causes polynya heat loss on interannual time scales.
We further explore the role of large-scale zonal wind variability by examining its relationship with approximated high-latitude mixed-layer salinity (SALT) anomalies, averaged over the top 100 m of the water column. Our interest in the role of salinity is motivated by observations (Campbell et al. 2019) and model simulations (Cheon et al. 2014) that demonstrate the triggering of deep convection by high SALT anomalies concurrent with the opposite, positive phase of the SAM (poleward shifted westerlies). First, we find that positive SALT anomalies both within and upstream of the polynya mask region predict polynya heat loss on lag-order time scales of 18–24 months ([SALT, Fp], Fig. 5a). The lag-order time scale and upstream predictive region are suggestive of westward-propagating, high-salinity anomalies that often precede polynya episodes in E3SMv0-HR, initially identified in Kurtakoti et al. (2018). Given this resemblance, we proceed to represent these upstream anomalies with a single time series of SALT, horizontally averaged over the eastern Weddell Sea (SALTEWS), in particular over the region 62°–68°S, 15°E−60°E (denoted by the dashed box in Fig. 5). Then, as with [Fp, U], we relate this single time series of mixed-layer salinity anomalies to the time series of zonal wind anomalies in each grid cell. In Fig. 5b, [SALTEWS → U] displays a similar pattern of Granger causality to [Fp → U] (Fig. 3a), although the regression coefficients are smaller. This similar wind pattern further supports the role of mixed layer salinity east of the polynya region in predicting polynya formation. By contrast, [U → SALTEWS] (Fig. 5c) displays a different pattern of Granger causality from [U → Fp] (Fig. 3b): positive salinity anomalies are predicted by an intensification of the circumpolar westerly winds between 70° and 50°S. As suggested by past observations and model experiments, the zonal wind pattern in Fig. 5c, which is representative of the SAM’s positive phase, imposes a negative wind stress curl over the Weddell Sea basin (Cheon et al. 2014, 2018; Campbell et al. 2019). It is also worth noting that this zonal wind pattern only Granger causes salinity anomalies in the upstream region (SALTEWS); U does not predict salinity anomalies when they are averaged over the polynya mask region (not shown). Again, this distinction implies that salinity anomalies are first introduced upstream of the polynya mask before being advected westward.
Several mechanisms could explain the generation of positive salinity anomalies in the eastern Weddell Sea. First, the negative wind stress curl spins up the Weddell Gyre circulation; and indeed, polynya heat loss is predicted by a stronger barotropic streamfunction in the Weddell Sea (not shown). This strengthening pulls in more salty water from the Antarctic Circumpolar Current in the southern branch of the Weddell Gyre. However, this causal pathway is somewhat complicated by the fact that the optimal lag order of the predicted anomalies in Fig. 5c is only about 2 months, hence much shorter than the 18–24 months in Fig. 5a. This lag order seems too short to be consistent with this advective mechanism. Alternatively, it is conceivable that the upstream salinity anomalies are caused by wind-driven upwelling. To investigate this possibility, we estimate Ekman suction velocities from wind stress curl. Excepting the coastal region, eastern Weddell Sea Ekman vertical velocities are approximately 0.5 m day−1 in the annual mean and exhibit small positive anomalies prior to polynya years (not shown). However, these small anomalies (0.05 m day−1), if sustained over several months, would only produce approximately 5 m of additional vertical displacement, which seems to be insufficient to elevate the core of the Weddell Deep Water to the surface. We conclude that other wind-driven processes may play an additional role in introducing the upstream salinity anomalies, such as anomalous brine rejection and ice production along the coast.
In summary, we show that large-scale zonal wind variability Granger causes both polynya heat loss and positive upper-ocean salinity anomalies in the Weddell Sea. The distinct time scales and spatial patterns that we isolate suggest a complex relationship between SAM variability and open-ocean convection. First, a prolonged negative phase of the SAM (Fig. 3b) preconditions the polynya mask region for open-ocean convection by inhibiting upward surface heat flux during non-polynya years (Fig. 4). Then, a poleward intensification of the westerlies (Fig. 5c) acts as a short-term trigger for deep convection by spinning up the Weddell Gyre circulation and introducing positive salinity anomalies that destratify the water column (Fig. 5a). This relationship between open-ocean convection and large-scale zonal wind patterns was not initially evident in the composite difference between polynya years and non-polynya years (Fig. 1d), likely because a composite difference does not distinguish between the zonal wind variability that predicts polynya heat loss (Fig. 3b) and the more localized zonal winds that respond to polynya heat loss (Fig. 3a).
c. Changes in poleward heat transport and ocean heat content
Last, we move outside the Weddell Sea to analyze larger-scale changes in meridional heat transport and ocean heat content during polynya years, as well as the relative partitioning of transport between the atmosphere and the ocean. In our initial analysis, we find that AHT and OHT are negatively correlated poleward of 60°S and positively correlated equatorward of 60°S (Fig. 6). The correlation values increase when we use decadal, rather than annual averages, reflecting the decreasing importance of oceanic heat storage on longer time scales. Short-term variability in OHT is thus likely buffered by heat storage, while long-term variability is compensated by surface fluxes. This distinction motivates a switch from monthly time series resolution to annual resolution in the Granger causality analysis for heat transport; a VAR model with monthly time steps would require an excessive lag order to capture multidecadal variability in the ocean. We accordingly conduct Granger–VAR analysis for AHT, OHT, and Fp using yearly time step intervals, with separate VAR models being constructed at each latitude. The composite difference in heat transport anomalies between polynya and non-polynya years is shown in Fig. 7, along with the significant Granger causality results and VAR model lag orders.
First, no significant Granger causality is found in the [AHT → Fp] direction or the [OHT → Fp] direction (not shown). In other words, polynya heat loss is not predicted by changes in meridional heat transport in the atmosphere or ocean. Second, we find a decrease in poleward atmospheric heat transport predicted by polynya heat loss, shown with the [Fp → AHT] causality test (Fig. 7a). Polynya heat loss predicts a positive anomaly in AHT (poleward weakening) between 70° and 40°S, on a time scale of 1 year. The sign of the predicted atmospheric response, as well as the maximum composite difference around 60°S, suggests that the changes in atmospheric heat transport are connected to the reduced meridional temperature gradient described in section 3b, and an associated reduction in dry static energy transport. Reduced atmospheric heat transport in response to polar-amplified warming has also been seen in the Arctic (Hwang et al. 2011). A small subset of VAR models for [Fp → AHT] have a larger lag order between 50° and 40°S, a notable visual irregularity. Given that a relatively fast response time is expected in the atmosphere, we interpret this different lag order as a statistical outlier associated with the AIC selection algorithm, rather than a distinct geophysical process.
Finally, we find that significant Granger causality is also found with the [Fp → OHT] test, which predicts a negative anomaly (poleward strengthening) of ocean heat transport south of 60°S, and a weakening of ocean heat transport north of 60°S (Fig. 7b). The changing sign of the predicted response with latitude is consistent with a poleward shift in the ocean heat flux convergence. The predicted behavior of ocean heat transport also resembles the zonal-mean turbulent heat flux anomalies (Fig. 1a), which display a similar poleward shift. The lag order of the OHT VAR models is spatially distinct: the window of predictability is 5–6 years in the high latitudes, and 1 year in the midlatitudes. This maximum lag order difference may reflect separate mechanisms for the predicted ocean heat transport change in each region. In the high latitudes, it likely takes several years for polynyas to deplete the deep ocean heat reservoir associated with open-ocean convection in the Weddell Sea, and for the ocean to replenish the reservoir. Accordingly, there would be a relatively long response time for ocean heat transport to compensate polynya heat loss in the high latitudes. By contrast, the faster midlatitude OHT changes predicted by polynya heat loss are likely mediated by the atmosphere. As seen in our analysis of Fig. 3a, the localized zonal wind response to polynya heat loss can promote anomalous ocean heat gain north of the climatological sea ice edge, despite the concurrent reduction in poleward atmospheric heat transport (Fig. 7a). This atmospherically driven heat gain is consistent with the midlatitude OHT response. In summary, poleward ocean heat transport reduction in the midlatitudes represents the fast response to polynya heat loss, whereas ocean heat transport enhancement in the high latitudes represents the slow response to polynya heat loss.
Our finding that meridional heat transport anomalies do not drive polynya heat loss, but are rather a response to it, is counterintuitive; prior modeling studies of Weddell Sea deep convection suggest that a buildup of the high-latitude ocean heat reservoir is an important precondition for the occurrence of a polynya (Latif et al. 2013; Cabré et al. 2017). Therefore, to clarify our understanding, we conduct a similar Granger causality test for the ocean heat content south of 65°S ([Fp, OHC]), where OHC is calculated as the cumulative temporal integral of advective and surface heat fluxes according to Eq. (4). The results for these causality tests are shown in Fig. 8. As expected, ocean heat content features a negative anomaly in the high latitudes during polynya years, which is driven by polynya heat loss ([Fp → OHC]; Fig. 8a). Most importantly, however, ocean heat content is also a predictor of polynya heat loss in the high latitudes ([OHC → Fp]; Fig. 8b), with a decadal buildup of heat predicting stronger polynya heat loss. Since OHC is inferred here from the difference between OHT and the meridional integral of F, we can assess the relative contribution of these two terms to the total OHC. Figure 8c shows that the buildup of the heat reservoir south of 65°S is driven by a persistent reduction in ocean heat loss (positive OHC anomalies) before polynya years, when more sea ice can insulate the ocean and prevent heat escape. This heat buildup is counteracted by a reduction in poleward ocean heat transport (negative OHC anomalies). Hence, surface heat flux changes dominate the total trend in ocean heat content.
4. Summary and discussion
We apply statistical causal inference techniques to a high-resolution, fully coupled climate model, drawing robust connections between polynya heat loss and several metrics of atmosphere–ocean variability in the high-latitude climate system. Under preindustrial conditions, we find that upward turbulent heat fluxes shift poleward during polynya years, promoting enhanced ocean heat loss from the Weddell Sea, a reduction in the meridional temperature gradient, and an intensification of the Weddell Sea cyclonic wind circulation. Polynya heat loss also has remote impacts on meridional energy transport: poleward heat transport decreases in the atmosphere, while ocean heat flux convergence anomalies closely track the upward turbulent heat flux anomalies. Finally, our statistical analyses disentangle the local wind response to polynya heat loss from the larger-scale atmospheric variability that controls open-ocean convection. In agreement with prior model simulations as well as observations, our predictive causality models show that polynya formation is influenced by zonal wind variability representative of the SAM (Campbell et al. 2019; Cheon et al. 2014, 2018). Open-ocean polynyas in E3SMv0-HR show the capability to both drive and respond to atmospheric variability in the Weddell Sea.
As with the 1974–76 Weddell polynyas, the Maud Rise seamount appears to be a critical bathymetric feature for triggering open-ocean polynyas and setting the site of deep convection in our simulation. (Holland 2001). In E3SMv0-HR, the dynamical impact of Maud Rise is a direct consequence of its accurate representation by the fine model resolution (Kurtakoti et al. 2018). The resulting spatial structure of sea ice reduction and increased surface temperature sets a zonally asymmetric pattern of the atmospheric response (i.e., strengthened zonal winds over the Weddell Sea during polynya years; Fig. 3a).
We find bidirectional Granger causality between polynya heat loss and zonal wind, showing that changes in Southern Hemisphere atmospheric circulation both drive and respond to open-ocean convection, with the coupled relationship being mediated by upper-ocean salinity and sea ice. This result suggests that future investigations of atmosphere–ocean–ice variability in the Southern Ocean, whether through targeted modeling experiments or observational analysis, should acknowledge the potential bidirectional causality between these climate processes. For instance, idealized experiments perturbing the westerly wind field over the Southern Ocean have produced polynyas through a wind-driven spinup of the Weddell Gyre, promoting Ekman divergence of the surface waters in the Weddell Sea (Cheon et al. 2014, 2018). While our analysis is consistent with the results of these experiments, we also find that the Weddell Sea subpolar westerlies are enhanced in response to polynya heat loss. In our study, detailed statistical analysis allows us to reveal this coupled relationship.
The meridional heat transport response to polynya heat loss (Fig. 7) can offer insight into the interannual, episodic nature of open-ocean convection in E3SMv0-HR. We see that, north of 60°S, the ocean does not compensate the decrease in poleward atmospheric heat transport, a behavior one might have expected from the compensation theory of Bjerknes (1964). In the high latitudes, our analysis of ocean heat content illustrates that each episode of polynya heat loss is preceded by periods of reduced surface heat flux, larger sea ice coverage in the Weddell Sea, and a buildup of the deep ocean heat reservoir. In this model, it thus appears that sea ice cover alone can lead to heat buildup and a periodic occurrence of polynyas, despite counteracting anomalies in ocean heat transport. Future examinations of this process will be key for constraining the intervals between deep convection events, which vary by two orders of magnitude (1–100 years) among the 25 convecting CMIP5 models (De Lavergne et al. 2014).
While it is unknown what hydrographic changes preceded the 1974–76 Weddell polynyas, Campbell et al. (2019) argue that the warming trend of Weddell Deep Water in recent decades is insufficient to trigger the small open-ocean polynyas observed in 2016 and 2017. They diagnose a 0.03 K decade−1 warming trend of Weddell Deep Water (WDW) from 2002 through 2016, consistent with the 0.032 K decade−1 mean trend observed between 1977–2001 by Smedsrud (2005). In E3SMv0-HR, this rate is on the order of 0.1 K decade−1 for non-polynya years. It is certainly possible that our model unrealistically limits surface heat loss, either due to biases in freshwater forcing (Stössel et al. 2015), or unrealistically weak mixing of the Southern Ocean under sea ice (Heuzé et al. 2013). Even so, any warming trend must reflect an imbalance between surface heat loss and lateral heat supply; and any process that further limits surface heat loss, like anthropogenic freshwater forcing (De Lavergne et al. 2014), would tend to increase this imbalance. Therefore, unless the ocean advective heat supply will adjust and neutralize the heat content buildup, an episodic release of the stored heat will happen eventually. Future simulations of open-ocean polynyas will need to include both anthropogenic forcing and a longer simulation time to fully resolve these ambiguities.
Fortunately, the improved representation of open-ocean polynyas in E3SMv0-HR appears to be clearly tied to the resolution of the model components (Kurtakoti et al. 2018). As resolution continues to improve among other fully coupled models, so too should the representation of deep convection in the Weddell Sea. New opportunities will become available to study how polynya heat loss influences high-latitude climate variability, as well as how this variability imprints onto the regionally distinct responses to anthropogenic climate change. Our statistical analysis of a preindustrial-control simulation can be used to motivate and inform future model experiments studying this important aspect of the climate system.
This research was supported by the Regional and Global Model Analysis (RGMA) component of the Earth and Environmental System Modeling (EESM) program of the U.S. Department of Energy’s Office of Science, as contribution to the HiLAT-RASM project. Additional support for Z.S.K. was provided by the National Science Foundation Graduate Research Fellowship Program under Grant NSF DGE-1842400. We thank Prajvala Kurtakoti for useful discussions and for providing some of the post-processed data used in this study.
Denotes content that is immediately available upon publication as open access.