Long-Term Arctic Snow/Ice Interface Temperature from Special Sensor for Microwave Imager Measurements

: The Arctic sea ice region is the most visible area experiencing global warming-induced climate change. However, long-term measurements of climate-related variables have been limited to a small number of variables such as the sea ice concentration, extent, and area. In this study, we attempt to produce a long-term temperature record for the Arctic sea ice region using Special Sensor for Microwave Imager (SSM/I) Fundamental Climate Data Record (FCDR) data. For that, we developed an algorithm to retrieve the wintertime snow/ice interface temperature (SIIT) over the Arctic Ocean by counting the effect of the snow/ice volume scattering and ice surface roughness on the apparent emissivity (the total effect is referred to as the correction factor). A regression equation was devised to predict the correction factor from SSM/I brightness temperatures (TBs) only and then applied to SSM/I 19.4 GHz TB to estimate the SIIT. The obtained temperatures were validated against collocated Cold Regions Research and Engineering Laboratory (CRREL) ice mass balance (IMB) drifting buoy-measured temperatures at zero ice depth. It is shown that the SSM/I retrievals are in good agreement with the drifting buoy measurements, with a correlation coefﬁcient of 0.95, bias of 0.1 K, and root-mean-square error of 1.48 K on a daily time scale. By applying the algorithm to 24-year (1988–2011) SSM/I FCDR data, we were able to produce the winter-time temperature at the sea ice surface for the 24-year period. SIITs are in good agreement with drifting buoy measurements of the SIIT (temperature at the zero ice depth), with a correlation coefﬁcient of 0.95, bias of 0.1 K, and RMS error of 1.48 K on a daily time scale. These results suggest that a long-term ice temperature record can be produced by applying the developed algorithm to well-calibrated long-term FCDR data composed of SSM/I and SSMIS data spanning over thirty years.


