Sea Ice Thickness Retrieval Based on GOCI Remote Sensing Data: A Case Study

: The accurate monitoring and measurement of sea ice thickness (SIT) is crucial for understanding climate change and preventing economic losses caused by sea ice disasters near coastal regions. In this study, a new method is developed to retrieve the SIT in Liaodong Bay (LDB) based on the Rayleigh-corrected reﬂectance from Geostationary Ocean Color Imager (GOCI) images in the winters of 2012 and 2013. Compared with previously developed SIT retrieval methods (e.g., the method based on the thermodynamic principle of sea ice) using remote sensing data, our method has signiﬁcant advantages with respect to the inversion accuracy (achieving retrieval skill scores as high as 0.86) and spatiotemporal resolution. Moreover, there is no signiﬁcant increase in the computational cost with this method, which makes the method suitable for operational SIT retrieval in the global ocean.


Introduction
Sea ice plays crucial roles in the earth's climate system and ecosystems [1][2][3][4][5], although it covers only approximately 7% of the global ocean [6]. It has a paramount position in thermodynamics by modulating sea ice ablation and accretion processes [7]. In addition, it can affect the vertical distribution of marine hydrological elements and the movement of sea water. However, once a sea ice disaster occurs, a series of problems for offshore operations, ship navigation, port transportation, and marine fisheries follow [5,7,8]. The Shengli and Liaohe oil fields distributed in Liaodong Bay (LDB) and the oil and gas drilling platforms in this region often suffer much damage due to sea ice disasters [9]. For example, the heavy sea ice disaster in LDB in 2010 led to an economic loss of approximately 292 million dollars [8,10]. Since the severe ice cover event in the Bohai Sea in 1969, more attention has been given to monitoring and investigating sea ice to provide scientific support for the prevention and mitigation of sea ice disasters. Therefore, research on the sea ice of the Bohai Sea is required.
Sea ice thickness (SIT) is an important index for describing sea ice conditions. It represents an integrated measure of total sea ice production, the surface salinity and freshwater fluxes to the ocean, and the total ocean-atmosphere heat loss [11]. Generally, the thicker the ice the higher its resistance and, thus, the lower its maximum speed [12]. Therefore, SIT data are also the basis for studying sea ice drift. Although more accurate SIT measurements can be obtained by relying on empirical macroscopic observations and in situ measurements from on-ship surveys than through other approaches, timely SIT measurements are difficult due to the highly variable distribution and morphology of sea ice. Numerical simulation is an effective way to study SIT because it can be used to calculate and simulate SIT time series data. In addition, coupled sea ice-ocean models can provide insight into the heat flux exchanges and dynamic processes of sea ice growth and melting cycles [13]. However, model results may deviate from observed data because many uncertainties exist in model implementation, such as the initial conditions, boundary conditions, forcing, and model parameterization. Fortunately, satellite remote sensing presents a solution to these obstacles to obtaining sea ice observations.
In the last four decades, first-and multiyear ice has been monitored with remote sensing techniques [14][15][16][17][18][19]. The thickness of multiyear ice and thick first-year ice has been estimated by radar and laser altimeters on board CryoSat-2 and ICESat [20,21]. However, altimeter-based ice measurements have large uncertainties for ice thicknesses less than 1 m [20]. The L-band radiometer on board ESA's Soil Moisture and Ocean Salinity (SMOS) mission has proven effective for detecting ice thicknesses less than 1 m during the winter period [22,23]. Passive microwave radiometer data from the Special Sensor Microwave Imager (SSM/I) and Advanced Microwave Scanning Radiometer-Earth Observing System (AMSR-E and AMSR-2) have been used to estimate very thin ice, i.e., ice with a thickness of 10-20 cm [24,25]. However, these active and passive remote sensing data cannot satisfy the need to map the hourly SIT at fine scale in local regions such as the Bohai Sea. For example, although synthetic aperture radar (SAR) has the capability to penetrate clouds and fog to realize all-weather observations throughout the day, it is unable to monitor thick sea ice due to the shortcomings of radar technology and the accuracy of ice thickness inversion needs to be verified [26,27]. Pan [28] analyzed the conditions in the formation, development and melting of sea ice in the Bohai Sea based on Moderate-Resolution Imaging Spectroradiometer (MODIS) data. Using MODIS data and meteorological data, Zhang [13] found that the sea ice area and its spatial distribution in LDB are very sensitive to cumulative frozen days, wind, and temperature. Ouyang [29] calculated the SIT of the Bohai Sea using a thermodynamic heat balance model based on MODIS data and atmospheric reanalysis data. However, it is difficult to monitor the thickness of flowing ice using remote sensing data because the replay cycle is very long. To address the problem of coarse spatiotemporal resolution, Liu [9] proposed a method to estimate SIT using recent Geostationary Ocean Color Imager (GOCI) data and found that GOCI data have the potential to yield extensive, accurate, and near-real-time (with a time resolution of 1 h) monitoring of sea ice changes. Recently, Yan [30] monitored the characteristics of Bohai Sea ice using GOCI data, and Landsat Thematic Mapper (TM) data with higher spatial resolution than the GOCI data were used to validate the data. However, few SIT retrieval models based on Rayleigh scattering of GOCI images exist. Rayleigh scattering has a large contribution to the total radiation of the remote sensing receiver [31], and has the greatest influence on the blue-green band, especially in the blue band up to 90% [32]. Ice has a strong reflection effect on the blue-green band (400-500 nm), and this part of the spectrum is more sensitive to different SITs [33], so we use the reflectance corrected by Rayleigh scattering to invert the SIT.
In this paper, a novel multiband model is provided that retrieves SIT based on GOCI images. We focus on the improvement of inversion accuracy and the optimization of the retrieval process. The paper is organized as follows: in Section 2, the study area is introduced. The information and processing of GOCI data are presented in Section 3, and the retrieval methods are elaborated in Section 4. Section 5 describes the results regarding the temporal and spatial variation in SIT. Section 6 discusses the advantages of our method and the effect of snow on sea ice in the retrieval algorithm. Finally, the study is summarized in Section 7.

