Discrimination Algorithm and Procedure of Snow Depth and Sea Ice Thickness Determination Using Measurements of the Vertical Ice Temperature Profile by the Ice-Tethered Buoys

Snow depth and sea ice thickness in the Polar Regions are significant indicators of climate change and have been measured over several decades by ice-tethered buoys. However, sea ice temperature profiles measured by ice-tethered buoys are rarely used to infer snow depth and sea ice thickness owing to the lack of automatic discrimination algorithms, restricting the use of the data for sea ice thermodynamics studies. In this study, snow depth and sea ice thickness were retrieved through the measurements of sea ice temperature profiles using discrimination algorithms of the change point and the maximum likelihood detection methods. The data measured by 50 ice-tethered buoys were used to evaluate the accuracy of the results determined by the algorithm. Influences on the seasonal sea ice thermodynamic state, vertical interval of temperature sensors on the buoys, and initial ice thickness on the estimation errors were also evaluated. The performance of the discrimination algorithm for the data from the Arctic and Antarctic regions was also compared. There were no identifiable differences between the estimation errors from the Arctic and Antarctica. Increases in both the interval of the temperature sensors and the initial ice thickness enlarged the error for the estimation of ice thickness. A procedure developed in this study strengthens the potential application of measurements from the ice-tethered buoys only with the measurements of the vertical temperature profile of the layer of snow-covered ice, but not the measurements of ice basal and surface positions using acoustic sounding.


Introduction
Sea ice is a sensitive indicator of climate change in the Polar Regions. The thickness of Arctic sea ice has declined dramatically in recent decades [1-3] and its accelerated melting of the sea ice in Arctic has led to changes in the ice-albedo feedback [4][5][6]. Observations of sea ice mass balance can improve the fundamental understanding of the role and sensitivity of sea ice in global climate change. In the study of sea ice mass balance, ice thickness is widely defined as the most integrated parameter for describing the sea ice conditions [7,8]. The snow on the sea ice can affect the surface albedo and heat exchange between the atmosphere and ice. The depth of snow is crucial for evaluating surface energy equilibrium and growth of sea ice, and for the retrieval algorithm of sea ice thickness using the data of satellite altimeter.    The installations and schematic diagrams of three ice-tethered buoys are shown in Figure 2.

IMB Data
The IMB can be considered as one of the most reliable and widely deployed observation instruments for sea ice mass balance. A potential deployment scheme for the IMB buoy can be described as follows. The IMBs are composed of a data logger, which is in charge of data acquisition, processing and remote transmission, a sea ice mass balance observation system and a thermistor string to measure sea ice temperature profiles. The IMB uses two acoustic sensors to measure sea ice

IMB Data
The IMB can be considered as one of the most reliable and widely deployed observation instruments for sea ice mass balance. A potential deployment scheme for the IMB buoy can be described as follows. The IMBs are composed of a data logger, which is in charge of data acquisition, processing and remote transmission, a sea ice mass balance observation system and a thermistor string to measure sea ice temperature profiles. The IMB uses two acoustic sensors to measure sea ice thickness and snow depth. These two acoustic rangefinder sounders are fixed on a supporting structure, which is designed to be deployed through a 0.10 m diameter ice hole. The upward looking sonar whose accuracy is 0.01 m is fixed extending from the ice hole at 3.8 m under the snow-ice interface to measure the distance from the sonar to the bottom of the ice to invert ice bottom growth and ablation. Another sensor with an accuracy of 0.001 m is installed at 1.2 m from the snow surface to measure the distance from the sensor to the top of the snow to invert snow depth. The thermistor string with a length of 4.5 m consists of three sub-strings, each of which has a length of 1.5 m. There are 45 temperature probes evenly arranged on the temperature string at intervals of 0.01 m. And the IMB data are obtained from Ice Mass Balance (IMB) Buoy Program. Here, we used 47 IMB buoys deployed from 2002 to 2014. There were 11 buoys completing stable operation for more than one year and a buoy running for more than two years (25 months).

TUT Buoy Data
The Polar Research Institute of China and Taiyuan University of Technology have started to develop and deploy a new type of ice-tethered buoy (TUT buoy) since CHINARE-2014 in the Arctic Ocean. The TUT buoy uses digital temperature sensor chips with a resolution of 0.0625 • C and an accuracy of 0.1 • C. A 4.5-m-long thermistor string is linked by 15 flexible printed circuit boards (FPCBs) with a length of 0.3 m. Two acoustic rangefinder sounders (Campbell SR50A and Teledyne-Benthos PSA916, respectively) were used to measure snow depth and ice thickness. Based on a 32-bit low-power microcontroller, we designed the data acquisition instrument and the iridium satellite data transmission system, which allowed data to be sent back to the lab computer in real time. The TUT buoy (TUT-2016A) used in this paper was deployed over multiyear ice in the Arctic on 9 August 2016 during the CHINARE-2016. The measurements, including snow depth, ice thickness, and sea ice temperature profiles of TUT-2016A, lasted from 9 August 2016 to 11 July 2017, with the data sampling interval of one hour.

