Transiting consolidated ice strongly influenced polynya area during a shrink event in Terra Nova Bay in 2013

Coastal polynyas in Antarctica are a window of air-sea energy exchange and an important source of Antarctic Bottom Water production. However, the relationship between the polynya area variation and the surrounding marine environment is yet to be fully understood. Here we quantify the influence of the volume of transiting consolidated ice on the Terra Nova Bay Polynya area with ice thickness data. Changes in transiting consolidated ice volume are shown to dominate the evolution and variation of the polynya during a typical polynya shrinking event that occurred between 19 June to 03 July, 2013, rather than katabatic winds or air temperature, which are commonly assumed to be the main drivers. Over the cold seasons from 2013 to 2020, the Terra Nova Bay Polynya area is highly correlated to the transiting consolidated ice volume. We demonstrate that thick transiting ice limits the polynya area by blocking the newly-formed sea ice from leaving. Thick transiting ice blocked newly formed sea ice from exiting the Terra Nova Bay polynya during a shrinking event in 2013 and had a stronger influence on polynya area than katabatic winds or air temperature changes, according to an analysis of sea ice and meteorological data.

P olynyas, the nonlinearly shaped open water and/or thin ice surrounded by sea ice 1 , are a frequent phenomenon in the polar region. In Antarctica, by one estimate, there are 37 typical coastal polynya areas 2 . They play a key role in the polar and global climate and marine system 3 . When warm water is exposed to the cold air, there is a large air-sea heat exchange, which results in enhanced sea ice formation 4,5 accompanied by brine rejection 6,7 . Despite their small areas in Antarctica (~1% of the maximum Antarctic ice area) 8 , coastal polynyas generate over 10% of the Antarctic sea ice volume 9 . Associated with the formation of sea ice, the brine is rejected, and in some polynyas, it produces Dense Shelf Water as an ingredient for Antarctic Bottom Water formation 10,11 , which is exported in the lower limb of the global overturning circulation 12,13 . Moreover, high salinity and convection in polynyas lead to high primary production, which is important for the marine ecosystem 2,14,15 .
In general, the overall surface area is the most important characteristic of polynyas. Its variation is thought to result from differing offshore winds speeds, which export the newly formed sea ice [16][17][18][19] . Air temperature is another factor, which controls ice production [20][21][22] . Moreover, with a case study of Mertz Glacier Polynya, Massom et al. 23 qualitatively described that the thicker fast-ice passing the outflowing zone of the polynya can also reduce its area. On interannual time scales, the polynya area is related to climate indices such as the Southern Annular Mode and El Niño Southern Oscillation [24][25][26] .
In this study, we focus on a coastal polynya: the Terra Nova Bay Polynya (TNBP) in the western Ross Sea (Fig. 1). The TNBP is well known for a great heat flux per area and high ice formation rates 4,9,11 , although the area of TNBP is small (annual averagẽ 4200 km 2 ) 27 . Numerical modeling suggests that 33% of Ross Sea high-salinity shelf water, a major component of Antarctic Bottom Water, forms in TNBP 28 . The TNBP is thought to be caused by the characteristics of the strong offshore katabatic winds from Reeves Glacier driving the sea ice and the presence of the Drygalski Ice Tongue (DIT) blocking sea ice carried by the winds and currents from the south [29][30][31][32][33] . There are 9 of the 13 major coastal polynyas that are influenced by both offshore winds and land-fast sea ice/ice tongue in Antarctic 13 . TNBP can represent them to a certain extent.
As with other coastal polynyas 34 , previous studies of TNBP areas focused only on the influence of winds and atmospheric temperature (e.g., refs. [35][36][37][38] .). Massom et al. 23 qualitatively described that the upstream fast-ice can limit the extent of Mertz Glacier Polynya, but to the best of our knowledge, no previous study has quantified the effect of the ice drift from the south and transiting the TNBP. Novel Antarctic sea-ice thickness data derived from the brightness temperatures measured by the European Space Agency's Soil Moisture and Ocean Salinity (SMOS) mission were recently made available 39,40 , providing a new opportunity to study such events and examine the impact of the upstream sea ice on the TNBP area variations.
In this study, we analyze the influence of the thickness of transiting consolidated ice (TCI) from the south on the TNBP area in addition to the air temperature and wind speed. To demonstrate the importance of TCI to the TNBP area, we show first a case that clearly shows that the TCI dominates the extent change of the polynya, then a statistical analysis in the cold seasons of 2013-2020. Finally, we confirm the mechanism of TCI affecting polynya area by the analysis of sea ice in the TNBP potentially affected region.