Study Area
The Bohai Sea ranges from 117 • 35'E to 122 • 15'E longitude and 37 • 07'N to 41 • 00'N latitude and includes LDB, Bohai Bay, Laizhou Bay, the Central Bohai Sea, and the Bohai Strait. LDB is in the northern temperate zone, the long-term average temperature in this region is 10.7 • C, and the average annual rainfall is approximately 500 mm. The average annual wind speed is 6.2 m/s and the predominant wind direction is NNW [34].
The Bohai Sea, an important economic development zone in China due to its rich resources (such as fishing and tourism resources, a harbor, oil, and sea salt) and valued environmental conditions, is the most southern area in the Northern Hemisphere with sea ice in winter [35]. As one of the bays encompassed by a semi-enclosed sea, LDB is the northernmost water body with ice and the source of sea ice in China [36]. It has an area of 30,597.35 km 2 and its maximum water depth is 32 m [13]. Sea ice in the bay generally forms in early winter and advances along the western coast to the south with cold air from Siberia, and it retreats at the end of February [13,35,37]. According to the growth process of sea ice, the sea ice of the Bohai Sea is first-year ice; it can be further classified as fast ice or drift ice according to its movement state. For a long time, fast ice was consistently observed near stations in LDB; it is thick, and the fast ice close to the northern coastline of LDB has been separated from other ices. The sea ice extent reaches its maximum in early February [38], and the annual ice age is approximately 3 to 4 months. In LDB, sea ice typically appears as early as December and lasts approximately three months, to the following February. Its extent is approximately 70 km from the shoreline. The typical ice thickness in this region is approximately 10-25 cm for the drifting ice and approximately 20-40 cm for the fast ice [39].