Determination of the Sea Ice Thickness/Snow Depth and Interfaces Among Air, Snow, Ice, and Water
The thickness of sea ice and the depth of snow can be calculated by determination of the media interfaces: top interface of air-snow or air-ice (when the snow has melted completely), snow-ice interface and the lower interface of ice-water. The interfaces can be derived from the temperature profiles observed by the ice-tethered buoys. Through identifying changes in the positions of interfaces, snow depth and sea ice thickness can be obtained.
The snow-ice interface is considered to be constant except for slush formation associated with flooding [22]. This phenomenon can be observed by the heating mode of the SIMBA. In this paper, IMB data with a temperature sensor interval of 0.10 m was mainly analyzed. The interval of TUT buoy is 0.03 m, but similar to the IMB, the TUT buoy does not have a heating mode as SIMBA to identify the interface of snow-ice and observe the evolution of snow to ice. Thus, the algorithm to identify the snow-ice interface could not be thoroughly considered and we did not develop the algorithm to solve this problem in this paper. Hereafter, we simply assumed the position of the snow-ice interface was constant.
The interfaces between air and snow or ice, and between ice and water were identified using the temperature profiles, the daily amplitude, which are equal to the differences between the maximum and minimum values of the daily temperature profiles of each sensor, and the vertical temperature gradients. The temperature profile, daily amplitude of temperature profiles and the vertical gradients of the temperature were widely used to estimate the interface between air and snow (or ice) [16,21,22]. The vertical gradient in air was different from in snow and ice. The daily amplitude of temperature profile was larger in air than in snow. The temperature profiles were measured by the buoys, and the daily amplitude of temperature profiles and the vertical temperature gradients were calculated from the temperature profile.

Interface between Air and Snow or Sea Ice
The evolution of snow depth over sea ice may be affected by synoptic processes phenomena such as storms, snowfall, and sleet, or due to melting caused by solar radiation. Thus, the temporal fluctuations of snow depth are much larger than for sea ice thickness. The daily amplitude of temperature profiles and the vertical gradients of temperature were examined because, in some cases, it was difficult to accurately distinguish the top interface by using the measurements of the temperature alone. From the data of the daily amplitude of temperature profiles and the vertical gradients of temperature, we checked the corresponding values one by one from the top of the temperature profiles to find the change points. Figure 3 shows algorithm to solve this problem in this paper. Hereafter, we simply assumed the position of the snowice interface was constant. The interfaces between air and snow or ice, and between ice and water were identified using the temperature profiles, the daily amplitude, which are equal to the differences between the maximum and minimum values of the daily temperature profiles of each sensor, and the vertical temperature gradients. The temperature profile, daily amplitude of temperature profiles and the vertical gradients of the temperature were widely used to estimate the interface between air and snow (or ice) [16,21,22]. The vertical gradient in air was different from in snow and ice. The daily amplitude of temperature profile was larger in air than in snow. The temperature profiles were measured by the buoys, and the daily amplitude of temperature profiles and the vertical temperature gradients were calculated from the temperature profile.

Interface between Air and Snow or Sea Ice
The evolution of snow depth over sea ice may be affected by synoptic processes phenomena such as storms, snowfall, and sleet, or due to melting caused by solar radiation. Thus, the temporal fluctuations of snow depth are much larger than for sea ice thickness. The daily amplitude of temperature profiles and the vertical gradients of temperature were examined because, in some cases, it was difficult to accurately distinguish the top interface by using the measurements of the temperature alone. From the data of the daily amplitude of temperature profiles and the vertical gradients of temperature, we checked the corresponding values one by one from the top of the temperature profiles to find the change points.    In winter or sea ice growth period, sea ice was covered by deep snow, so the air-snow interface was considered as the top interface. As shown in Figure 3a, the obvious inflection point of the temperature profiles occurred at the snow-ice interface, with a near constant temperature gradient above this interface. During this period, the internal temperature profile through the sea ice was basically linear (Figure 3b). In the warming period, the snow cover may sublimate or melt. When the ice cover reached a thermodynamic balance, the ice thickness increased to the maximum value for the entire ice season (Figure 3c).
The temperature profile of sea ice at the end of the warming period was C-shaped. Correspondingly, the temperatures measured by the sensors exposed to the air were above 0 °C. The intersection of the top temperature profile with 0 °C was identified as air-ice interface (Figure 3d). The daily amplitude of temperature profile of sea ice can be described as follows: where D(i) is the daily amplitude of a certain point on the temperature profile; Tmax(i) and Tmin(i) are the maximum temperature and minimum temperature at a certain time of a day, respectively. The measurements at the top of temperature profiles had significant daily amplitudes of about 4.5 °C, as shown in Figure 4. However, the bottom of temperature profiles had minor daily changes in our In winter or sea ice growth period, sea ice was covered by deep snow, so the air-snow interface was considered as the top interface. As shown in Figure 3a, the obvious inflection point of the temperature profiles occurred at the snow-ice interface, with a near constant temperature gradient above this interface. During this period, the internal temperature profile through the sea ice was basically linear (Figure 3b). In the warming period, the snow cover may sublimate or melt. When the ice cover reached a thermodynamic balance, the ice thickness increased to the maximum value for the entire ice season ( Figure 3c).
The temperature profile of sea ice at the end of the warming period was C-shaped. Correspondingly, the temperatures measured by the sensors exposed to the air were above 0 • C. The intersection of the top temperature profile with 0 • C was identified as air-ice interface (Figure 3d). The daily amplitude of temperature profile of sea ice can be described as follows: where D(i) is the daily amplitude of a certain point on the temperature profile; T max (i) and T min (i) are the maximum temperature and minimum temperature at a certain time of a day, respectively. The measurements at the top of temperature profiles had significant daily amplitudes of about 4.5 • C, as shown in Figure 4. However, the bottom of temperature profiles had minor daily changes in our study. By examining the sharp change point of daily temperature profiles, the top interfaces can easily be determined. The temperature profile of snow remained linear in winter, but its vertical gradients of temperature differ from the sea ice because the different coefficients of heat conduction among air, snow, and ice [16,21,22]. study. By examining the sharp change point of daily temperature profiles, the top interfaces can easily be determined. The temperature profile of snow remained linear in winter, but its vertical gradients of temperature differ from the sea ice because the different coefficients of heat conduction among air, snow, and ice [16,21,22]. The temperature gradient of temperature profile can be described as follows [16,21,22]: where G is the temperature gradient of the temperature profile; T is the temperature measured by buoys; z is the vertical coordinate. As shown in Figure 5, the interface between air and snow can be estimated from the vertical gradients of temperature by seeking the sharp change point of temperature gradients.

