Development of ANN-Based Algorithm to Estimate Wintertime Sea Ice Temperature Profile Over the Arctic Ocean

The thermal structure of the Arctic sea ice is a critical indicator in the atmosphere–sea ice–ocean energy budget and, thus, for understanding Arctic warming and associated climate change. Therefore, understanding this thermal structure and its monitoring should be vital. However, it is challenging to obtain a 3-D view of the thermal structure of the sea ice (such as the temperature profile) through satellite measurements because of the lack of understanding of the nonlinear relationship between sea ice emission and measured radiance at the top of the atmosphere. In this study, a model was developed to estimate the temperature profile within the Arctic sea ice during winter using satellite-borne passive microwave measurements. An artificial neural network (ANN) technique based on deep learning was introduced, and the nonlinear relationship between satellite-measured brightness temperatures and buoy-measured sea ice temperature profiles was learned. The ANN model was mapped and verified using the tenfold cross-validation technique. The developed ANN model was able to restore the sea ice temperatures at all specified levels with correlation coefficients > 0.95, absolute biases < 0.1 K, and root mean square errors < 1.6 K. The retrieved temperature results well represent expected thermal structures, in addition to the snow–sea ice interface temperature similar to that in the published literature. Besides the data for validating climate model simulations, the results also promise applications for improving the sea ice growth model performance by tightly constraining the vertical thermal structure in the sea ice growth model.


I. INTRODUCTION
M ONITORING the thermal state of Arctic sea ice is important for understanding Arctic climate change associated with global warming [1], [2]. In particular, the credibility of climate model results regarding the atmosphere-sea ice-ocean energy budget at the Arctic basin scale should be increased when the observed sea ice thermal features support the model results [3], [4]. However, conventional measurements, including buoy or other in situ measurements, are far from available for the diagnosis of large-scale and long-term thermal structures of sea ice, owing to sparse data coverage.
Polar-orbiting satellite measurements often serve as an alternative means to monitor meteorological parameters over Arctic sea ice regions. In particular, spatial and temporal coverage is wide and frequent enough to continuously monitor sea ice parameters over the entire Arctic Ocean [5]. For this purpose, microwave measurements have been effectively used [6]. However, so far, microwave measurements have been applied mainly to monitor the geographical coverage of sea ice, such as concentration, area, type, and edge, based on different features of radiative characteristics between ocean and sea ice and between first-year (FY) and multiyear (MY) sea ice [5], [6].
More recently, however, there has been some success in inferring the thermal conditions of sea ice using satellite-borne microwave measurements. Tonboe et al. [7] and Kilic et al. [8] noted a close relationship between Advanced Microwave Scanning Radiometer-E/2 (AMSR-E/2) vertically polarized (V-pol) 6.9-and 10.7-GHz brightness temperatures (TBs) and temperature at the snow-sea ice interface, suggesting that the snow-sea ice interface temperature could be successfully determined from 6.9-or 10.7-GHz TB measurements. Lee and Sohn [9] demonstrated that the snow-sea ice interface temperature could be estimated from 6.9-GHz TB measurements based on the application of the so-called "combined Fresnel equations" [10]. The basic theory behind these successes is that microwaves can penetrate snow and reach the upper layer of sea ice, whose penetration depths depend on microwave frequencies. Lower frequency tends to reach deeper layers, and at the same time, absorption takes place at varying depths depending on frequencies as well [11], [12]. Thus, microwave radiation emitted from sea ice may embed the thermal information of a certain absorbing sea ice layer. If this is the case, space-borne microwave multichannel measurements of TBs may provide information on the vertical temperature profile within the sea ice column. However, the nonlinear relationship between the emission layer temperature and measured TBs, caused by the optical thickness of snow and the emissivity of sea ice, should be resolved to extract temperature information from satellite measurements.
Machine learning, a branch of artificial intelligence that implements relationships between input and output data, is widely used in various geoscience and remote sensing studies [13], [14]. In particular, advances in computational power and optimization techniques, and improvements in the overfitting problem have led to innovation in deep learning technology that allows interpreting and generalizing the core relationships of high-level abstract data into a network of a specific structure. Therefore, deep learning has the distinct advantage of being more easily able to handle difficult-toprogram applications due to nonparametric, nonlinear, and multivariate problems [15]. Considering that the main problem in retrieving vertical temperature profiles from microwave measurements is to resolve the nonlinear relationship between sea ice temperature and satellite-measured radiance, deep learning may help solve the problem.
In this study, we estimated the winter sea ice temperature profile over the Arctic Ocean from AMSR2 measurements by applying a deep learning-based artificial neural network (ANN) algorithm, which was developed based on the relationship between the buoy-measured temperature profiles and satellite-measured microwave window channel TBs. A detailed account of data used, along with its preprocessing procedure, is provided in Sections II and III, respectively. The physical basis and ANN model development are described in Sections IV and V, while Section VI outlines the factors influencing model performance. Subsequently, Sections VII and VIII present the 3-D structure and representativeness of Arctic sea ice temperature over the past decade (2013-2022), followed by a summary and discussion in Section IX.

