Long-term coastal-polynya dynamics in the southern Weddell Sea from MODIS thermal-infrared imagery

. Based upon thermal-infrared satellite imagery in combination with ERA-Interim atmospheric reanalysis data, we derive long-term polynya characteristics such as polynya area, thin-ice thickness distribution, and ice-production rates for a 13-year investigation period (2002–2014) for the austral winter (1 April to 30 September) in the Antarctic southern Weddell Sea. All polynya parameters are derived from daily cloud-cover corrected thin-ice thickness composites. The focus lies on coastal polynyas which are important hot spots for new-ice formation, bottom-water formation, and heat/moisture release into the atmosphere. MODIS has the capability to resolve even very narrow coastal polynyas. Its major disadvantage is the sensor limitation due to cloud cover. We make use of a newly developed and adapted spatial feature reconstruction scheme to account for cloud-covered areas. We ﬁnd the sea-ice areas in front of the Ronne and Brunt ice shelves to be the most active with an annual average polynya area of 3018 ± 1298 and 3516 ± 1420 km 2 as well as an accumulated volume ice production of 31 ± 13 and 31 ± 12 km 3 , respectively. For the remaining four regions, estimates amount to 421 ± 294 km 2 and 4 ± 3 km 3 (Antarctic Peninsula), 1148 ± 432 km 2 and 12 ± 5 km 3 (iceberg A23A), 901 ± 703 km 2 and 10 ± 8 km 3 (Filchner Ice Shelf),