GOCI Images
The GOCI is a geostationary orbit satellite image sensor used to observe and monitor ocean colors around the Korean Peninsula [40]. It was loaded on the Communication, Ocean, and Meteorological Satellite (COMS) of South Korea, which was launched in June 2010, and its band information is shown in Table 1. Its spatial resolution is approximately 500 m, and the range of the target area is approximately 2500 × 2500 km, centered on the Korean Peninsula (with the center of the observation area located at 36 • N and 130 • E). The coverage area includes the Bohai Sea, the Yellow Sea, and part of the East China Sea [9]. The sensor has a high temporal resolution of 1 scene per hour and a spatial resolution of 500 m [41], which makes wide spatial coverage and near-real-time SIT monitoring in LDB possible.

SIT Observation
The in situ sea ice observations were carried out at two offshore oil platforms (JZ9-3 and JZ20-2) in LDB ( Figure 1) and included ice type, ice concentration, and ice drift speed and direction [42]. The measurements were carried out hourly during daytime (07:00-18:00 local time). The ice thickness range around each oil platform was also observed and is used in this paper for comparison with the GOCI-retrieved ice thickness data ( Table 2).

Data Processing to Eliminate Atmospheric Effects
Many studies have focused on cloud cover in remote sensing data. Although various methods have been proposed to address the lack of data caused by cloud coverage while removing the influence of clouds, the loss of spectral information of ground objects inevitably occurs. Most of these methods employ algorithms that remove clouds by homomorphic filtering, the replacement method, or time-sequence filtering [43]. However, changes in sea ice are extremely spatiotemporally sensitive and different types of ice may be distributed over a few minutes or a few square meters. The existence of ice water mixtures and crevasses makes interpolation algorithms inappropriate. Other studies using GOCI images have chosen clear and cloudless images as data sources and have not focused on cloud cover processing [9,44]. In this study, to eliminate atmospheric effects, we similarly selected clear and cloudless GOCI images of LDB. These images, from 2012 to 2013, were processed for geometric rectification and atmospheric correction, and images of the region of interest were extracted by the GOCI Data Processing System (GDPS) platform. After atmospheric correction the effect of Rayleigh scattering was eliminated and Rayleigh-corrected remote sensing reflectance (R rc ) was obtained [45]: where L * t is the radiance of the sensor after radiation calibration, F 0 is the vertically incident solar irradiance outside the atmosphere, θ 0 is the solar zenith angle, and R r is the Rayleigh scattering reflectance estimated by the 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) model [45].
Rayleigh scattering (molecular scattering) occurs when the particle size is much smaller than (less than one tenth) the wavelength of the incident light. The intensity of scattered light is inversely proportional to the fourth power wavelength of the incident light [46]. The Rayleigh scattering spectrum is directly related to the velocity distribution of molecules and atoms and thus contains information about the thermodynamic parameters of the flow, including density, temperature, and bulk velocity [47]. The sea ice surface albedo is much greater than the albedo of sea water and other coastal features. It was possible to ignore the influence of clouds because cloud-free GOCI images of the study area were chosen. Thus, the Rayleigh-corrected reflectance could be used for retrieval to distinguish SITs.
To determine the Rayleigh-corrected reflectance of sea ice of different thicknesses, we chose 8 positions in LDB and calculated the Rayleigh-corrected reflectance in 8 bands at these positions ( Figure 2). We found that band 1 can well distinguish sea ice and open water.

Sea Ice Extent Extraction
To extract and classify sea ice, visual interpretation, thresholding [48], and analysis of simple spectral characteristics of satellite images [10,49] were used. Accurate determination of sea ice extent is very significant for understanding the considerable decrease in sea ice [50]. The sea ice extent and its variability were defined following the method of Shapiro [50]. The extent was extracted by visual interpretation, assisted by the thresholding method based on the features of the Rayleigh-corrected reflectance spectral curve (Figure 2). For the convenience of description, the Rayleigh-corrected reflectance of the x-th band was denoted by R rcx . The threshold is defined repeatedly in the spectrum extraction experiment and is obtained by numerical approximation. The threshold based on R rc1 > 0.055 was used to extract a first estimate of sea ice extent to separate it from open water. To avoid the misjudgment of sea ice caused by offshore suspended sediment and broken edges of sea ice in band thresholding, the extraction of sea ice extent was corrected by visual interpretation. Some of the extracted results are shown in Figure 3.