Results
A shrinking event in late June 2013. Between 19 June to 03 July, 2013, the TNBP shrank rapidly (Figs. 2b and 3). In the beginning, the TNBP was~5200 km 2 expanding to the east of DIT and sustained for about 5 days. On 24 June the area of TNBP decreased to~10 km 2 . During the last 9 days, the mean area of TNBP was about 47 km 2 .
Previous studies suggested that the offshore katabatic wind speed and air temperature force the change of TNBP area 29,35 . However, during this shrink event, the mean air temperature of TNBP increased from 20 June to 25 June, then decreased. 3 days before 24 June, the offshore katabatic wind speed reduced from 17.1 m s −1 to 4.3 m s −1 , and on 03 July, although the offshore winds were as strong as that of 19 June, while the TNBP was still far smaller than that of 19-June (Fig. 2b). This indicates that the trend of rapid area decrease from 23 June to 24 June and steadiness for the rest of the period cannot simply be the result of the offshore katabatic wind and air temperature. Moreover, we also calculated the cumulative offshore wind anomaly and sea ice production (see Supplementary Discussion 2), which are also found to not have dominant roles in this case ( Supplementary  Fig. 5). Figure 2b shows that the trend of the thickness of TCI was highly correlated with the TNBP area. From 23 June to 25 June, the TCI thickness observed upstream of the current (the black dashed box in Fig. 3) rose from~0.5 m to~0.8 m, corresponding to a rapid increase in thickness. The correlation coefficient between the TCI thickness and TNBP area is −0.93.
The TCI thickening is related to the ice drift caused by wind. The large-scale wind/pressure field (Fig. 4) shows that there was a cyclone in the north of the Oates Land until 22 June. Due to a barrier of land, it had less influence on the TNBP. On 23 June, the cyclone moved to the east of Cape Adare, and on 24 June, the winds over the region of TCI abruptly turned from westerly to southerly, which corresponds to the abrupt TNBP area change.
Following the wind direction variation caused by the cyclone, the upstream sea ice (the black dashed box in Fig. 3) drifted to the south of DIT rather than bypassing it. Meanwhile, due to the blockage of DIT, the TCI increased gradually. With wind and ocean currents 41 , sea ice accumulated from south to north and formed a blocking in the east of TNBP, which led to this shrink event (Fig. 3). After 24 June, the northward TCI flow and accumulation of sea ice has not stopped, which contributed to the increase of ice thickness in the dashed box ( Fig. 3), while at this time, the TNBP cannot shrink any further (the area is 0 km 2 ).
The connection between the TNBP area and the transiting consolidated ice thickness. Figure 5b shows the seasonal evolution of the TNBP area (See Supplementary Methods). The TNBP always connects to the open sea directly in summer (Fig. 5c). During the cold season, the mean polynya area is 2361 km 2 (in winter (JJA), the mean area is 1740 km 2 ) and the mean ice thickness in the Terra Nova Bay is 0.25 m (Fig. 5a). The mean thickness of TCI in the cold season is 0.6 m, which is within the SMOS ice thickness measurement range. More statistical results are obtained to support the deduction in the case. Similar to the shrink event above, the relationship between the TNBP area and the TCI thickness is also important in the cold seasons from 2013 to 2020 (Table 1). The partial correlation analysis of these factors shows that the TNBP area has a large partial correlation with the TCI thickness (R 2 = 0.22, p < 0.005), which is comparable to the air temperature (R 2 = 0.19, p < 0.005) and offshore wind speed (R 2 = 0.09, p < 0.005), and much larger than the winds parallel to the coast (R 2 = 0.01). This implies that the upstream TCI flow (Fig. 1c) plays a primary role in the extent of TNBP.
In Table 1, we also show the correlation coefficients in winter (June-August, JJA). Although there is a lower correlation in winter, due to the expansion/shrinkage of the polynya and the thinning/thickening of TCI, it is still significant (p < 0.005). This indicates that our conclusion is still valid.
Another interesting result is that the relationship between the TNBP area and the TCI thickness may be nonlinear, which can be found in the scattering of the two (Fig. 6a). Taking the logarithm of the transiting ice thickness, the partial R-squared (partial R 2 ) between the area and logarithmic TCI increases by 14% (from 0.22 to 0.25, both of them p < 0.005), and the R 2 also rises from 0.30 to 0.33 (Fig. 6b, p < 0.005), which suggests a nonlinear relationship. In general, the area-to-wind-speed; and area-to-air-temperature relationship are often considered to be linear. For the offshore wind speed, the assumption is reasonable in a simplified one-dimensional model; 19,42 for the air temperature, through scale analysis (parameters from refs. 43,44 ), the effect can also be simplified to a linear model from a polynomial function 45 . However, the extreme deformation measured in the outlet zone of the TNBP strongly suggests a nonlinear process of TCI blocking 46 . More details will be given in the sections below.  Thomas 46 . The grey line is the mean extent of the TNBP. The red box is the study area of the TNBP area, the blue box is where the TCI thickness is obtained, and the yellow box is the potentially affected region. The grey land mask comes from the GSHHG (a Global Self-consistent Hierarchical, High-resolution Geography database) coastline resources. d As in (c) but for the probability of polynya occurrence over TNBP.
Furthermore, with the rolling correlation, there are some periods when the correlation coefficient exceeds 0.8 (dotted line in Supplementary Fig. 6), in which the TCI thickness is more crucial than other factors. This high correlation may also relate to the large-scale wind field. Under strong southerly wind, the ice drift is closer to the TNBP rather than deflecting eastward, and the influence of TCI on the area can be enhanced. More frequent cyclones leading to stronger southerly winds in the Ross Sea may result in a greater correlation between the area and the TCI thickness, as shown in the Supplementary Discussion 3.
Relationships between the TNBP area and the transiting consolidated ice in different situations. Some other characteristics of TCI are shown by the correlation coefficient between the TNBP area and the TCI thickness in different intervals of air temperature and offshore wind speed in Fig. 7. The correlation is insignificant when both the temperature and wind speed are low. Under this condition, the TNBP is small and tight against DIT (such as after 25 June in 2013 in Fig. 3), which is far away from the transiting ice. With the increase in air temperature or offshore wind speed, the TNBP expands and approaches the TCI flow; the correlation between the TCI thickness and TNBP area increases. The correlation reaches a maximum (R 2 = 0.45) when the air temperature is between −11.5°C and −17.5°C, and the offshore wind speed is between 10 m s −1 and 20 m s −1 . Figure 7 also shows that when the TNBP extent is not beyond the shelter of DIT (>3500 km 2 , the area of Terra Nova Bay), i.e., the TNBP is still at the western part of the bay, the influence of TCI, which is always east of DIT, is very significant. This indicates that the influence of TCI flow is so strong that its impact on the area can be observed when it is at a distance from TNBP. However, if the offshore wind is strong enough it can force the ice-flow direction further east, reduce its effect and lower the correlation. The partial correlation coefficient in different intervals also shows similar characteristics ( Supplementary  Fig. 8).
The sea ice thickness in the potentially affected region. To explain the nonlinear relationship between the TNBP area and TCI, we do a decomposition on the normalized thickness of sea ice in the potentially affected region (the yellow box in Fig. 1c) by factor analysis (Fig. 8). The factor analysis depends on the Rotated Empirical Orthogonal Function (REOF) decomposition and reflects the latent common source of the change of variables (see Data and Methods) 47 . Figure 8 shows that the first 3 dominant modes explain 33%, 29%, and 22% of the signal respectively. The second mode decreases gradually with the distance from TNBP (the dark grey line). The temporal coefficient of mode 2 (PC2) is correlated with the TNBP area significantly (R 2 = 0.36, Fig. 8c). Moreover, PC2 also closely relates to the offshore wind speed (R 2 = 0.15), air temperature (R 2 = 0.27), and TCI thickness (R 2 = 0.17). These three factors are the key roles governing the TNBP area. Therefore, we assume that mode 2 represents the effect of polynya controlled by offshore wind, temperature, and TCI. The spatial pattern also indicates the direct influence of TNBP is within 25 km of the polynya.
The third mode has a structure of increasing amplitude with distance from the TNBP (Fig. 8a). The PC3 is correlated with the sea ice thickness in the meridional extension region where mode 3 is of high amplitude (Fig. 8e). This implies association with the ice flow in Fig. 1c, and therefore mode 3 reflects the effect of flowing sea ice in the open Ross Sea.
The first mode only has strong amplitude northeast of the TNB, where the newly-formed sea ice in TNBP outflows and the TCI passes through (Fig. 1c). It leads to a significant correlation between the PC1 and TCI thickness (R 2 = 0.16). The collision of newly-formed sea ice and TCI may be a major source of the signal in PC1. The convergence of sea ice flux estimated by winds from ECMWF Reanalysis v5 (ERA-5) with the winds-ice drift relationship in the potentially affected region can verify it partially, although there are some uncertainties in ice drift (see Data and Methods). After a 30-day smoothing, the convergence of sea ice flux exhibits a correlation coefficient of 0.57 with PC1 (Fig. 8b), suggesting that mode 1 represents the sea ice deformation, which is driven by ice convergence. Similar results were obtained using higher resolution wind fields from Antarctic Mesoscale Prediction System (Supplementary Fig. 10). Hollands and Dierking (2016) 48 which used sea ice flow data from Synthetic Aperture Radar to estimate sea ice divergence also suggests that sea ice accumulates in this region and deforms. Figure 7 shows a nonlinear relationship that when the TCI is thick enough, the thickness effect that limits the polynya size becomes "saturated" and further increases in ice thickness have less effect on the polynya area. The internal ice stress is the key factor and reason that TCI is nonlinearly related to polynya area. As the ice gets thicker it is able to transmit the internal stress generated by the ice in the open ocean into the Terra Nova Bay region. The thinner ice formed in the polynya pushes against the older ice, but the internal stress prevents it from moving, ice deformation occurs, and the polynya growth becomes limited.