Introduction
Coastal polynyas are recurring areas of thin ice and open water which are generally formed by divergent ice motion, e.g., by strong offshore winds or ocean currents. During winter, polynyas are important hot spots for ice production (IP), deep-water formation and gas ventilation of the ocean (e.g., Morales Maqueda et al., 2004).
Due to the dependency of ocean-to-atmosphere heat flux on the ice thickness (Maykut, 1978) and the nature of thermodynamic ice growth by release of latent heat at the ice/ocean interface, high rates of ice production occur predominately under thin-ice/frazil-ice areas as well as in openwater areas during austral winter (e.g., Tamura et al., 2011Tamura et al., , 2007Renfrew et al., 2002;Thorndike et al., 1975). The knowledge of ice production in polynyas is of large interest. Tamura et al. (2008) estimate that about 10 % of the Antarctic sea ice is produced within the major polynyas.
In this study, the geographical focus lies on the southern Antarctic Weddell Sea region between 69-78 • S and 65-20 • W (Fig. 1). With an average maximum winter extension of about 4.75 × 10 6 km 2 , the Weddell Sea has the largest fraction of the Antarctic sea-ice cover (Drucker et al., 2011). While the Weddell Sea features polynya activity along almost its complete coastal area, the majority of studies focus only on selected subregions. The investigation area was limited spatially to the southern Weddell Sea, since this area is considered to be a key region for the formation of high-salinity shelf water, and observational reference data of polynya dynamics are needed for sea-ice/ocean models (Haid et al., 2015).
The majority of recent studies dealing with polynya dynamics in this region make use of passive-microwave sensors. In comparison to thermal-infrared sensors, passive- microwave sensors retrieve data with no limitations due to day-/nighttime or clouds. Drucker et al. (2011) use SSM/I and AMSR-E data to derive ice-production rates and estimate ice export rates for the Weddell and Ross seas in the years from 1992 to 2008. Kern (2009) derives pan-Antarctic coastal polynya areas (POLA) from SSM/I data for the same time period. Tamura et al. (2011) use SSM/I data to estimate surface heat and salt fluxes that are associated with sea-ice growth/melt in the Southern Ocean. Nihashi and Ohshima (2015) conduct a combined circumpolar mapping of polynya area and fast ice as well as ice production based on passivemicrowave as well as thermal-infrared remote sensing.
In this study, we present long-term results of coastalpolynya dynamics in the southern Weddell Sea that are derived from Moderate-Resolution Imaging Spectroradiometer (MODIS) thermal-infrared imagery. Remote sensing of sea ice using thermal-infrared data yields the opportunity to monitor thin-ice thicknesses (TITs) and distribution on a regular basis (Adams et al., 2013;Willmes et al., 2010;Yu and Rothrock, 1996). While the MODIS cloud mask (MOD35) shows in general good results for polar daytime (e.g., Frey et al., 2008;Liu et al., 2004), its performance decreases during polar nighttime (e.g., Holz et al., 2008;Frey et al., 2008;Riggs et al., 2006;Hall et al., 2004). This results from the lack of visible channels, as well as the low thermal contrast between clouds and the underlying snow/ice (Liu and Key, 2014;Frey et al., 2008;Liu et al., 2004). We therefore use additional criteria to identify cloud-covered areas and use a new approach for gap filling with respect to thin-ice thickness (Paul et al., 2015a). While the basic thin-ice retrieval was used before in Arctic regions (e.g., Preußer et al., 2015;Adams et al., 2013;Willmes et al., 2011Willmes et al., , 2010, we present now the first continuous and cloud-cover corrected time series of polynya dynamics during the austral winter period (April to September) and for the complete coastal area in the southern Weddell Sea for years from 2002 to 2014.

MODIS data (MOD/MYD29)
In this study we make use of the MODIS Sea Ice product (MOD/MYD29, Riggs et al., 2006;Hall et al., 2004), which is a swath-based product that is derived from both MODIS sensors on-board the NASA Aqua and Terra satellites and distributed by the US National Snow and Ice Data Center (NSIDC). Our work is based on the ice-surface temperature data set (IST), which is one part of the sea-ice product. The provided data have a spatial resolution of 1 km × 1 km at nadir and each swath covers an area of 1354 km (across track) × 2030 km (along track). The spatial resolution decreases off nadir to about 2 km × 5 km at the swath edge.
The overall IST accuracy of the MOD/MYD29 sea-ice product under ideal (i.e., clear-sky) conditions is 1-3 K (Riggs et al., 2006;Hall et al., 2004) and is derived based on a constant emissivity for snow/ice (Hall et al., 2015). The MOD/MYD29 data product was already corrected for cloud cover using the MODIS cloud mask. However, the performance of the MODIS cloud mask is drastically reduced during polar nighttime due to the lack of visible channels and the low thermal contrast between clouds and the underlying snow/ice (Liu and Key, 2014;Frey et al., 2008;Liu et al., 2004). We used these data for the austral winter period from 1 April to 30 September for the years from 2002 to 2014 that comprise a total of 78 696 MODIS swaths which cover the southern Weddell Sea area. This equals an average of 33 single swaths per day.

AMSR2 data -ARTIST sea ice (ASI)
For comparison, we use daily sea-ice concentration data derived from the Advanced Microwave Scanning Radiometer 2 (AMSR2) provided by the University of Bremen (Beitsch et al., 2013;Spreen et al., 2008). The ASI algorithm estimates sea-ice concentration from polarization differences in brightness temperatures obtained from the AMSR2 89 GHz channel. Since 1 August 2012, ASI seaice concentrations are available with a spatial resolution of 3.125 km × 3.125 km.

European Centre of Medium-Range Weather Forecast (ECMWF) ERA-Interim reanalysis
The set of input data is complemented by the ECMWF ERA-Interim atmospheric reanalysis data (Dee et al., 2011). The ERA-Interim data are distributed at a spatial resolution of 0.75 • and a temporal resolution of 6 h (at 00:00, 06:00, 12:00, and 18:00 UTC). For our analysis we used the 2 m air temperature, the 10 m wind-speed components, the mean sea-level pressure, and the 2 m dew-point temperature. These data sets supply all necessary input to calculate the energy balance for the thin-ice retrieval.
In addition to these meteorological variables, cloud-cover information was also taken from the ERA-Interim reanalysis data. In our case, we used the medium-level cloud-cover fraction data, which were found to show the in general best agreement with the MODIS cloud mask (Liu and Key, 2014). This data set is used to complement the MODIS cloud mask during nighttime and was also spatially and temporarily interpolated to fit each MODIS swath.

Thin-ice retrieval
The derived TIT of up to 0.5 m is calculated by using a surface energy balance model (Adams et al., 2013). This model utilizes the fact that the ice-surface temperature of thin ice is related to ice thickness via the energy balance equation (Drucker et al., 2003;Yu and Rothrock, 1996).
Since the total atmospheric flux (i.e., the energy loss to the atmosphere Q atm , Eq. 1) is balanced by the conductive heat flux through the ice (Q i , Eq. 2), the ice thickness (h i ) can be calculated from Eq. (3). Here, Q 0 is the net radiation balance, H 0 and E 0 are the turbulent fluxes of sensible and latent heat, κ i is the thermal conductivity of sea ice (2.03 W(m K) −1 ) as well as T s and T fp which are the ice-surface temperature and the ice/ocean-interface temperature, respectively.
Necessary assumptions comprise thin and snow-free sea ice, a linear temperature gradient throughout the ice, and the temperature at the ice/ocean interface to be constant and at freezing point of sea water (T fp ). In addition, all processing is carried out for nighttime only, to avoid dealing with shortwave radiation and snow/ice-albedo effects. We therefore removed daytime pixels based on their solar-incidence angle. This was done on a per-pixel basis for all available MODIS swaths between 1 April and 30 September for the years from 2002 to 2014. Thereby, we also simplify the net radiation to longwave radiation terms. The turbulent fluxes of sensible and latent heat as well as the long-wave radiation balance are calculated following Adams et al. (2013) using the ERA-Interim atmospheric reanalysis data sets. For the flux calculation, the iterative stability-dependent bulk approach of Launiainen and Vihma (1990) is used. The long-wave radiation balance is calculated using an improved parametrization by Jin et al. (2006) for the downward long-wave radiation and the Stefan-Boltzman law for the upward long-wave radiation.
Studies of the ERA-Interim downward long-wave radiation data revealed that ERA-Interim overestimates the cloudcover fraction of low clouds (Zygmuntowska et al., 2012). We found this effect also to appear in clear-sky situations diagnosed from visual screening of the corresponding MODIS imagery for the Antarctic. Overestimated cloud-cover fractions lead in turn to too-high downward long-wave radiation values (Tastula et al., 2013). Tastula et al. (2013) found a large positive bias in incoming long-wave radiation for ERA-Interim in their analysis of different atmospheric reanalysis data sets over Antarctic sea ice. Hence, we decided against the use of long-wave radiation data from Era-Interim reanalysis.
The exact procedure of the thin-ice thickness retrieval as well as all equations are thoroughly described in Adams et al. (2013). For ice thicknesses between 0.0 and 0.2 m, Adams et al. (2013) state an average uncertainty of ±4.7 cm. Icethickness data above a 0.2 m threshold yield substantially higher uncertainties (Adams et al., 2013), which is why we limit our investigation to a maximum ice thickness of 0.2 m.

Area of investigation and general procedure
For a better comparison to other studies, the southern Weddell Sea coastal area was divided into six subregions ( Fig. 1), namely from west to east: Antarctic Peninsula (AP), Ronne Ice Shelf (RO), the area around the grounded iceberg A23A (IB), Filchner Ice Shelf (FI), Coats Land (CL), and the Brunt Ice Shelf (BR).
All MODIS swaths were processed in several steps as visualized in Fig. 2.
First, the swaths were projected onto a common equirectangular grid with an average spatial resolution of 2 km × 2 km using a nearest-neighbor approach. The spatial resolution of this type of grid decreases slightly with decreasing latitude and vice versa. No interpolation between projected MODIS pixels was applied (i.e., data gaps are possible). Subsequently, ERA-Interim data, which exhibit a much coarser spatial (0.75 • ) and temporal (6-hourly) resolution than the MODIS data, were linearly interpolated to spatially and temporally fit the MODIS data on the common reference grid.
Based on these adjusted swath-based data sets, TIT is calculated pixel-wise and daily TIT composites are calculated based on the median thin-ice thickness and corresponding ice-surface temperature of all available swaths per pixel. The resulting composites comprise TIT and IST data together with the daily swath-based median energy-balance compo- Work flow based on the described input data and thinice retrieval. Established processing steps that were already used in other studies are marked in blue, while new additions are marked in red. These additions comprise a persistence index (PIX), the newly defined combined cloud-cover information of MODIS cloud mask and ERA-Interim medium-level cloud cover. Missing data due to cloud cover are accounted for by a two-step procedure that comprises spatial feature reconstruction (Paul et al., 2015a) and proportional extrapolation . The combined cloudcover information separates between confident clear-sky (ccs), definite cloud-cover (dcc), mixed-covered pixels (mcp), and uncovered pixels (ucp).
nents of each thin-ice pixel. In the next step, cloud-covered data were identified and flagged as will be described in the following subsections.

Identification of cloud-cover affected data
As already mentioned, the MODIS cloud mask shows deficiencies in the nighttime cloud-cover detection. We therefore use ERA-Interim medium-level cloud-cover information to reduce the influence of cloud-affected pixels in the aggregation process of the daily TIT composites. We apply a simple cloud-cover threshold (75 %, Fig. 3) to identify cloudcovered regions from ERA-Interim data and thereby provide an additional indicator to the MODIS cloud mask for an identification of potentially cloud-affected pixels. The effect of different thresholds is shown in Fig. 3a-d. The binary cloud-cover information (i.e., depending on whether a pixel is cloud influenced or not) of both the MODIS cloud mask as well as the ERA-Interim medium-level cloud-cover fractions is pixel-wise aggregated into daily composites and separated into four different cloud-cover dependent classes (Fig. 2). Pixels that are considered as clear sky by means of the MODIS cloud mask as well as the ERA-Interim cloudcover information for all available swaths per day are assigned to the confident clear-sky class (ccs). In contrast, pixels that feature cloud cover in both data sets for all available swaths are assigned to the definite cloud-cover class. All pixels that show a mix of cloud-covered and clear-sky conditions are assigned to the mixed-covered pixel class (mcp) with a clear-sky index based upon the ratio of the number of clearsky input swaths to the total number of input swaths. Pixels that are not covered by any swath during 1 day are assigned to the uncovered pixel class.
On average, about five swaths cover each pixel in a daily composite. This number of swaths can vary mainly due to cloud cover with up to 20 different swaths covering a single pixel or region.
In order to reduce the inherent misclassification problems of the MODIS cloud mask with clouds taken for either thin or thick ice, we introduced an additional procedure that assigns a persistence index (PIX; Figs. 2 and 3f). Hereby, we make use of the temporal and spatial persistence of the temperature signal in polynyas in contrast to mobile cloud patterns. The persistence index is the ratio of the number of swaths per pixel per day that contain thin ice (TIT ≤ 0.2 m) to the total number of clear-sky acquisitions covering that pixel.
Both the number of swaths and the persistence index are then used to further reduce the effect of undetected cloud cover on later processing steps (Fig. 3g).

Correction scheme for cloud-covered data
Employing the above mentioned two procedures of cloudcover dependent classification and persistence index calculation, we are able to identify cloud-contaminated MODIS data in each daily TIT composite. These are then complemented by a two-step procedure utilizing a combination of spatial feature reconstruction (SFR; Paul et al., 2015a) and proportional extrapolation (PE; Preußer et al., 2015).
In the SFR approach, the information of a 7-day interval (doi −3...−1 and doi 1...3 ) centered around the day of interest (doi) is weighted according to its temporal proximity to the initial day of interest. This yields a probability of thin-ice occurrence for the day of interest based on the surrounding 6 days (p TIT , Eq. 4).
Information about polynya area is on average significantly correlated within at least 3 days and > 90 % per-pixel gaps are shorter than 4 days (Paul et al., 2015a). We use the set of weights (w 3 = 0.02, w 2 = 0.16, w 1 = 0.32) and probability threshold (th = 0.34) of this study, which featured the highest spatial correlation in the analysis of Paul et al. (2015a). A detailed description and analysis of the SFR approach and its setup is given in Paul et al. (2015a).
The PE approach on the other hand assigns thin ice to cloud-covered areas in the same proportion as it is detected in the cloud-free area. For example, if a region is 80 % cloud free and 50 % of the cloud-free area features a thin-ice signal, then 50 % of the cloud-covered region is considered as thin ice.
In the two-step procedure, we first apply the SFR approach to all cloud-free pixels (i.e., pixels in the ccs class and mcp class where the majority of pixels show clear-sky conditions) that also feature a PIX value greater than 0.5 (i.e., pixels with a thin-ice thickness ≤ 0.2 m present in more than 50 % of the swaths covering that pixel) in the 7-day interval and are covered by at least three swaths. Pixels that do not match these two criteria are considered of minor quality and are discarded from the SFR approach.
Based on Eq. (4), the binary information (i.e., thin ice or no thin ice) of the 6 days surrounding the day of interest is matched with their corresponding weights and the polynya probability estimated. A probability value above 0.34 is classified as polynya area in a resulting binary image.
Subsequently, cloud-contaminated pixels with a probability above the threshold (0.34) are assigned a pixel-wise weighted average ice thickness and ice-surface temperature value based on the 6 days surrounding the initial doi (TIT add /IST add , Fig. 2). Weights are applied in the same configuration as for the SFR approach itself. Additionally, pixels categorized as minor quality previously but feature a probability above the threshold of 0.34 in the SFR output are "upvalued" from minor to high quality and are assigned their original thin-ice thickness and ice-surface temperature values.
The remaining coverage gaps that could not be corrected for by this approach, e.g., due to temporal gaps longer than 3 consecutive days, are filled by the PE scheme .
In the case that after the application of the SFR approach more than 50 % of the investigated subregion is cloudcontaminated, daily estimates of polynya area and ice production will be interpolated between neighboring days with sufficient (i.e., above 50 %) cloud-free coverage.

Derived output polynya parameters
From our cloud-cover corrected daily thin-ice thickness composites, we then derive daily polynya area (defined as area with open water and thin-ice between 0.0 and 0.2 m thickness) as well as the accumulated wintertime ice production from heat loss for each POLA pixel (e.g., Tamura et al., 2011;Willmes et al., 2011).
Assuming that the complete heat loss through the ice (Q i , Eq. 5) over 1 day contributes to new ice formation ( ∂h ∂t ), the volume ice-production rate ∂V ∂t can be calculated by multiplying ∂h ∂t with the area of each POLA pixel. In Eq. (5), we assume a sea-ice density (ρ i ) of 910 kg m −3 and the latent heat of fusion (L f ) for sea ice to be 3.44 × 10 5 J kg −1 . All estimates of volume ice production are integrated over the period of 24 h, assuming constant ice-thickness and ice-surface temperature. Hence, we are not able to account for sub-daily variations with this procedure.
Furthermore, thin-ice thickness distributions of daily POLA as well as frequency distributions of thin-ice occurrence are calculated for each subregion. The results are then put in context with other recent remote sensing and model studies.

Inter-comparison MODIS/AMSR2
A threshold of 70 % was used for the definition of a polynya pixel from AMSR2. The comparison of frequency of polynya pixel occurrence from AMSR2 (Fig. 4a) and MODIS (thinice thicknesses equal or below 0.2 m, Fig. 4b) for the years 2013 and 2014 reveals an overall good agreement. Both data sets were projected onto a common equirectangular reference grid, similar to the one used in the general data processing procedure but with a coarser average resolution of approximately 3.5 km × 3.5 km.
The difference in frequencies is shown in Fig. 4c. The majority of pixels lie in a ±10 % range, and there are positive as well as negative differences. Along the coastline of Filchner Ice Shelf and Coats Land, MODIS shows overall higher frequencies, whereas especially close to the ice-shelf edge of Ronne Ice Shelf and Brunt Ice Shelf AMSR2 shows higher values of polynya occurrences. However, effects of a different land/ice-shelf mask in AMSR2 data compared to MODIS are visible (e.g., at Brunt Ice Shelf). As we will discuss later in the study, there are difficulties in discriminating between fast ice, adjacent icebergs, ice shelves, and thin ice in passivemicrowave data (Tamura et al., 2008). Additionally, this effect of a per-pixel shift between both data sets may partly result from the applied nearest-neighbor forward projection when transferring the data to the common reference grid. The small underestimation by MODIS near the Ronne Ice Shelf may also be caused by persistent cloud cover in that area, particularly during the presence of a polynya.

Cloud-cover corrections
The spatial coverage of cloud-free MODIS data without a cloud-cover correction is roughly between 40 and 100 % due to the large amount of daily swaths that at least partially cover a large region like the southern Weddell Sea. The use of the spatial feature reconstruction approach (Paul et al., 2015a) can mitigate the effect of cloud cover to a great extent yet not completely. At the same time, this is done with high confidence that no thin ice is randomly added in areas where a polynya is unlikely to appear, which can potentially happen when solely using proportional extrapolation . However, to achieve full coverage the additional use of PE is necessary.
An average of 78 % in all MODIS composites per day over all six subregions is classified as cloud free. After the application of the SFR approach, we achieve an average coverage of 97 % per day, where the remaining 3 % get corrected for by PE.
When comparing our two-step cloud-cover correction procedure to the sole use of proportional extrapolation (Table 1) we find that the latter yields generally higher correction estimates for the three smaller regions (IB, FI, and CL) as well as for the Antarctic Peninsula region. However, this observation is reversed for the highly active regions in front of the RO and the BR. Here, the sole use of PE results in lower estimates. Especially in spatially large regions with low polynya activity like the Antarctic Peninsula (Fig. 5) sole use of PE reveals its disadvantage by more than doubling the uncorrected value (Table 1). Whenever no coverage with cloud-free information above 50 % after the use of SFR can be achieved, the POLA information of that day will be interpolated between days with sufficient (i.e., above 50 %) coverage. Under the same premise, the PE approach is used for comparison to the twostep procedure (SFR + PE; Table 1). However, when solely applying the PE approach, interpolation is necessary on average on 34 occasions per year and region due to present cloud coverage in the MODIS composites. With the use of the SFR approach, interpolation is on average only necessary on 2 days per year and region. The use of SFR potentially leads to a decrease uncertainties, especially when investigating single days.

Thin-ice frequency and thickness distribution
Almost the complete coastal area in the south and east of the investigated southern Weddell Sea features a recurrent thinice signal for the years from 2002 to 2014 (around 25 % and above with extremes of up to 80 % for some years, Fig. 5). The overall high recurrence of thin ice is a very important finding when considering that the primary focus of many studies lies solely on the RO region when investigating the Weddell Sea (e.g., Nihashi and Ohshima, 2015). This neglects the importance of e.g., the BR region to the overall ice production and ocean-atmosphere heat exchange. It also underestimates the interannually highly variable contribution of the area around the grounded IB, FI, and CL to the bottomwater formation due to salt release during ice formation.
The relative frequency of thin-ice occurrences in the AP region is spatially focused around smaller grounded icebergs and rather low compared to the other subregions. The interannual contribution also decreases during our investigation period in the years 2007 to 2010 (Fig. 5f-i), when the group of icebergs detach from the ground.  Preußer et al., 2015) and a combined step-wise approach using both, spatial feature reconstruction (SFR; Paul et al., 2015a), and proportional extrapolation. Values are presented for each of the six subregions (Fig. 1)  The very light blue areas in the overall thin-ice frequency distribution (Fig. 5n) found in the northeast of the investigation area correspond to a low sea-ice concentration area due to the far south lying marginal ice zone in April in the years 2005, 2006, and to a smaller extent in 2007 near the Brunt Ice Shelf region (at around 75 • S, taken from AMSR-E observations, not shown). This is also reflected in our thin-ice frequency distributions (Fig. 5d-f).
There is a sharp separation between a zone with present activity in the north and northeast to a zone with almost no activity in the south and southwest closer to the Ronne and Filchner ice shelves in the total frequency distribution (Fig. 5n, A). This line of separation (Fig. 5n, A) coincides very well with the course of the continental slope ( Fig. 1; Arndt, 2013). However, no conclusive explanation can be given. Markus (1996) studied the effect of the grounded iceberg A23A on the ice production in front of the Filchner Ice Shelf and found a drastic average increase in sea-ice concentration (i.e., a decrease in the amount of thin-ice area and hence ice production) during the freezing period. While the recurring formation of a fast-ice bridge is also visible in our results (a low frequency area between the coast and the grounded iceberg A23A; Fig. 5n, B), we still find high polynya activity in front of CL and FI. Especially the area around the grounded IB and the westward side of the ice bridge show high thin-ice occurrences throughout all years (Fig. 5a-m). It is also possible to see differences in extent and persistence of this ice bridge from our data set via the absence of thin-ice occurrences.
In comparison to results found by Nihashi and Ohshima (2015), our frequency estimates for the RO region between 30 and 40 % are in good agreement.
The overall distribution of thin-ice thicknesses for the southern Weddell Sea (Fig. 6) does not show strong regional nor interannual differences as can be seen from the low standard deviations. However, the standard deviations are larger for the smallest and largest thickness classes. In general, the thin-ice thickness distribution features a steep increase of occurrences from smaller to larger thickness classes. Thin−Ice Thickness Class [cm] Figure 6. Averaged distribution of thin-ice thicknesses in the southern Weddell Sea over the whole 13-year time interval between April and September. The class width is 2 cm; i.e., the 1 cm class includes thicknesses between 0 and 2 cm. Error bars represent ± 1 standard deviation.