II. USED DATA
For the ANN training, this study uses TBs measured by the AMSR2 as inputs and sea ice temperature profiles measured by the Cold Regions Research and Engineering Laboratory (CRREL) and the Alfred Wegener Institute (AWI) Ice Mass Balance (IMB) buoys as outputs. A detailed description of the input and output data is provided in the following.

A. IMB Buoy Sea Ice Temperature Profile
The output training data are constructed with two sources: CRREL and AWI IMB buoys that are deployed over Arctic sea ice twice a year in the spring and fall. The CRREL buoy, whose vertical extent covers the surface air layer to the water below sea ice, measures the near-surface air temperature, snow-sea ice temperature profile with a 10-cm interval, snow depth, and sea ice thickness every 1-4 h [16]. The AWI buoy measures the near-surface air temperature and the snow-sea ice temperature profile with a 2-cm interval every 6 h [17].
It is generally known that deployed drifting buoys experience obstacles such as melting, wind, and deformation, which can cause slant alignment or movement of the reference point indicating a snow-sea ice interface. As time progresses, the chances of encountering such disturbances are more likely higher, causing erroneous buoy observations. To reduce such chances and thus strengthen the data credibility, we only use the first winter data amongst measurements taken by buoys deployed in the fall only. The ten CRREL and 17 AWI IMB buoy measurements showing no substantial data gap during the winter (November-March) were taken during the 2012-2020 period. Fig. 1 shows the trajectories of the total 27 CRREL and AWI IMB buoys. Detailed information on the geographical location, date, snow depth, and sea ice thickness and type at the initial deployment is provided in Table I [18]. The AMSR2 level-1R TB data spatiotemporally collocated with measurements from 27 selected buoys with every 6-h interval, specifically at 00, 06, 12, and 18 UTC, were used as inputs for the ANN training. The spatial difference of 10 km and the temporal difference of 3 h are considered to be criteria for constructing the collocated input training data. A total of 7459 points were available spanning the years 2012-2022 during the November-March period for only the area where the sea ice concentration (SIC) was greater than 95%. In addition, to produce long-term sea ice temperatures over the Arctic Ocean from 2013 to 2022, AMSR2 level-3 daily TB data with a spatial resolution of 10 × 10 km were used as input data. Table II lists the detailed  characteristics of the AMSR2 sensors. ANN training was performed only over the sea ice area, whose SIC was over 95%. We used AMSR2-derived SIC data from the NASA Team 2 algorithm based on the use of 36.5-GHz V-pol and 18.7-and 89-GHz V-/H-pol TBs [19]. The daily SIC data in a 10 × 10 km spatial resolution were available for the AMSR2 observation periods.
III. DATA PREPROCESSING CRREL offers the acoustic/thermal mass balance buoy data in preprocessed form. As a preprocessing step, it determines boundaries of air, snow, sea ice, and ocean in the temperature profile by combining the acoustic sensor measurements and the vertical temperature gradient. Specifically, the acoustic sensors are used to measure the snow-top and sea ice-bottom, while the vertical temperature gradient is used to estimate the snow-sea ice interface. The preprocessed CRREL data include the temperature profile rearranged to have 0 cm at the determined snow-sea ice interface and, accordingly, adjusted snow depth and sea ice thickness.
Different from CRREL preprocessed data, AWI provides the thermal mass balance buoy data in raw form (i.e., without  Table I. CRREL-like preprocessing). In addition, the AWI data do not include acoustic sensor measurements. Thus, to make AWI data consistent with CRREL data, we took the preprocessing procedures using the method proposed by Shi et al. [20]. This method estimates the snow-sea ice interface in the temperature profile by utilizing the monthly mean vertical temperature gradient. The temperature profile is then rearranged to be 0 cm at the determined snow-sea ice interface.
In addition to the estimation of the snow-sea ice interface, it was needed to implement a preprocessing for the GPS geolocations of the AWI buoys. This is due to the fact that a significant number of the buoys displayed anomalous drifting trajectories. Fig. 2(a) provides an example of the drifting trajectories throughout the observation periods for two buoys, i.e., AWI 2014T33 and AWI 2018T46. This figure clearly depicts an issue concerning the geolocation data, i.e., irregular trajectories. Therefore, we introduced additional preprocessing to eliminate these irregularities in the geolocation data. Specifically, the mean and 3-standard deviation for the 6-h drifting distance for each buoy were computed, and if the distance between two adjacent buoy locations within 6 h exceeded a threshold (i.e., 3-standard deviation), it was considered abnormal and removed from the training data. Fig. 2(b) depicts the trajectories of two buoys after the preprocessing, demonstrating the successful removal of data showing abnormal locations. In this study, two preprocessing steps are implemented for AWI buoy data that consist of determining the snow-sea ice interface and eliminating the irregular geolocation.
IV. PHYSICAL BACKGROUND To provide a physical background for the temperature profiling over sea ice from microwave window channel measurements, we constructed a conceptual model for the snow/sea ice (referred to as "model sea ice") from which   Fig. S1 in the Supplementary Material). Over those regions, vertical distributions of temperature, salinity, and density are obtained from the simulation of sea ice for January 1, 2021 [21], and representing profiles for each type of sea ice are obtained by taking area averages. Mean snow depth and sea ice thickness are found to be 26 cm and 2.01 m for the MY sea ice region and 18 cm and 1.39 m for FY sea ice, respectively. Resultant profiles are provided in Fig. 3(a).
WFs for the AMSR2 channels at given frequencies of 6.9, 10.7, 18.7, 23.8, 36.5, and 89.0 GHz were calculated using the Snow-Sea Ice Emission Model (SSIEM [22]) with inputs given in Fig. 3(a). The SSIEM is a model integrating a sea ice radiative transfer module with the Microwave Emission Model of Layered Snowpacks (MEMLS [23]), which enables the radiative transfer calculation in the snow/sea ice layers.
Calculated WFs for MY and FY sea ice are given in Fig. 3(b) and (c), respectively. It is clear that the emission is contributed mostly from progressively deeper layers as the channel frequency decreases for both sea ice types. This is because the imaginary part of the refractive index decreases with decreasing frequency [24], and thus, the penetration becomes deeper with decreasing frequency. It is also apparent that the WFs overlap to some extent, allowing for finite channel data to define the temperature profile but only in the upper part of the sea ice layer. However, compared to WFs for MY sea ice, WFs for FY sea ice are located mostly in the layers shallower than 0.5 m except for the 6.9-GHz channel, whose WF extends relatively deeper layers up to 1 m. The relatively shallower emission layer for FY ice is due to the higher salinity of the FY ice. Snow layers are also capable of contributing to the emission, but the 89.0-GHz channel appears to carry more information on the snow emission, possibly providing the snow information to the ANN model. Hence, the frequency-dependent emission layers and overlapping WFs suggest that thermal information on sea ice can be extracted from satellite-based microwave window channel measurements. However, here, instead of developing a pure theoretical physical retrieval algorithm, we employ the machine learning approach to link the microwave radiances to the temperature profile.
It should be noted that this well-established physical understanding is only applicable during the cold season because the water contamination in the microwave signal during the melting season obscures the emitted temperature signal. Even slight wetness inside a snowpack can significantly alter the emission WF, such as significantly skewing WFs toward a wet layer [25], [26], [27]. Therefore, we focus on the cold season of the November-March period, in which dry snow conditions are prevalent [28].