Introduction
The increasing emission of man-made greenhouse gases has caused global warming of 0.8 • C since 1900 [1][2][3], but the temperature increase is not the same everywhere over the globe. Especially, over the Arctic polar region, where snow and ice cover the ocean, the increase in the near-surface is nearly twice that in the tropics-a phenomenon called Arctic amplification [4,5]. Such a rapidly rising surface temperature over the Arctic is manifested by changes in the sea ice extent over last three decades. The sea ice extent has been declining at a rate of ≈3.8%/decade [6]. The largest decline is observed in September with an alarming rate of ≈13%/decade, corresponding to more than 30% of sea ice loss compared with the late 1970s [7,8].
The impact of ice loss appears to be profound, not only influencing the local climate, hydrological cycle, ecosystems, and natural habitats, but also affecting the global weather/climate system. With the declining sea ice, more ocean water will be exposed to cold and windy atmospheric environments, probably inducing enhanced evaporation to the cold atmosphere, which can lead to more water vapor, clouds, and more precipitation [9]. The precipitation change in response to the sea ice melt is particularly interesting because the precipitation itself (regardless of snow or rain) has the potential to feedback the climate [10]. As expected, increased precipitation is noted in the Arctic [11], mainly due to the increased local evaporation caused by the sea ice loss [12]. Furthermore, a higher percentage of precipitation throughout the Arctic summer time will be in liquid form, namely rain, due to the warmer low atmosphere, forcing the snow cover to further decline [10]. It is projected that rain will be the dominant form of precipitation in the Arctic by the end of the 21st century [13]. The increased rain, a hydrological response to dramatically decreasing sea ice (or to the warming atmosphere) should lead to a positive feedback to ice melting, through the substantially decreasing surface albedo, establishes a self-feeding cycle where warming temperatures will lead to more rain and ice melt [10] and cause more warming due to even more rain and ice melt.
With respect to the influence of the decline of Arctic sea ice on the weather/climate elsewhere, there is cumulating evidence that Arctic warming has been identified as an important element of the variability of mid-latitude, large-scale atmospheric circulation including storm tracks, planetary waves, jet stream, and energy transport [14,15]. It is reported that the atmospheric variability due to Arctic warming is closely linked to more frequent events of atmospheric blocking patterns [16] and increased cold surges with snowy winters in East Asia, Europe, and Northern America [14,17,18]. Central and northern Europe have experienced six years (2007-2012) of extraordinary wet summer, likely due to the southward shift of the summer jet forced by Arctic sea ice loss [19]. More extreme weather may be expected worldwide in the future due to continuous Arctic warming and the associated sea ice loss.
Although current monitoring of sea ice extent and understanding of subsequent impacts on climate components are in high confidence, dramatic changes of the Arctic climate system and expected continuous change in the future make it difficult to assess the future evolution of Arctic environments and/or predict global climate change. To facilitate our understanding of Arctic climate change and thus to reduce the uncertainty of future climate change assessment, observations of climate variables, such as the sea ice extent and its depth, snow cover, and depth are needed, in addition to observations of thermal structures throughout the atmosphere-snow-ice column. At least, much needed thermodynamic variables may be the surface air temperature, snow top temperature, and snow-ice interface temperature (SIIT), from which the thermal structure of the snow-ice system can be studied and related to the sea ice melt or sea ice depth change. The SIIT is intriguing because, if ice depth is known, the heat content of an ice column can be calculated by assuming a nearly constant lapse rate within the ice layer, at least during winter. The amount of heat content held by the ice should be the main factor controlling upward heat transfer from the ocean water to the atmosphere. It builds up a preconditioning of how fast the ice melts in spring and summer, or how fast the ice grows in fall and winter. Thus, the SIIT is valuable information, which can lead to better climate monitoring over the Arctic. In fact, the aforementioned sea ice decline, changes in the Arctic hydrological cycle, and influences on the mid-latitude weather are all related to the SIIT because the SIIT is the variable that can give information on the thermal state of the ice.
In this study, we propose an algorithm to produce the long-term SIIT from SSM/I microwave (MW) measurements, which can penetrate the snow layer without much interference. A recent study showed that 6.9 GHz brightness temperature (TB) measurements are theoretically related to the sea ice emissivity and temperature [20]. Although the retrieved ice temperature from 6.9 GHz TBs is in good agreement with the buoy-measured SIIT, the acquisition of long-term data may not be possible, because 6.9 GHz TB measurements before the launch of the Advanced Microwave Scanning Radiometer-EOS (AMSR-E) in August 2002 are not available.
In reality, long-term global MW observations are available, as shown in the intercalibrated Fundamental Climate Data Record (FCDR) of TBs [21], covering more than 30 years, but without the 6.9 GHz channel. In the higher frequency region, such as 19.4 GHz, it is difficult to apply the Fresnel assumption for the emissivity retrieval due to surface and volume scatterings. However, if the influences of surface and volume scatterings are known and included in the radiative transfer calculation within the sea ice, it may be possible to retrieve the temperature even using higher frequencies.
Aiming at producing long-term SIIT over the Arctic sea ice, we develop an algorithm for retrieving the SIIT from SSM/I TB data, by counting surface and volume scattering effects on the measured TBs, which were examined in the previous study [22]. For this objective, we introduce a correction factor to remove surface and volume scattering influences from measured TBs. Subsequently, the radiative transfer equation is solved for the physical emission temperature from SSM/I (or SSMIS alike) TBs only using the correction factor and combined Fresnel equation. Based on this methodology, we attempt to produce long-term sea ice temperature over the Arctic. The obtained long-term SIIT will certainly help to understand how the global warming has influenced the Arctic climate change in the past decades.