Optimization of the Ice Thickness Inversion Algorithm
The SIT is strongly related to the sea ice surface albedo [51][52][53]; therefore, empirical retrieval based on remote sensing reflectance has been attempted in numerous stud-ies [9,10,49,52]. Liu [9] calculated SIT in the Bohai Sea using the relationship of shortwave broadband albedo and SIT proposed by Grenfell and Perovich [54]. GOCI-derived shortwave albedo was obtained through linear regression of the shortwave broadband albedo estimated by MODIS data and the GOCI reflectance of each band [9]. This method was established under only the best light conditions each day (10:00-12:00). Therefore, it cannot be used under other light conditions. We attempted to avoid this shortcoming in our method by using in situ SIT measurement data and multiband reflectance for the regression analysis.

Retrieval of SIT Using Several Bands
We conducted multiple linear regression on random combinations of three bands to select the most representative of the 8 bands; there were 336 combinations in all. The root mean square error (RMSE) and skill score were calculated for each combination, and the results are shown in Figure 4, We found that the permutation "band 2 (443 nm, blue), band 4 (555 nm, green), and band 8 (865 nm, near-infrared)" yielded the minimum RMSE and the highest skill score. (The formulas of RMSE and skill score are shown in Equations (3) and (4), respectively.) Therefore, this combination was considered the most suitable band combination for ice thickness inversion. The retrieved results, using Equation (2) and derived from these three bands, were compared with the in situ data in Figure 5, and the correlation coefficient was 0.7.

Retrieval of SIT Using All Bands
The setting of the GOCI bands does not specifically target sea ice detection. Therefore, considering only individual bands could result in the loss of some important information. As a next step, we considered all the components in the ice thickness inversion algorithm. The ice thickness was fitted and simulated by multiple regression methods and multiple linear regression, respectively. The multiple linear regression equation of all the bands is shown in Equation (5).
We also used principal component analysis (PCA) to obtain the SIT. PCA was originally developed as a tool to reduce a large data set to a small set and extract the main information from a given data sample. This method has been used to identify distinct spatial and spectral patterns in the case of multiband remote sensing data [55]. In this study, we performed PCA over all bands of GOCI, and the results are listed in Table 3. Principal component 1 (PC1) and principal component 2 (PC2) were chosen to build the SIT retrieval method. PC1 accounted for 78.23% of the total variance, and PC2 accounted for 7.26% of the total variance. We used these two principal components to retrieve SIT by multiple linear regression in Equation (6). This paper focuses on the development of the SIT inversion method, and ice thickness changes and predictions will be developed in our future work. We compared the RMSE values and skill scores of our retrieval results with those obtained using Ouyang's method [29]. We found that multiple linear regression of 8 bands was best suited for retrieving the SIT in LDB ( Figure 5).  Table 2.

Temporal and Spatial Variation of the SIT in LDB
The SIT during 2012 and 2013 was retrieved by multiple linear regression (Equation (5)). Its distribution is mapped in Figure 6. The SIT is largest in the northeast part of LDB, especially in the region close to the northern shoreline, and it is consistent with the photos from Figure 6v,w. The sea ice becomes increasingly thin from north to south, and the same trend appears on the southwestern "tail". In the east-west dimension, the SIT in the bay near the west and east coasts is small relative to that elsewhere, with the largest SIT usually occurring in the middle of the bay. This phenomenon can be described as follows: thicker ice is always situated between the two "corners" of the region of Liaoning Province.
Our method based on GOCI images can capture SIT changes from the interannual scale to the hourly scale. Comparing the SIT between 2012 and 2013, we can see that the SIT in February 2013 is larger than that in February 2012, as shown in Figure 6a-l,r-u. This finding is consistent with Yan's results [30]. Focusing on January and February 2013, SIT tends to decrease over time, and the findings confirm the statistical results: the sea ice in LDB began to melt gradually at the end of February and early March. Moreover, the change in SIT is obvious over a week (Figure 6m-o) and even over a day (Figure 6c-f). The inversion method established can clearly capture the diurnal variation in sea ice, which was most apparent in the thin ice area. This phenomenon can be explained as follows: thinner ice melts more easily during the daytime than at night due to the solar radiation during the day, especially near noon; furthermore, thinner ice drifts more easily than thicker ice. However, the changes over months or seasons are complex and depend on the wind, currents, temperature, precipitation, and the ice's internal structure. The SIT is consistently larger in January than in February, but its distribution varies between the different years.

