Multidecadal anomalies of Bohai Sea ice cover and potential climate driving factors during 1988–2015

Despite the backdrop of continuous global warming, sea ice extent has been found not to consistently decrease across the globe, and instead exhibit heterogeneous variability at middle to high latitudes. However, the existing studies are focused primarily on high latitude frozen seas, while studies on the long-term variability of sea ice cover at middle latitudes are generally lacking. Afforded by continuous satellite imagery, evolution of sea ice cover over nearly three decades from 1988 to 2015 in the Bohai Sea as a peculiar mid-latitude frozen sea area is reported for the first time. An anomalous trend of slight overall increase of 1.38 ± 1.00% yr–1 (R = 1.38, i.e. at a statistical significance of 80%) in Bohai Sea ice extent was observed over the 28 year period. The detrended annual average ice area (AAIA) was further found to correlate with a slight decreasing mean ice-period average temperature (IAT, r = –0.58, p < 0.01) of 11 meteorological stations around the Bohai Sea as well as a mild increasing cumulative freezing degree days (CFDD, r = 0.65, p < 0.01). Correlation with decreasing Arctic Oscillation (AO) index (r = –0.60, p < 0.01) and North Atlantic Oscillation (NAO) index (r = –0.69, p < 0.01) over the study period suggested AO and NAO as the primary large-scale climate factors for Bohai Sea ice. In addition, the seasonal cycle of ice cover showed a single peak with longer freezing phase than melting phase, due to the different temperature change rate during the freezing and melting phases. The results can provide important references for monitoring the recent climate change in the region and beyond.