Polynya area and ice production
Results for POLA and IP for the years from 2002 to 2014 between 1 April and 30 September are presented in Figs. 7 and 8, respectively. The investigated six subregions can be divided into two classes: the highly active Ronne Ice Shelf (Fig. 7b) and Brunt Ice Shelf (Fig. 7f) regions with also a large contribution in ice production ( Fig. 8b and f)  Generally, by means of multi-year trends over the course of 13 years in POLA and IP (Figs. 7 and 8), all regions except CL exhibit a negative trend for both POLA and IP. However, not all of those trends are statistically significant. AP shows a significant trend (α = 5 %) with approx. −57 km 2 a −1 and −0.54 km 3 a −1 , respectively. The region around the grounded IB also features a multi-year negative significant trend in POLA with approx. −66 km 2 a −1 .
Another important fact that is shown by our data set is the importance of the FI and CL region. Together with the area around the grounded IB, each shows a higher average POLA and IP than the AP region. However, the interannual variability is very high, especially in the IB and FI region. The combined average POLA of IB, FI, and CL is in the same order of magnitude as the larger contributors of the Ronne and Brunt ice shelves (Figs. 7 and 8). This also holds for the ice production. Without the IB region, FI and Cl show an average polynya area of 1400 km 2 and accumulated wintertime ice production of 15 km 3 .
The frequency of days exceeding a certain polynya area threshold is shown in Fig. 9. This underlines again the high activity in the Ronne (Fig. 9b) and Brunt (Fig. 9f)  1000 km 2 is present in front of the Brunt Ice Shelf. The same is true for 50 % of the investigated days in the Ronne Ice Shelf region and more than 25 % of the investigated days for the area around iceberg A23A. A threshold of 8000 and 9000 km 2 is exceeded in 10 % of all observed days from 2002 to 2014 (April to September) in the Ronne Ice Shelf and Brunt Ice Shelf regions, respectively.
Annual average polynya days (Fig. 10) reveal additional information on the temporal distribution of polynya dynamics (Fig. 9). For example, appearances of POLA values above 10 000 km 2 for the Filchner Ice Shelf (Fig. 10d) region are mainly focused on the year 2003. Here, about one-third of found polynyas feature a polynya area above 10 000 km 2 . The interannual variability in polynya-size distribution per day is also quite high in other regions, especially in the RO and BR regions as well as in the area surrounding the grounded IB. Years with much higher than average annual POLA values (Fig. 7) also show a high amount of days with polynya days in the > 10 000 km 2 class.
The interannual winter average POLA and the wintertime accumulated IP are compared to the results from Timmermann (2013) in Figs. 11 and12. Haid andTimmermann (2013) use the Finite-Element Sea-ice-Ocean Model (FESOM) coupled with a dynamic-thermodynamic sea-ice model to derive polynya area and ice-production rates. Their data set covers the austral winter period from May to September for the years from 1990 to 2009. Haid andTimmermann (2013) used NCEP reanalysis data (Kalnay et al., 1996) with a spatial resolution of 1.875 • as their atmospheric forcing and defined areas with either a sea-ice concentration below 70 % or an ice thickness below 0.2 m as POLA. Haid and Timmermann (2013) calculate sea-ice production from the release of latent heat but also take the oceanic heat flux into account. The model resolution is 3-5 km near the Ronne Ice Shelf. Due to different investigation intervals and subregions, we excluded the April data and focused on the AP, RO, and the BR for the period from May to September in compliance with the model estimates from Haid and Timmermann (2013).
For the AP region, the model estimates are consistently higher than our results (Figs. 11 and 12). A general reason for this could be the use of the NCEP reanalysis data (Kalnay et al., 1996). Especially the narrow mountain range of the Antarctic Peninsula is not sufficiently represented in the atmospheric reanalysis data. Compared to high-resolution simulations, NCEP data show much higher wind speeds at the Antarctic Peninsula (Haid et al., 2015). Topographic effects on the wind such as katabatic winds and barrier winds influence a broader region due to the smoothed topography in NCEP. Hence, wind-induced coastal polynya formation cannot be reproduced, particularly in the Antarctic Peninsula region. This presumably leads to overall unrealistically high POLA estimates and hence ice-production rates in the study from Haid and Timmermann (2013) for the AP region.
Focusing on polynya-area comparison for the two subregions RO and BR (Fig. 11b and c), the presented comparison shows that the our results exceed the model estimates every year with the exception of 2003 for the Ronne Ice Shelf region (Fig. 11b). However, the extent of that difference can change a lot between each year and subregion. For the years from 2004 to 2009 in the Brunt Ice Shelf region, while our results again exceed the model estimates, the pattern of years with relative high/low POLA is equal between our results and the model estimates from Haid and Timmermann (2013). For the remaining years, no apparent similarities are present.
By means of ice production, estimates for the RO as well as the BR region are in the same order of magnitude and oftentimes very similar (with the exception of the years 2002 and 2003 for RO; Fig. 12). Generally, the model estimates are expected to be smaller than comparable satellite-derived results due to the neglect of the oceanic heat flux in our method. A positive oceanic heat flux reduces sea-ice production and at the same time potentially increases polynya area due to an increased amount of thinner ice. Hence, neglecting the oceanic heat flux increases ice production and potentially reduces polynya area in our satellite observations. However, Haid and Timmermann (2013) show that the oceanic heat flux is only of the order of 10 % of the total atmospheric heat flux.
One can only speculate about the reasons for the increased average annual POLA and IP in the years 2002 and 2003 as there is no statement in the study of Haid and Timmermann (2013). Potentially, a mix of below-average air temperatures in combination with above average wind speeds caused this increase in ice production with at the same time increased polynya size/persistence.
Assuming that the satellite-based estimates represent the upper limit for ice production, the relatively high model estimates (Haid and Timmermann, 2013) result from a larger polynya area and a higher ice production per polynya area compared to our results, which likely depends on the computation of the atmospheric heat fluxes.
In the study of Tastula et al. (2013) ERA-Interim and NCEP data are compared to observational data of a drifting ice station in the Weddell Sea for 3 months in 1992. Their results indicate a warm bias of about 4 K for NCEP and a warm bias of about 2 K for ERA-Interim. With respect to wind speed, NCEP shows a positive bias (about 1 ms −1 ) while ERA-Interim features a negative bias (about −0.75 ms −1 ). Different atmospheric forcing can be one of the reasons for the increased POLA and IP in the years 2002 and 2003 for the Ronne polynya in the study of Haid and Timmermann (2013). As shown by Haid et al. (2015), the atmospheric forcing has a strong impact on ice production not only with respect to bias but also with respect to regional distribution. In addition, the method for the computation of turbulent atmospheric fluxes is different between our method and Haid and Timmermann (2013) ues particularly of the sensible heat flux for the sea-ice/ocean model. A comparison between different studies based on satellite observations as well as models is difficult due to different investigation periods, diversity in spatial extent of used POLA masks as well as varying methods and changes in atmospheric forcing. For the following comparison of our results with several recent studies and investigations of the southern Weddell Sea (Tables 2 and 3), we tried to match the used investigation period and regional extent in the best possible way.
Based on average wintertime polynya area (Table 2), we find the best agreement with POLA estimates of Nihashi and  (5665)  Ohshima (2015). The estimates are based on the period from May to August for the years 2003 to 2011 and are acquired from a combination of AMSR-E 89 and 36.5 GHz data. A general problem with passive-microwave data is that it can be problematic to distinguish between the signals from adjacent icebergs, fast ice or ice shelves, and the thin-ice signal of coastal polynyas (Tamura et al., 2008). The effect of mixed pixels causes passive-microwave approaches to potentially misinterpret polynya signals. Due to the generally rather coarse resolution, algorithms based on passive-microwave sensor data can only resolve relatively large polynyas while fairly narrow coastal polynyas, like those in the Weddell Sea, may remain partly undetected (Tamura et al., 2007).
Kern (2009) derived polynya areas from SSM/I data for an investigation period from April to September for the years 1992-2008. While our results in Table 2 only present 4-year averages compared to the much longer averaging period in Kern (2009), results for the BR area are in good agreement. For the RO area, our results are lower than the estimates by Kern (2009). However, Kern (2009) summarizes our subregions of RO and FI as well as parts of IB, CL, and AP, which makes comparisons in this region very difficult compared to the BR region.
Our multi-year averages exceed the results from Haid and Timmermann (2013) in the Ronne Ice Shelf region and especially in the Brunt Ice Shelf region (Table 2)  Antarctic Peninsula region our results are much lower than the model estimates. The good agreement in results between our study and Nihashi and Ohshima (2015) and the differences to Haid and Timmermann (2013) are also potentially related to different atmospheric forcing (Haid et al., 2015). The comparison of ice production between ours and other recent studies in the southern Weddell Sea (Table 3) reveals that the model results (Haid and Timmermann, 2013) for the Antarctic Peninsula are again much higher than our es-timates. In the remaining two regions, the MODIS results exceed the model estimates in the BR region but are lower in the RO region, which results from the differences in IP for the years 2002 and 2003. Drucker et al. (2011) found in their study based on AMSR-E 12.5 km × 12.5 km 36 GHz data for the years 2003 to 2008 and for a slightly longer investigation period from April to October an average cumulative wintertime ice production of 99 km 3 for the Ronne Ice Shelf, 112 km 3 for the Brunt Ice Shelf, and 30 km 3 for the region around the grounded iceberg A23A. These estimates are about 3 times larger than our results in ice production for each of the three investigated regions (Table 3). This effect is very unlikely to be solely based on the additional month used in the investigation period. However, in the study by Nihashi and Ohshima (2015) a comparison between ice production from AMSR-E and SSM/I with regard to the spatial resolution of the sensors showed that lower spatial resolution oftentimes goes along with higher estimates of ice production. Estimates from SSM/I are likely to exceed the results based on AMSR-E data with a much higher resolution (12.5 km 2 vs. 6.25 km 2 , respectively). In the Ronne Ice Shelf region, the difference between SSM/I and AMSR-E is especially large (42 km 3 for AMSR-E compared to 71 km 3 for SSM/I; Ohshima, 2015, based upon Tamura et al., 2011). Given an increase in resolution due to the here-used MODIS data compared to the low resolution of the AMSR-E 36 GHz channel used by Drucker et al. (2011), this might explain our overall lower results to some extent. Additionally, Drucker et al. (2011) stated an error of their results of ±30 %.
Similar to the POLA results, the most recent study by Nihashi and Ohshima (2015) for the years 2003 to 2010 again shows the best agreement with our results (Table 3). While their study is also based on a slightly longer wintertime investigation period from April to October, the increase in ice production for October is more in line with expectations compared to Drucker et al. (2011). This is potentially related to the combined use of AMSR-E 89 and 36 GHz data, utilizing the increased spatial resolution of the 89 GHz data.