Discussion and conclusions
In this study, we show that a rapid TNBP shrinking event from 19 June to 03 July, 2013 is explained by TCI changes, but not the commonly referenced factors of katabatic wind and air temperature. The upstream TCI thickness, which explains the TNBP area change well, is in turn driven by large-scale winds transporting more sea ice to the outlet zone. The leading role played by the TCI thickness can be also seen in the cold seasons from 2013 to 2020, especially when the background southerly winds are strong. Further analysis of sea ice thickness in the potentially affected region shows that it may relate to nonlinear processes associated with the deformation with TCI.
The study area in this article is close to the coast, meaning that the ERA-5 winds have larger uncertainties 49,50 , and there may be   Fig. 3 The sea ice thickness surrounding the TNBP from 19 June to 03 July, 2013. The color map responds to the sea ice thickness with the land mask from NASA's MOD44W (in white pixels). The vector shows the daily mean winds from ERA-5 data. The black box is the TNBP study area, and the black dashed box is where we obtain the TCI thickness. In the red dashed area, the edge of the TNBP is marked with a magenta line. The grey land mask is from the GSHHG coastline as in Fig. 1 51 , while it can reflect the real variation of TCI SIT (see Sea ice data). However, our study provides strong evidence that the TCI blocks outflowing ice and can govern the TNBP area change. It implies that the sea ice production and even the formation of dense shelf water can be affected by TCI indirectly. The impact of TCI on the polynya area is unlikely to be unique to the TNBP. For example, in Mertz Glacier polynya, which is the 4th largest polynya in Antarctica and located to the west of the Mertz Glacier Tongue, the inhibitory effect of TCI on the polynya extent was also qualitatively obtained 23 . Nihashi and Ohshima (2015) point out that the upstream land-fast sea ice, which is common for Antarctic polynyas, can also play a similar role as the ice tongues 13 . There are 9 of 13 major coastal polynyas in Antarctica where ice tongues/land-fast sea ice promote their formation 52 , including the Cape Darnley polynya, which is considered an important source of Antarctic Bottom Water 10 . We suspect that these polynyas are also influenced by TCI, though it is hard to verify at this time due to the thick ice surrounding them being beyond the SIT measurement range of SMOS. Moreover, in these polynyas, associated with the effect of thicker TCI tending to be saturated, TCI's influence may decrease, although it lacks a quantification verification. Our result suggests the pressing need to investigate the role of TCI blocking in governing the area of these polynyas. reanalysis (pink box in Fig. 1) from the European Centre for Medium-Range Weather Forecasts with a horizontal resolution of 0.25°× 0.25°. The ERA-5 data is also used to analyze the large-scale atmospheric circulation from 19 June to 03 July, 2013. While the reanalysis dataset can be used for synoptic weather situations over space and time, the meteorology data observed by AWS is more accurate. The average correlation between ERA-5 and Manuela AWS is 0.95 for air temperature, 0.6 for wind speed, and 0.4 for wind direction. Due to the larger error of the wind field in ERA-5 at Manuela, we use the data from Manuela AWS in the nearshore areas. Here we define the prevailing wind direction in the cold season perpendicular to the coastline (263.7°) as the offshore direction.