Interface between Ice and Ocean
To identify the lower ice interface (ice-ocean interface), temperature profiles for the lower ice layer were obtained from some thermistors near the bottom of the sea ice. The seawater temperature was determined using the lower five thermistors, which generally had a negligible temperature gradient from the bottom of the sea ice. The points where the temperature profile of the lower ice layer intersected the ocean temperature were regarded as the ice-ocean interface (Figure 3b). The iceocean interface determined by the method of seeking described above had a good accuracy in winter The temperature gradient of temperature profile can be described as follows [16,21,22]: where G is the temperature gradient of the temperature profile; T is the temperature measured by buoys; z is the vertical coordinate. As shown in Figure 5, the interface between air and snow can be estimated from the vertical gradients of temperature by seeking the sharp change point of temperature gradients.
study. By examining the sharp change point of daily temperature profiles, the top interfaces can easily be determined. The temperature profile of snow remained linear in winter, but its vertical gradients of temperature differ from the sea ice because the different coefficients of heat conduction among air, snow, and ice [16,21,22]. The temperature gradient of temperature profile can be described as follows [16,21,22]: where G is the temperature gradient of the temperature profile; T is the temperature measured by buoys; z is the vertical coordinate. As shown in Figure 5, the interface between air and snow can be estimated from the vertical gradients of temperature by seeking the sharp change point of temperature gradients.

Interface between Ice and Ocean
To identify the lower ice interface (ice-ocean interface), temperature profiles for the lower ice layer were obtained from some thermistors near the bottom of the sea ice. The seawater temperature was determined using the lower five thermistors, which generally had a negligible temperature gradient from the bottom of the sea ice. The points where the temperature profile of the lower ice layer intersected the ocean temperature were regarded as the ice-ocean interface ( Figure 3b). The iceocean interface determined by the method of seeking described above had a good accuracy in winter

Interface between Ice and Ocean
To identify the lower ice interface (ice-ocean interface), temperature profiles for the lower ice layer were obtained from some thermistors near the bottom of the sea ice. The seawater temperature was determined using the lower five thermistors, which generally had a negligible temperature gradient from the bottom of the sea ice. The points where the temperature profile of the lower ice layer intersected the ocean temperature were regarded as the ice-ocean interface ( Figure 3b). The ice-ocean interface determined by the method of seeking described above had a good accuracy in winter or sea ice growth period. This method became unreliable in summer, especially in ice melting period when the temperature gradient across the lower interface weakened. In summer, the temperature profile of sea ice became nonlinear with a C-shaped curve. Then the lower ice interface was determined from the obvious inflection point in the vertical C-shaped ice temperature profile (Figure 3d). In winter, temperature profile of sea ice remained linear and temperature of the basal ice layer was colder than the upper ocean. There will be a sharp inflection point occurring at the interface. Thus, ice-ocean interface can be estimated from the vertical gradients of sea ice temperature profile.

Change Point Detection Method
In Section 3.1, we introduced the principle of distinguishing the top and lower interfaces based on the sea ice temperature profile, the daily amplitude of temperature profiles, and the vertical temperature gradients. The calculation accuracy of sea ice thickness and snow depth depends on the accuracy and consistency of interface discrimination.
An off-line change point detection method was used to obtain the optimal estimation of inflection points derived from data of ice-tethered buoys, which were used to identify the top and lower ice interfaces. We assumed that there were N points in the sea ice temperature profiles, so that the sea ice temperature profile is described as follows: where T 1 is the temperature of air or snow; T 2 is the temperature of sea ice; T 3 is the temperature of ocean. In this study, the vertical gradients of temperature and daily amplitude of temperature profiles were also used to distinguish snow depth and ice thickness. For these estimations, T 1 , T 2 , T 3 represented the vertical gradients of temperature or daily amplitude of temperature profiles, respectively. S k was the temperature value from the temperature profile. t 1 was the location of the ice-snow interface, t 2 was the location of the ice-ocean interface, and e k was the noise in the data of sea ice temperature profile. The detection of multiple change points can be accomplished with the single change point detection method, the intent of which is to estimate of the change points.
The principle of the least-squares change point estimation is to take a change point to minimize the sum of the squares of the segment errors between the two subgroups where the change point is the demarcation point. Thus, this method can be applied to the discrimination of the top and lower interfaces. The target formula can be described as follows: where I top is the estimate of the top interface and I lower is the estimate of the lower interface. C stands for a point in the sequence number of data cluster. E 1 (m), E 2 (m), E 3 (m) and E 4 (m) are both segmented mean estimations of the data cluster: We introduced the least-absolute deviations in the estimation of the change point. The principle of this estimation method is to define a change point to minimize the absolute value of the segment error, for which the least-absolute deviations is more robust to outliers than the least-squares method on the handling of errors. The least-squares change point estimation brings in larger errors, which are often the outliers in the measurements and may lead to false results of change point estimation. The least-absolute deviations method can effectively reduce the effect. The least-absolute deviations estimation of change point can be described as follows: In the change point detection method, we propose a selection algorithm for change point detection to improve the accuracy of inflection point detection. The selection algorithm of two change points is described as follows: • Input: The sea ice temperature profile or the daily amplitude of temperature profiles or the vertical temperature gradients S k (k = 1, 2, 3, . . . , N).
The error of change point estimation is E 1 .