Advantages and Limitations of the Retrieval Methods
As shown in Figure 5 and Table 4, our method has advantages when retrieving the SIT in LDB. The "lotus leaf ice" close to the north shoreline in Figure 6 is more consistent with the interpretation masks than is that in Figure 7. Another advantage is the high resolution of the images. Figure 7 presents the results of Ouyang's algorithm. Ouyang's method was developed based on MODIS data with 1-km resolution, and clear jagged edges with similar mosaic shapes are apparent in Figure 7c; it is almost impossible for such regular block-like sea ice to exist far from the shore in nature. Therefore, the results retrieved by the GOCI data appear more realistic. Table 4. The RMSE and skill score values of different SIT retrieval methods. Ouyang denotes Ouyang's method [29]; 8 bands denotes the eight-band regression method, and 3 bands denotes the 3-band regression method.  Previous studies on ice thickness have often ignored the precise extraction of the ice extent [56], but in our calculation process a clear ice extent was extracted first, which ensured that the sea ice area with thinner ice could be included. Although the accuracy of the inversion needs to be improved, the mask processing of the ice-free area contributes greatly to the improvement in calculation efficiency. In this sense, our method is suitable for producing operational SIT data.

Station
Another important factor affecting the thickness variability is ice drift. Ice is affected by precipitation and air temperature during its growth and development. The wind acts as an external force that causes ice drift, which blows away the surrounding ice. The dispersed ice can receive more solar radiation, which facilitates melting and ablation. In our next work, the diurnal variation in the SIT and extent of LDB will be retrieved to study the relationships between ice thickness variation and wind, precipitation, temperature, etc.
In general, three bands can be used to synthesize false-color images for the visual interpretation of remote sensing data, which is more suitable for supervised classification with artificial intervention. Equation (5) is similar to unsupervised classification, and the method considers that all bands contribute to sea ice identification. Bands 5, 7 and 8 are related to atmospheric correction, bands 2 and 3 are related to chlorophyll, and bands 1, 4 and 6 are related to yellow substances, suspended sediment, and fluorescence signals, respectively (Table 1). Suspended sediment and thin ice may sometimes be misjudged due to the phenomenon of "foreign body with spectrum" [56], while suspended sediment has a certain correlation with the spectrum of fluorescent substances [57]. In addition, suspended sediment will affect the water transparency, thus affecting the distribution of phytoplankton and chlorophyll [58]. Therefore, this relationship can also be used to assist in the retrieval of suspended sediment. Once the problem of the retrieval of suspended sediment is solved, the accuracy of thin ice identification and thickness inversion will be improved. Therefore, the information of all bands is needed for ice thickness inversion. In Yan's paper [30], the information of all the bands is also used to carry out ice thickness inversion. Even with combinations of 3 bands, the retrieval results have slightly lower skill scores for JZ20-2 than for JZ9-3 because the in situ SIT falls within a range, and the variation range for JZ20-2 is larger than that for JZ9-3. We cannot obtain accurate and reliable in situ observation data [59]. In addition, the uncertainty of the SIT inversion is related to the oil platform location. The SIT is more susceptible to the influence of tides, waves, and ice drift [60].
The atmospheric correction includes Rayleigh correction and aerosol correction. The former is dominant and accurate because Rayleigh scattering makes a large contribution to the total radiation of the remote sensing receiver. The latter is difficult to have accurate correction for, because aerosol scattering changes with time and space. The aerosol correction is accurate only in clear water areas, but it is difficult to conduct in turbid water and sea ice [61]. Atmospheric correction is necessary for ocean remote sensing because the ocean has a weak signal. However, it is not necessary for our study, because the ice signal is stronger than that of the atmosphere. After Rayleigh correction removes the major atmospheric contribution, the residual effect of aerosols is very small. The watercolor retrieval based on Rayleigh-corrected reflectance is already widely applied for turbid water with high concentrations of sediment [62] or algae [63]. Rayleigh correction is more feasible on the sea ice since the sea-ice signal is stronger than that of the turbid water. The results show that the accuracy of the inversion results after Rayleigh scattering correction is higher than that of previous results. However, it is better to have an aerosol correction on the reflectance of GOCI observations and we will apply this suggestion in the follow-up research work.