V. DEVELOPMENT OF THE ANN MODEL
In this study, to estimate the temperature profile within sea ice from satellite-measured microwave multichannel TBs, we adopted an ANN technique based on deep learning and implemented a nonlinear relationship between two parameter datasets (i.e., IMB buoy-measured sea ice temperature profiles and AMSR2-measured TBs). The following subsection describes the procedures of ANN model development and validation.

A. k-Fold Cross-Validation Method
Holdout cross-validation is one of the most widely used techniques for mapping an ANN structure and verifying its performance. However, when the dataset is small, it is susceptible to data leakage and has low reliability. Therefore, because limited buoy data were used for training in this study, we applied a k-fold cross-validation method that improves the shortcomings of the holdout method (see Fig. 4 for the outline of our employed method). The k-fold method randomly divides the training dataset of the input/output into k clusters (or folds) based on random sampling without replacement. Kohavi [29] and Rodriguez et al. [30] recommended the optimal cluster number to be 10 (k = 10) in terms of the bias-variance tradeoff and computational cost. Accordingly, this study sets k to 10 (i.e., tenfold cross-validation) and built the ANN structure with the highest accuracy by testing various model configurations (e.g., the number of hidden layers, number of neurons, activation function, optimizer, inputs, and outputs) during the verification process. As shown in Fig. 4, one of the ten clusters was set aside as a cluster for verification, whereas the remaining nine clusters were used as clusters for learning. This procedure was repeated until each of the ten clusters could be used as a verifying cluster.
Thus, the k-fold method can utilize all the training data while repeating the learning and verification processes. Briefly, the k-fold method iterates the holdout method k times independently and estimates a general prediction error.