Passive MW Data
In this study, we use FCDR SSM/I 19.4 and 37.0 GHz vertically and horizontally polarized TBs to retrieve the sea ice temperature. The SSM/I sensors have measured MW TBs with a viewing angle near 53 • since 1987. The sensors aboard different platforms are physically consistent but depend on varying orbital and navigational quantities such as the incident angle, resolution, and observation time. The objective of SSM/I FCDR data is to generate consistent and well-calibrated MW TBs from 1987 onwards by using multiple techniques for intercalibration [21,23]. In this study, these SSM/I FCDR orbital data were converted to 25 km Equal-Area Scalable Earth (EASE)-grid data, on a daily basis, to collocate SSM/I data with gridded AMSR-E TBs and SSM/I bootstrap sea ice concentration data [24]. More information on the SSM/I FCDR dataset can be found in References [21,23] and on the webpage of the Precipitation Research Group of the Department of Atmospheric Science of the Colorado State University (http:/rain.atmos.colostate.edu/FCDR/).
The AMSR-E 6.9 GHz TBs are used for retrieving the adjusted real refractive index and emitting layer temperature of Arctic sea ice, adapting a method described in Reference [20]. Level 3 daily TBs at 6.9 GHz were obtained from the Global Change Observation Mission (GCOM) data center site (http://gcom-w1.jaxa.jp/). In addition, bootstrap sea ice concentration data are applied to minimize the influences of water contamination. When sea ice coexists with open water within the same pixel, a heterogeneous problem arises in the interpretation of the MW brightness temperature, making the retrieval difficult. In this study, we minimized the water contamination effect in the retrieval by imposing 98% of sea ice concentration as a criterion. Sea ice concentration data were downloaded from the National Snow and Ice Data Center (NSIDC) website (http://nsidc.org/) [24].

CRREL Buoy Measurements
To validate the developed algorithm, we compared the pixel-based retrieved sea ice temperature with in situ measurements from U.S. Navy Cold Regions Research and Engineering Laboratory (CRREL) ice mass balance (IMB) buoy measurements at a daily time scale [25]. The buoy incorporates a thermometer string instrument to measure the temperature profile from the near-surface atmosphere, through the sea ice, and to sea water. In addition, the snow depth and sea ice thickness are measured with an acoustic sounder installed at the buoy. Daily averages were taken from 2-hourly measured data, which were downloaded from the CRREL website (http://imb-crrel-dartmouth.org/imb.crrel/). Temperature at zero ice depth (which is also referred to as SIIT but from buoy measurement) is used for the comparison with collocated SSM/I-derived SIIT.
Two types of CRREL IMB buoy data are available, i.e., "Provisional Buoy" data and "Archived Buoy" data. The "Provisional Buoy" data were archived without quality control, and thus often abnormal features are found in the "Provisional Buoy" data, i.e., (1) the air temperature is substantially higher than snow/ice interface (or zero ice depth) temperature even during the winter time; (2) a strong inversion layer exists within the sea ice; (3) the air temperature variation is smaller than that of temperature at zero ice depth; and (4) the air/snow interface or at snow/ice interface show no vertical temperature gradient change. The "Archived Buoy" data are quality-controlled, but, even some of "Archived Buoy" data show abnormal features similar to that found in the "Provisional Buoy" data, although there are relatively few cases. In this study, we only used "Archived Buoy" data for the comparison because of the abnormal features found in "Provisional Buoy" data. However, more control checks are necessary even for "Archived Buoy" data, and we apply quality assurance checks (1)-(4) to "Archived Buoy" data before the comparison. Details of the quality checks are provided in Section 4.

Theoretical Background and Methodology
The theoretical background for the algorithm development was developed and introduced in a previous study ( [22], this reference, (Lee, Sohn, and Shi, 2018), is now referred to as LSS18). Here, we provide a summary of the theoretical backgrounds used in that study, and an algorithm development for SSM/I measurements. Detailed explanations are provided in Appendix A.

Simple Snow/Ice Model for Radiative Transfer
In this study, sea ice used for the radiative transfer is composed of three layers, i.e., a snow layer at the top, surface ice layer (or freeboard sea ice layer) in the middle below the snow layer, and congelation ice layer as the lowest layer [26]. The snow layer (top layer) and upper part of the surface ice layer (middle layer) are presumed to be the volume scattering layer that attenuates the radiation emanating from the emission layer below. The model has one emission layer with temperature and related emissivity (T E and ε, respectively). The emission layer is located below the volume scattering layer and composed of the lower part of the surface ice layer and/or upper part of the congelation ice layer (see Figure 1 in Reference [22]). For the radiative transfer, surface and volume scatterings occurring above the emission layer are treated independent of the emission layer, not allowing interactions between the two layers.
If TB sfc is defined as an upwelling TB at the snow top, then the radiative transfer through the snow-ice layer can be expressed as follows: where τ sca is the optical thickness associated with the snow/ice volume scattering.

Apparent Emissivity
We defined the apparent emissivity (ε A ) by taking the ratio between T E and TB sfc : Apparent emissivities at the SSM/I 19.4 GHz channels are obtained by taking the ratio of SSM/I 19.4 GHz TB sfc to the sea ice emission temperature (T E ), according to Equation (2). TB sfc here is from the measured TB after removing atmospheric effects, which are calculated using the line-by-line radiative transfer model of the Satellite Data Simulator Unit (SDSU)-version 2.1 [27] with inputs of atmospheric temperature and humidity profiles from European Centre for Medium-Range Weather Forecasts Reanalysis (ERA)-Interim data [28]. The collocated T E is obtained from the AMSR 6.9 GHz retrieved SIIT, based on the finding that its emission weighting function is similar to that for 19.4 GHz. Figure 1 shows the spatial distributions of the apparent emissivities for both polarizations at 19.4 GHz on 1 January 2004. The obtained apparent emissivities are used for calculating the so-called "correction factor" described in Section 3.3. The apparent emissivity can also be related to scatterings, as shown in Equation (1): By using known ε A distributions in Figure 1, we try to solve the unknowns on the right-hand side of Equation (3).
Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 15 By using known εA distributions in Figure 1, we try to solve the unknowns on the right-hand side of Equation (3).

Correction Factor
Since ε in Equation (1) represents the emissivity for a rough surface, we needed to express ε with an emissivity for the smooth surface (εs). This is because the emissivity characteristics of a smooth surface based on the Fresnel equations are well understood, and thus we intended to use the combined version of the Fresnel equation [29,30] to solve the unknowns. In LSS18, a concept regarding the surface roughness effect was developed by relating the smooth surface emissivity to the rough surface emissivity: where F is the roughness parameter. Based on the concept in Equations (3) and (4), F and τsca at SSM/I 19.4 GHz frequency can be determined by introducing a two-dimensional surface roughness model and using the adjusted real refractive index obtained together with TE from AMSR-E 6.9 GHz TB measurements. The detailed method including the two-dimensional surface roughness model to estimate F and τsca for a given εA can be found in LSS18.
The correction factor (CF) can be defined as the ratio between the apparent emissivity and smooth surface emissivity. Since the smooth surface emissivity is altered by the surface roughness and volume scattering extinction in the constructed model, CF can be written as: After solving F and τsca for 19.4 GHz with given εA, the correction factor CF for 19.4 GHz is calculated using Equation (5). The CF data for 19.4 GHz are constructed over the Arctic sea ice region for the December-January-February (DJF) period of 2003-2011 (AMSR-E period). We present the CF distribution on 1 January 2004 in Figure 2a.
The Equation (5) can be written for two polarized components (vertical (V) and horizontal (H) components for CF, εA, εs, and F). Thus, with the help of the combined Fresnel equation for emissivity given in Equation (6), εs can be obtained by solving Equation (5) [30].

Correction Factor
Since ε in Equation (1) represents the emissivity for a rough surface, we needed to express ε with an emissivity for the smooth surface (ε s ). This is because the emissivity characteristics of a smooth surface based on the Fresnel equations are well understood, and thus we intended to use the combined version of the Fresnel equation [29,30] to solve the unknowns. In LSS18, a concept regarding the surface roughness effect was developed by relating the smooth surface emissivity to the rough surface emissivity: where F is the roughness parameter. Based on the concept in Equations (3) and (4), F and τ sca at SSM/I 19.4 GHz frequency can be determined by introducing a two-dimensional surface roughness model and using the adjusted real refractive index obtained together with T E from AMSR-E 6.9 GHz TB measurements. The detailed method including the two-dimensional surface roughness model to estimate F and τ sca for a given ε A can be found in LSS18. The correction factor (CF) can be defined as the ratio between the apparent emissivity and smooth surface emissivity. Since the smooth surface emissivity is altered by the surface roughness and volume scattering extinction in the constructed model, CF can be written as: After solving F and τ sca for 19.4 GHz with given ε A , the correction factor CF for 19.4 GHz is calculated using Equation (5) The Equation (5) can be written for two polarized components (vertical (V) and horizontal (H) components for CF, ε A , ε s , and F). Thus, with the help of the combined Fresnel equation for emissivity given in Equation (6), ε s can be obtained by solving Equation (5) [30].
So far, ε A and CF distributions are constructed for the winter (DJF) of the 2003-2011 period during which AMSR-E and SSM/I data coexist, to develop an SSM/I-only algorithm.

Development of the SSM/I-Only Algorithm for the Correction Factor
The retrieval of long-term sea ice TE data from SSM/I FCDR data beyond the AMSR-E period depends on how we can obtain the CF from SSM/I measurements alone. For that purpose, we utilize a finding of strong linear relationship between the CF and spectral gradient ratio (GR in Equation (7)). Here, we attempted to develop a regression algorithm determining CF at 19.4 GHz using SSM/I TB measurements. Now we relate the estimated CF (in Figure 2a) to collocated SSM/I TBs and GR between 19.4 and 37.0 GHz at vertical polarization: 37 19 19 ,37 37 19 where ai (i = 0, 1, 2, and 3) are regression coefficients (Table 1), which are obtained using daily CF and TBs data for DJF months of the 2003-2011 period, except for 2004, over the area with a sea ice concentration >98%.  It is noted that the obtained coefficients are slightly different between polarizations, but spatial distributions and magnitudes of CF R at both polarizations are nearly identical (not shown). It is thought that the surface scattering influence on the apparent emissivity appears to be minor for sea ice in the case of the 19.4 GHz frequency, in spite of surface scattering capable of affecting the

Development of the SSM/I-Only Algorithm for the Correction Factor
The retrieval of long-term sea ice T E data from SSM/I FCDR data beyond the AMSR-E period depends on how we can obtain the CF from SSM/I measurements alone. For that purpose, we utilize a finding of strong linear relationship between the CF and spectral gradient ratio (GR in Equation (7)). Here, we attempted to develop a regression algorithm determining CF at 19.4 GHz using SSM/I TB measurements. Now we relate the estimated CF (in Figure 2a) to collocated SSM/I TBs and GR between 19.4 and 37.0 GHz at vertical polarization: where a i (i = 0, 1, 2, and 3) are regression coefficients (Table 1), which are obtained using daily CF and TBs data for DJF months of the 2003-2011 period, except for 2004, over the area with a sea ice concentration >98%. It is noted that the obtained coefficients are slightly different between polarizations, but spatial distributions and magnitudes of CF R at both polarizations are nearly identical (not shown). It is thought that the surface scattering influence on the apparent emissivity appears to be minor for sea ice in the case of the 19.4 GHz frequency, in spite of surface scattering capable of affecting the polarization state [22]. By applying the regression equation, geographical distributions of CF R are obtained for 1 January 2004 (Figure 2b), which can be compared with the AMSR-derived CF in Figure 2a. The figure shows that they are in good agreement, with small regional differences within 0.05 over the Arctic Ocean (Figure 3a). The scatterplot between CF and CF R for one month (January 2004) also clearly shows that regression-based SSM/I CF R retrievals are comparable to the direct CF, showing a correlation coefficient of 0.98. The bias and RMSE in this comparison are −0.0004 and 0.0053, respectively. In order to examine the robustness of the regression Equation (8) Table 2. This suggests that the regression coefficients obtained in this study work well for other years. Because of this, we apply the coefficients in Table 1 to the whole SSM/I period.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 15 polarization state [22]. By applying the regression equation, geographical distributions of CF R are obtained for 1 January 2004 (Figure 2b), which can be compared with the AMSR-derived CF in Figure  2a. The figure shows that they are in good agreement, with small regional differences within 0.05 over the Arctic Ocean (Figure 3a). The scatterplot between CF and CF R for one month (January 2004) also clearly shows that regression-based SSM/I CF R retrievals are comparable to the direct CF, showing a correlation coefficient of 0.98. The bias and RMSE in this comparison are −0.0004 and 0.0053, respectively. In order to examine the robustness of the regression Equation (8) Table 2. This suggests that the regression coefficients obtained in this study work well for other years. Because of this, we apply the coefficients in Table 1 to the whole SSM/I period.
Consequently, the results demonstrate that the CF can be predicted from Equation (8) with SSM/I TBs only. Therefore, by incorporating CF R into Equation (5), εA can be calculated by applying CF to εs. Once εA is available, it is straightforward to estimate the emitting layer temperature TE from TBsfc using Equation (2).

SSM/I Snow/Ice Interface Temperature (SIIT)
Since the emissivities for the smooth surface and regressed CF R values are now available for SSM/I 19.4 GHz, the apparent emissivities can be calculated from Equation (5), i.e., The spatial distributions of regression-based apparent emissivities are given in Figure 4. The distributions are almost the same as that for the directly calculated apparent emissivity in Figure 1, demonstrating that apparent emissivity can be estimated by solely using SSM/I measurements. Because

R A
 is now available for the given SSM/I measurements, the emission temperature TE can be estimated. However, because TE is in fact the temperature between snow and ice layers, as expressed in Reference [11], the obtained emission temperature is also referred to as SIIT.  Consequently, the results demonstrate that the CF can be predicted from Equation (8) with SSM/I TBs only. Therefore, by incorporating CF R into Equation (5), ε A can be calculated by applying CF to ε s . Once ε A is available, it is straightforward to estimate the emitting layer temperature T E from TB sfc using Equation (2).

SSM/I Snow/Ice Interface Temperature (SIIT)
Since the emissivities for the smooth surface and regressed CF R values are now available for SSM/I 19.4 GHz, the apparent emissivities can be calculated from Equation (5), i.e., ε R A = CF R ε s . The spatial distributions of regression-based apparent emissivities are given in Figure 4. The distributions are almost the same as that for the directly calculated apparent emissivity in Figure 1, demonstrating that apparent emissivity can be estimated by solely using SSM/I measurements. Because ε R A is now available for the given SSM/I measurements, the emission temperature T E can be estimated. However, because T E is in fact the temperature between snow and ice layers, as expressed in Reference [11], the obtained emission temperature is also referred to as SIIT.    Figure 6a for 1 January 2004. The temperature difference distributions may be attributed to the difference distributions between CF and CF R because the unexplained CF variance in the regression should be propagated into the temperature difference. The two difference maps show very similar features (i.e., Figures 3a and 6a). It is noted that the temperature difference is generally within 1.5 K regardless of the ice type. A scatterplot of the temperature difference between two for January month of 2004 is provided in Figure 6b. Most of the pairs are near the one-to-one line, implying that both temperatures have a strong   Figure 6a for 1 January 2004. The temperature difference distributions may be attributed to the difference distributions between CF and CF R because the unexplained CF variance in the regression should be propagated into the temperature difference. The two difference maps show very similar features (i.e., Figures 3 and 6). It is noted that the temperature difference is generally within 1.5 K regardless of the ice type.   Figures 3a and 6a). It is noted that the temperature difference is generally within 1.5 K regardless of the ice type. A scatterplot of the temperature difference between two for January month of 2004 is provided in Figure 6b. Most of the pairs are near the one-to-one line, implying that both temperatures have a strong

Validation against the Buoy Snow/Ice Interface Temperature
Validating the developed algorithm, retrieved SSM/I sea ice temperatures are compared with CRREL IMB drifting buoy measurements over the Arctic sea ice area [19]. For the comparison, gridbased SSM/I-retrieved temperatures are collocated with buoy data by applying the following steps: (1) Is the center point of any given SSM/I target with 25 km × 25 km resolution target located within a 12.5 km distance from the buoy location at a given time? (2) Are the corresponding sea ice concentrations greater than 98%? (3) Are the selected satellite pixels entirely covered by sea ice? (4) Select a closest target value if multiple SSM/I targets satisfy conditions (1), (2), and (3). Although the buoy data pass the above-mentioned criteria, the data are rejected if snow/ice interface information is not provided or if the data failed the quality assurance criteria given in Section 2.2. The remaining measurements along 13 deployed tracks are taken for the validation. The observation periods and start and end latitudes and longitudes are given in Table 3. Among the 13 buoy tracks, 8 were found in the northern Beaufort Sea and 5 were found in the Central Arctic Ocean, as shown in Figure S1 in the Supplementary Material. Figure 7a shows a scatterplot of the buoy-measured SIITs and corresponding SSM/I values along all 13 buoy tracks, with a correlation coefficient of 0.86, bias of -0.36 K, and RMS error of 3.42 K. The plot shows a scattered pattern, in which the daily variation is greater than 7 K. Except for the blue dots, a large bundle of black points is located along the 1:1 line. It is noted that the scattered blue dots are associated with three buoys (#2, #6, and #8 in Table 3) in which temperature variations at the snow/ice interface are greater than 7 K for some data points. Close examination of those data reveals that such large variations are not due to the shallow ice depth or ice type (not shown). It is thought that they are likely due to problems related to the buoy instruments themselves. Rationally thinking, such large temperature variations within the ice may not be acceptable such that 7 K is subjectively applied as a criterion to remove such buoy measurements.
Because of those reasons, we exclude the blue dots in Figure 7a from the paired data set for the comparison. Comparison results obtained after applying quality assurance checks indicate a correlation coefficient of 0.95, bias of 0.15 K, and RMS error of 1.48 K (Figure 7b). It is important to note that the SSM/I-derived physical temperature of the sea ice surface shows a good 1:1 correspondence to buoy-measured SIIT, indicating that the temperature obtained for SSM/I in this

Validation against the Buoy Snow/Ice Interface Temperature
Validating the developed algorithm, retrieved SSM/I sea ice temperatures are compared with CRREL IMB drifting buoy measurements over the Arctic sea ice area [19]. For the comparison, grid-based SSM/I-retrieved temperatures are collocated with buoy data by applying the following steps: (1) Is the center point of any given SSM/I target with 25 km × 25 km resolution target located within a 12.5 km distance from the buoy location at a given time? (2) Are the corresponding sea ice concentrations greater than 98%? (3) Are the selected satellite pixels entirely covered by sea ice? (4) Select a closest target value if multiple SSM/I targets satisfy conditions (1), (2), and (3). Although the buoy data pass the above-mentioned criteria, the data are rejected if snow/ice interface information is not provided or if the data failed the quality assurance criteria given in Section 2.2. The remaining measurements along 13 deployed tracks are taken for the validation. The observation periods and start and end latitudes and longitudes are given in Table 3. Among the 13 buoy tracks, 8 were found in the northern Beaufort Sea and 5 were found in the Central Arctic Ocean, as shown in Figure S1 in the Supplementary Material. Figure 7a shows a scatterplot of the buoy-measured SIITs and corresponding SSM/I values along all 13 buoy tracks, with a correlation coefficient of 0.86, bias of -0.36 K, and RMS error of 3.42 K. The plot shows a scattered pattern, in which the daily variation is greater than 7 K. Except for the blue dots, a large bundle of black points is located along the 1:1 line. It is noted that the scattered blue dots are associated with three buoys (#2, #6, and #8 in Table 3) in which temperature variations at the snow/ice interface are greater than 7 K for some data points. Close examination of those data reveals that such large variations are not due to the shallow ice depth or ice type (not shown). It is thought that they are likely due to problems related to the buoy instruments themselves. Rationally thinking, such large temperature variations within the ice may not be acceptable such that 7 K is subjectively applied as a criterion to remove such buoy measurements. study represents the temperature near the snow/ice interface defined using buoy measurements in which the typical emitting layer lay in MW frequencies.  As previously discussed, we intended to produce long-term SSM/I-derived temperature information over the Arctic sea ice area. Since we developed an algorithm to retrieve the SIIT, we give an attempt to produce 24 years of SSM/I SIIT. For demonstration purposes, we generated the mean December SIIT from 1988 to 2011 and the results are given in Figure 8. The figure clearly indicates that the temperature is significantly colder in the early years of the period, suggesting that the Arctic sea ice had experienced a rapid temperature rise during the 1988-2011 period. This type of data provides important information about the evolution of climate change over the Arctic Ocean in the last decades, if we extend this effort to include the SSMIS, which effectively replaces the SSM/I. Because of those reasons, we exclude the blue dots in Figure 7a from the paired data set for the comparison. Comparison results obtained after applying quality assurance checks indicate a correlation coefficient of 0.95, bias of 0.15 K, and RMS error of 1.48 K (Figure 7b). It is important to note that the SSM/I-derived physical temperature of the sea ice surface shows a good 1:1 correspondence to buoy-measured SIIT, indicating that the temperature obtained for SSM/I in this study represents the temperature near the snow/ice interface defined using buoy measurements in which the typical emitting layer lay in MW frequencies.
As previously discussed, we intended to produce long-term SSM/I-derived temperature information over the Arctic sea ice area. Since we developed an algorithm to retrieve the SIIT, we give an attempt to produce 24 years of SSM/I SIIT. For demonstration purposes, we generated the mean December SIIT from 1988 to 2011 and the results are given in Figure 8. The figure clearly indicates that the temperature is significantly colder in the early years of the period, suggesting that the Arctic sea ice had experienced a rapid temperature rise during the 1988-2011 period. This type of data provides important information about the evolution of climate change over the Arctic Ocean in the last decades, if we extend this effort to include the SSMIS, which effectively replaces the SSM/I. Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 15

Discussion and Conclusions
Aiming at producing a long-term temperature record over the Arctic sea ice area, we developed an algorithm to retrieve the wintertime SIIT over the Arctic Ocean from SSM/I FCDR data. In so doing, a simple radiative transfer model was set up, in which the snow-ice is composed of three layers (i.e., snow layer, volume scattering ice layer, and emitting ice layer, from the top). In the model, the emitting radiation from the emission layer is attenuated via volume scatterings from snow/ice layers above as well as by ice surface scattering. By using the collocated AMSR-E 6.9 GHz derived SIIT as an emission temperature for SSM/I 19.4 GHz, the apparent emissivity was determined from SSM/I 19.4 GHz TBs. The obtained apparent emissivity was then related to volume scattering and surface roughness to examine how these parameters can be related to the correction factor (i.e., the ratio between apparent emissivity at the surface and specular emissivity at the emission layer). Then the correction factor was defined and predicted from a combination of SSM/I TBs (i.e., TB19V and TB37V). The correction factor from SSM/I measurements can be used to predict the corresponding apparent emissivity from the combined Fresnel equation. It is straightforward to estimate the SIIT from SSM/I 19.4 GHz measurements once the correction factor is determined from SSM/I measurements.
The developed algorithm has been successfully implemented for SSM/I TB measurements available from July 1987. The retrieved temperatures were validated against the temperatures at the snow-ice interface (or the temperature at the top of the ice layer) from CRREL IMB drifting buoy measurements. The validation results indicated that the SSM/I-derived SIITs are in good agreement with drifting buoy measurements of the SIIT (temperature at the zero ice depth), with a correlation coefficient of 0.95, bias of 0.1 K, and RMS error of 1.48 K on a daily time scale. These results suggest that a long-term ice temperature record can be produced by applying the developed algorithm to well-calibrated long-term FCDR data composed of SSM/I and SSMIS data spanning over thirty years.

Discussion and Conclusions
Aiming at producing a long-term temperature record over the Arctic sea ice area, we developed an algorithm to retrieve the wintertime SIIT over the Arctic Ocean from SSM/I FCDR data. In so doing, a simple radiative transfer model was set up, in which the snow-ice is composed of three layers (i.e., snow layer, volume scattering ice layer, and emitting ice layer, from the top). In the model, the emitting radiation from the emission layer is attenuated via volume scatterings from snow/ice layers above as well as by ice surface scattering. By using the collocated AMSR-E 6.9 GHz derived SIIT as an emission temperature for SSM/I 19.4 GHz, the apparent emissivity was determined from SSM/I 19.4 GHz TBs. The obtained apparent emissivity was then related to volume scattering and surface roughness to examine how these parameters can be related to the correction factor (i.e., the ratio between apparent emissivity at the surface and specular emissivity at the emission layer). Then the correction factor was defined and predicted from a combination of SSM/I TBs (i.e., TB19V and TB37V). The correction factor from SSM/I measurements can be used to predict the corresponding apparent emissivity from the combined Fresnel equation. It is straightforward to estimate the SIIT from SSM/I 19.4 GHz measurements once the correction factor is determined from SSM/I measurements.
The developed algorithm has been successfully implemented for SSM/I TB measurements available from July 1987. The retrieved temperatures were validated against the temperatures at the snow-ice interface (or the temperature at the top of the ice layer) from CRREL IMB drifting buoy measurements. The validation results indicated that the SSM/I-derived SIITs are in good agreement with drifting buoy measurements of the SIIT (temperature at the zero ice depth), with a correlation coefficient of 0.95, bias of 0.1 K, and RMS error of 1.48 K on a daily time scale. These results suggest that a long-term ice temperature record can be produced by applying the developed algorithm to well-calibrated long-term FCDR data composed of SSM/I and SSMIS data spanning over thirty years.
It should be mentioned that caution should be exercised if the developed algorithm is applied for thin ice because the penetration depth of the MW spectrum of interest in this study can penetrate deeper to reach the ocean water below the ice layer if the ice depth is less than approximately 16 cm, which is twice the penetration depth for first-year sea ice at 6.9 GHz with typical dielectric constant of 3.3-i0.15 [31]. In other words, a portion of the signal may be from underlying water although the magnitude is expected to be small. Additionally, the algorithm is only valid for high-concentrated sea ice pixels because the regression coefficients for the correction factor are based on TB data over the ice area with an ice concentration >98%. Therefore, when the sea ice concentration (SIC) decreases with melt ponds on the sea ice during the summer ice melt, the current approach may produce erroneous results. Furthermore, the regression coefficients for CF R may depend on the SIC products because we tried to remove possible water-contaminated pixels. We used 98% of the SIC obtained from National Snow Ice Data Center (NSIDC) product based on the National Aeronautics and Space Administration (NASA)'s bootstrap algorithm as a criterion to remove such pixels. It is known that different algorithms for the SIC (e.g., Ocean and Sea Ice Satellite Application Facility (OSI SAF), NASA Team, and Bristol algorithms among others) result in different SIC distributions with respect to seasons and regions [32][33][34][35]. However, considering that the water contamination should be largely confined within the marginal sea ice area during the winter time which we focus on in this study, the use of different data sets may provide results that do not significantly deviate from the current one. Nonetheless, it may be worthwhile to examine how different SIC data give different results, especially over the marginal sea ice area.