Data and methods
Sea ice data. The polynya area in the cold season of 2013-2020 is obtained by the polynya signature simulation method algorithm based on the Advanced Microwave Scanning Radiometer 2 brightness temperature data, with a resolution of 6.25 km × 6.25 km. Sea-ice thickness data derived from SMOS brightness temperatures are provided by the Alfred Wegener Institute for the period of 2010-2020 for Antarctica. Table 1 Correlation coefficients among the TNBP area and its inflecting factors in the cold season (15-April~15-October) and winter (June-August). The first number in each cell is R 2 , and the number in the brackets is r. The first row shows the partial correlations. The asterisk (*) means that it has not passed the significant test (p = 0.005). The parameter numbers in the first row match parameters in the first column. Offshore Winds Air Temperature TCI Corr. Partial Corr. Fig. 6 The scatter between the TNBP area and the TCI thickness, and the correlation between the factors and the area. a Scatter and trend line (red) of negative exponent. b Offshore wind speed, air temperature, and TCI thickness. The left bars represent the results from the raw ice thickness and the right bars represent the results that the ice thickness is taken the logarithm. The outside bars show the correlation coefficient and the inside bars show the partial correlation coefficient.
Its resolution is 12.5 km × 12.5 km. SMOS ice thickness is retrieved under the assumption of 100% ice coverage, which leads to considerable underestimation of ice thickness when strong ice drift and ice divergence exist 40,54 . We define the thickness of TCI as the mean sea-ice thickness over the area southwest of the DIT and east of the Franklin Island (dark grey box in Fig. 1). The sea ice >1 m is not frequently observed here and is within SMOS ice thickness measurements range. Moreover, we verify that the SIT obtained by SMOS here with the underway observation and ICESat-2/CroySat-2 data (see Supplementary Discussion 1). With the observation by Rack et al. 51 , the SMOS SIT data can reflect the variation of ice thickness here, while there may be an underestimate. The potentially affected region is defined as the area including the TNBP and the region to its east (yellow box in Fig. 1). The ice drift is estimated by an empirical equation which is that we take 2% of the wind speed as the sea-ice drift speed, and the drift direction is 15°to the left of the wind direction 55 . It should be noted that due to the influence of topography and current, there should be some uncertainties in the sea ice drift predicted by this method in the near-coastlines region 56 .
Partial Correlation. In this study, we use the linear partial correlation coefficient to describe the relationship between one meteorological factor y and the TNBP area, excluding effects correlated to other factors. With the multiple linear regression, we first remove the signal of other factors from the TNBP area and y. The correlation of two residuals resulting from the regression is defined as the partial correlation coefficient. For example, to address the relationship between the TNBP area (A TNBP ) and the TCI thickness (TCI), and remove the effect of air temperature (T) and offshore wind speed (U), we first calculate the residual e Area and e TCI from linear regression of A TNBP with T, U and of TCI with T, U: where X is the matrix of controlled variables (i.e., [T U] T ) with the size n × 2, where n is the length of time series; W is the regression coefficients which can be estimated by the least-squares method. It is a 1 × 2 vector. The correlation coefficient between e Area and e TCI is the linear partial correlation coefficient, and the variance explained R 2 Area TCIÁT W is where Cov is the covariance; Var is the variance.
Factor analysis. Following Fischer et al. 47 , we employ factor analysis to decompose the sea ice thickness in the potentially affected region. The general step of factor analysis is to first decompose the variables by EOF and then rotate the major modes to maximize the sum of the variance of the factor loadings in dominating each mode, i.e., make them close to 0 or ±1. This rotation is called varimax rotation. Lian and Chen 57 discussed the advantages of the rotation applied to EOF. In EOF, a single mode may contain multiple signals, which leads to a too large weight of this mode and is difficult to explain. For example, in this study ( Supplementary  Fig. 9), the first mode of EOF accounts for over 65% of the total variance and the amplitude difference between different regions is small. We expect that the modes can highlight regions with large signals while the variances explained by each mode remain close. Therefore, we do a linear transformation on the primary modes, which redistributes the signals concentrated on a single mode. Supplementary Fig. 11 shows the progress of transformation. Compared with the raw modes (black), the factor loadings of the new mode 1 (blue) are polarized. Spatially, it means some regions are highlighted (e.g., mode 1 in Fig. 8a and Supplementary Fig. 9). The decrease of most factor loadings in mode 1 also leads to the decrease of the variance occupied by the first mode, which indicates that mode 1 may not contain too many signals. This helps us to understand the physical meaning of each mode. Supplementary Fig. 11 also shows that the linear transformation is actually a rotation. Moreover, the rotation of factor analysis relaxes the orthogonality constraint in EOF, which may avoid the generation of some artificial signals.
Here are the steps of factor analysis in this study. First, the matrix of normalized sea ice thickness (ASIT) is decomposed by EOF ASIT n;m ¼ X n;p T p:m ð3Þ where rows (n) of ASIV represent spatial positions and columns (m) represent time; T is principal components and X is factor loadings of modes. Then we take the first k principal components (here k = 3) and maximize the sum of the variance of each column of X by an orthogonal rotation matrix Γ, i.e., under the condition that Eq. 4b is satisfied, we need to obtain a matrix of Γ, which maximize S: where the I is an identity matrix. With the calculated Γ, the ASIT can be written as ASIT n;m ¼ X 0 n:k T 0 k:m þ ε n;m ð5aÞ X 0 n:k ¼ X n;k Γ k ð5bÞ -2 9~-2 3 -2 7~-2 1 -2 5~-1 9 -2 3~-1 7 -2 1~-1 5 -1 9~-  Fig. 7 The correlation coefficient between the TCI thickness and the TNBP area in different situations. The filled contour map shows the correlation between the TCI volume and the TNBP area in different ranges of offshore wind speed and air temperature; the contour map with labels shows the mean area in this range. The range span of air temperature is 6°C; that of offshore wind is 10 m s −1 . Dots on the map mark the insignificant results (p > 0.01). The partial correlation/correlation coefficient in winter is shown in Supplementary Fig. 8.
where ε is noise, which is called specific factors, X′ is the rotated modes, which is called the loading matrix, T′ is rotated principal components, which is also called the matrix of latent common factors.

Data availability
All data used in this study are publicly accessible from these websites:  c.

d. e.
Potentially Affected region