B. Artificial Neural Network
We used the ANN approach, which enables the formation of complex networks by combining many neurons and layers, and then handles nonlinear relationships by optimizing the combination strength (weight) between networks through learning. A neural network consisting of two or more hidden layers is generally referred to as a deep neural network and has the advantage of generalizing more complex data efficiently with fewer nodes, as against a single neural network [31], [32], [33].
The ANN was set up using TensorFlow [34] and Keras [35]. The input training data were the AMSR2-measured TBs (in K) at six V-pol channels. Only V-pol TBs are used because the sea ice emissivity of the V-pol component is approximately 10%-30% higher than that of the H-pol component [12], [36], thereby facilitating the extraction and analysis of sea ice information from measured radiances at the top of the atmosphere (TOA). The output training data comprised 11 levels of sea ice temperatures (in K) from the snow-sea ice interface (=0 m) to the 1-m level in sea ice depth, with a 10-cm interval, which is the designated observation interval for the IMB buoy. Regarding the maximum retrieval depth, we tested various output levels and found that the ANN consistently achieved high accuracy across all levels, provided that the maximum depth was less than 1 m. Based on these results, the maximum sea ice depth of 1 m was selected. This is acceptable because the microwave emission deeper than 1 m of sea ice depth is small, as shown in the distributions of AMSR2 channel WFs in Fig. 3. Here, we assume that the temperature between two adjacent levels is linear. Each input/output training data vector was scaled to a range between 0 and 1 using min-max normalization as a preprocessing step.
At the same time, the number and connections of hidden layers and neurons were also determined by testing various configurations to ensure that the ANN results had the highest accuracy. In this study, a neural network was implemented by deploying two hidden layers that perform calculations between the input and output layers. The first and second hidden layers were composed of 14 and 11 neurons, respectively (see Fig. 5). All neurons were fully connected between the hidden layers. In the hidden layer, each neuron (x) was defined as expressed in (1), with the respective weights (W ) corresponding to the upper layer values (a) along with the bias (b). Subsequently, the bias-corrected weighted sum was delivered to the next layer (h) by passing it through the activation function ( ), as defined in (2) where "m" denotes the index of the hidden layers (m = 1 and 2 as the first and second hidden layers, respectively), "i" is the index of the connected upper layer values (n = 6 and 14 for the first and second hidden layers, respectively), and " j" indicates the index of neurons in the hidden layer ( j = 1-14 and 1-11 for the first and second hidden layers, respectively). A rectified linear unit (ReLU) was adopted as the activation function in the following equation: The cost is defined as the mean square error (mse) between the predicted and actual output values. Learning the ANN model means adjusting the weights and biases to achieve the minimum cost. First, the network's weights and biases are initialized to "1." They are then iteratively updated using the Adam optimizer [37] until the cost meets the prescribed convergence criteria from early stopping.

C. Verification of the ANN Model
As shown in Fig. 4, ten learning cycles were conducted, and each cycle contained a verifying cluster of data. The learned ANN model at each cycle was applied to the AMSR2-measured TBs to estimate sea ice temperatures at 11 different sea ice depths. The verification results, in which the estimated sea ice temperatures from all verifying clusters during the ten learning cycles were compared with buoy measurements, are shown in Fig. 6 as scatter plots along with the associated correlation coefficient (r ), bias, and root mse (RMSE) values. At all 11 layers from 0 to 100 cm, the ANN model was able to produce temperatures with r -values > 0.95, absolute biases < 0.1 K, and RMSEs < 1.6 K, strongly indicating that the vertical distribution of the temperature within the 1-m depth can be well produced from satellite-borne microwave measurements. In addition, the RMSEs are found to decrease as the depth increases, which can be attributed to the relatively larger temperature variation in the upper layer of sea ice compared to the lower layer-supporting material for this can be found in Fig. S2 in the Supplementary Material.
Encouraged by the good validation results, relearning was performed using all training data to finalize the ANN model. This relearning procedure is depicted at the bottom of Fig. 4. Note that the estimated temperature was set to 272 K if the estimated value was higher than 272 K (i.e., a postprocessing step was performed). This is because the maximum sea ice temperature of the buoy observations was approximately −1.15 • C.