Energy-balance components
The mean wintertime area-weighted atmospheric fluxes are shown in Fig. 13 based upon the here-used three components of turbulent fluxes of sensible (H ) and latent heat (E) as well as the long-wave radiation balance (L * ). Shortwave radiation was omitted due to the stated limitation to only use nighttime MODIS swaths for the thin-ice retrieval. The estimates from Haid and Timmermann (2013) were adjusted accordingly and also only show these three components.
As expected, we find the sensible heat flux to be the largest contributor to the total atmospheric flux (e.g., Maykut, 1978). The interannual variability is rather low for the BR and the RO region, while it is comparatively pronounced for the AP. Regional differences are present, but they are rather low in comparison to the model estimates from Haid and Timmermann (2013). Overall, the results of Haid and Timmermann (2013) exceed ours, especially in the RO region. The discrepancy in the Antarctic Peninsula region is most likely again related to the use of NCEP atmospheric reanalysis data (Kalnay et al., 1996).
We stated earlier that the model estimates seem to achieve a higher ice production per polynya area compared to the MODIS-derived results. This is potentially related to large differences between the used atmospheric reanalysis data of NCEP (Kalnay et al., 1996) compared to ERA-Interim (Dee et al., 2011). While the found differences in ERA-Interim and NCEP might be the primary cause for the differences in the energy-balance components (as stated in the previous chapter based on Tastula et al., 2013), a clear explanation cannot be given. Renfrew et al. (2002) found in their study for the years 1992-1998 and using an adaptive, in general much longer freezing period (roughly from February to November) a contribution of 63, 22, and 15 % for the sensible, latent, and radiative components of the atmospheric heat flux, respectively. The average contribution to the total atmospheric flux amounts to 56, 7, and 37 % for the sensible, latent, and longwave radiative components in our results. In our data as well as in the results from Haid and Timmermann (2013), the long-wave radiation shows an on average higher contribution to the total atmospheric flux than in the study of Renfrew et al. (2002). Our average freezing season area-integrated energy exchange with the atmosphere amounts to 7.04 × 10 18 J compared to 3.48 × 10 19 J in Renfrew et al. (2002). However, Renfrew et al. (2002) used a much longer freezing period. Additionally, there is no temporal overlap between our study and the work of Renfrew et al. (2002), which makes an intercomparison very difficult.