Introduction
Sea ice cover has long been recognized as a sensitive and important indicator of the climate system in both global and regional observation and modelling studies (Vinnikov et al 1999, Bai et al 2011, Notz and Stroeve 2016. Due to the continuous global warming in recent decades (Vihma 2014), there has been a decline in ice extent in most polar and sub-polar frozen seas, including the Baltic Sea (Omstedt and Chen 2001), the sea of Okhotsk (Harada et al 2014), the Hudson Bay and the Canadian Archipelago . The Bohai Sea is a semi-enclosed sea located in North China (figure 1), which is the southernmost frozen sea in the Northern Hemisphere. A better understanding of the variability of ice cover in this region is important for monitoring regional climate change (Gong et al 2007, Bai et al 2011. In addition, sea ice hazards in the Bohai Sea can severely affect aquaculture, marine transportation, offshore oil field operation and other economic activities in the Bohai Rim as a national economic hub of China , Zhang et al 2016. And yet when managed properly, sea ice can also be desalinated to produce freshwater for agriculture or industry, and thus has the potential to help alleviate the severe fresh water shortage in the region as one of the most water-scarce areas in China . Moreover, knowledge of the ice cover is necessary to decipher the sea ice contaminant pathways (Beattie et al 2014, Tovar-Sánchez et al 2010. Therefore, research on the longterm variability of sea ice cover can provide important references for climate change and disaster monitoring and risk management, assessing potential pollutants release, and estimating sea ice resources. Despite the backdrop of continuous global warming, sea ice extent has been found not to consistently decrease across the globe, and instead exhibit heterogeneous variability at middle to high latitudes. While a near-zero trend for Arctic sea ice extent from 2007 to 2013 was observed (Swart et al 2015), increasing trends were reported in other regions including the Gulf of St. Lawrence , the Bering Sea , and the Antarctic (Gagne et al 2015). However, the existing studies including those mentioned above are focused on high latitude frozen seas while studies on the long-term (e.g. multidecadal) variability of sea ice cover at middle latitudes including the Bohai Sea are generally lacking. For those addressing the annual or multiannual variability of sea ice cover in the Bohai Sea, the variability has been attributed to climate factors, such as local climate forcing (Su et al 2012) and Arctic Oscillation (AO) (Bai et al 2011) in the literature. Specifically, Su et al (2012) found that local climate forcing had a very direct and immediate impact on ice cover, and the rapid formation of sea ice in the Bohai Sea during the 2009/10 winter was caused by continuous cold snaps coincident with anomalously cold weather. Bai et al (2011) and Gong et al (2007) reported that AO played a prominent role in mid-high latitude climates via teleconnections, consequently affecting the ice conditions in the Bohai Sea.
In the present study, the evolution of sea ice cover in the Bohai Sea during the winters of 1988-2015 was analyzed using satellite remote-sensing imagery, to elucidate the multidecadal variability in this peculiar mid-latitude frozen sea area. 640 satellite images of sea ice cover during the study period were analysed. Statistical methods such as correlation analysis and linear regression were used to explore the potential climate driving factors during this period.
The paper is organized as follows: section 2 introduces the study area, remote sensing and relevant climate factor data sources as well as the methods used to extract and process the data. The intersensor calibration of extracted sea ice area, the primary results on the multidecadal and seasonal variability of the sea ice cover from 1988 to 2015, as well as the relationships between ice extent and climate factors are presented in section 3, together with the relevant discussion. The paper is concluded with a summary of the major findings. These data were subsequently used to calculate the cumulative freezing degree days (CFDD). The monthly AO index, North Atlantic Oscillation (NAO) index and the Niño 3.4 region sea surface temperature anomaly index, which was used as a proxy of the El Niño-Southern Oscillation (ENSO) variability, from 1988 to 2015 was obtained from the NOAA/ Climate Prediction Center (www.cpc.ncep.noaa.gov/). The North Pacific (NP) index were downloaded from National Center for Atmospheric Research (https:// climatedataguide.ucar.edu/climate-data/north-pacif ic-np-index-trenberth-and-hurrell-monthly-and-win ter). The monthly Pacific Decadal Oscillation (PDO) data were downloaded from the University of Washington/JISAO (http://research.jisao.washington. edu/pdo/). Correspondingly, annual circulation data were obtained by averaging the monthly data during the winter months (December-March).

Methodology 2.3.1. Extraction of sea ice area
The AVHRR images were pre-processed using PaNDA (Package for NOAA Data Analysis) software package (Shimoda et al 1998, which included format conversion, radiometric correction, atmospheric correction and geometric correction. In the format conversion, the hierarchical data format was converted to meta data format for subsequent correction. Then the radiometric correction was undertaken to eliminate the image distortion caused by radiation errors, which was followed by the atmospheric correction to remove the effects of the atmosphere on the reflectance values of the satellite images. The geometric correction was then performed to remove potential geometric distortions from a distorted image. The pre-processed images were used to extract the sea ice area. The extraction method is primarily based on the differences in reflectance and surface temperature among the sea ice, sea water and land. The first step was demarcation of sea ice and land. The near-infrared portion of the spectrum (wavelength greater than 0.8 mm), corresponding to band 2 of NOAA/AVHRR, was used to delineate land from sea areas. The sea ice and seawater were then differentiated by the zonal threshold method combining bands 2 (reflectance) and 4 (brightness temperature). The following criteria were used here for icewater seperation, Ice : Band 2 > T 2 and Band 4 < T 4 ð1Þ where the thresholds T 2 and T 4 were set in the range of 0.08∼0.14 and 266.0∼269.0 for band 2 and 4 for different zones (see figure 1 for the zone division), respectively. Detailed description of the method can be found in Xie et al (2003) and Liu et al (2013). The zonal threshold method was found to extract the sea ice area with an accuracy of ∼80% (Guo et al 2008).
The MODIS images were initially processed for geometric correction, format conversion and atmospheric correction in ENVI. The desired sea ice features were then extracted from these processed images using an object-based segmentation and feature extraction method, which simplified and compressed the image to objects that were more identifiable (Zhao et al 2012, Hussain et al 2013, Miao et al 2015. Specifically, the images were first segmented into multiple objects or groups of pixels. After image segmentation, supervised classification was performed to sort the objects. Four types of examples, namely, sea ice, land, cloud and seawater, were chosen for the example-based feature extraction workflow in ENVI. The classification considered not only the spectral values but also the spatial patterns including the texture, shape, and contextual properties of each individual object. Finally, the k-nearest neighbour algorithm (k = 5) was used to generate the shapefile, which was further corrected through Environ. Res. Lett. 12 (2017) 094014 visual interpretation in ArcGIS, i.e. removing any polygon that was incorrectly classified as sea ice and vice versa. Zhao et al (2012) reported that the accuracy of extracting sea ice area using the object-based method can reach ∼90%. Figure 2 shows an example of sea ice area extraction.
A total of 640 remote-sensing images, each of which corresponds to a separate date during the study period, were selected after removing those covered by clouds or with inferior quality from more than 2500 downloaded images. The annual and monthly image distributions over the entire study period are shown in figure 3. Notably, there are fewer images in the 1991/92 (n = 6) and 1999/2000 (n = 13) winters due to severe cloud cover. As far as monthly distribution is concerned, the images in January and February (period with maximum ice cover) account for 39.8% and 31.7% of the total images over the entire study period, respectively, whereas those in December (freezing period) and March (melting period) account for 21.9% and 6.6%, respectively. The 640 images were analyzed to determine daily ice area (DIA). Annual average ice area (AAIA) and weekly average ice area (WAIA) were subsequently derived from the aggregates of DIAs. Annual maximum ice area (AMIA) was also identified from DIAs for each winter.

Calculation of CFDD
The formation of sea ice in a specific region is closely related to how long the air temperature has been below the freezing point of the seawater, which is quantified by the CFDD (Tamura-Wicks et al 2015). The formula of the CFDD is as follows: where the CFDD is the cumulative freezing degree days; T f is the freezing point of the seawater, T a is the mean of the daily air temperature of 11 meteorological stations surrounding the Bohai Sea; t is time measured in days over which CFDD is calculated; and d s and d e are the start and end date of the period considered, respectively. The critical freezing temperature T f is set at -2.0°C, which is consistent with the salinity of 30 PSU in the Bohai Sea (Su et al 2012).

Trend and correlation analyses
The trends presented in this study were calculated by linear least squares regression. Slopes and standard deviations of fitted lines were calculated following Taylor (1997). A rough indicator of statistical significance as non-0 was determined by calculating the ratio (R) of the magnitude of the slope to its standard deviation (Parkinson 2014, Parkinson and. A trend was significant at a 90% or 95% confidence levels if R exceeded 1.70 or 2.06 using the Student's t-test with 26 degrees of freedom (the number of years minus two). To assess the correlations between sea ice area and potential climate driving factors, the Pearson correlation method was adopted using detrended data. The detrended time series of sea ice area and potential climate driving factors were calculated by subtracting the value determined by the linear least squares trend from the corresponding original value for each year (Serreze et al 2016).

Results and discussion
3.1. Intersensor calibration of extracted sea ice area Since Guo et al (2008) found that the extracted sea ice area in the Bohai Sea from MODIS data appeared to be more accurate than that from AVHRR data, intersensor calibration against MODIS results was performed to enhance the accuracy of AVHRR results. A dataset (n = 91) were selected from two sensors on the same dates from 2000 to 2010 and compared in figure 4(a). Despite largely consistent trends, significant deviation is noted between MODIS and AVHRR results, and the rootmean-square error (RMSE) of sea ice area extracted from these two sensors is 26.6%. As such, linear regression was used to match the extracted sea ice area as in Guo et al (2008), and the resultant agreement is satisfactory (R 2 = 0.97, p < 0.01). The obtained linear regression model shown in figure 4(b) was used to calibrate the AVHRR results, and the calibrated sea ice area is also shown in figure 4(a) (solid line and dots). The calibrated result is much closer to the MODIS result, and the RMSE is reduced to 11.6%.

Seasonal cycle of sea ice cover in the Bohai Sea
The seasonal cycle of ice extent in the Bohai Sea was derived from the WAIA averaged over the 1988-2015 period ( figure 6). Overall, the seasonal cycle of average WAIA shows a single peak. The typical seasonal cycle of sea ice area in the Bohai Sea kicked off on an ice freeze-up date in early December. Afterwards, the ice growing phase lasted approximately 8 to 9 weeks until late January. The ice cover in the Bohai Sea reached the seasonal maximum around late January to early February. From mid-February to mid-March, the WAIA dropped rapidly, and the ice break-up phase lasted approximately 5 weeks. Notably, the ice growing phase (8 weeks) was considerably longer than the breakup phase (5 weeks) in the Bohai Sea, suggesting that the processes of sea ice formation and melting were not symmetrical. The different freezing and melting rates were closely related to the local temperature change rate, which were -0.08°C day À1 and 0.19°C day À1 during the sea ice freezing and melting phases, respectively, after averaged over the 28 years study period. Notably, Wang et al (2012) also reported the asymmetry of the seasonal ice cycle in the Great Lakes.
The standard deviation associated with the average WAIA is also shown in figure 6, and the large interannual variability indicates that the local synoptic forcing such as surface air temperature and wind may play an essential role in ice conditions in the Bohai Sea. Wang et al (2012) found that the natural variability of ice cover in the Great Lakes was largely affected by the internal climate forcing. The predictability of sea ice cover by statistical or numerical models was generally poor due to the large STD, especially at the multidecadal time scales. Figure 7 shows that the sea ice coverage, i.e. the percentage of the total Bohai Sea area covered by ice, exhibits significant monthly changes over the study period. In December, the sea ice gradually formed and attained low coverage (84% frequency < 10%). Whereas in January and February, the frequency of the sea ice coverage in the range of 10%-20% increased to 48% and 43%, respectively, and high coverage (>25%) occurred exclusively in these two months as well. As the sea ice melted in March, the ice coverage decreased accordingly, resulting in its domination in the lowest coverage (<5%) category. From the frequency distribution of the entire ice season, the ice coverage was most concentrated in the low and medium proportions (0%-15%) over the study period. Environ. Res. Lett. 12 (2017) 094014

Correlations between sea ice cover and climate factors
The Bohai Sea is a shallow sea with an average depth of 18 m, and the local synoptic forcing presumably plays an essential role in the ice conditions as mentioned above. Since our primary focus is the dynamics of the sea ice area of the entire Bohai Sea, and further in viewing the influence of the various large-scale climate factors such as AO, ENSO, PDO and NPon the local synoptic forcing in the Bohai Sea, only climate factors and derivatives at regional scale and beyond were considered in this study.
To determine the role of climate forcing on the evolution of the sea ice cover in the Bohai Sea, correlative relationships between the sea ice cover and the potential climate factors were examined.

The multidecadal trends of variation of sea ice cover and climate factors
The AAIA, the mean of the ice-period average temperature (IAT, defined as the average temperature for the ice period from December 1 to March 15 in the following year) of 11 meteorological stations around the Bohai Sea as well as the associated CFDD, and the AO index for the 1988-2015 period are shown in figure 8. The straight trend lines presented in the figure were obtained using the least squares fitting. The linear regression returns a slight overall increase of 1.38% ± 1.00% yr -1 (R = 1.38, i.e. at a statistical significance of 80%), which is in parallel with the slight decrease in IAT (R = 0.74), the mild increase in CFDD (R = 1.04) and decreasing AO index (R = 1.74, i.e. at a statistical significance of 90%) (see figure 8).
Despite the backdrop of continuous global warming, sea ice area has been found not to consistently decrease across the globe. Ji et al (2014) reported that global warming was spatiotemporally heterogeneous, i.e. the change in local temperature could exhibit different trend (rise or fall). The sea ice extent in the Antarctic had increased slightly from 1979 to 2012, in contrast to the dramatic decline in Arctic sea ice extent over the same period (Gagne et al 2015). The observed increase in sea ice area in the Antarctic was attributed to the natural climate variability. In particular, the increased cyclonic atmospheric flow over the Amundsen Sea in the Antarctic was identified as a key factor, which brought cold air from the south to the north . In the Northern Hemisphere, an increase, although statistically not significant with R = 0.43, in the average ice area in the Bering Sea was also reported. The anomalous trend was attributed to the weakening of the Aleutian Low, resulting in less warm air from the south being adverted into the Bering Sea region, thus favouring northerly winds that promoted ice formation Parkinson 2012, Stroeve et al 2012).

Correlation with local temperature
A moderate negative correlation between the detrended time series of AAIA and the mean value of the IAT (r = -0.58, p < 0.01) of 11 meteorological stations around the Bohai Sea during the study period was shown in figure 9(a), suggesting that the average sea ice condition in each year was closely related to the IAT. A moderate positive correlation (r = 0.65, p < 0.01) between the detrended time series of AAIA and CFDD is shown in figure 9(c), which suggests that CFDD is another important climate factor for sea ice formation. Similar correlations were also found between AMIA and IAT/CFDD (see figures 9(b) and 9(d)). This suggests that the local temperature is a controlling factor for the sea ice conditions in the Bohai Sea, which is consistent with the findings of Zhang et al (2016).

Correlation with AO, NAO, ENSO, PDO and NP
The AO is the principal empirical orthogonal function (EOF) mode of wintertime sea level pressure anomalies of the extra-tropic Northern Hemisphere (Ogi et al 2016). The NAO refers to a redistribution of atmospheric mass between the Arctic and the subtropical Atlantic, which is measured by the Environ. Res. Lett. 12 (2017) 094014 difference between the sea level pressure over Stykkisholmur, Iceland and Lisbon, Portugal (Hurrell 1995). As shown in figure 10, the Bohai Sea ice extent tended to be smaller in high-AO and NAO phases and larger in low-AO and NAO phases. The detrended Bohai Sea ice area exhibits moderate negative correlation with the AO index (r = -0.60, p < 0.01) and NAO index (r = -0.69, p < 0.01). During the negative phase of the AO, enhanced atmospheric pressure over the Arctic induced weaker westerly winds in the upper atmosphere, which allowed the cold air to reach further southern regions, resulting in colder winter in these regions (Bai et al 2011). At the same time, when the NAO index was in negative phase, the pressure difference was small, which decreased the storminess in the westerly wind belt and increased the zonal circulation and pushed large amounts of cold Arctic air to the south, resulting in significant cooling in Northeast Asia including the Bohai Sea (Bai et al 2011. Presumably due to the above mechanisms, the Bohai Sea usually experienced higher-than-normal ice cover during low AO and NAO winter. All harsh winters (1997/98, 2000/ 01, 2009/10, 2010/11 and 2012/13) occurred during the negative phase of AO and NAO. In the 2009/10 winter in particular, very strong negative AO and NAO occurred, resulting in the maximum ice extent over the entire study period.
To estimate the impact of modes other than the AO and NAO on the evolution of sea ice area during 1988-2015, we assessed the correlation of sea ice extent to relevant large-scale climate factors such as ENSO, PDO and NP again using detrended time series (The relevant plots are not shown here for brevity). The detrended Niño 3.4 index exhibits a weak and insignificant correlation (r = 0.10, p = 0.62) with Bohai

Conclusions
The evolution of ice cover in the Bohai Sea from 1988 to 2015 was analyzed using satellite remote-sensing imagery. The correlation between the ice cover and the potential climate driving factors that influenced the ice conditions were also explored. The major findings from the present study are as follows: 1. Despite the backdrop of continuous global warming, the ice extent in the Bohai Sea exhibited a small overall increasing trend of 1.38% ± 1.00% yr -1 during the study period (R = 1.38, i.e. at a statistical significance of 80%).
2. The seasonal cycle of the ice extent in the Bohai Sea shows a single peak with longer freezing phase than melting phase. The different freezing and melting rates were closely related to the local temperature change rate, which were -0.08°C day À1 and 0.19°C day À1 during the freezing and melting phase, respectively.  Environ. Res. Lett. 12 (2017) 094014 3. The detrended AAIA was further found to correlate with a slight decreasing mean ice-period average temperature (IAT, r = -0.58, p < 0.01) of 11 meteorological stations around the Bohai Sea as well as a mild increasing cumulative freezing degree days (CFDD, r = 0.65, p < 0.01). Correlation with decreasing Arctic Oscillation (AO) index (r = -0.60, p < 0.01) and North Atlantic Oscillation (NAO) index (r = -0.69, p < 0.01) over the study period suggested AO and NAO as the primary large-scale climate factors for Bohai Sea ice.
Due to the inherent constraint of the remotesensing analysis, we only considered sea ice extent in the present study, and follow-up study incorporating sea ice thickness would be of interest in the future.