The Effect of Snow on Sea Ice in the Retrieval Algorithm
The effect of snow on sea ice should be considered. The snow surface has high reflectivity to solar radiation. Generally, the surface reflectance of fresh or clean snow can reach 86-95%, whereas that of gaps or unclean snow reaches only approximately 45%. The surface reflectance of sea ice is approximately 40-65% [64,65]. The weather in the Bohai Sea is dry and rainless in winter. Incidental snowfall will increase the surface reflectance of sea ice. However, how to remove the effect of snow on sea ice from remote sensing images remains a difficult problem. In addition, ice and snow begin to melt away at the end of winter and early spring. The albedo decreases with the increases in grain size and water content of granular snow for the thick snow layer [66][67][68], leading to underestimates in our inversion results in the thin ice area with snow melting.
The remote sensing image data we selected were mostly clear and cloudless. Therefore, we believe that the weather conditions were good and that there was no precipitation. Of course, a more rigorous approach is needed to take into account the cumulative precipitation time and the duration that snow can persist. The presence of snow will undoubtedly increase the spectral signal received by the sensor. Moreover, it can be seen from the comparison between our results and the measured data that the ice thickness values of the remote sensing inversion are consistent with the measured data. Therefore, we believe that this algorithm is suitable for ice thickness inversion in LDB with clear skies.

Conclusions
In this study, a new method was developed to retrieve the SIT in LDB based on the Rayleigh-corrected reflectance from GOCI images in the winter of 2012 and 2013. The main results are as follows: (1) The Rayleigh-corrected reflectance can be used for SIT inversion based on cloudfree GOCI images. The reflectance is high in thick sea ice areas and low in thin ice areas.
(2) We compared 4 GOCI inversion methods for SIT retrieval and found that multiple linear regression of 8 bands is best suited for retrieving the SIT in LDB. The skill score reached 0.86, which is higher than that of previous studies.
( 3) The inversion method based on GOCI images can capture SIT changes from the interannual scale to the hourly scale. Diurnal variation in the sea ice was apparent in the thin ice. This phenomenon can be explained as follows: thinner ice melts more easily during the daytime than at night due to the solar radiation, especially near noon; furthermore, thinner ice drifts more easily than thicker ice.
GOCI data have the potential to achieve extensive, accurate, and near-real-time monitoring of sea ice and thus can provide a suitable index to guide maritime operations and productions. This method makes sea ice evaluation in the Bohai Sea possible, guaranteeing the infrastructure safety of the region and promoting the exploration of resources and development of this region.
Author Contributions: Each of the authors were actively involved in the product generation or preparation of the paper. B.H. and Q.Y. conceived the idea for the study. F.G. and R.Z. performed the retrieval methods and prepared the paper. X.T.-K., B.H., L.Z., T.C. and Q.Y. processed the review and editing of paper. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: GOCI products was obtained from the Korea Ocean Color Satellite Center (KOSC) and are available from http://kosc.kiost.ac.kr/ (accessed on 22 January 2021). The measured data presented in this study are available on request from the corresponding author. The data are not publicly available due to 3rd party data policy.