•
Step 3: For the profile S k (k = C 0 + 1, C 0 + 2, C 0 + 3, . . . , N), we use the Equation (11)  Step 4: If E 1 < E 2 , C 11 and C 12 are two change points representing the top and lower interfaces. Otherwise, C 21 and C 22 are two change points. • Output: Two change points and calculated ice thickness and snow depth.

Maximum Likelihood Detection Method
Here, we used an off-line mean change point detection method for a normal distribution based on the maximum likelihood method for the optimal estimation of the change point position. First, we discuss the situation of single change points, and then expand to the detection method of multiple points. S k was the temperature value from the temperature profile and t was the only change point of the profile: The noise of this profile can be de described as follows: If the noise is Gaussian white noise, the likelihood function is described as follows: Sensors 2018, 18, 4162

of 22
The parameters E 1 , E 2 , σ 1 and σ 2 are both uncertain. To acquire the likelihood maximum, we calculated the derivative of E 1 and E 2 , respectively and made the derivative equals zero. Then we get a set of equations: We plugged the results of Equations (15)- (18) to simplify the Equation (14), and we can get a simplified equation as follows: The maximum of the Equation (19) is equivalent to the minimum value of the following log-likelihood, which is the maximum likelihood estimation of the change point of the normal distribution: After applying a logarithm to Equation (20), a new equation was obtained: The t that makes the minimum value of the upper form is the maximum likelihood estimation of the position of the change point. However, in the calculation of sea ice thickness and snow depth, two change points are used. Here, S k are the temperature values from the temperature profile and have two change points of t 1 and t 2 . A set of 4D vectors T (0, t 1 , t 2 , N) is introduced. The criterion of change point estimation by maximum likelihood method is that seeking out a group of T (0, t 1 , t 2 , N) to make the next equation have the minimum value of solution:

The Results Calculated Using Ice Temperature Profile
We combined the calculated top and lower interfaces to obtain the ice thickness and compared with "actual" ice thickness. The "actual" ice thickness of IMBs and TUT was measured by acoustic sounders and the "actual" ice thickness of SIMBAs was measured by drilled holes. The maximum of average error of sea ice thickness was 0.0948 m in IMB-2009D, whereas the minimum of average error of sea ice thickness was 0.0205 m in SIMBA-2013A. The estimation error of the entire samples of sea ice buoys is distributed in a relatively narrow range. We set the estimation error within 0-0.06 m as a small error, denoted as SE, which were 92% of the total samples. In detail, the errors in SE were mainly within 0.02-0.03 m (54%). There were 19 groups with errors ranging from 0.03-0.06 m, and there were five groups with sea ice thickness errors of 0.03-0.04 m. The number of buoys with errors of 0.04-0.05 m and 0.05-0.06 m were seven groups, respectively. The errors of sea ice thickness within 0.07-0.09 m were described as a medium error affected 6% of the 50 buoys, denoted as ME. The sea ice thickness errors of more than 0.09 m can be considered as a large error and was 6%, denoted as LE.

The Results Calculated Using Daily Amplitude of Temperature Profile
The maximum error of snow depth using daily amplitude of temperature profiles was 0.093 m obtained for the IMB-2011C, and the minimum error of snow depth was 0.019 m obtained from two buoys of IMB-2003C in Arctic and SIMBA-2013A in the Antarctica. The range of estimation error of snow depth is small, which there were 46 buoys in SE, three buoys in ME, and one buoy in LE.

The Results Calculated Using Vertical Gradients of Temperature
The maximum error of sea ice thickness was 0.094 m in IMB-2009D. The minimum error of sea ice thickness was 0.014 m in IMB-2012J. There were two of the total buoys in ME, and one buoy in LE. The maximum error of snow depth was 0.097 m in IMB-2011C. The minimum error of snow depth was 0.016 m in SIMBA-2013A.

