Sea ice in the polar regions plays an important role in global climate by interacting with the ocean and atmosphere [1]. Snow-covered sea ice reflects a significant amount of incoming solar radiation back to space [2-5], thus working as an insulator between the ocean and atmosphere [6]. The annual variations of sea ice cover can influence ocean and atmospheric circulation [7,8], surface temperature [9], and ecology [10,11]. According to previous studies, sea ice extent and thickness in the Arctic has decreased since the 1960-1970s, mainly due to the dramatic increase in atmospheric CO2 emissions and its associated warming, which has a stronger signature in high northern latitudes [4,[12][13][14][15][16]. In contrast with the Arctic, several Antarctic studies have reported that sea ice extent increased for last 40 years, but started to decrease dramatically after a record high in 2014 [17,18]. Sea ice thickness and volume have also slightly increased, but a definitive trend is not available [19,20].
In particular, the Ross Sea is one of the most noticeable regions showing a significant increase of sea ice extent in the Southern Ocean [17,21,22]. There, massive amounts of sea ice are produced in the Ross Sea polynya which is the largest recurring polynya area in Antarctica [23][24][25]. Thus, sea ice dynamics in the Ross Sea has been the subject of a number of scientific studies [1,[26][27][28][29][30]. Due to the harsh weather and long polar nights, spaceborne remote sensing has been the most common and efficient approach for these studies; spaceborne remote sensing allows obtaining sea ice observations over a large area with regular time spans without the need of physical visits to the study region.
In recent years, the National Aeronautics and Space Administration (NASA)'s Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) has been providing precise sea ice freeboard and thickness estimates over the polar regions. ICESat-2 was launched in September 2018 with the objectives of measuring changes in ice sheets, glaciers, and sea ice thickness [31]. The Advanced Topographic Laser Altimeter System (ATLAS) on ICESat-2 measures the surface heights by using 532 nm wavelength laser photons. ATLAS has six beams of three pairs of strong and weak beams with 11 m footprints spaced by 0.7 m [32], providing a better spatial resolution compared to the previous ICESat mission [33]. Various studies have utilized ICESat-2 for the characterization of sea ice freeboard and thickness in the polar regions, taking advantage of this outstanding spatial resolution [34][35][36].
However, despite the good precision and resolution of ICESat-2, the use of ICESat-2 for studying sea ice is substantially limited by the spatial sparsity of the instrument's measurements (Figure 1). Since ICESat-2 data is collected along tracks, there are large sections of area for which there is no data because they are not under a track. In addition, ICESat-2 cannot measure freeboard in cloudy conditions because the laser photons do not penetrate clouds. Thus, ICESat-2 data have limitations in sampling sea ice elevation over a large area during a short sampling period, such as one or two weeks. Due to this limitation, ICESat-2 has been so far only used to generate monthly sea ice freeboard or thickness maps over the polar regions based on the near-monthly sub-cycle of its orbit [31,[35][36][37]. Even though these monthly maps are useful in understanding the overall annual trend of sea ice thickness, they cannot capture short-term formation and dynamics of sea ice, especially in regions where sea ice moves very rapidly, such as near the polynya ice production areas. Since the wind speed and wind direction in the Ross Sea change very rapidly even on a daily scale [38][39][40], the spatial distribution of sea ice freeboard over the Ross Sea can fluctuate significantly in a short period. Thus, weekly maps of the sea ice freeboard over the Ross Sea should provide an improved characterization of sea ice production and drift in the polynya areas when compared to monthly maps.
Remote Sens. 2021, 13, x FOR PEER REVIEW 2 of 18
In particular, the Ross Sea is one of the most noticeable regions showing a significant increase of sea ice extent in the Southern Ocean [17,21,22]. There, massive amounts of sea ice are produced in the Ross Sea polynya which is the largest recurring polynya area in Antarctica [23][24][25]. Thus, sea ice dynamics in the Ross Sea has been the subject of a number of scientific studies [1,[26][27][28][29][30]. Due to the harsh weather and long polar nights, spaceborne remote sensing has been the most common and efficient approach for these studies; spaceborne remote sensing allows obtaining sea ice observations over a large area with regular time spans without the need of physical visits to the study region.
In recent years, the National Aeronautics and Space Administration (NASA)'s Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) has been providing precise sea ice freeboard and thickness estimates over the polar regions. ICESat-2 was launched in September 2018 with the objectives of measuring changes in ice sheets, glaciers, and sea ice thickness [31]. The Advanced Topographic Laser Altimeter System (ATLAS) on ICESat-2 measures the surface heights by using 532 nm wavelength laser photons. ATLAS has six beams of three pairs of strong and weak beams with 11 m footprints spaced by 0.7 m [32], providing a better spatial resolution compared to the previous ICESat mission [33]. Various studies have utilized ICESat-2 for the characterization of sea ice freeboard and thickness in the polar regions, taking advantage of this outstanding spatial resolution [34][35][36].
However, despite the good precision and resolution of ICESat-2, the use of ICESat-2 for studying sea ice is substantially limited by the spatial sparsity of the instrument's measurements (Figure 1). Since ICESat-2 data is collected along tracks, there are large sections of area for which there is no data because they are not under a track. In addition, ICESat-2 cannot measure freeboard in cloudy conditions because the laser photons do not penetrate clouds. Thus, ICESat-2 data have limitations in sampling sea ice elevation over a large area during a short sampling period, such as one or two weeks. Due to this limitation, ICESat-2 has been so far only used to generate monthly sea ice freeboard or thickness maps over the polar regions based on the near-monthly sub-cycle of its orbit [31,[35][36][37]. Even though these monthly maps are useful in understanding the overall annual trend of sea ice thickness, they cannot capture short-term formation and dynamics of sea ice, especially in regions where sea ice moves very rapidly, such as near the polynya ice production areas. Since the wind speed and wind direction in the Ross Sea change very rapidly even on a daily scale [38][39][40], the spatial distribution of sea ice freeboard over the Ross Sea can fluctuate significantly in a short period. Thus, weekly maps of the sea ice freeboard over the Ross Sea should provide an improved characterization of sea ice production and drift in the polynya areas when compared to monthly maps. In order to address the spatial sparsity of ICESat-2 and produce weekly freeboard maps over the Ross Sea, we employ a geostatistical kriging approach. Kriging is a popular and well-established spatial interpolation methods for estimating unknown values [41]. In addition, one of the advantages of using kriging is that it is useful to examine the spatial autocorrelation of data (e.g., statistically representative range of ICESat-2 data) [42]. Since a mining engineer, Krige, used statistics to estimate ore reserves [43], geostatistical analysis and the kriging method have been used for the mapping of many environmental fields. In terms of the cryosphere, Herzfeld et al. [44] and Iacozza and Barber [45] characterized the snow depth on sea ice from geostatistical variograms, and Huang et al. [46] estimated the spatial distribution of snow depth by using kriging. Kriging is a common method to interpolate ice shelf thickness [47,48] or ice sheet mass balance [49][50][51]. Lindsay et al. [52] and Kurtz et al. [53] used the ordinary kriging to interpolate sea surface heights for the estimation of sea ice thickness, and Wang and Hou [54] used kriging to determine the spatial distribution of surface temperature in the Antarctic. Herzfeld et al. [55] used variograms for the morphological characterization of ice surface types from glacier-roughness sensor (GRS) data.
Although geostatistics is common for various applications, there have been few attempts to apply this geostatistical kriging to short-term mapping of sea ice freeboard or thickness. In this study, we characterize the spatial autocorrelation of the ICESat-2 freeboard measurements in the Ross Sea and apply kriging interpolation to estimate sea ice freeboards in locations where the ICESat-2 track is not available. Finally, based on this geostatistical approach, we generate weekly freeboard maps and analyze the short-term sea ice variability in the Ross Sea in 2019.
The area of interest for this study is the Ross Sea, which is the southernmost sea of Antarctica located in a deep embayment of the Southern Ocean between Victoria Land to the west and Marie Byrd Land to the east (Figure 2). To generate the weekly maps, the Ross Sea is subjectively limited to the region 78 • S-70 • S latitude, and between the Victoria Land/Oates Land coast and the 150 • W longitude. In addition, considering the general distribution of polynyas and ocean currents in this region, we divide the entire Ross Sea area into four sub-sectors: northwest (NW), northeast (NE), southwest (SW), and southeast (SE) as shown in Figure 2.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 18
In order to address the spatial sparsity of ICESat-2 and produce weekly freeboard maps over the Ross Sea, we employ a geostatistical kriging approach. Kriging is a popular and well-established spatial interpolation methods for estimating unknown values [41]. In addition, one of the advantages of using kriging is that it is useful to examine the spatial autocorrelation of data (e.g., statistically representative range of ICESat-2 data) [42]. Since a mining engineer, Krige, used statistics to estimate ore reserves [43], geostatistical analysis and the kriging method have been used for the mapping of many environmental fields. In terms of the cryosphere, Herzfeld et al. [44] and Iacozza and Barber [45] characterized the snow depth on sea ice from geostatistical variograms, and Huang et al. [46] estimated the spatial distribution of snow depth by using kriging. Kriging is a common method to interpolate ice shelf thickness [47,48] or ice sheet mass balance [49][50][51]. Lindsay et al. [52] and Kurtz et al. [53] used the ordinary kriging to interpolate sea surface heights for the estimation of sea ice thickness, and Wang and Hou [54] used kriging to determine the spatial distribution of surface temperature in the Antarctic. Herzfeld et al. [55] used variograms for the morphological characterization of ice surface types from glacier-roughness sensor (GRS) data.
Although geostatistics is common for various applications, there have been few attempts to apply this geostatistical kriging to short-term mapping of sea ice freeboard or thickness. In this study, we characterize the spatial autocorrelation of the ICESat-2 freeboard measurements in the Ross Sea and apply kriging interpolation to estimate sea ice freeboards in locations where the ICESat-2 track is not available. Finally, based on this geostatistical approach, we generate weekly freeboard maps and analyze the short-term sea ice variability in the Ross Sea in 2019.
The area of interest for this study is the Ross Sea, which is the southernmost sea of Antarctica located in a deep embayment of the Southern Ocean between Victoria Land to the west and Marie Byrd Land to the east (Figure 2). To generate the weekly maps, the Ross Sea is subjectively limited to the region 78°S-70°S latitude, and between the Victoria Land/Oates Land coast and the 150°W longitude. In addition, considering the general distribution of polynyas and ocean currents in this region, we divide the entire Ross Sea area into four sub-sectors: northwest (NW), northeast (NE), southwest (SW), and southeast (SE) as shown in Figure 2. One of the important characteristics of the Ross Sea is its large polynyas, mostly in the SW sector [23,24]. Polynyas can be defined as isolated areas of open water and thin ice within the ice pack in a polar region [56]. The Ross Sea has three large persistent polynyas (Figure 2): the Ross Ice Shelf polynya (RISP), Terra Nova Bay polynya (TNBP), and McMurdo Sound polynya (MSP) [26,40]. The formation of these polynyas from March to November is largely related to katabatic winds from the Ross Ice Shelf and Victoria Land [57,58].
In this study, sea ice freeboard data over the Ross Sea region from the ICESat-2 ATL10 (release 003) were obtained from NASA Earthdata (earthdata.nasa.gov). The time period selected for this study was 1 March to 30 November 2019 which covered one annual icegrowth cycle from the start of the freezing season to the start of the melting season. We note that ICESat-2 data from 26 June 2019 to 9 July 2019 were not available because ATLAS was off during this time in the safe-hold mode.
The ATL10 freeboard product is derived from the sea ice surface heights product (ATL07) which in turn was derived by the aggregation of 150 photons from the ATLAS global geolocated photon data product (ATL03). The surface heights of ATL07 are calculated by assuming a Gaussian distribution of the aggregated 150 photons and the precision of this surface height is known as less than 2 cm for flat surfaces [59]. Based on this surface height h s , the sea ice freeboard h f of the ATL10 product is calculated as follows:
where h re f is a sea surface reference height estimated from a collection of sea ice leads within a 10 km segment. Herein, this freeboard represents the total freeboard (i.e., the ice freeboard and snow depth), and freeboard is calculated only for where daily sea ice concentration (SIC) is greater than 15% [60]. In the ATL10 product, sea ice floes and leads are classified by three photon factors: photon rate, photon distribution, and background rate [37]. Since the surface heights are derived from the aggregations of 150 geolocated photons, the spatial resolution (i.e., height segment length of 150 photon aggregations) depends on the surface reflectance or roughness. Generally, the segment lengths are 10-200 m for strong beams, and 40-800 m for weak beams [37]. Because the strong beams have better resolutions than weak beams, only strong beams are used in this study.
The purpose of this study is to characterize and estimate the overall distribution of sea ice freeboard over the entire Ross Sea by using ICESat-2 observations and the kriging approach. However, for some sea ice floes, ridged ice or thin ice in leads can distort the global freeboard values for level sea ice (Figure 3a). Thus, the ATL10 freeboard data are first resampled every 3 km to reduce the impacts of sea ice deformation (ridged ice or sea ice leads) on the level ice freeboard values (Figure 3b). Geostatistical approaches can have difficulty predicting these local and linear sea ice deformation features because they are based on the assumption of data stationarity [61].
It is known that the sea ice freeboard generally follows a lognormal distribution [62] and the modal freeboard represents the freeboard for level ice [63,64]. Therefore, in order to reduce the effect of the sea ice deformation, the freeboard distributions are fitted to lognormal distributions for every 3 km cell and the modal freeboards for the 3 km cells are considered as representative freeboards of the level sea ice (Figure 3b-d). Since three strong beams of ICESat-2 are separated with 3.3 km intervals, we make the ICESat-2 track data look like ~3 × 3 km 2 pseudo-square cells (Figure 3b).
After the ATL10 data are aggregated as modal freeboards every 3 km, we build geostatistical models for interpolating freeboard values to select locations where ICESat-2 observations are not available. A geostatistical model first measures the spatial dependence or autocorrelation from a semivariogram that represents the relationship between the semivariance γ(h) and the lag distance h. The average semivariance at a lag distance h is defined by:
where ( ) is the measured sample value (i.e., ATL10 freeboard in this study) at the point , ( + ℎ) is the measured sample value at another point displaced from the point by a lag distance h, and N is the number of observations for a lag distance h. After an experimental semivariogram is derived from the measurements, the experimental semivariogram must be fitted with a mathematical model to estimate the semivariance at any given distance [42]. This fitted semivariogram has three elements: nugget, effective range, and sill. The nugget ( ) is the semivariance at the distance 0, which means measurement error. The effective range is the distance at which the semivariogram starts to be flat, which means the data have no spatial correlation beyond this range. The sill is the semivariance at this leveling point. Among various types of semivariogram models, three representative variogram models are calculated and compared in this study: exponential, Gaussian, and spherical. These models are defined by the equations:
After the ATL10 data are aggregated as modal freeboards every 3 km, we build geostatistical models for interpolating freeboard values to select locations where ICESat-2 observations are not available. A geostatistical model first measures the spatial dependence or autocorrelation from a semivariogram that represents the relationship between the semivariance γ(h) and the lag distance h. The average semivariance at a lag distance h is defined by:
where z(x i ) is the measured sample value (i.e., ATL10 freeboard in this study) at the point x i , z(x i + h) is the measured sample value at another point displaced from the point x i by a lag distance h, and N is the number of observations for a lag distance h. After an experimental semivariogram is derived from the measurements, the experimental semivariogram must be fitted with a mathematical model to estimate the semivariance at any given distance [42]. This fitted semivariogram has three elements: nugget, effective range, and sill. The nugget (c 0 ) is the semivariance at the distance 0, which means measurement error. The effective range is the distance at which the semivariogram starts to be flat, which means the data have no spatial correlation beyond this range. The sill is the semivariance at this leveling point. Among various types of semivariogram models, three representative variogram models are calculated and compared in this study: exponential, Gaussian, and spherical. These models are defined by the equations:
Based on these variogram models, ordinary kriging is employed to estimate sea ice freeboards. The general equation for estimating the z value at a point is:
where Z 0 is the estimated value, Z x is the known value at point x, W x is the weight on the point x, and s is the number of sample points for estimation. The weights are derived from solving a set of simultaneous equations [42]. In this study, Esri's ArcGIS software is used for the derivation of semivariogram models and ordinary kriging.
In order to cross-validate these kriging models, we divide the weekly data into a training dataset and validation dataset. For every week, we randomly select two tracks as a validation dataset and build kriging models using the other available tracks (Figure 4). This model is then compared to the validation dataset to check if this geostatistical approach is reliable for the estimation of the unknown values. In this process, we compare the three different variogram models to determine which model is the most accurate. After the accuracy of the geostatistical model is assessed, we map sea ice freeboard over the Ross Sea weekly. These freeboard maps have a grid size of 25 km which is comparable to previous sea ice mapping results [35,36,65] and the ICESat-2 Level 3 product of monthly gridded sea ice freeboard (ATL20) [60]. The freeboard for the center of each 25 km grid cell is estimated by applying kriging with the 3 km mode sampled ATL10 data. We note that the mapped freeboard in this study is the total freeboard for level sea ice because the impacts of sea ice deformation (i.e., ridges and leads) are minimized by the previous 3 km mode sampling (Section 4.1). The number of ATL10 tracks used for each weekly map is generally 15-25; the number of used tracks depends on the time and location of ground tracks and cloud cover.
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 18
Gaussian:
Based on these variogram models, ordinary kriging is employed to estimate sea ice freeboards. The general equation for estimating the z value at a point is:
where is the estimated value, is the known value at point x, is the weight on the point x, and s is the number of sample points for estimation. The weights are derived from solving a set of simultaneous equations [42]. In this study, Esri's ArcGIS software is used for the derivation of semivariogram models and ordinary kriging.
In order to cross-validate these kriging models, we divide the weekly data into a training dataset and validation dataset. For every week, we randomly select two tracks as a validation dataset and build kriging models using the other available tracks (Figure 4). This model is then compared to the validation dataset to check if this geostatistical approach is reliable for the estimation of the unknown values. In this process, we compare the three different variogram models to determine which model is the most accurate. After the accuracy of the geostatistical model is assessed, we map sea ice freeboard over the Ross Sea weekly. These freeboard maps have a grid size of 25 km which is comparable to previous sea ice mapping results [35,36,65] and the ICESat-2 Level 3 product of monthly gridded sea ice freeboard (ATL20) [60]. The freeboard for the center of each 25 km grid cell is estimated by applying kriging with the 3 km mode sampled ATL10 data. We note that the mapped freeboard in this study is the total freeboard for level sea ice because the impacts of sea ice deformation (i.e., ridges and leads) are minimized by the previous 3 km mode sampling (Section 4.1). The number of ATL10 tracks used for each weekly map is generally 15-25; the number of used tracks depends on the time and location of ground tracks and cloud cover.
Figure 5 describes the cross-validation results by time. In general, the correlation coefficient (R) is greater than 0.4, which indicates a moderate correlation [66]. R value is relatively greater in the thicker-ice season from July to November (R > 0.5), but it is relatively low in the thinner-ice season from March to June (R ~0.4-0.6) (Figure 5a). Similarly, while the root mean square error (RMSE) and mean absolute error (MAE) are relatively low in March to June, they are relatively high in July to November (Figure 5b,c). The higher errors in July to November may be attributed to higher freeboard values (i.e., thicker ice) during these months.
Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 18
Figure 5 describes the cross-validation results by time. In general, the correlation coefficient (R) is greater than 0.4, which indicates a moderate correlation [66]. R value is relatively greater in the thicker-ice season from July to November (R > 0.5), but it is relatively low in the thinner-ice season from March to June (R ~ 0.4-0.6) (Figure 5a). Similarly, while the root mean square error (RMSE) and mean absolute error (MAE) are relatively low in March to June, they are relatively high in July to November (Figure 5b,c). The higher errors in July to November may be attributed to higher freeboard values (i.e., thicker ice) during these months. The overall cross-validation results of the kriging freeboard estimation for three variogram models are depicted in Figure 6. The exponential, Gaussian, and spherical variogram models show similar results to each other: an R of 0.56-0.57, RMSE of ~0.12 m, MAE of 0.07 m. Although this level of R value indicates only moderate correlation between kriged freeboard and actual freeboard, the p-value (< 0.001) shows that this correlation is significant. Additionally, this cross-validation result implies that the geostatistical approach interpolates and estimates of freeboards in the Ross Sea with an expected RMSE of approximately 0.12 m, and the differences between the three variograms are negligible. The overall cross-validation results of the kriging freeboard estimation for three variogram models are depicted in Figure 6. The exponential, Gaussian, and spherical variogram models show similar results to each other: an R of 0.56-0.57, RMSE of ~0.12 m, MAE of 0.07 m. Although this level of R value indicates only moderate correlation between kriged freeboard and actual freeboard, the p-value (< 0.001) shows that this correlation is significant. Additionally, this cross-validation result implies that the geostatistical approach interpolates and estimates of freeboards in the Ross Sea with an expected RMSE of approximately 0.12 m, and the differences between the three variograms are negligible. Furthermore, we compare the cross-validation parameters for each subsector: NW, SW, NE, and SE (Figure 7). As shown in Figure 7a, the SW and SE sectors show higher R values than the NW and NE sectors. This is because the ICESat-2 tracks are more densely distributed in the southern regions (Figure 1), which consequently guarantees higher reliability of kriging interpolation for the southern sectors. However, since the RMSE and MAE depend on the freeboard value as well, the SE sector, where the thickest ice is located, shows a higher RMSE and MAE compared to the other sectors (Figure 7b,c). The SW sector generally shows the lowest RMSE and MAE among the four sectors, which should be attributed lower freeboard of the SW sector due to the existence of large polynyas there. More details about ice freeboard and thickness for these sectors and the reason for it are discussed in the following sections. Furthermore, we compare the cross-validation parameters for each subsector: NW, SW, NE, and SE (Figure 7). As shown in Figure 7a, the SW and SE sectors show higher R values than the NW and NE sectors. This is because the ICESat-2 tracks are more densely distributed in the southern regions (Figure 1), which consequently guarantees higher reliability of kriging interpolation for the southern sectors. However, since the RMSE and MAE depend on the freeboard value as well, the SE sector, where the thickest ice is located, shows a higher RMSE and MAE compared to the other sectors (Figure 7b,c). The SW sector generally shows the lowest RMSE and MAE among the four sectors, which should be attributed lower freeboard of the SW sector due to the existence of large polynyas there. More details about ice freeboard and thickness for these sectors and the reason for it are discussed in the following sections. Furthermore, we compare the cross-validation parameters for each subsector: NW, SW, NE, and SE (Figure 7). As shown in Figure 7a, the SW and SE sectors show higher R values than the NW and NE sectors. This is because the ICESat-2 tracks are more densely distributed in the southern regions (Figure 1), which consequently guarantees higher reliability of kriging interpolation for the southern sectors. However, since the RMSE and MAE depend on the freeboard value as well, the SE sector, where the thickest ice is located, shows a higher RMSE and MAE compared to the other sectors (Figure 7b,c). The SW sector generally shows the lowest RMSE and MAE among the four sectors, which should be attributed lower freeboard of the SW sector due to the existence of large polynyas there. More details about ice freeboard and thickness for these sectors and the reason for it are discussed in the following sections.
We also compare the temporal variations of the variogram factors (i.e., range and sill). As shown in Figure 8a, the effective ranges of the variograms generally range in 5-10 km, but the exponential variogram model shows greater range than the other variogram models because of its mathematical approach (Equation ( 3)). The exponential model takes a longer distance to reach the level sill points than the other models due to the mathematical characteristics, even though they have all the same experimental semivariograms. Since the ranges of the variograms represent the spatially autocorrelated distance, it can be said that the ICESat-2 ATL10 freeboards are representative of freeboards on scales of 5-10 km surrounding sea ice floes in the Ross Sea.
takes a longer distance to reach the level sill points than the other models due to the mathematical characteristics, even though they have all the same experimental semivariograms. Since the ranges of the variograms represent the spatially autocorrelated distance, it can be said that the ICESat-2 ATL10 freeboards are representative of freeboards on scales of 5-10 km surrounding sea ice floes in the Ross Sea.
The sills increase from March to November for all three variogram models (Figure 8b). The sills are less than 0.01 m 2 in April-June, but it increases from July, so reaches up to 0.02-0.03 m 2 in November. Since sills can be regarded as the maximum variations of the kriging interpolations, this result indicates that the variations of the kriging estimation are larger in winter to spring than in autumn. Given that sea ice is thicker in winter to spring, it is reasonable that winter and spring show large variations. Moreover, although we attempt to reduce the impacts of the sea ice deformation by applying 3 km mode sampling, some sea ice deformation (ridges) can still introduce higher variations during the winter to spring period because a significant portion of sea ice is heavily deformed in the Ross Sea [67]. Meanwhile, it is also noted that late October to November shows the largest sill values. Since late October to November is the start of melting season, we can deduce that the melting of sea ice and reduction of the sea ice covers can lead to higher variations of this kriging interpolation. The sills increase from March to November for all three variogram models (Figure 8b). The sills are less than 0.01 m 2 in April-June, but it increases from July, so reaches up to 0.02-0.03 m 2 in November. Since sills can be regarded as the maximum variations of the kriging interpolations, this result indicates that the variations of the kriging estimation are larger in winter to spring than in autumn. Given that sea ice is thicker in winter to spring, it is reasonable that winter and spring show large variations. Moreover, although we attempt to reduce the impacts of the sea ice deformation by applying 3 km mode sampling, some sea ice deformation (ridges) can still introduce higher variations during the winter to spring period because a significant portion of sea ice is heavily deformed in the Ross Sea [67]. Meanwhile, it is also noted that late October to November shows the largest sill values. Since late October to November is the start of melting season, we can deduce that the melting of sea ice and reduction of the sea ice covers can lead to higher variations of this kriging interpolation.
Based on the cross-validation results, we generate weekly freeboard maps over the Ross Sea by using the spherical variogram model. Figure 9 shows the weekly variations of the sea ice freeboard for each sector retrieved from these freeboard maps. In April, the average freeboard for the entire Ross Sea is approximately 0.15 m, and it increases by October. October shows the maximum average sea ice freeboards up to 0.25 m. After October, freeboard started to decrease and returned to under 0.2 m in November.
It is notable that the four sectors in the Ross Sea show their own distinctive freeboard variations. The SE sector shows the highest and most variable freeboard values compared to the other sectors. In April to May, the average of the freeboard for the SE sector is approximately 0.12-0.15 m, but it started to increase from late May. Consequently, the averaged freeboard for the SE sector reached up to 0.5 m in the second week of October, while the other sectors show freeboard less than 0.25 m. Contrary to the SE sector, the SW sector shows the lowest freeboard values for the whole year: less than 0.15 m. The freeboard in the SW sector shows relatively constant values of freeboard, without large increases or decreases. In April to June, the four sectors show only small differences in their freeboard. During June to November, however, the SE sector shows the highest increases in freeboard values, followed by NE, NW, and SW. In order to better understand the spatial and temporal variations of freeboard, we compare weekly freeboard maps and their distributions over the Ross Sea for each sector and for three seasons: the slow-freezing season (April to June), fast-freezing season (July to September), and melting season (October to November). Freeboard maps and distributions in April to June are shown in Figure 10. From April 12 to May 16, the freeboard distributions and the mean freeboard values look similar without large differences between them. Nevertheless, it is noted that the SW sector generally shows bimodal distributions such that the lower peak freeboard (~0.05 m) is lower than the other sectors (~0.1-0.15 m). This is attributed to the existence of large polynyas in the Ross Sea, especially the RISP. Clearly, the thin ice region near the RISP remains during the entire slow-freezing season from April to June (Figure 10). Two other polynyas (the TNBP and MSP) should be occasionally observed in the Ross Sea, but they are hardly overlapped with ICESat-2 tracks because of their relatively small spatial scales compared to the RISP (Figure 2). Thus, these relatively small polynyas can be only observed if the ICESat-2 tracks are exactly coincident with the occurrence of these polynya events. There was no coincidence between the ICESat-2 tracks and the TNBP or MSP polynya events It is notable that the four sectors in the Ross Sea show their own distinctive freeboard variations. The SE sector shows the highest and most variable freeboard values compared to the other sectors. In April to May, the average of the freeboard for the SE sector is approximately 0.12-0.15 m, but it started to increase from late May. Consequently, the averaged freeboard for the SE sector reached up to 0.5 m in the second week of October, while the other sectors show freeboard less than 0.25 m. Contrary to the SE sector, the SW sector shows the lowest freeboard values for the whole year: less than 0.15 m. The freeboard in the SW sector shows relatively constant values of freeboard, without large increases or decreases. In April to June, the four sectors show only small differences in their freeboard. During June to November, however, the SE sector shows the highest increases in freeboard values, followed by NE, NW, and SW. In order to better understand the spatial and temporal variations of freeboard, we compare weekly freeboard maps and their distributions over the Ross Sea for each sector and for three seasons: the slowfreezing season (April to June), fast-freezing season (July to September), and melting season (October to November).
Freeboard maps and distributions in April to June are shown in Figure 10. From April 12 to May 16, the freeboard distributions and the mean freeboard values look similar without large differences between them. Nevertheless, it is noted that the SW sector generally shows bimodal distributions such that the lower peak freeboard (~0.05 m) is lower than the other sectors (~0.1-0.15 m). This is attributed to the existence of large polynyas in the Ross Sea, especially the RISP. Clearly, the thin ice region near the RISP remains during the entire slow-freezing season from April to June (Figure 10). Two other polynyas (the TNBP and MSP) should be occasionally observed in the Ross Sea, but they are hardly overlapped with ICESat-2 tracks because of their relatively small spatial scales compared to the RISP (Figure 2). Thus, these relatively small polynyas can be only observed if the ICESat-2 tracks are exactly coincident with the occurrence of these polynya events. There was no coincidence between the ICESat-2 and the TNBP or MSP polynya events during our study period. Although there is no large increase or decrease in the freeboard for all four sectors in this slow-freezing season, freeboard starts to increase in the eastern sectors (NE and SE) from the week of 24 May. In particular, the freeboard starts to increase mainly from the east side of the SE sector, which indicates the inflow of thick multiyear ice from the Amundsen Sea, east of the Ross Sea area [68][69][70].
Remote Sens. 2021, 13, x FOR PEER REVIEW 11 of 18 during our study period. Although there is no large increase or decrease in the freeboard for all four sectors in this slow-freezing season, freeboard starts to increase in the eastern sectors (NE and SE) from the week of 24 May. In particular, the freeboard starts to increase mainly from the east side of the SE sector, which indicates the inflow of thick multiyear ice from the Amundsen Sea, east of the Ross Sea area [68][69][70].
Figure 11 shows the weekly freeboard maps and distribution changes during the fastfreezing season from July to September. The SW sector still shows a bimodal distribution during this season, and the first lower mode has approximately 0.05 m of freeboard, which indicates the thin ice around the RISP in this sector. The freeboard in the SE sector has dramatically increased during this period until the first week of September. The SE sector also shows a bimodal distribution, but both modal freeboards continuously shift to the right, unlike the consistent bimodal freeboard values for the SW sector. Additionally, while the right mode increases much faster, the left mode moves with similar increase rates to the NW and NE sector. It is logical that increase in freeboard in the SE sector is mainly led by two different sources: (1) inflow of thick ice from the Amundsen Sea (greater mode), and (2) internal ice growth in the SE sector (lower mode). Thus, the SE sector shows a very wide range of freeboard compared to the other sectors. The NE and NW sectors usually show unimodal distributions, but sometimes they also show bimodal distributions. Specifically, the thin ice from the SW sector is sometimes advected into the NW sector, shown as the small mode in the distribution, and thicker ice from the SE sector is sometimes advected into the NE sector, shown as the greater mode in the distribution in these weeks.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 18 also shows a bimodal distribution, but both modal freeboards continuously shift to the right, unlike the consistent bimodal freeboard values for the SW sector. Additionally, while the right mode increases much faster, the left mode moves with similar increase rates to the NW and NE sector. It is logical that the increase in freeboard in the SE sector is mainly led by two different sources: (1) inflow of thick ice from the Amundsen Sea (greater mode), and (2) internal ice growth in the SE sector (lower mode). Thus, the SE sector shows a very wide range of freeboard compared to the other sectors. The NE and NW sectors usually show unimodal distributions, but sometimes they also show bimodal distributions. Specifically, the thin ice from the SW sector is sometimes advected into the NW sector, shown as the small mode in the distribution, and thicker ice from the SE sector is sometimes advected into the NE sector, shown as the greater mode in the distribution in these weeks. Freeboards reached a maximum in the week of 11 October and began to decrease after that (Figure 9). During this melting season, the SW sector contains the lowest freeboard, and we can still identify the spatial distribution of the polynya-related thin ice in this sector (Figure 12). The NW sector shows the second lowest freeboard due to the impacts of thin ice advected from the SW sector. The SE sector shows the highest freeboard and the widest distribution of the freeboard, but the freeboard decreases from the week of 24 November. An interesting finding here is the highest freeboard region shifts northward in the SE sector. During the weeks of 11-31 October, the highest freeboard regions are located in the southeast coastal area of the SE sector. As time passes, however, we can observe that these highest freeboard regions are shifted to the northern side of the SE sector (i.e., border between the SE and NE sectors) in the weeks from 1 November to 28 November. This can be associated with the movement of thick sea ice towards the north of the Ross Sea. According to previous studies [68,69,71], sea ice moves northward in this sector, so thick ice was pushed away to the north side. Freeboards reached a maximum in the week of 11 October and began to decrease after that (Figure 9). During this melting season, the SW sector contains the lowest freeboard, and we can still identify the spatial distribution of the polynya-related thin ice in this sector (Figure 12). The NW sector shows the second lowest freeboard due to the impacts of thin ice advected from the SW sector. The SE sector shows the highest freeboard and the widest distribution of the freeboard, but the freeboard decreases from the week of 24 November. An interesting finding here is that the highest freeboard region shifts northward in the SE sector. During the weeks of 11-31 October, the highest freeboard regions are located in the southeast coastal area of the SE sector. As time passes, however, we can observe that these highest freeboard regions are shifted to the northern side of the SE sector (i.e., border between the SE and NE sectors) in the weeks from 1 November to 28 November. This can be associated with the movement of thick sea ice towards the north of the Ross Sea. According to previous studies [68,69,71], sea ice moves northward in this sector, so thick ice was pushed away to the north side.
We compare the weekly freeboard estimated by kriging with the National Ice Center (NIC) weekly ice chart (Figure 13). The NIC ice charts are produced through the manual interpretation of various data sources including in situ, remote sensing (e.g., visible, infrared, active/passive microwave), and model outputs [72]. Although they do not provide accurate thickness measurements, they consist of multiple polygons showing the categorical stage of development of sea ice (e.g., nilas, grey ice, young ice, first-year thin ice, first-year thick ice, old ice, etc.). This ice chart data has been used to examine the spatiotemporal variations of ice thickness and volume in the Ross Sea area [69]. In this study, we convert the categorical ice type of each polygon into numerical ice thickness values using the method of DeLiberty et al. [69] Thus, since each polygon of the ice chart has a single thickness value, we compare the converted ice thickness for each polygon with the ICESat-2 freeboard averaged over that polygon.
thickness value, we compare the converted ice thickness for each polygon with the ICESat-2 freeboard averaged over that polygon.
As shown in Figure 13, the kriged freeboard and the NIC ice thickness show a significant correlation (R = 0.587, p-value < 0.001) to each other. This implies that our weekly freeboard estimation agrees with the existent sea ice thickness products. Furthermore, the slope of this linear fit (2.61) is close to that from the literature that converted sea ice freeboard into ice thickness for the Ross Sea area in early spring (2.45) [73].
However, it should be noted that they do not have a perfect linear relationship because of the major differences between them. First, while our freeboard maps represent the total freeboard, the NIC ice charts represent the ice thickness. Considering the variations of snow/ice density or snow depth over the Ross Sea [74], the freeboard and ice thickness can be deviated from the complete linear fit. Second, while our freeboard maps are based on the quantitative measures from ICESat-2, the NIC ice charts are based on the manual interpretation of various data. Thus, the thickness from the NIC ice chart is not an accurate reference but a rough categorical estimation of sea ice thickness. Third, although we try to reduce the impacts of sea ice deformation using the 3 km mode sampling, our freeboard maps can still involve these dynamic impacts. In contrast, the thickness from the NIC ice chart is merely based on the thermodynamic ice growth [69]. Therefore, herein, we only check a rough coherency between our estimation and the NIC ice chart, rather than the accurate validation of our freeboard maps. We would suggest the weekly ICESat-2 freeboard maps, if they can be produced in near real-time, to be incorporated into a realtime weekly NIC ice chart to improve the potential accuracy of the chart. As shown in Figure 13, the kriged freeboard and the NIC ice thickness show a significant correlation (R = 0.587, p-value < 0.001) to each other. This implies that our weekly freeboard estimation agrees with the existent sea ice thickness products. Furthermore, the slope of this linear fit (2.61) is close to that from the literature that converted sea ice freeboard into ice thickness for the Ross Sea area in early spring (2.45) [73].
However, it should be noted that they do not have a perfect linear relationship because of the major differences between them. First, while our freeboard maps represent the total freeboard, the NIC ice charts represent the ice thickness. Considering the variations of snow/ice density or snow depth over the Ross Sea [74], the freeboard and ice thickness can be deviated from the complete linear fit. Second, while our freeboard maps are based on the quantitative measures from ICESat-2, the NIC ice charts are based on the manual interpretation of various data. Thus, the thickness from the NIC ice chart is not an accurate reference but a rough categorical estimation of sea ice thickness. Third, although we try to reduce the impacts of sea ice deformation using the 3 km mode sampling, our freeboard maps can still involve these dynamic impacts. In contrast, the thickness from the NIC ice chart is merely based on the thermodynamic ice growth [69]. Therefore, herein, we only check a rough coherency between our estimation and the NIC ice chart, rather than the accurate validation of our freeboard maps. We would suggest the weekly ICESat-2 freeboard maps, if they can be produced in near real-time, to be incorporated into a real-time weekly NIC ice chart to improve the potential accuracy of the chart.
This study presents a geostatistical characterization and weekly mapping of total freeboard in the Ross Sea by using ICESat-2 ATL10 products. The spatial autocorrelation ranges of ATL10 freeboard are quantified from the geostatistical variograms, and freeboards are interpolated by ordinary kriging based on these variogram models. The resulting weekly freeboard maps successfully identify the spatial distributions of polynya areas in the Ross Sea.
However, the geostatistical approach used here can have several limitations. From the statistical point of view, this method can introduce significant uncertainties if the number of ICESat-2 observations are not sufficient in a given region [75,76]. In particular, since ICESat-2 laser photons cannot penetrate into clouds, overcast weather conditions can cause a lack of measurements. This is evidenced by a somewhat discontinuous appearance in some of the freeboard maps (Figures 101112).
Another issue of the geostatistical approach used in this study is that it assumes that the sea ice across the Ross Sea remains static over a week's period (i.e., stationarity). However, the sea ice distribution over the Ross Sea is very dynamic and can change significantly even within one week, especially in the SW sector where frequent katabatic wind events can drive significant sea ice production. Since our geostatistical kriging model is based on the assumption of static conditions for the training dataset, freeboard changes occurring with a shorter than one-week scale cannot be described by this approach.
In addition, it is also necessary to consider the assessment or validation of the ATL10 total freeboard product for the Ross Sea. Although some studies assessed the ATL10 freeboard data in comparison with other air-borne missions, such as Operation Ice Bridge (OIB) for the Arctic [59], the accuracy of ATL10 freeboard over the Ross Sea has not been evaluated. Due to the lack of validation data in the Ross Sea, the cross-validation of this study is conducted with the ATL10 freeboard product itself, rather than independent sea ice freeboard measurements. Therefore, it is necessary to assess the accuracy of the ATL10 freeboard product by using sea ice field data in the near future. In addition, herein we only analyze the freeboard values, not the sea ice thickness. Although there are several methods to convert freeboard into thickness for the Ross Sea area [28,73], it is not clear what method would best describe the sea ice conditions over Ross Sea for the entire year. Therefore, it is also critical to examine how to associate the freeboard changes with the thickness variations over the Ross Sea.
In this study, we implement a geostatistical analysis to examine the spatial autocorrelation characteristics of ICESat-2 ATL10 freeboards over the Ross Sea. Based on this spatial autocorrelation, we generate weekly freeboard maps by applying ordinary kriging to detect spatiotemporal variations of sea ice over the Ross Sea. Before applying kriging interpolation for mapping, we conduct a cross validation to check the reliability of this kriging method for freeboard mapping. According to the cross-validation results, the kriging-estimated freeboards show significant correlation (R > 0.5), RMSE ~0.12 m, and MAE ~0.7 m with the ATL10 freeboard measurements. Three different variogram models (exponential, Gaussian, and spherical) show insignificant differences in their accuracy. This cross-validation result shows that the geostatistical approach can be useful for estimating freeboards in missing spaces where ICESat-2 data are not directly available. The spatial correlation range of the variogam models indicates that the ATL10 freeboard product can statistically represent the freeboard on scales of 5-10 km from the source data track. For the area farther away from the ATL10 measurements than this effective range, the variograms show the sill values less than 0.01 m 2 for April to July, but greater than that for August to November. The uncertainty of this kriging method increases in winter to spring as the ice grows, but the largest uncertainty is expected for the melting season.
The resulting weekly maps are used to examine the spatiotemporal changes of sea ice freeboard over the Ross Sea. We divided the Ross Sea into four sectors (SW, SE, NW, and NE) to better characterize short-term regional sea ice changes. Our weekly freeboard maps capture the spatial distributions of large polynya areas (i.e., RISP) over the Ross Sea. The polynya areas are mainly located in the SW sector, but sometimes they extend to the SE or NW sectors. The SW sector shows the lowest and less variable freeboard changes with bimodal freeboard distributions due to the large polynya areas in this sector. In contrast, the SE sector shows the largest and more variable changes of the freeboard from the season to the melting season, which is partially attributed to the thick multiyear ice advected from the Amundsen Sea. Meanwhile, we also find that the thick ice in the SE sector drifts northward in the melting season. If our short-term freeboard mapping is combined with other multispectral or radar satellite data (e.g., Landsat, Sentinel-1, Sentinel-2), we are confident that we can efficiently estimate sea ice production in the polynya areas over the Ross Sea.