VI. FACTORS INFLUENCING THE ANN MODEL
This section investigates the factors likely to influence the ANN model's performance. Three factors are noted to be possibly affecting the model results: overfitting, generalization ability, and sensitivity to snow depth. The relevant countermeasures and tests for these three factors are presented in the following.

A. Overfitting
In this study, to reduce the potential risk concerning overfitting issues, we take the following four measures during the course of developing the ANN model.
1) The model capacity refers to the degree of complexity in the network structure, which increases as the hidden layers are stacked deeper or the number of neurons per layer is increased. The model with high capacity bears a high probability of encountering the overfitting problem by memorizing too many features, even noise, of the learning data [38]. Thus, it is imperative to achieve a balance between model capacity and training data, but there is no general way to determine the appropriate model capacity for the training data [39].
Regarding the model capacity, we evaluated the model performance by gradually reducing the number of hidden layers and neurons. Through this process, we could determine a concise network structure that maintains the model's high performance. In that way, we made an attempt to design an efficient network architecture for the training data to mitigate overfitting. 2) The k-fold method that we adopted here is considered to be one of the preferred cross-validation methods to diagnose overfitting. This is because the k-fold method repeatedly performs model learning and validation k times by randomly shuffling the learning and validation datasets without duplication. Consequently, if no overfitting occurs k times, it can be deemed that the model's robustness and reliability are ensured [40], [41], [42]. 3) Monitoring the relative behavior of learning and validation losses is widely used as a means of diagnosing overfitting in machine learning. If the number of epochs substantially increases, the loss for learning data would decrease, but this may lead to overfitting (i.e., the loss for validation data increases). Thus, there is a need to stop learning at the point when the minimum value of the validation loss maintains for a specific duration to avoid overfitting [43]. In this study, we employed the early stopping technique terminating the learning process at an appropriate point to prevent overfitting. 4) Normalization of the training data to improve the optimization also serves as a means of alleviating the potential risk of overfitting [32].

B. Generalization
To examine the issue of generalization in the developed ANN model, we conducted an additional experiment across learning-verifying-testing. In this experiment, the training database B was partitioned into three distinct random datasets (40% for learning, 10% for validation, and 50% for testing), i.e., B = B learn + B val + B test . The learning dataset (B learn ) was used to learn the ANN model, the validation dataset (B val ) was used to validate the ANN model and monitor the Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. overfitting problem, and the testing dataset (B test ) was used to assess the model's generalization ability. It is worth noting that the point of difference between B val and B test is that B val is not learned but still involved in the training process, whereas B test is a completely independent dataset that does not participate in the training at all. Hence, the testing with B test ensures a more reliable evaluation of the model's performance for unseen data.
The ANN model was successfully learned and validated using B learn and B val (albeit not presented herein). Subsequently, we tested the ANN model using B test and plotted the evaluation results for each layer in Fig. 7. All layers exhibited r -values > 0.95, absolute biases < 0.1 K, and RMSEs < 1.6 K along the one-to-one corresponding line. These testing results were in excellent agreement with those obtained from tenfold cross-validation. Therefore, the proposed model architecture is capable of estimating the sea ice temperatures from unseen data with a level of accuracy comparable to that of its validation performance. This suggests that the ANN model has good generalization ability.