Statistical Errors of the Two Methods
In general, all the calculated average errors compared to the observed values of ice thickness obtained by the methods of change point and maximum likelihood are within a range of 0-0.11 m. After we used two methods to process the data of temperature profile, daily amplitude of temperature profile and vertical gradients of temperature profile, the maximum average errors of ice thickness occurred in IMB-2009D. The minimum average errors of ice thickness occurred in SIMBA-2013A (which appeared five times) and IMB-2012J (it appeared once). Correspondingly, all the calculated average errors of snow depth were 0.016-0.105 m. The maximum average errors of snow depth occurred in IMB-2011C. The minimum average errors of snow depth occurred in SIMBA-2013A (it appeared five times) and IMB-2003C (it appeared two times).
In Figure 8 we show the statistical results of error of snow depth (top interface) and ice thickness (lower interface) using the methods of change point detection method and maximum likelihood detection method for the data measured by 50 ice-tethered buoys.    In Figure 8a, the histogram has peaks of error of 0.02-0.03 m and 0.03-0.04 m. Change point theory had better performance with 14% more groups appearing in the smaller error range (<0.03 m) compared to the maximum likelihood detection method. However, the distribution trends of error in the top interface presented by the two methods are basically the same.

Factors Influencing on the Estimation Error
The results of calculating ice thickness using the temperature profiles by two methods are shown in Figure 8b. The difference only occurred in the intervals of 0.02-0.03 m and 0.03-0.04 m. Other sets of results remained highly consistent.
In Figure 8c, maximum likelihood detection method can reduce estimated errors of snow depth using daily amplitude of temperature compared to the change point method. There were 10% errors of 50 buoys in SE using maximum likelihood detection method more than the other method. As results of change point theory, estimated errors of snow depth was larger than 0.02 m.
The calculation of top interfaces using maximum likelihood method was closer to the observed values compared to that obtained using the change point method, with average errors of 0.038 m to 0.043 m, respectively, as shown in Figure 8d.
The results of calculating ice thickness using the vertical gradients of temperature by two methods are shown in Figure 8e. Two methods had comparable average error of ice thickness (0.038 m and 0.037 m). Errors calculated by the maximum likelihood detection method were 43% appearing in 0.02-0.03 m, while the probability was 28% for the change point detection method.
The maximum and minimum errors of ice thickness and snow depth calculated by the two methods based on three parameters are shown in Table 2. It indicates that the performance of the discrimination of ice thickness by the maximum likelihood theory based on temperature profile and the discrimination of snow depth by the change point theory using the same parameter are better than other combinations.

Seasonal Dependence of Estimation Error
As shown in Figure 9, IMB-2011C was an IMB buoy deployed in the Arctic in late April and the maximum (minimum) value of the average deviation of sea ice thickness was 0.038 m (0.018 m) obtained in September (April) 2011.

Seasonal Dependence of Estimation Error
As shown in Figure 9, IMB-2011C was an IMB buoy deployed in the Arctic in late April and the maximum (minimum) value of the average deviation of sea ice thickness was 0.038 m (0.018 m) obtained in September (April) 2011.  We showed 48 buoys with monthly ice thickness error in the Arctic and two buoys in Antarctica. The maximum value of average error was 0.026 m obtained in July and the minimum value of average error was 0.022 m obtained in January in the Arctic (Figure 10a). For SIMBA-2013A, the maximum value of average error was 0.023 m in November and the minimum value of average error was 0.019 m in April (Figure 10b). For SIMBA-2014A, the maximum value of average error was 0.024 m in December and the minimum value of average error was 0.019 m in June (Figure 10c). In summer, sea ice was at the melting period, the temperature of the upper ocean was colder than the ice cover, and a large number of bubbles and pores appeared in the middle layer of sea ice. At this time, sea ice temperatures may be in an isothermal state. And in such state, the calculated sea ice thickness will have larger error. In winter, when sea ice was in its freezing period, the sea ice internal temperature gradient was established very well, and the difference of temperature between ice and upper ocean was more obvious. Thus, the calculated sea ice thickness had a smaller error.

The Influence of the Vertical Interval of Temperature Sensors
When we used the discrimination algorithm to calculate the sea ice thickness, the temperature sensor intervals on the thermistor strings influenced the accuracy of sea ice thickness. The thermistor strings of IMBs have 45 temperature sensors equidistantly distributed with spacing of 0.10 m. The SIMBAs have 240 temperature sensors with 0.02 m interval on a flexible chain. TUT buoy has 150 temperature sensors and its temperature sensors interval is 0.03 m. We selected a TUT and two SIMBA buoys to analyze the effect of vertical spacing of temperature sensors on the calculation of ice thickness ( Figure 11) and resampled the temperature profile data from 0.02 m to 0.10 m. In summer, sea ice was at the melting period, the temperature of the upper ocean was colder than the ice cover, and a large number of bubbles and pores appeared in the middle layer of sea ice. At this time, sea ice temperatures may be in an isothermal state. And in such state, the calculated sea ice thickness will have larger error. In winter, when sea ice was in its freezing period, the sea ice internal temperature gradient was established very well, and the difference of temperature between ice and upper ocean was more obvious. Thus, the calculated sea ice thickness had a smaller error.

