A New Approach for Monitoring the Terra Nova Bay Polynya through MODIS Ice Surface Temperature Imagery and Its Validation during 2010 and 2011 Winter Seasons

Polynyas are dynamic stretches of open water surrounded by ice. They typically occur in remote regions of the Arctic and Antarctic, thus remote sensing is essential for monitoring their dynamics. On regional scales, daily passive microwave radiometers provide useful information about their extent because of their independence from cloud coverage and daylight; nonetheless, their coarse resolution often does not allow an accurate discrimination between sea ice and open water. Despite its sensitivity to the presence of clouds, thermal infrared (TIR) Moderate Resolution Imaging Spectroradiometer (MODIS) provides higher-resolution information (typically 1 km) at large swath widths, several times per day, proving to be useful for the retrieval of the size of polynyas. In this study, we deal with Aqua satellite MODIS observations of a frequently occurring coastal polynya in the Terra Nova Bay (TNB), Ross Sea (Antarctica). The potential of a new methodology for estimating the variability of this polynya through MODIS TIR during the 2010 and 2011 freezing season (April to October) is presented and discussed. The polynya is observed in more than 1600 radiance scenes, after a preliminary filter evaluates and discards cloudy and fog-contaminated scenes. This reduces the useful MODIS swaths to about 50% of the available acquisitions, but a revisit time of less than 24 h is kept for about 90% of the study period. As expected, results show a high interannual variability with an opening/closing fluctuation clearly depending on the regime of the katabatic winds recorded by the automatic weather stations Rita and Eneide along the TNB coast. Retrievals are also validated through a comparison with a set of 196 co-located high-resolution ENVISAT ASAR images. Although our estimations slightly underestimate the ASAR derived extents, a good agreement is found, the linear correlation reaching 0.75 and the average relative error being about 6%. Finally, a sensitivity test on the applied thermal thresholds supports the effectiveness of our setting.