C. Snow Depth
The snow lying over sea ice is capable of absorbing and scattering microwave radiation; thus, it is desirable to include snow information explicitly in the training procedures. However, since this study does not use explicit snow information as input, it is worthwhile to examine how the ANN model performance is sensitive to the snow depth (or snow equivalent water) overlying the sea ice when buoy observations are made. Instead of using direct snow observations, we conducted sensitivity tests comparing the results of two extreme snow depth groups. We identified nine CRREL IMB buoys that had available data on snow depth throughout the training period, except for the 2015J buoy whose snow depth measurements are missing. For each buoy, we calculated the mean snow depth during the training period, providing a representative value. The mean snow depths of all nine buoy observations were arranged in ascending order (see Fig. 8). Two extreme groups were selected; snow depths of 20 and 30 cm were used as criteria for determining one shallow group and one thick group. Buoys showing the three thinnest snowpacks (<20 cm) and buoys with the two thickest snowpacks (>30 cm) were found. Three buoys (2021L, 2014F, and 2013I) showing the three thinnest snow depths were referred to as Group 1, whereas 2012G and 2013F buoys showing the two thickest snow depths were referred to as Group 2. Then, the ANN model was learned using a dataset without Group 1 and verified with Group 1 and vice versa. Fig. 9 shows the scatter plots for the sensitivity test results with the corresponding r , bias, and RMSE values. In the figure, scatter plots are the verification results of Groups 1 and 2 when Groups 1 and 2 are excluded from learning the ANN model. All layers had r -values > 0.99, absolute biases < 0.2 K, and RMSEs < 1.0 K, which were similar to the tenfold cross-validation results (i.e., r -values > 0.95, biases < 0.1 K, and RMSE < 1.6 K). Notably, the model performance was nearly the same for the two extreme snow depth cases, suggesting that the ANN model is capable of resolving the influence of snow depth on TBs. This can be attributed to the fact that TB changes in 6.9-, 18.7-, and 36.5-GHz V-pol channels, typically used in the snow depth retrieval algorithms, can carry the information on the snow depth [44], [45]. Furthermore, the use of 89 GHz likely provided snow information to the ANN model, as shown by distributions of AMSR2 channel WFs in Fig. 4 and [46] for inferring the snow scattering type from 89-GHz measurements.

VII. 3-D VIEW OF SEA ICE TEMPERATURE
OVER THE ARCTIC OCEAN With the developed ANN model, we aimed to produce a 3-D distribution of sea ice temperature up to the 1-m depth Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. from AMSR2 TB measurements over the Arctic Ocean during the winter. In this study, to avoid microwave signals being contaminated by ocean water, the retrieval was limited to an area whose SIC was larger than 95% over the Arctic Ocean. The January-March (JFM) composite was produced instead of the December-March period because the sea ice thickness is generally deeper than 1 m during the JFM period (see Fig. S3 in the Supplementary Material); hence, the ANN model can be applied without additional measures.
The developed ANN model provides sea ice temperatures at 11 levels from 0 cm at the snow-sea ice interface to 100 cm in depth, with a vertical resolution of 10 cm. As a case, ten-year (2013-2022) daily sea ice temperature profiles within 1 m of deep sea ice over the Arctic Ocean were produced by applying the developed ANN model to the AMSR2 daily V-pol TB data. Then, we determined the ten-year JFM mean temperatures to examine whether the obtained results depicted known features related to the Arctic sea ice. To assess the quality of the ANN-estimated sea ice temperature profile, we also provided the Kilic-estimated temperature at the 0-cm depth. In the study by Kilic et al. [8], snow depths were first retrieved from microwave measurements based on the regression equation relating AMSR2 6.9-, 10.7-, and 36.5-GHz V-pol TBs to buoy-measured snow depth. Subsequently, the retrieval equation for the sea ice temperature at 0 cm was obtained by regressing the AMSR2 6.9-GHz V-pol TB (or 10-GHz V-pol TB) to the buoy-measured snow-sea ice interface temperature with regression-based snow depth. Therefore, it can be presumed that snow-sea ice interface temperature in the study by Kilic et al. [8] is from three AMSR2 V-pol TBs (6.9, 10.7, and 36.5 GHz). Thus, the physical basis of this study is similar to that of the study by Kilic et al. [8], despite the use of different approaches (multiple regression versus deep neural network).
The comparison of the JFM mean temperature at 0 cm between the study by Kilic et al. [8] [see Fig. 10(a)] and this study [see Fig. 10(b)] shows good agreement in geographical pattern although the temperature in this study appears to be slightly colder. The two temperatures at 0 cm are in broad agreement with r -value = 0.82, mean deviation = −2.24 K, and root mean square deviation (RMSD) = 2.75 K. This result suggests that the temperature estimated in this study is at least comparable to the regression-based Kilic-estimated temperature, and the deep learning approach successfully resolves the nonlinear characteristics inherent between sea ice temperatures and TOA TBs.
Furthermore, our approach produced temperatures at ten additional levels (from 10 to 100 cm). Regarding accuracy, we claim that accuracies at other sea ice levels should be of the same order of accuracy at 0 cm, as demonstrated against buoy observations and against the result by Kilic et al. [8]. This is because the verification results against buoy observations (see Fig. 6) are nearly identical to the 0-cm results. With this notion in mind, we present the geographical distributions of temperature at different depths in Fig. 10(b)-(f) in the deepening order. The oceanic heat flux to the bottom of sea ice is small and nearly constant during winter [47]; thus, the temperature of sea ice is inevitably dominated by atmospheric variation [21]. From such facts, it is expected that the geographical contrasts shown near the surface level will tend to weaken at deeper levels and become more homogenous. Consistent with this expectation, it is clearly shown that the sea ice temperature gradually increases with depth, and the geographical contrast weakens with depth. Such a 3-D view appears to reflect well-known features influencing the thermal structure of sea ice.