The Influence of the Vertical Interval of Temperature Sensors
When we used the discrimination algorithm to calculate the sea ice thickness, the temperature sensor intervals on the thermistor strings influenced the accuracy of sea ice thickness. The thermistor strings of IMBs have 45 temperature sensors equidistantly distributed with spacing of 0.10 m. The SIMBAs have 240 temperature sensors with 0.02 m interval on a flexible chain. TUT buoy has 150 temperature sensors and its temperature sensors interval is 0.03 m. We selected a TUT and two SIMBA buoys to analyze the effect of vertical spacing of temperature sensors on the calculation of ice thickness ( Figure 11) and resampled the temperature profile data from 0.02 m to 0.10 m. As shown in Figure 11, the estimation error of sea ice thickness increases as the vertical interval of ice temperature increases and the maximum error is at interval of 0.10 m. This is because the larger vertical interval will bring in larger uncertainties for the measurements of sea ice temperature and for the internal temperature gradient of sea ice.

The Influence of Initial Ice Thickness
Typically, ice floes in the Arctic and Antarctic with thicker initial ice thickness would have smaller growth rate and earlier onset of ice melt. Figure 12 shows the relationship between the initial sea ice thickness and the average sea ice thickness error, suggesting that the buoys with the smaller initial ice thickness would have the smaller estimation error of sea ice thickness, which is likely to attributed to the fact that the vertical gradient of temperature at the basal layer is hard to be established and still ambiguous even in the early freezing season for thicker ice thickness. The correlation between average error of sea ice thickness and initial ice thickness is clear. In the future, as the thinning of Arctic sea ice cover, our algorithm will become more reliable. As shown in Figure 11, the estimation error of sea ice thickness increases as the vertical interval of ice temperature increases and the maximum error is at interval of 0.10 m. This is because the larger vertical interval will bring in larger uncertainties for the measurements of sea ice temperature and for the internal temperature gradient of sea ice.

The Influence of Initial Ice Thickness
Typically, ice floes in the Arctic and Antarctic with thicker initial ice thickness would have smaller growth rate and earlier onset of ice melt. Figure 12 shows the relationship between the initial sea ice thickness and the average sea ice thickness error, suggesting that the buoys with the smaller initial ice thickness would have the smaller estimation error of sea ice thickness, which is likely to attributed to the fact that the vertical gradient of temperature at the basal layer is hard to be established and still ambiguous even in the early freezing season for thicker ice thickness. As shown in Figure 11, the estimation error of sea ice thickness increases as the vertical interval of ice temperature increases and the maximum error is at interval of 0.10 m. This is because the larger vertical interval will bring in larger uncertainties for the measurements of sea ice temperature and for the internal temperature gradient of sea ice.

The Influence of Initial Ice Thickness
Typically, ice floes in the Arctic and Antarctic with thicker initial ice thickness would have smaller growth rate and earlier onset of ice melt. Figure 12 shows the relationship between the initial sea ice thickness and the average sea ice thickness error, suggesting that the buoys with the smaller initial ice thickness would have the smaller estimation error of sea ice thickness, which is likely to attributed to the fact that the vertical gradient of temperature at the basal layer is hard to be established and still ambiguous even in the early freezing season for thicker ice thickness. The correlation between average error of sea ice thickness and initial ice thickness is clear. In the future, as the thinning of Arctic sea ice cover, our algorithm will become more reliable. The correlation between average error of sea ice thickness and initial ice thickness is clear. In the future, as the thinning of Arctic sea ice cover, our algorithm will become more reliable.

Comparison between the Results Observed by the Buoys
Data of temperature profiles measured on two sets of SIMBA deployed in Antarctica near Zhongshan Station in 2013, 2014 were used to discern ice thickness with algorithms. Because of regional variations in sea ice physical properties or variations in the external environment, such as air temperature, snowfall, and salinity, the applicability of the discrimination algorithm needs to be reassessed. In the same ice season of freezing period, for example, vertical temperature gradients of sea ice were different between Arctic and Antarctica. In Antarctica, the sea ice where SIMBA-2013A and SIMBA-2014A were deployed was first year ice with very thin snow depth (<0.10 m). Thus, the internal temperature of sea ice may be easily affected by air temperature changes, resulting in larger estimation errors. More measurements from temperature sensors or interpolated temperature were fitted to linear profile to calculate the intersection to improve the accuracy of calculating the ice thickness. As shown in Figure 13, we analyzed the monthly average sea ice thickness error of the buoys in the Arctic and Antarctica. This phenomenon may be attributed to the difference in vertical temperature gradients of sea ice.

Comparison between the Results Observed by the Buoys
Data of temperature profiles measured on two sets of SIMBA deployed in Antarctica near Zhongshan Station in 2013, 2014 were used to discern ice thickness with algorithms. Because of regional variations in sea ice physical properties or variations in the external environment, such as air temperature, snowfall, and salinity, the applicability of the discrimination algorithm needs to be reassessed. In the same ice season of freezing period, for example, vertical temperature gradients of sea ice were different between Arctic and Antarctica. In Antarctica, the sea ice where SIMBA-2013A and SIMBA-2014A were deployed was first year ice with very thin snow depth (<0.10 m). Thus, the internal temperature of sea ice may be easily affected by air temperature changes, resulting in larger estimation errors. More measurements from temperature sensors or interpolated temperature were fitted to linear profile to calculate the intersection to improve the accuracy of calculating the ice thickness. As shown in Figures 13 we analyzed the monthly average sea ice thickness error of the buoys in the Arctic and Antarctica. This phenomenon may be attributed to the difference in vertical temperature gradients of sea ice.