Summary and conclusion
In this study, we present a data set of MODIS-derived polynya-area estimates and ice-production rates as well as thin-ice frequency distribution, thin-ice thickness distribution, and energy-balance components for the southern Weddell Sea, Antarctica. This was done for a 13-year time interval (2002 to 2014) during the austral winter period from April to September for the complete coastal area separated www.the-cryosphere.net/9/2027/2015/ The Cryosphere, 9, 2027-2041, 2015 into six subregions. For that, we utilized the higher spatial resolution of MODIS compared to passive-microwave sensors such as SSM/I and AMSR-E/AMSR2. The addition of a more strict exclusion of cloud-covered data using ERA-Interim data to the established thin-ice retrieval as well as the adaptation of a new approach to compensate for cloudcovered areas in daily MODIS composites is presented and discussed.
The data set is unique in a way that it is the first longterm investigation of polynya dynamics based on cloudcover corrected thermal-infrared data that covers the complete southern Weddell Sea coastal area. The results were discussed in comparison to recently published studies using a variety of different methods and approaches (satellite sensors and models). On average over 13 years, we find the sea-ice areas in front of the Ronne and Brunt ice shelves to be the most active with an annual average polynya area of 3018 ± 1298 and 3516 ± 1420 km 2 as well as an accumulated volume ice production of 31 ± 13 and 31 ± 12 km 3 , respectively. For the remaining four regions, our estimates amount to 421 ± 294 km 2 and 4 ± 3 km 3 (Antarctic Peninsula), 1148 ± 432 km 2 and 12 ± 5 km 3 (iceberg A23A), 901 ± 703 km 2 and 10 ± 8 km 3 (Filchner Ice Shelf), as well as 499 ± 277 km 2 and 5 ± 2 km 3 (Coats Land).
Given these results and also the presented thin-ice frequency distribution, the neglect of certain regions in other studies, namely the area around the grounded iceberg A23A as well as the area off the coast of the Filchner Ice Shelf and Coats Land, showed up to drastically underestimate the total average polynya area and ice production in the southern Weddell Sea. Together, all three regions contribute comparably to the most active regions in front of Ronne and Brunt ice shelves. These regions should be further investigated by upcoming studies as they also contribute to the bottom-water formation. Supplementary data are available at doi:10.1594/PANGAEA.848612 (Paul et al., 2015b).