Introduction
In the polar regions, the exchange of heat, energy, mass, and momentum between the ocean and the atmosphere is strongly affected by the presence and the evolution of the sea ice cover [1] continuous formation of new ice by leaving the relatively warm open water exposed to a very cold atmosphere. The persistence of the TNBP is also due to the blocking effect of the Drygalski Ice Tongue, which shelters the bay from the ice drift from the south directed to NE [31] and controls, through its length, the polynya extent [33,37]. Observational and modeling studies indicate that the concentration of the pack ice outside the polynya also plays a role in the dynamics of TNBP [15,38]. It was estimated that this ice factory produces approximately 10% of the sea ice formed annually in the Ross Sea [15,39], and that it represents the largest producer of High Salinity Shelf Water (HSSW) [2,40,41], playing a crucial role in the formation of Antarctic Bottom Waters (AABW). Therefore, a better knowledge of the spatio-temporal variability of TNBP is needed to better understand the involved ocean-ice-atmosphere interactions and to provide a quantitative prerequisite to improve (i) the estimation of the associated ice production e.g., [23,42,43]; (ii) the analysis of the mass transport necessary to compensate the surface heat budget over the Ross Sea continental shelf [2,44] and (iii) the tuning of coupled ocean-atmospheric models and numerical simulations of the thermohaline circulation induced by polynyas [33,[45][46][47]. To this end, the most useful variable that can be retrieved from satellite imagery is the polynya area. Section 2 summarizes the previous efforts in the remote sensing of this variable in the Terra Nova Bay (hereafter TNB) through different sensors, while Section 3 presents a detailed description of satellite and Automatic Weather Stations (AWS) dataset included in this study as well as the developed methodology. Results are discussed in Section 4 where we also present: (i) an analysis of the katabatic events identified by AWS located along the TNB coast; (ii) a validation of our methodology through a comparison with several co-located satellite radar images; (iii) a comparison with results achieved through the application of a different procedure [27]; and (iv) an evaluation of the procedure sensitivity, in the prospective of its future application for studying the winter season long-term (2003-2017) variability of TNBP by using MODIS TIR observations. Conclusions are reported in Section 5.

Remote Sensing of the Terra Nova Bay Polynya
The coastal TNBP is generated and maintained by persistent katabatic winds (from Larsen, Reeves, Priestley, and David Glaciers) in the area included between the Drygalski Ice Tongue (at south) and the Campbell Ice Tongue (at north) ( Figure 1). During winter time, the open water mean extent usually ranges from 1000 to 1300 km 2 [15], but maxima up to 8500 km 2 were observed [48].
Several methods have been developed for estimating this extent using both PMW sensors and TIR imagery, even though the values obtained are slightly different. PMW sensors, such as SSM/I (/IS) and AMSR-E (-2) are routinely used for ice monitoring and proved to be excellent tools to study the location and the extent of the polynya that is typically extracted from maps of ice concentration derived from these data e.g., [16,28]. These techniques allow us to estimate the TNBP extent once or twice per day with a spatial resolution of several kilometers. However, Kern et al. [10] demonstrated that TNBP extent, location and orientation can change very fast, suggesting that openings and closings occur on timescales shorter than 24 h, usually in response to the regime of the katabatic events [36]. Ciappa and Pietranera [30] also suggested that fluctuations of the open water area can occur within very short time intervals, in the order of hours, and peaks of areal growing rate may exceed 300 km 2 per hour during extreme katabatic events.
These results pointed out the necessity of a more intensive sub-daily sampling frequency than the one allowed by PMW sensors and explains the recent efforts to improve the use of TIR data for routine monitoring of the TNBP activity [27,48]. Polynyas are clearly distinguished from the surrounding ice by TIR sensors, such as MODIS and the Advanced Very High-Resolution Radiometer (AVHRR). Several studies estimating TNBP area from AVHRR [39,49] exist, but they are severely affected by the limited number of observations due to the few satellite thermal sensors operational at that time. Including the two sensors operational on Aqua and Terra satellites, MODIS ice surface temperatures (IST) can be retrieved with a revisit time of few hours over polar regions, a spatial resolution of about 1 km, and a limited root mean square error widely sufficient to discriminate between the warm polynya core and the surrounding cold ice pack. Despite this shorter revisit time and the finer spatial resolution, up to 2010 TIR imagery has not been commonly used for monitoring TNBP; atmospheric disturbance is a basic limit and processing cloud corrections over many scenes requires a considerable amount of time and computing resources [11,27].
In recent years, some methodological studies presented the use of SAR active microwave systems, mostly operating at C-and X-bands, for the estimation of the TNBP extent e.g., [6,50]. These methods are based on the differences in backscattered radar intensity and image texture which allow the separation of open water, thin ice, and offshore pack ice. Spatial resolution (ranging from a few meters to about 100 m per pixel) is much finer than the PMW and TIR but the revisit time of individual SAR systems is often in the order of days, being SAR swaths much narrower than TIR or PMW sensor swaths; thus, SAR imagery is not typically used for routine monitoring activities. Moreover, even though more high-resolution SAR products are available today (e.g., those provided through the Sentinel constellation), only few images have been acquired over the TNBP during the last two decades.

Data and Methods
Several products were used to fulfill the objective of this study. Several MODIS datasets were processed for retrieving ice surface temperature (IST) and acquiring atmospheric information useful to identify scenes partially or completed covered by clouds and/or fog. Then, the AWS observations and the products used for validation purposes are introduced and described. Finally, the procedure that we implemented to select the informative MODIS scenes and to extract the TNBP extent is illustrated.

MODIS Thermal-Infrared Dataset
Hollands and Dierking [6] argue that the correct selection of specific satellite sensors for studying polynyas strongly depends on the scientific objectives. As we stated above, in this paper we are interested in developing a methodology for the long-term reconstruction of the TNBP variability during the last fifteen years. To this aim we opted for using TIR data acquired by MODIS sensors operational aboard the polar orbiting satellite Aqua of the NASA Earth Observing System. In particular, we used level-1b granules provided by the MODIS Atmosphere Archive and Distribution System NASA website (LAADS) to separately analyze all the swath-based scenes acquired during a day. Several products have been included in our procedure, as follows: MYD02 products have been used to retrieve the IST involved in the procedure estimating the polynya extent, as explained later in this section; geolocation parameters have been extracted from MYD03 products which include geodetic latitude and longitude, surface height above geoid, solar zenith and azimuth angles, satellite zenith and azimuth angles, land/sea mask for each 1 km sample, and additional information to enable the calculation of the approximate location of the center of the detectors of any of the 36 MODIS bands; MYD35 level-2 data have been used to collect swath-based cloud information to identify, and eventually discard, scenes partially or completely covered by clouds. All the observations available for the austral winter period (1 April to 31 October) from 2010 to 2011 have been analyzed.
This comprises a total of 31,030 images per year, i.e., 145 single scenes per day, having a spatial resolution of 1 km × 1 km at nadir and covering an area of about 1354 km (across track) × 2030 km (along track). Among them we extracted all the scenes observing the TNB study area (see Figure 1). Cloudy scenes have been completely discarded from our analyses, while both clear-sky and fog-contaminated scenes have been included, thus getting a revisit time of few hours in absence of thick clouds (see details in Section 3.4). In fact, even though the actual accuracy of the IST retrieval remains substantially unknown [51], it has been previously assessed that it is possible to map the sea ice also through fog [11]. In total, 1630 valid MODIS scenes, 868 in 2010 and 762 scenes in 2011, have been used for the wintertime TNBP extent estimation (Table 1).

AWS Wind Measurements
MODIS scenes allow us to identify the traces of the wind field too, since the katabatic drainage has a warm signature detected by IST. This effect, due to the adiabatic compression and to the strong vertical mixing within the wind layer [34], make the katabatic events visible in the IST maps by their sharply defined wakes as the air stream leaves the glacier valleys [11]. However, these traces do not allow the calculation of the true wind speed but only suggest the main katabatic pathways and the presence or absence of wind forcing at a given time. Therefore, in this study we included the wind information collected by the AWS located in the proximity of the Terra Nova Bay coast (see Figure 1) and maintained by the Italian Programma Nazionale di Ricerca in Antartide (PNRA). Rita AWS (74.72 • S, 164.03 • E) and Eneide AWS (74.69 • S, 164.12 • E) have been installed within the Meteo-Climatological Observatory of the PNRA near the Italian base "Mario Zucchelli", downstream of the Priestley Glacier. Even though these AWS are not exactly located within the main katabatic flow, previous studies demonstrated their usefulness in the analysis of the TNBP opening and dynamics [2,9,33,36]. A third AWS, Manuela (74.95 • S, 163.69 • E), is located on Inexpressible Island downstream of the Reeves Glacier, one of the main routes for the katabatic flows from the interior of Antarctica; unfortunately, no useful continuous data are available from this AWS during the 2010 and 2011 winter periods.
All the AWS are equipped with sensors to measure air temperature, relative humidity, atmospheric pressure, wind speed, wind direction and solar radiation. In this study, we analyzed only one-hourly intervals wind direction ( • N) and wind speed (m·s −1 ) data to characterize the katabatic wind regime in the Terra Nova Bay area during the 2010 and 2011 winter seasons.

Validation Products
C-band ENVISAT Advanced Synthetic Aperture Radar (ASAR) products were included in this paper for validating MODIS derived TNBP extent estimations. We analyzed a huge set of ENVISAT ASAR images acquired on (i) Wide Swath (WS); (ii) Image, Wave and Alternating Polarization (IM); and (iii) Global Monitoring (GM) modes. WS images have a pixel size of 75 m, and represent about 72% of our ASAR dataset, while IM (3%) and GM (25%) scenes have a resolution of about 30 m and 100 m, respectively. Furthermore, it is important to point out that WS and GM images have 400 km or larger wide swaths, while IM ones extends up to 100 km; all of them offer a perfect view of TNBP conditions. All scenes were retrieved from the EOLi-SA European Space Agency (ESA) online catalogue and processed with the NEXT SAR Toolbox. Unfortunately, only a limited number of scenes were acquired over the TNB area during the study period (an average of two scenes per month). Nonetheless, a total of 196 Envisat ASAR images were selected for our validation aim. Both backscatter and textural criteria were used to identify polynya extension; the typical ice streaks produced on open water by katabatic winds were used for a visual check of the polynya area after its computation [52]. Then, TNBP area was manually computed by tracing the polygon encircling the polynya and multiplying the number of pixels inside the polygon by the area of one pixel [28,29]. As expected, a great variability in the polynya size was observed in different ASAR scenes, but it was not possible to follow this variability daily due to the limited number of available scenes. The same occurs for other SAR products, e.g., those acquired through ERS, Cosmo-SkyMed (CSK), and Radarsat-2.
Finally, TNBP extent estimations retrieved from MODIS through the daily polynya area (POLA) procedure provided in Paul et al. [27] were used for a direct intercomparison with our estimations. Their technique derives daily polynya area (defined as area with open water and thin-ice between 0.0 and 0.2 m thickness) from cloud-cover corrected daily thin-ice thickness composites retrieved through the combination of all available MODIS TIR imagery and ERA-Interim atmospheric reanalysis data. To this aim, they make use of a spatial feature reconstruction scheme to account for cloud-covered areas [53].

The Polynya Extent Estimation Technique
In a clear-sky MODIS scene, the open water region is rather easy to estimate because it is visible as a warm and homogeneous area, with temperatures above the saline water freezing point. Nonetheless, we usually consider as polynya also the ocean surface that is partially covered by newly formed thin ice (i.e., pancake, frazil, and grease ice) around the open water. Thus, the area of polynyas is defined as the sum of the surfaces of open water and thin ice, where the thickness of the sea ice ranges from zero to a given threshold. The temperatures of these sea ice pixels depend on several factors (i.e., the ice thickness and concentration, and the air temperature) and can change from scene to scene in MODIS imagery. Therefore, a given single threshold is not enough to give a correct estimation of the TNBP extent. On the other hand, the presence of strong gradients of temperature, associated to the transition between open water and ice pack, allow us to detect its edge with a good accuracy, independently from the exact IST value that we are observing. Here we propose an automated procedure made up of several iterative steps for computing the extent of the TNBP area from MODIS TIR images ( Figure 2). First, MYD03 geolocation data are used to select all the MYD02 swath-based scenes covering completely the TNB study area ( Figure 1). The procedure inspects if every single swath covers the region and if the swath width is large enough to encompass entirely the study area. Land-and ocean-pixels are discriminated by using the NOAA National Geophysical Data Center coastline data for the Ross Sea.
Then, for each selected MYD02 scene, MYD35 data are analyzed in order to discard the scenes affected by atmospheric disturbance in which cloud cover compromises a confident evaluation of the ocean surface temperature [54]. This cloud mask provides four level of confidence regarding whether a pixel is thought to be seen through clear-sky or not. For each scene the percentage of pixels labeled confident clear (>0.99%), probably clear (>0.95), uncertain (>0.66) and not clear (≤0.66) over the ocean surface of the survey area is calculated. All the scenes, within which the percentage of the oceanic pixels defined by low confidence level (≤0.66) is less than 50%, are selected and processed. All the remaining scenes are definitively discarded. It is worthy to remark that the use of the MYD35 products, as discussed in Hall et al. [51], provides a conservative estimate of cloud cover but often does not map fog as cloud. Moreover, the accuracy of this cloud mask is reduced during the polar night due to the lack of visible channels, as well as the low thermal contrast between clouds and the underlying snow/ice [55,56]. Ciappa et al. [48] suggest selecting usable scenes by visual inspection, preferred to the use of the MYD35. Visual inspection allows the best screening for discarding scenes where the shape of the TNBP is not clearly detectable behind the fog. On the other hand, this approach involves several problems, i.e., the fact that unquantifiable errors might be introduced by the human inspector's ability, experience and attention, the absence of a standard parameter, and the enormous amount of time needed for an accurate inspection of thousands of scenes (especially when using level-1b calibrated radiances in place of IST MYD29 products).
Recently, Paul et al. [27] complemented the MODIS cloud mask information with ERA-Interim re-analysis of medium-level cloud-cover fraction data which were spatially and temporally interpolated to fit MODIS swaths. In this way they were able to introduce additional criteria to identify cloud-covered areas and to apply gap filling procedures, with respect to thin-ice thickness information, where cloud cover is identified. This approach is very promising, but the creation of a combined adjusted dataset includes several additional computing passages, i.e., the temporal and spatial interpolation of MODIS and ERA-Interim data on a common grid, a spatial feature reconstruction and the consequent proportional extrapolation [57], with the relative propagation of uncertainties. Still, this approach forces the user to move to daily mean values, thus losing the information about the TNBP extension variability within each single day.
In this study, we aim to develop an automatic procedure which allows retrieving as much reliable estimations of the TNBP as possible for each single day, with less help possible from any external user. Hence, we skipped visual inspection and limited the further analysis to the usable scenes, selected as described above.
In the third step, calibrated radiance data from MODIS IR channels 31 and 32, centered at approximately 11.0 µm and 12.0 µm, are converted to brightness temperatures by inversion of the Planck's function [58]. Then, IST values are calculated for the clear-sky scenes using the split-window technique originally developed by [59] and including the regression coefficient set for the Antarctic region [60]. The accuracy of this technique during winter conditions is discussed in Hall et al. [51] and Scambos et al. [61], which report no significant offset over the sea ice and an RMS error of about 1 K. The overall accuracy of retrieved IST was comprised between 0.3 K and 1.3 K when compared with standard MYD29 IST products.
Finally, IST is used as a discriminatory variable for the discrimination of the TNBP from the surrounding sea ice cover in each MODIS scene. Following statistics previously reported by [6,48], upper and lower threshold values are set to 265 K and 255 K respectively for distinguishing sea ice from open water. Still, a visual analysis of MODIS scenes in which TNBP was clearly open was carried out; this analysis confirmed that the 265 K value correctly locate the edge of the polynya. Hence, the algorithm classifies all the pixels characterized by an IST lower than/equal to 255 K as sea ice; those with an IST value equal to/greater than 265 K are defined as open water. Using this method, a large part of each single MODIS image is classified into a binary sea ice map with values of 1 (sea ice) and 0 (open water) at a 1 km resolution. This binary map is then used for a better selection of pixels characterized by IST values ranging between 255 K and 265 K.
An automated iterative procedure is then applied for classifying these intermediate pixels based on their neighborhood. For each of them, the procedure looks at the surrounding pixels, thus observing an area of 3 km × 3 km, and checks whether they have been classified as sea ice or open water. Accordingly, the undefined pixel is assigned to the most frequently occurring category. This operation is iterated moving from the coast to offshore.
Therefore, the sequence of the selected IST scenes (up to 13 per day, at time intervals between 15 min and 4 h) is used to estimate, for each selected scene, the extent of the TNBP by counting the number of 1 km × 1 km pixels classified as open water. Finally, TNBP retrieved extents are automatically compared with the values estimated for the previous and the following scenes to verify their consistency. Suspect IST scenes are pointed out in order to allow a visual check and the possibility of discarding them (less than 7%). As we expected, in most cases this happens during the closing phase when the TNBP is characterized by the presence of large areas of thin ice produced during subsequent katabatic events [30,42].

Terra Nova Bay Polynya Variability during Wintertime 2010 and 2011
According to the methodology described above, the TNBP area has been estimated for the 2010 and 2011 winter periods. Results are presented in Figure 3 where both single MODIS scene and daily average values are plotted. These values also include the open water disclosures at the extremity of the Drygalski Ice Tongue when present; according to previous studies, open water areas were found east of the glacier tongue tip in 21% of analyzed scenes. A high interannual variability is reported for both years as a result of successive opening and closing phases due to katabatic winds regime. As expected, short polynya events (2-5 days) are largely more frequent than long-lasting ones (10-15 days) during both 2010 and 2011 winter periods; the same happens for small openings (500-1000 km 2 ) which clearly represent the highest percentage of retrieved polynya events. In particular, TNBP openings larger than 1000 km 2 are present in less than 25% of the investigated scenes, while polynya areas larger than 4000 km 2 are identified only in 1% of the analyzed images. The relative frequency of MODIS scenes exceeding certain area thresholds is shown in Figure 4. An average polynya area of 1008 km 2 (Figure 3a) and 878 km 2 (Figure 3b) was retrieved for winters 2010 and 2011 respectively; estimates below 100 km 2 were not included in this analysis. However, we observed open water areas larger than 100 km 2 in the proximity of the Nansen Ice Shelf for more than 77% of the study period; this information is consistent with previous analyses [23,48]. In 2010, a major polynya event up to 5000 km 2 occurred during 19-21 August (days 231 to 233), followed by at least four large openings spanning between 3000 and 4000 km 2 . Large open water areas detected during the last week of October (days 277 to 304) seem to be associated to a premature beginning of the sea ice melting season more than to a typical polynya event. In 2011 three major TNBP openings reaching about 4000 km 2 can be detected during April, May and first half of October, with a polynya extent which is lower on average than in the previous year. Most surprising, only a single large opening (around 3000 km 2 ) was detected through MODIS imagery during mid-winter (August).
Even though previous studies estimated TNBP extents up to 7100 km 2 [28] and 8600 km 2 [62], using ENVISAT/ASAR and ERS/SAR products respectively, the observed fluctuations of the open water area within the TNBP are typical, with annual maxima which can exceed 4000-5000 km 2 , e.g., [32,48]. These results also confirm the important effect of the Drygalski Ice Tongue, as it blocks the sea ice advection from the south and frequently represents an upper limit to the TNBP offshore widening (the bay itself measures about 3300 km 2 when considering the glacier tongue length as its longitudinal side).
Generally, the retrieved mean winter area of about 900-1000 km 2 agrees with previous estimates based on thermal infrared data acquired by AVHRR ranging between 1000 and 1300 km 2 for different years [15,32,39]. Still, it is not far also from the average area of 1300 km 2 in May and June 2009, with a maximum extent of about 3000 km 2 , found by [11] using MODIS IST products. As for 2010, our results completely agree with corresponding scene by scene observations reported by figure 6 in [11] for the same winter season but using MODIS-derived IST products through a different processing scheme. The two TNBP area time series differ only on 20 August when they estimated a larger polynya extent of about 8000 km 2 , almost twice our estimates. Nevertheless, the analysis of a co-located ENVISAT ASAR image clearly confirms the magnitude reported in Figure 3a. This case study is described in detail in Section 4.2. As for 2011 (Figure 3b), a general agreement is found with TNBP extent fluctuations on 17-27 July (days 198 to 208) and 9-21 September (days 252 to 264) provided by Figure 4 in [30], through MODIS IST products and X-band VV-polarized CSK SAR imagery. In particular, position and magnitude of the relative peak detected on 23 July (around 2100 Km 2 ) perfectly correspond to our estimations; the peak on 15-16 September also corresponds but its extent is not comparable, probably due to the limited number of MODIS clear-sky scenes available for our procedure during those days.
On the other hand, our results contrast with the mean winter extent of about 3000 km 2 estimated, for different years, through passive microwaves by several authors [10,23,29,63]. Of course, these large discrepancies are due to the coarse spatial resolution of the passive microwave products and the adopted thresholds for ice concentration and thickness, rather than to the interannual variability of the TNBP. Different systems are observed when analyzing TIR and PMW dataset and a clear comparison is difficult: TIR observations provide a measure of the skin layer of the water/sea ice, while PMW derived products represent an integrated measure of their upper layers [51]. This means that MODIS IST enables to better separate thin ice from mixtures of open water and sea ice.  Furthermore, it is worthy to remark that the technique presented in this paper does not allow retrieving the TNBP area during persistent cloudy conditions. Nonetheless, it was possible to detect its extent with a revisit time shorter than 24 h for more than 89% of the study period, and with a revisit time shorter than 12 h for more than 74% of the time; only in 4% of the time it was impossible to get valid information about the state of the TNB region for periods longer than 72 h. Furthermore, the great majority of missing data in the time series, i.e., the discarded MODIS scenes, concerns the closing phase of TNBP when there is a higher probability of clouds and fog obscuring irreparably the study area in TIR imagery. This is consistent with the fact that, despite they typically form during the polynya growing phase [48,64], low clouds and fog are usually pushed offshore by strong winds during the katabatic events.
The retrieved polynya extent estimations have been compared with katabatic events detected by AWS Rita and Eneide wind speed records. Previous studies confirmed that the airflows from the different glaciers around the bay influence different shapes and magnitude of the TNBP opening events [6,47,48]. Airflows from the Reeves Glacier ( Figure 1) are always active when the wide polynya is open, while wind currents from Priestley and Larsen Glacier influence polynya extent in the northern and southern side of the bay respectively; the katabatic flows from the David Glacier, instead, cause the open water disclosure around the Drygalski Ice Tongue. This suggests that the AWS Manuela, which is located within the pathway of the Reeves Glacier airflows, usually represents the best source of katabatic wind observations. Unfortunately, its measurements were irreparably corrupted during wintertime 2010 and 2011 and could not be used in this study. On the other hand, AWS Eneide and Rita are very useful when analyzing strong polynya events caused by katabatic winds descending from the Antarctic plateau to the Nansen Ice Shelf through multiple pathways. These measurements are available for the whole sea ice season during the study period as reported in Figure 5. According to the strong directional constancy of the AWS wind observations, only speed measurements are reported. The linear correlations between the wind records collected by AWS Eneide and Rita and the TNBP area extent are summarized in Figure 5c. The estimated R-values are rather low but consistent with previous analysis coupling wind records acquired by the two AWS with the open water area in the whole TNB [48]. Of course, these results are influenced by several factors, i.e., the marginal location of the included AWS which does not allow to capture all the airflows responsible of polynya openings, and the fact that a minimum wind flow (for a certain amount of time) is necessary to maintain the polynya open. However, the correlations obtained for different time lags of the wind dataset confirms that (i) the TNBP does not respond immediately to the observed wind forcing and (ii) its opening is delayed by several hours. This delay seems to change from event to event, so that the average time of response deducted for 2010 and 2011 winter periods is considerably different, with peaks of correlation at 9 h and 14 h respectively. Results are also different from what observed by [48] for winter season 2009 (5 h delay). Furthermore, higher R-values (ranging between 0.29 and 0.51) have been derived when limiting the regression analysis to intense wind forcing (greater than 20 m·s −1 ) and long-lasting (more than 8 h) events. This is consistent with results previously presented by [36], who argued that the open water percentage increases in correspondence with katabatic events of long duration blowing with speed greater than 25 m·s −1 .

Terra Nova Bay Polynya Daily Variability
Despite their sensitivity to clouds, MODIS thermal infrared products give us the opportunity to acquire more TNBP images per day. Hence, it is often possible to monitor the complete evolution of its growing and closing phases within one or more days. Fluctuations of the open water area can occur within very short time intervals of the order of few hours [6,10]. Here we present two case studies about TNBP events showing a significant hourly variability dated to 18-23 August 2010 ( Figure 6) and 26-29 August 2010 (Figure 7). In between the two events (24-25 August) the bay was kept closed by sea ice cover due to the absence of intense airflows blowing from the continent. Both figures show some peculiar aspects of the growing and closing phases of a typical coastal polynya through the sequence of selected MODIS IST scenes and the corresponding TNBP extent estimations. Furthermore, wind regimes as measured by AWS Eneide and Rita are also reported, as well as high-resolution ENVISAT ASAR Wide Swath images co-located in space and time. The sequence of MODIS scenes reported in Figure 6 shows that the TNBP area increases sum a few hundred km 2 (19 August 03:10) to more than 4700 km 2 (20 August 12:00) in one day and a half, then remains stable for about 24 h (21 August 11:05) and finally disappears during the next 36 h (23 August 6:00). Figure 6 also shows that a few hours later (23 August 14:10) the situation is changing again, with small openings appearing along the Nansen Ice Sheet and the northern side of the Drygalski Ice Tongue. Unfortunately, MODIS TIR scenes on 22 August are critically affected by atmospheric disturbance. Those reported in Figure 7 describe an even higher level of variability, with TNBP opening in less than 8 h during 26 August, then being stable at about 1000-1300 km 2 for few hours (from 26 August at 13:00 to 27 August at 04:00), disappearing after less than 24 h (27 August at 13:45), and finally growing quickly again up to about 800-1000 km 2 (29 August). As expected, the TNBP area variability seems to be closely related to the wind speed regime recorded by AWS Eneide and Rita (Figures 6b and 7b) during both case studies. In particular, it is interesting to notice that, even though several times wind speed exceeded 20-25 m·s −1 , only when this condition persisted for some time (12 h and more) the TNBP opened extensively. This is clearly visible for the TNBP openings identified after long-lasting wind events on 26 August and mostly on 20 August, when the polynya grew up to about 4800 km 2 due to persistent katabatic winds blowing at more than 30 m·s −1 for about 24 h. Furthermore, the MODIS IST scenes point out that the TNBP area is clearly characterized by a sharp thermal contrast between the warm open polynya and the cold sea ice offshore on 20 August, the same as during its growing phase represented through the previous scenes. Conversely, the TNBP presents a different shape, characterized by less homogeneous thermal pattern on 21 August; this allows the polynya to extend further in distance from the coastline but results in a smaller estimated extent. The main difference from the previous day scene is the decrease of the IST observed in several portions of the former homogeneous open polynya area. Possibly, this is due to the formation of thin ice, whose temperature is much closer to the ambient temperature than to the freezing point one. This difference suggests that the local maximum of the polynya area probably occurred between these two scenes.
The aspects identified for the 21 August scene are typical of the closing phase; they are also recognizable in the MODIS thermal images of 27 August (Figure 7). In fact, the closing phase of coastal polynyas is usually associated with a rapid slowdown of the katabatic wind (as visible in Figures 6b and 7b), or with its change of direction, which allows the mixture of open water and frazil ice to progressively refreeze. This process starts from the marginal zone of the polynya offshore, where the smoother transition between the warm polynya core values and the cold ice pack temperatures is observed, toward the coast. This explains why an open polynya area survives near the coast also when katabatic winds stop, i.e., on the 28 August (Figure 7), and suggests that the estimate of the TNBP area from thermal data is more difficult during the closing phase than during the growing one [42]. For this reason, the methodology applied in this study tries to minimize potential errors due to the choice of a thermal threshold that could lead to the inclusion or exclusion of large areas of thin ice, through the pseudo-ship neighborhood analysis of all the pixels characterized by critical IST values (255-265 K), as described in Section 3.4.
The results derived from the analysis of the MODIS TIR products have been validated through the comparison with co-located high-resolution ENVISAT ASAR images of the TNBP. In the SAR imagery polynyas are clearly visible and well defined thanks to an exceptional spatial resolution (up to 5 m). This resolution decreases to 150 m for the ASAR WS images reported in this study. Nevertheless, they allow us to define the shape of the TNBP and to estimate its area with a remarkable precision as for ASAR images reported in Figures 6c and 7c. We can easily distinguish the brighter wind-generated Langmuir streaks, which develop parallel to the wind direction, associated with frazil and pancake ice formation [6,65]. In Figure 6, the evident separation of these structures along the downwind border is an indication that the katabatic wind was very strong when the image was acquired by ENVISAT. This means that the airflow was able to maintain the sea ice parallel bands up to the margin of the polynya, even though its strength obviously decreased rapidly due to surface friction as the wind moved away from the Nansen Ice Sheet. For the same reason, the offshore border of the polynya, which separates this newly formed ice from the thicker compact sea ice, is really sharp on this date (as seen in IST scenes). The nearly parallel bands of homogeneous sea ice that we can see outside the polynya area represent the old polynya borders which were originated by the previous repeated polynya activity, supported by quasi-periodic katabatic events, and then shifted offshore [6,42].
The comparison of TNBP extents retrieved through corresponding MODIS and ASAR products provided very interesting results for both case studies in Figures 6 and 7, with differences of 69 km 2 , 63 km 2 and 125 km 2 , for 20, 28 and 29 August respectively. A more extensive validation of MODIS retrieved estimates for the whole 2010 and 2011 winter season is provided in the following section.

Technical Validation of the MODIS Derived Polynya Extent Estimation
The TNBP area estimated using the methodology presented in Section 3.4 has been verified through comparison with 196 coincident ENVISAT ASAR images acquired during the study periods. To this aim, the polynya area on each ASAR scene, manually digitized according to the contour of the high backscattering area and the presence of the typical Langmuir structures as described above, is assumed to be the true extent in account of its best spatial resolution. Of course, large leads occasionally present in the MODIS images could also be included in the total open polynya area computed by our automated procedure. Nonetheless, we are confident that their limited extents (few tenths of km 2 ) do not alter significantly the efficiency of the estimations.
Some meaningful examples of ASAR-MODIS TNBP extent intercomparison are shown in Figure 8a,b. The comparison is really excellent for two of the selected 2010 examples, those dated to 3 May and 2 August, while a difference of about 600 km 2 emerges for the 16 May. Time delay between MODIS scenes (earlier) and ASAR acquisitions is very short (60 to 90 min) and comparable between images, thus synopticity seems not to be responsible to the differences on 16 May. Similarly, the phase of the observed polynya is not a possible cause because only the 2 May event can be ascribed to a closing polynya event, as suggested by the smoothed open water edges (differently from the sharp IST gradients visible on the other scenes). Figure 8b presents three examples from the 2011 time series. Two of them, 10 and 12 May, give an excellent agreement between MODIS and ASAR estimations during the early and late closing phase of a significant polynya event. Conversely, the third example from 13 October shows a lower correspondence between the two methodologies, with a difference of about 300 km 2 . Again, it is not easy to explain this underestimation but, probably, it is just a result of the polynya growing rapidly during the four hours which separate MODIS acquisition from ASAR one. This hypothesis is supported by MODIS estimations time series (Figure 3b) which reports wider TNBP areas on the next days (14-15 October-days 287-288). Actually, a general underestimation of the TNBP area derived through ASAR imagery emerged from time series in Figure 3, with MODIS retrievals never exceeding significantly the coincident ASAR extents. On the other hand, they seem to provide a reliable approximation of the TNBP, ensuring the high revisit time and the hourly monitoring that SAR and optical products cannot guarantee despite their excellent spatial resolution. Furthermore, as discussed for Figures 6 and 7, it is very interesting to remark that MODIS IST imagery can easily distinguish the actual polynyas, made of open water and recently formed thin ice, from the downwind relatively thicker bands of sea ice, probably created during the previous openings and then made compact by temporary decreases of the airflows from the coast. In this sense, the high capability of this automated methodology is confirmed by the absence of MODIS overestimations of the TNBP area, in comparison to ASAR imagery, due to misinterpreted accumulated sea ice in the outlet zone. Generally, this comparison produced an average relative error of around 6%, with a maximum value of 29% found on 16 May 2010; these values, estimated assuming ASAR as reference, are consistent with previous results reported by [48]. Figure 9a shows a scatterplot of all the near coincident ASAR manual and MODIS automated estimations, which time lags were shorter than 24 h. The MODIS retrievals' occasional underestimation described above is confirmed and a linear correlation coefficient of 0.75 is derived. No significant differences emerge when considering the different ASAR products (i.e., WS, GM and IM). Finally, MODIS daily retrievals for the same dates obtained through the POLA algorithm [27] were also used for an additional comparison. These estimations showed a higher correlation than ASAR to our retrievals (R ≈ 0.83, as shown in Figure 9b); this is not surprising since both MODIS algorithms retrieve TNBP areas through IST products derived by the same dataset (thus having the same spatial resolution). Conversely, a lower correlation is found between POLA and ASAR products (R ≈ 0.69, not shown). Probably, this is due to the fact that POLA estimations represent daily average values. Moreover, POLA includes an automated cloud cover reconstruction procedure [53], thus including several scenes partially affected by atmospheric disturbance which are completely discarded during the preliminary steps of the new procedure presented through this study. Furthermore, no additional manual input or filtering is applied during the POLA computations. In general, POLA overestimations of our retrievals do not match with polynya openings which extension is underestimated by our methodology, as confirmed by Figure 9b and the lower correlation coefficient. Differences seem rather to be associated with occasional offshore bands of thin sea ice or larger leads occurring downwind in the TNBP outlet zone and with intrinsic TNBP variability within each single day (averaged by the POLA algorithm).

Sensitivity Analysis
A sensitivity analysis of the value used as upper threshold to distinguish sea ice from open water is reported in Figure 10. The observed slight underestimation showed by MODIS retrievals in comparison with ASAR imagery (Figure 9a), let the reader suppose that a different setting of the algorithm thermal thresholds could improve the accuracy of our estimates. For this reason, a sensitivity test was performed and TNBP area was estimated for winter season 2010 using higher and lower values for the upper threshold, usually set to 265 K. These test values ranged between 260 K and 266 K, at steps of 1 K. Figure 10. Thermal threshold sensitivity test for winter season 2010. The study period is divided in two intervals: (a) 1 April-4 July; (b) 5 July-31 October. TNBP extent retrievals are represented through different colors when using different thresholds, as described in the legend. Yellow circles report the ASAR Envisat co-located retrievals. Black rectangles present a zoom on days (a) 92-93 (2-3 April); and (b) days 232-233 (20-21 August), the case study fully analyzed in Figure 6.
As expected, results show that higher values, i.e., 266 K (red values in Figure 10), implicate a further underestimation of reference ASAR data. On the other hand, TNBP extents derived through the application of 260 K (black line) lead to a large overestimation, often providing unreliable values. By far, the best choice lies between 264 K (green line) and 265 K (blue line), whose linear correlation is nearly the same (R = 0.74 and 0.75 respectively). Among these we opted for the second one, following what previously suggested by [6]. A clear example of the typical variability linked to the threshold sensitivity is provided by the zooms in Figure 10, where we highlight the information concerning the 20-21 August case study, described in detail in Figure 6.

Conclusions
The detection of polynya events, and the magnitude of their extents, are mandatory for studying the physical processes connected to their openings, i.e., heat flux exchange and ice production [2]. Moreover, reliable satellite polynya area estimations are also necessary for validating and/or initializing ocean-ice-atmosphere models [23] and for better understanding the ocean-atmosphere coupling which is recognized as a leading factor influencing sea ice variability [66]. Hence, the goal of this study was to demonstrate the potential of a new procedure for estimating the wintertime polynya area at a high spatial resolution (1 km) through MODIS TIR imagery. To this aim, the frequently occurring TNBP was observed in more than 1600 MODIS IST clear-sky scenes acquired during the winter seasons (April to October) 2010 and 2011, and its area was computed following the processing scheme summarized in Figure 2. Despite the well-known atmospheric disturbance on TIR data, results confirmed that MODIS acquisitions allow a continuous monitoring of the TNBP showing a revisit time shorter than 24 h (12 h) for more than 89% (74%) of the study period. This also allowed observing a high TNBP area daily variability as shown through several case studies.
A high interannual variability was found for the analyzed winter seasons, with an opening/closing fluctuation clearly depending on the regime of the katabatic winds funneled by TNB coastal glaciers. As expected, short polynya events (2-5 days) were more frequent than long-lasting ones (10-15 days); similarly, small openings up to 1000 km 2 (75%) were more frequent than larger ones, with areas over 4000 km 2 identified only in 1% of the analyzed images. Recurring open water disclosures at the extremity of the Drygalski Ice Tongue were also retrieved in 21% of analyzed scenes. The mean annual extent of the open water area was at lower side of the known range (around 900-1000 km 2 ) but in good agreement with previous visible and infrared satellite observations. A very good agreement was also found with MODIS derived TNBP estimates reported by previous studies [30,48] for sea ice seasons 2010 and 2011.
The reliability of our MODIS estimates was supported by the comparison with a set of 196 contemporaneous high-resolution ENVISAT ASAR images. This dataset was particularly useful for retrieving accurate values of the TNBP area, thanks to the evident contrast between the bright backscattering polynya, typically characterized by the Langmuir structures, and the uniform sea ice pack offshore. Results showed a linear correlation of 0.75 between MODIS and ASAR retrievals, with an average relative error of about 6% computed assuming ASAR as reference. Several scene-to-scene comparisons confirmed the excellent agreement between the ASAR bright backscatter regions and the IST warm cores. Nonetheless, MODIS TNBP retrievals generally underestimate co-located ASAR values during winter 2010 and 2011. Sensitivity tests proved that this does not depend on the setting of the thermal thresholds included in the main step of our procedure. To minimize such bias, the combined use of MODIS derived products, active microwave and optical data would be more advisable, as well as the implementation and the use of specific segmentation schemes to separate different ice conditions; this could lead to a better determination of the extent of polynyas and their outlet zone [6]. However, even though the ESA Sentinel constellation missions (i.e., Sentinel-1 and -2) reduced the temporal gaps between successive data acquisitions, SAR imagery is not yet able to provide the needed revisit time for a complete monitoring of the daily variability of coastal polynyas.
Moreover, MODIS products are available since 2003 allowing the study of polynya variability at 1 km resolution during the last 15 years. In a forthcoming study, we will apply the methodology presented here to analyze the TNBP evolution during the 2003-2017 period in order to identify possible variations in its behavior and eventually discuss them in the framework of global climate changes [66][67][68]. To this aim, a preliminary analysis will be carried out to evaluate the opportunity of using also TIR observations acquired by the MODIS sensor onboard the Terra satellite, not included in this study. Then, we intend to adapt and apply this methodology to other Antarctic coastal polynyas (e.g., Mertz, Ross Sea, Cape Darnley), being conscious that automated algorithms developed for a specific location cannot be directly used globally since sea ice physical properties may differ significantly from region to region [6], as would the oceanography of the region in question (e.g., freshness and stratification of the surface and mixed layers of the ocean). We also plan to use these polynya retrievals to estimate their heat budgets and ice production rates [2,69] and eventually infer thin ice thickness [27] and AABW formation rates [70].