Procedure to Determine Snow Depth and Ice Thickness
To efficiently use the snow depth and ice thickness discrimination algorithm, a procedure of discrimination ( Figure 14) was built. All sea ice temperature data must first go through the pretreatment process. We found that there were a lot of incomplete, inconsistent, and abnormal values in the original temperature profile data from the ice-tethered buoys, seriously affecting the efficiency of the algorithm. During data mining, the original data preprocessing process was also viewed as data cleaning. In general, data cleaning is mainly to delete irrelevant data or duplicate data, and to smooth the noise data in the original data set and filter out irrelevant data. In various thermistor strings of ice-tethered buoys, if the temperature-sensitive element is a thermal resistance type, such as platinum resistances or thermocouples, data loss in the temperature profile is more likely to occur than with semiconductor temperature sensors. We used interpolation to handle missing values by creating an interpolation function using a known point near the missing value to replace an unknown value; if the temperature value represents the temperature of the upper ocean, the fluctuations in the upper seawater temperatures are basically at a very small level, we used the average value of the neighboring values near the missing value as a substitute for this point.
The procedure to determine snow depth and ice thickness using the morphological temperature profile (linear or C-type) was according to the discussed determination modes. Information of ice season (growth period, warming period and melting period) were also considered. The top and lower

Procedure to Determine Snow Depth and Ice Thickness
To efficiently use the snow depth and ice thickness discrimination algorithm, a procedure of discrimination ( Figure 14) was built. All sea ice temperature data must first go through the pretreatment process. We found that there were a lot of incomplete, inconsistent, and abnormal values in the original temperature profile data from the ice-tethered buoys, seriously affecting the efficiency of the algorithm. During data mining, the original data preprocessing process was also viewed as data cleaning. In general, data cleaning is mainly to delete irrelevant data or duplicate data, and to smooth the noise data in the original data set and filter out irrelevant data. In various thermistor strings of ice-tethered buoys, if the temperature-sensitive element is a thermal resistance type, such as platinum resistances or thermocouples, data loss in the temperature profile is more likely to occur than with semiconductor temperature sensors. We used interpolation to handle missing values by creating an interpolation function using a known point near the missing value to replace an unknown value; if the temperature value represents the temperature of the upper ocean, the fluctuations in the upper seawater temperatures are basically at a very small level, we used the average value of the neighboring values near the missing value as a substitute for this point. interfaces were represented by change points to calculate snow depth and ice thickness. The discrimination of ice thickness was selected by maximum likelihood theory using the temperature profile. For the calculation of snow depth, change point theory was used to process temperature profile of sea ice. Figure 14. Schematic of the procedure for snow depth and ice thickness calculation during any period of ice season (ice growth, ice warming and ice melting) and deployment position (Arctic or Antarctica) by ice-tethered buoys. The bold numbers indicate major step in this approach, as noted in the paper.
For any given temperature profiles of sea ice measured by ice-tethered buoys, the following discrimination algorithm procedure can be followed: 1. Data preprocessing. Firstly, measured results of temperature profile of sea ice are processed by threshold filtering and removing the outliers, such as 30 °C or −60 °C. In the second place, smoothing filter is used to filter out irrelevant data. In the third place, we use interpolation to handle missing values by creating an interpolation function to replace the null values of the temperature profile of sea ice.
2. Based on the temperature profiles, the vertical gradients of temperature profiles are calculated. Then the period of sea ice is determined.
(1) Calculate the intersection of the temperature profile and 0 °C. If the intersection is lower than the initial snow-ice interface, the period of sea ice is the melting period. (2) If the intersection is higher than the initial snow-ice interface, the temperature profiles and the vertical gradients of temperature profiles are used to identify the period of sea ice. If the vertical temperature profile is linear, it is the ice growth period. If the temperature profile is C-type, it is the warming period. And in the growth and warming period, the vertical gradients of temperature profiles are different, thus the period can be identified as shown in Figure 15.
3. Use the change point theory and maximum likelihood theory to generate the output of the change points.
4. Use the change points to obtain the top and lower interfaces as described in Section 3.1. The snow depths and ice thicknesses are then obtained. Figure 14. Schematic of the procedure for snow depth and ice thickness calculation during any period of ice season (ice growth, ice warming and ice melting) and deployment position (Arctic or Antarctica) by ice-tethered buoys. The bold numbers indicate major step in this approach, as noted in the paper.
The procedure to determine snow depth and ice thickness using the morphological temperature profile (linear or C-type) was according to the discussed determination modes. Information of ice season (growth period, warming period and melting period) were also considered. The top and lower interfaces were represented by change points to calculate snow depth and ice thickness. The discrimination of ice thickness was selected by maximum likelihood theory using the temperature profile. For the calculation of snow depth, change point theory was used to process temperature profile of sea ice.
For any given temperature profiles of sea ice measured by ice-tethered buoys, the following discrimination algorithm procedure can be followed: 1.
Data preprocessing. Firstly, measured results of temperature profile of sea ice are processed by threshold filtering and removing the outliers, such as 30 • C or −60 • C. In the second place, smoothing filter is used to filter out irrelevant data. In the third place, we use interpolation to handle missing values by creating an interpolation function to replace the null values of the temperature profile of sea ice.