VIII. REPRESENTATIVENESS OF SEA ICE TEMPERATURE
TO THE EURASIAN SIDE OF THE ARCTIC OCEAN In the polar region, the remoteness, harsh environmental conditions, and high operational costs give rise to substantial spatiotemporal constraints, thereby hindering in situ observations. Moreover, since IMB buoys are commonly deployed on sea ice during the autumn months of September and October, there is a limitation that buoys cannot be deployed without sea ice formation during that period. In this regard, the Eurasian side of the Arctic Ocean is the region where buoy deployment is not feasible because most of it is an open ocean in the autumn. Furthermore, the ocean current there experiences a directional flow from the Eurasian continent toward North America during the winter. Consequently, the absence of buoy observation data on the Eurasian side of the Arctic Ocean in the winter is one of the constraints faced by all Arctic sea ice-related studies and this study.
In this section, we took an indirect way to examine the representativeness of the Eurasian side of the Arctic Ocean as to what degree of accuracy might be held in the ANN-estimated sea ice temperature. As previously mentioned, in the Arctic sea ice region during winter, the atmospheric forcing exhibits large and high variations, while the oceanic forcing remains relatively small and nearly constant. Hence, as the sea ice temperature is primarily changed by atmospheric forcing, we expect high consistency in temperature variation between near-surface atmosphere and sea ice-top. For analyzing the relevant claim, we first obtained the near-surface air temperature at 2 m from ERA5 reanalysis, which was completely independent of the ANN model. Then, the ERA5 near-surface air temperature reanalysis was compared with the ANN snow-sea ice interface temperature estimate. It was worth noting that we employed temperature anomalies here instead of absolute temperatures. This is because the temperature anomalies provide a better representation of temperature variation over time rather than absolute temperatures and can also preclude the potential impact of any biases in both the ERA5 and ANN models [48]. Fig. 11(a) illustrates the time series of the three-month mean and its standard deviation for both temperature anomalies in the East Siberian Sea-Laptev Sea-Nansen Basin region [ESLN region in Fig. 11(b)] for ten years from 2013 to 2022. Here, temperature anomalies are quantified from the respective mean temperature for the entire analysis period, and the ESLN region is considered the representative region for the Eurasian side of the Arctic Ocean. The anomaly temperatures of ANN and ERA5 are in good agreement in terms of the general trend. However, the thermal variations of ANN values at the snow-sea ice interface are found to be smaller than ERA5 temperature at the near-surface atmosphere. Such a magnitude difference is highly reasonable when the insulation effect of the snow is considered. Based on these comparison results, we speculate that the ANN model has some degree of representativeness on the Eurasian side of the Arctic Ocean. However, further studies are necessary for this region to provide a more comprehensive understanding.