2.
Based on the temperature profiles, the vertical gradients of temperature profiles are calculated. Then the period of sea ice is determined.
(1) Calculate the intersection of the temperature profile and 0 • C. If the intersection is lower than the initial snow-ice interface, the period of sea ice is the melting period. (2) If the intersection is higher than the initial snow-ice interface, the temperature profiles and the vertical gradients of temperature profiles are used to identify the period of sea ice. If the vertical temperature profile is linear, it is the ice growth period. If the temperature profile is C-type, it is the warming period. And in the growth and warming period, the vertical gradients of temperature profiles are different, thus the period can be identified as shown in Figure 15.

3.
Use the change point theory and maximum likelihood theory to generate the output of the change points.

4.
Use the change points to obtain the top and lower interfaces as described in Section 3.1. The snow depths and ice thicknesses are then obtained.

Conclusions
This research demonstrated an approach for calculating snow depth and ice thickness in the Arctic Ocean and Antarctic from temperature profiles, daily amplitude and vertical gradients of temperature profile observed by the ice-tethered buoys. We transformed the problem of calculation of snow depth and ice thickness into a change point detection problem, and evaluated the applicability of theories of change point and maximum likelihood in the determination of snow depth and ice thickness. A procedure was developed to use data of temperature profiles of sea ice to calculate the top and lower interfaces of the layer of snow-covered ice applicable to various ice seasons, temperature sensors interval, initial ice thickness, deployment position. The results indicated that calculated snow depth and ice thickness maintain acceptable range of error. Estimation ice thickness errors in summer were significantly larger than in winter. The higher resolution of the temperature measurement can effectively reduce the error in actual calculation of ice thickness. By analyzing the effect of initial ice thickness, the algorithm and procedure are more adaptable to rapidly changing Arctic sea ice because the smaller the initial ice thickness could be easier to build sea ice temperature gradient thus have smaller estimation errors for ice thickness and snow depth. The practicality of this procedure was examined. Results also indicated that the application of discrimination algorithm and procedure can effectively improves the scientific value of the data of temperature profile through the snow-covered ice even without the joint deployment of acoustic sensors to measure the ice thickness and snow depth.
The ice-tethered buoys used in this study are IMBs, SIMBAs and one TUT buoy. The SIMBAs have a heating mode that effectively distinguishes the snow-ice interface. Therefore, the combination of the data of heating mode of SIMBAs and our discrimination algorithm together could be further used to distinguish the snow-ice interface in the next step. In the current study, we don't detail with the interface between snow and ice because the data measured by the SIMBAs is still limited. In addition, we plan to optimize observation instruments set up in the buoy. For instance, the TUT can

Conclusions
This research demonstrated an approach for calculating snow depth and ice thickness in the Arctic Ocean and Antarctic from temperature profiles, daily amplitude and vertical gradients of temperature profile observed by the ice-tethered buoys. We transformed the problem of calculation of snow depth and ice thickness into a change point detection problem, and evaluated the applicability of theories of change point and maximum likelihood in the determination of snow depth and ice thickness. A procedure was developed to use data of temperature profiles of sea ice to calculate the top and lower interfaces of the layer of snow-covered ice applicable to various ice seasons, temperature sensors interval, initial ice thickness, deployment position. The results indicated that calculated snow depth and ice thickness maintain acceptable range of error. Estimation ice thickness errors in summer were significantly larger than in winter. The higher resolution of the temperature measurement can effectively reduce the error in actual calculation of ice thickness. By analyzing the effect of initial ice thickness, the algorithm and procedure are more adaptable to rapidly changing Arctic sea ice because the smaller the initial ice thickness could be easier to build sea ice temperature gradient thus have smaller estimation errors for ice thickness and snow depth. The practicality of this procedure was examined. Results also indicated that the application of discrimination algorithm and procedure can effectively improves the scientific value of the data of temperature profile through the snow-covered ice even without the joint deployment of acoustic sensors to measure the ice thickness and snow depth.
The ice-tethered buoys used in this study are IMBs, SIMBAs and one TUT buoy. The SIMBAs have a heating mode that effectively distinguishes the snow-ice interface. Therefore, the combination of the data of heating mode of SIMBAs and our discrimination algorithm together could be further used to distinguish the snow-ice interface in the next step. In the current study, we don't detail with the interface between snow and ice because the data measured by the SIMBAs is still limited. In addition, we plan to optimize observation instruments set up in the buoy. For instance, the TUT can be assembled with a thermistor string of 10 m or longer, which a sensor spacing is 0.02 m or smaller. The accuracy of the thermistor string can be increased to reach 0.01 • C by replacing the semiconductor temperature sensor with a more accurate platinum resistor to improve the accuracy of snow depth and ice thickness estimation. The sensor of conductivity and temperature (CT) can be assembled to measure the salinity and temperature of seawater just under the ice cover. In this study, 48 buoys deployed in the Arctic Ocean and two SIMBA buoys deployed in Antarctica were used. Therefore, in future work, we will strengthen the analysis of the data of the ice-tethered buoys deployed in Antarctica, especially for the ice zones further offshore, to enhance the applicability of our algorithm. Different ice types should be considered and analyzed when using our algorithm. For instance, estimated ice thickness should be combined with ice types such as thick or thin multiyear ice, hummock, ponded ice, deformed ice, and ice ridges to enhance applicability and reduce the uncertainties of the discrimination algorithm under different ice types.