IX. CONCLUSION AND DISCUSSION
We estimated the vertical distribution of sea ice temperature by developing a deep learning-based ANN model using sea ice temperature profiles from CRREL and AWI IMB buoy observations and collocated AMSR2 microwave window channel TBs (i.e., 6.9-, 10.65-, 18.7-, 23.8-, 36.5-, and 89-GHz V-pol TBs). In the learning process, the model generalized the nonlinear relationship between satellite-measured microwave TBs and buoy-measured sea ice temperature profiles. The background of using passive microwave measurements to retrieve sea ice temperature is given as follows: 1) microwave radiation can penetrate deep into sea ice layers, and thus, emission layers at different microwave frequencies are different from each other and 2) the nonlinear relationship between temperature profiles and emission signals from different sea ice layers can be well resolved using deep learning techniques.
Machine learning algorithms learn not only the relationship between inputs and outputs but also the interrelationship between outputs themselves [49]. This mutually cooperative information between outputs serves to add one more stringent regulation, thereby enhancing the model's reliability and accuracy. Moreover, the temperatures of sea ice share a certain degree of correlation with each other due to being tied together as a single sea ice system. These correlations further strengthen this regulation. Consequently, even if the given TBs provide a relatively weak emission signal at a depth of 1 m in FY sea ice, the ANN model is capable of minimizing substantial errors in the output because the temperature estimate is made within a framework of the interrelationships learned from output training data.
Sea ice temperatures were estimated at 11 levels, from the snow-sea ice interface (0 m) to 1-m depth, with a vertical resolution of 10 cm during winter. The comparison of ANN-estimated sea ice temperature profiles with in situ measurements of drifting buoy showed good correspondence, with r -values > 0.95, absolute biases < 0.1 K, and RMSEs < 1.6 K for all 11 layers. This indicates that the ANN model can retrieve sea ice temperatures using satellite measurements.
Over the Arctic Ocean domain, daily sea ice temperature profiles for the JFM period from 2013 to 2022 were estimated by employing the developed ANN model to AMSR2 daily TBs, and the ten-year mean climatology was then obtained to examine whether the estimated temperature data depicted known features related to Arctic sea ice. The snow-sea ice interface temperature was compared with the Kilic-estimated temperature, and these two temperature values were consistent (r = 0.82, bias = −2.24 K, and RMSD = 2.75 K). From the perspective of the 3-D thermal structure, the ANN-estimated sea ice temperature distributions were consistent with known characteristics associated with the sea ice thermal structure. Over the Eurasian side of the Arctic Ocean domain, the time series of temperature anomalies of both the snow-sea ice interface temperature and ERA5 near-surface air temperature were in good agreement in terms of the general trend and expected variation. Overall, the ANN estimations are found to be coherent with ERA5 reanalysis within a reasonable bound. However, owing to the spatiotemporal constraints of in situ observations in the Arctic Ocean, the deployment of buoys in the training data is comparatively more concentrated in MY ice regions than in FY sea ice regions. Accordingly, we speculate that the uncertainty over the FY sea ice region is likely higher than that over the MY sea ice region.
3-D distributions of sea ice temperature should be valuable for studying the impact of Arctic sea ice on climate change. In particular, during winter, sea ice and overlying snow insulate the relatively warmer ocean from the colder atmosphere.
However, because of the recent thinning of sea ice, this insulating effect has become weaker, thereby causing considerable heat transfer through the sea ice to the atmosphere. Consistent with this notion, Landrum and Holland [50] demonstrated that atmospheric warming could be increased due to the increased wintertime conductive heat fluxes caused by thinning sea ice (and decreasing snow depth). As the conductive heat fluxes within sea ice can be determined by the vertical gradient of sea ice temperature, the heat fluxes from the ocean to the atmosphere can be quantitatively determined from the retrieved temperature profiles. These types of data can be used to examine the reliability of the climate model simulations.
Furthermore, the retrieved temperature data may be used to constrain the initial thermal state in the model simulations needed for the sea ice growth. For example, Kang et al. [21] demonstrated that Arctic sea ice thickness simulation can be considerably improved using the satelliteestimated snow/sea ice temperature data to nudge a sea ice growth model involving 1-D thermodynamic diffusion. In addition, Anheuser et al. [51] used the snow-sea ice interface temperature data from space to estimate the basin-wide thermodynamic sea ice thickness growth. This study offers a 3-D view of the upper sea ice layer temperature over the entire Arctic Ocean, which can further augment our ability to reproduce sea ice thickness across the Arctic by thermodynamically constraining the upper sea ice, although the current approach is applicable only during the winter season. His research interest is satellite remote sensing studying aviation hazard meteorology.
Eui-Jong Kang received the Ph.D. degree in Earth and environmental sciences from Seoul National University, Seoul, Republic of Korea, in 2021.
Since that time, he has been a Post-Doctoral Researcher with the Research Institute of Basic Sciences, Seoul National University. His research interests are the development and application of microwave data assimilation techniques, with a specific focus on the winter Arctic sea ice area. He is currently a Professor with the School of Earth and Environmental Sciences, Seoul National University. He is engaged in aerosol and cloud observations using in situ and remote sensing techniques. His research covers the radiative effects of light-absorbing aerosols and the microphysical properties of Arctic clouds. He was a Post-Doctoral Researcher with the Satellite Meteorology Laboratory, Seoul National University, from 2021 to 2023. He has been working on developing satellite retrieval algorithms for snow and sea ice parameters (e.g., snow depth, sea ice thickness, and emissivity) at the Danish Meteorological Institute, Copenhagen, Denmark, since 2023. His main scientific interest is to synergize multisatellite observations, including passive microwave, infrared, and active altimeter measurements, for improving current retrieval algorithms.