Snow parameter retrieval (SPR) algorithm for the GCOM-C/SGLI sensor: validation over the Greenland ice sheet

Abstract

This paper presents a new snow parameter retrieval (SPR) algorithm for the Global Change Observation Mission-Climate/Second Generation Global Imager (GCOM-C/SGLI) instrument (2018-present). This algorithm combines accurate radiative transfer model (RTM) simulations and Scientific Machine Learning (SciML) methods, Multi-Layer Neural-Network (MLNN) techniques in particular. It provides pixel-by-pixel optically equivalent snow grain size in two layers (i.e., a thin surface snow layer and a deep snow layer), snow impurity concentration and broadband blue- and black-sky albedo which constitute standard SGLI products. In addition, this RTM-SciML algorithm retrieves aerosol optical depth and provides an important retrieval error quality flag. This retrieval error flag, established by comparing reflectances obtained from RTM simulations using the retrieved snow and aerosol parameters as input with the measured reflectances, provides a pixel-by-pixel quality check of the retrieval parameters. Application of the RTM-SciML algorithm to SGLI images obtained over the Greenland Ice Sheet revealed a significant change in snow parameters from a cold July 2018 to a warm July 2019. The inferred blue-sky albedo was in general agreement with field measurements with RMSE = 0.0517 and MAPE = 4.64% for shortwave albedo at the SIGMA-A site, and the black-sky albedo, inferred from retrieved snow parameters, was found to be similar (within 5% relative difference) to the blue-sky values. Although developed specifically for application to data obtained by the SGLI imager, the SPR algorithm can easily be adapted for application to other similar multi-spectral sensors, such as MODIS (already done), VIIRS, and OLCI.

1 Introduction

Long-term global mapping of snow albedo and snow property parameters plays an important role in monitoring of the Earth climate system. Satellite remote sensing has offered a very valuable and powerful way to record the evolution of global snow extent and properties with high temporal and spatial resolution (Frei et al., 2012; Deltz et al., 2012; Tedesco, 2015). The visible and near infrared bands can be used to obtain snow coverage, broadband albedo, and snow physical parameters (Stamnes et al., 2007; Lyapustin et al., 2009b; Zege et al., 2011; Kokhanovsky et al., 2019). In particular, the retrieval of snow physical parameters provides a direct way to determine the spectral albedo and the snow surface radiation budget.

Remote sensing algorithms for snow and ice physical properties have been developed for different satellite sensors, such as the GLobal Imager (GLI) (Stamnes et al., 2007), the Moderate resolution imaging spectroradiometer (MODIS) (Lyapustin et al., 2009b; Painter et al., 2009; Zege et al., 2011), and the Ocean and Land Colour Instrument (OLCI) (Kokhanovsky et al., 2019). Some of these algorithms (Stamnes et al., 2007; Lyapustin et al., 2009b) were based on radiative transfer (RT) model generated lookup tables (LUTs), while others (Kokhanovsky et al., 2019; Zege et al., 2011) were based on analytical methods. A physically based, multiple endmember spectral mixture method (Painter et al., 2009) has also been applied to land surfaces consisting of a mixture of snow, vegetation, rocks, and soil. Sirguey et al. (2009) developed a comprehensive atmospheric correction and inversion scheme to retrieve snow parameters in mountainous areas. More recently, Bair et al. (2021) attempted to retrieve snow impurity concentration under complex sub-pixel mixing conditions using a generalized multispectral unmixing approach. However, these previous methods did not solve the radiative transfer equation for a coupled atmosphere-snow system with multiple layers in the atmosphere and a two-layer snowpack consisting of nonspherical snow grains.

The Second Generation Global Imager (SGLI) onboard the Global Change Observation Mission-Climate (GCOM-C) satellite operated by the Japan Aerospace Exploration Agency (JAXA) is aimed at global and long-term observations of the Earth environment. Launched in December 2017, SGLI has a wide spectral coverage from 380 nm to 12 m, a high spatial resolution of 250 m, and a field of view exceeding 1,000 km. In this paper, we present the radiative transfer model (RTM) simulations and Scientific Machine Learning (SciML) snow parameter retrieval (SPR) algorithm developed to retrieve standard SGLI snow products. The SPR algorithm shares a number of basic concepts and ideas with our algorithm designed for the previous GLobal Imager (GLI) sensor (Stamnes et al., 2007). It employs seven SGLI bands extending from 380 nm to 2,210 nm to retrieve snow products including optically equivalent snow grain size in two layers (a thin surface layer and the bulk of the snowpack), and snow broadband (blue- and black-sky) albedo on a pixel-by-pixel basis. The structure of this paper is as follows. Section 2 gives an overview and provides detailed information about the SPR algorithm. Section 3 describes the satellite and ground-based measurements used in this paper. Section 3.1 shows some case study results of the SPR algorithm where simultaneous measurements from the field and SGLI are available. Section 3.2 provides validation results based on SGLI images obtained over GrIS in July of 2018–2021. A summary and conclusions are provided in Section 4. Detailed information about the RTM and the snow and atmosphere Inherent Optical Properties (IOPs) used in the RT simulations is provided in the Supplementary Appendix (Section 5).

2 Snow parameter retrieval (SPR) algorithm

In contrast to the look-up-table approach adopted for GLI (Stamnes et al., 2007), several important improvements have been made in the RTM-SciML algorithm for SGLI. These improvements include (i) employing a comprehensive Machine Learning Classification Mask (MLCM) for simultaneous cloud screening and surface classification (Chen et al., 2018; Zhou et al., 2023); (ii) using a coupled atmosphere-surface RTM (Stamnes et al., 2018) to create a large synthetic dataset of simulated top-of-the-atmosphere (TOA) reflectances as a function of snow and aerosol physical parameters; (iii) using this synthetic dataset to train a multi-layer neural-network (MLNN) for the retrieval, which has led to significant improvements in both retrieval accuracy and speed; (iv) using a non-spherical (Voronoi) particle model, instead of a spherical snow particle model, which has been shown to be more realistic (Ishimoto et al., 2012; Ishimoto et al., 2018; Tanikawa et al., 2020); (v) using a pseudo-spherical correction to consider Earth curvature, which provides more accurate simulations in high latitude areas; (vi) using a “filtered” (instead of a random) distribution of snow and aerosol parameters to generate the synthetic dataset used for MLNN training, which mimics a more realistic snow situation, and leads to significantly improved retrievals. The “filtering” approach is described in Section 2.2.1.

The SPR algorithm is part of the RTM-SciML remote sensing algorithm framework as illustrated in Figure 1. In this paper we will discuss its application to snow-covered targets over land. The RTM-SciML framework consists of four key components: (i) a comprehensive Machine Learning Classification Mask (MLCM) for simultaneous cloud screening and surface classification as described in Section 2.1; (ii) a retrieval algorithm which employs multi-layer neural-network (MLNNs) trained by data generated using a forward RTM (Section 2.2.1) to retrieve snow/aerosol physical parameters as described in Sections 2.2.2 and 2.2.3, and (iii) an uncertainty estimation method as described in Section 2.2.5.

FIGURE 1

2.1 Machine learning classification mask (MLCM)

Cloud screening is an essential first step in the satellite retrieval chain. Therefore, a Machine Learning Classification Mask (MLCM) was first applied to SGLI images to identify clear-sky overlying snow-only pixels. The current version of the SPR algorithm will be applied only to cloud-screened (clear-sky), snow-only pixels. Hence, the quality of the MLCM will control the quality of the information inferred from the SPR retrieval algorithm. The MLCM algorithm employs Scientific Machine Learning (SciML) methods (Chen et al., 2018; Baker et al., 2019; Zhou et al., 2023) in conjunction with a large synthetic dataset generated by RTM simulations for a coupled atmosphere-surface system (Chen et al., 2018; Zhou et al., 2023). This large synthetic dataset includes inherent optical property (IOP) data for a variety of aerosol and cloud (liquid water and ice) types as well as bidirectional reflectance distribution function (BRDF) data for a variety of surface types. From this MLCM tool, we obtain not only cloud-free pixels, but also a pixel-by-pixel classification of the underlying surface into several categories such as snow-only, mixed snow/vegetation (Chen et al., 2018), sea-ice, and liquid water (Zhou et al., 2023). As shown by Stillinger et al. (2019) traditional threshold-based cloud mask methods have problems over complex snow-covered terrain due to the similarity of snow and clouds in most reflectance channels as well as the surface sub-pixel mixing conditions of snow and vegetation/soil/rock. We have shown that this approach can be significantly improved through the threshold-free approach provided by the MLCM algorithm (Chen et al., 2018; Zhou et al., 2023).

The MLCM algorithm is generic in nature and could in principle be applied to any sensor with a suitable set of channels (Chen et al., 2018; Zhou et al., 2023). However, for our purpose it is important to note that our MLCM approach has been validated by collocated CALIPSO and MODIS data (Chen et al., 2014; 2018), and that it has been tailored specifically for application to the SGLI sensor. Pixels from SGLI images, classified by this MLCM algorithm as clear-sky, snow-only pixels, were used to develop and test the SPR retrieval algorithm.

2.2 Snow parameter retrieval algorithm trained by RT simulations

2.2.1 Forward model simulations and synthetic dataset generation

A comprehensive forward RTM (Stamnes et al., 2018) was used to simulate upward radiances at the top of the atmosphere (TOA). We used the sub-arctic summer atmosphere profile (Anderson et al., 1986), aerosol properties described in Supplementary Section 5.1.1.1, snow IOPs described in Supplementary Section 5.1.1.2, and a two-layer snow configuration described in Table 1. The SPR algorithm was designed with two snow layers to take into account the snow vertical structure. We fixed the thickness of the upper snow layer to be 1.5 cm and the lower layer to be 98.5 cm. Since the penetration depth for the short wavelength channels can exceed 10 cm (see Section 3.1 for details), we chose a total snow thickness of 1 m to make sure there will be no signal coming from the underlying surface. The upper layer snow thickness is within the range of the light penetration depth at the m wavelength. The snow grain size was allowed to be different in the two layers, but the impurity concentration was assumed to be the same. The snow densities were assumed to be 0.15 and 0.25 for the upper and lower snow layers, respectively. The aerosol parameters used to generate the synthetic dataset required to develop the RTM-SciML algorithm are described in Supplementary Section 5.1.1.1.

The broadband albedo for each case in the synthetic dataset was also simulated in two ways. One is for blue-sky (cloud-free atmosphere) broadband albedo, and the other is for black-sky (no atmosphere) broadband albedo. Both broadband albedo values were simulated in the following 3 bands: visible (300–700 nm, VIS) albedo; near-infrared (700–3,000 nm, NIR) albedo; and shortwave (300–3,000 nm, SW) albedo. The broadband albedo can be configured to other wavelength ranges as needed.

We generated a synthetic dataset consisting of 2.7 million cases after filtering (see Section 2.2.4 for details) for the MLNN training. This dataset includes simulated TOA radiances and broadband albedo values (for blue-sky and black-sky conditions) in the VIS, NIR, and SW wavelength ranges, as well as the broadband albedo measured by instruments deployed at AWS stations, for snow and aerosol parameters as well as geometry angles appropriate for SGLI measurements. The ranges of the physical parameters and the geometry angles covered in the synthetic dataset are provided in Table 1. This table applies specifically to the conditions of the GrIS, but can easily be extended for global application.

It has been argued (Warren, 2013) that satellite remote sensing of snow impurity content (like black carbon) is unlikely to be successful, except in highly polluted industrial regions. Since the visible albedo depends on snow impurity content, we cannot ignore this effect when retrieving snow grain size (Warren and Wiscombe, 1980). Hence, we try to get information of impurity content as well (Aoki et al., 2014b), although the retrieved values are expected to have large uncertainties or may even be below the detection limit of currently available instrumentation in space for the low impurity contents encountered in pristine high-latitude regions like the Greenland Ice Sheet (GrIS) as discussed by Warren (2013).

2.2.2 Algorithm training, consistency check and output

In the SPR algorithm, we used the RTM-SciML method employing Multi-Layer Neural-Networks (MLNNs) (described below) to (i) retrieve snow and aerosol parameters; (ii) retrieve black-sky albedo and blue-sky albedo, and (iii) estimate the retrieval error. The MLNN or multilayer perceptron (MLP) is a feed-forward artificial neural network that is frequently used for pattern classification, recognition, prediction, and function approximation. It has been demonstrated that MLNNs with one or more hidden layers and a non-linear activation function can be used to approximate nonlinear functions very well (Chen et al., 1990; D’Alimonte and Zibordi, 2003; D’Alimonte et al., 2004; 2012). Therefore, they are suitable for solving our inverse problem.

The transfer (activation) function of the neurons for the MLNN was taken to be the hyperbolic tangent function:

In the output layer a linear transfer function was used to link the hidden layers to the output. The exact expression of this MLNN can be written as:where , is an element in the input layer. It is important to note that if the goal is to make a retrieval of state parameters directly from TOA reflectance measurements, then the input parameters in Equation 2 are the ground elevation, the measured TOA reflectances at the seven SGLI channels plus the solar/viewing geometry, and the output parameters are the desired retrieval (state) parameters. We shall refer to this MLNN as R2P (Radiance Parameter) for short. On the other hand, if the goal is to use the MLNN [Equation 2] as a fast interpolator to obtain the TOA radiances and associated Jacobians, then the input parameters should be the state parameters and the solar/viewing geometry, and the output parameters should be the TOA radiances (Stamnes and Stamnes, 2016). We shall refer to this MLNN as P2R (Parameter Radiance) for short.

We use a 3-hidden-layer structure with 15 neurons for all MLNNs. In Equation 2, , , and are the weights of the three hidden layers and the output layer, and , , , and are the biases of the three hidden layers and the output layer. The weights and biases are to be determined by the training. The parameter (see Equation 1) is the hyperbolic tangent function, and is the element in the output layer, which in the R2P case is one of the four retrieved snow and aerosol parameters (top and bottom snow grain sizes, snow impurity concentration, and aerosol optical depth); hence . Figure 2 shows the R2P training performance with the statistical analysis values , defined as the coefficient of determination, and the root mean square error (RMSE) value.

FIGURE 2

For all algorithm trainings (as listed in Figure 2), the synthetic dataset was randomly divided into two independent groups: a training group (85% of the dataset), and a validation group (15% of the dataset). The training group was used to optimize the weights and biases of the neural network and the validation group was used to validate the neural network after the training was finished. The MLNN was trained using the ‘adam’ stochastic gradient descent method (Kingma and Ba, 2017) in Scikit-Learn, which employs an adaptive learning rate (Amari et al., 2000) to reduce the training time. The training algorithm updates the weights and biases iteratively, based on the residual between the training target and the neural network output. In each iteration, the validation dataset was used to monitor the performance of the current neural network by computing the root mean square error (RMSE).

An increase in the RMSE implies that the neural network was over-trained. Therefore, the training was terminated if the RMSE increased for 10 consecutive iterations in order to avoid over-fitting. A L2 regularization scheme (Neumaier, 1998) was also applied to the training algorithm to minimize the possibility of over-fitting.

The training configurations of the R2P, P2R and albedo MLNNs are listed in Table 2. The training performance of the MLNNs is shown in Figure 2. In general, the performance is good, although the correlation coefficient for the aerosol optical depth is larger than for the other parameters, because aerosol retrieval over bright surfaces like snow is a challenging task.

2.2.3 Algorithm self-consistency and feature selection

To investigate the consistency of our trained MLNN, we apply an autoencoder neural network (aaNN) (see Section 2.2.4) between the training dataset and the satellite measurements. For this purpose, we trained a special neural network by using our simulated TOA reflectance dataset in such a way that it attempts to learn an “identity data set” (i.e., the input data are identical to the output data) with the least number of hidden layers required to create a “compressed” representation of the input.

We used satellite measured TOA reflectances as the input to this trained network. If the output reflectances differ significantly from the input reflectances, it would imply that the satellite measurements are inconsistent with the radiances available in the synthetic dataset. Hence, this neural network should help indicate if the algorithm is reliable, and identify other issues within the data such as instrument saturation, sensor degradation or calibration issues. A low value would indicate that the simulated TOA reflectances are consistent with satellite measured data. But a high value would imply that the simulated TOA reflectances are not representative of the satellite measurements, which could be an indication that the real snow is “out of the range” of (i.e., inconsistent with) our synthetic data. The reason could be that (1) the range of the synthetic dataset is not large enough to cover the different types of snow; (2) the sensor has calibration or degradation issues; (3) there are cloud screening or surface classification issues due to contamination by cloud, vegetation (or snow impurities like red algae, not considered in our forward RTM simulations), or melting snow/ice (we assumed a dry snowpack in our forward model simulations). This test is performed on a channel-by-channel basis. Hence, by carefully examining the values, certain channels (features) could be removed from the dataset to improve its robustness.

To generate the synthetic dataset for the test, we used input layer elements consisting of the three geometry angles, the ground elevation, and the synthetic TOA reflectances in eight SGLI channels (380 nm, 412 nm, 530 nm, 674 nm, 869 nm, 1,050 nm, 1,640 nm, and 2,100 nm). The output layer elements were the TOA reflectances in these eight SGLI channels. The number of neurons for the three hidden layers were 10, 6, and 10, respectively. The is defined as follows:where denotes any one of the eight SGLI channels. The values are the original eight SGLI satellite measurements, while the values are the eight output TOA reflectances produced by the trained autoencoder neural network.

Figure 3 shows detailed values in each SGLI channel for 10-day (July 1–10, 2018) average reflectances. For a good consistency, the value should be below about 5%. The largest values were found to occur in the 674 nm and 530 nm channels. We found that there is a constant bias in these two channels and in the same area. The large error in the 674 nm band is due at least partly to a signal saturation issue occurring over the southern area of the GrIS. Therefore, in our snow retrieval, we excluded this channel and used only the remaining seven SGLI channels. For the 530 nm SGLI band, we conjecture that the problem may be a calibration issue, since SGLI post-launch calibration shows a 4% higher value when compared with EGRIP data (https://suzaku.eorc.jaxa.jp/GCOM_C/data/prelaunch/index_cal.html). Hence, the high positive values for the 530 nm band could be due partly to the calibration bias. But, the values seem to exceed 10% over the Southern GrIS (Figure 3) which might be too high to attribute solely to a calibration issue. Besides the calibration issue, the maximum radiance for the 530 nm band is designed to be 350 which can lead to saturation over very bright targets leading to high values. Due to the high reflectance of snow, these calibration/saturation issues are quite noticeable, especially over the southern area of the GrIS. For now, we continue to include the 530 nm band in our algorithm until JAXA confirms or refutes our conjecture regarding its calibration/saturation.

FIGURE 3

2.2.4 Synthetic dataset filtering

To train the autoencoder, we constructed a dataset consisting of 500,000 data points of SW, VIS, and NIR albedo values randomly selected from 2 months of MCD43 snow albedo data over GrIS and used this dataset to train the autoencoder. The low dimensionality of this albedo dataset based on only four parameters (aerosol optical depth, top layer snow grain size, bottom layer snow grain size, and snow impurity concentration) implies that it is prone to overfitting during the training of the aaNN. Here we used a 5-layer network with 4 neurons to construct the aaNN. Identical input and output of solar zenith angles, as well as the SW, VIS, and NIR albedo values were used in the training. The histogram of SW, VIS, and NIR albedo values of the MCD43 training dataset and the training results of the aaNN are shown in Figure 4, which indicates excellent training accuracy. Figure 4 also shows the histogram of albedo distributions of the synthetic dataset (bottom panels), which are different from those of the MCD43 data. We then applied the trained aaNN to the albedo dataset generated by using the noncentral distribution of snow parameters discussed above.

FIGURE 4

FIGURE 5

One should note that the MODIS MCD43 SW albedo product pertains to the wavelength range 300–5,000 nm which is different from the wavelength range (300–3,000 nm) used in our albedo simulations. For wavelengths longer than 3,000 nm the snow albedo is extremely low (usually 0.1%) and the incoming solar irradiance is very small. Therefore, we consider the use of autoencoders trained by MCD43 albedo products justifiable in spite of the wavelength range difference.

2.2.5 Retrieval error estimation (convergence check)

In order to estimate the quality of SPR results generated by the RTM-SciML method, we first introduce a retrieval error flag on a pixel-by-pixel basis. We use the retrieved aerosol and snow parameters to calculate the TOA reflectances [ in Equation 4] at all satellite channels using a P2R MLNN as explained in Section 2.2.2. Then these computed reflectances are compared with the satellite measured ones to determine the retrieval error for channel as:The retrieval error flag is defined by averaging the over all channels:where is the number of satellite channels used in the retrieval. for a typical SGLI image is shown in Figure 6. It can be seen in Figure 6 that is in general below 2%–5% for most channels except for the shortwave infrared SW3 (1,630 nm) channel, which is notably wider and may have caused a larger discrepancy between simulated and measured reflectances (Chen et al., 2017). One can also find wide, slanted stripes from the composition of image from different times of the day. Although different in their absolute values, the for different channels show similar patterns as pixels with high values appear to be affected by inadequate cloud screening or melting snow. Hence we consider an averaged value of 10% to be the criterion for a “good” retrieval for SGLI.

FIGURE 6

3 Data, validation and application

The main focus of this paper is the development and validation of an SPR algorithm for the SGLI sensor, based on the RTM-SciML method and applied to snow-covered land areas. The MLCM and retrieval MLNNs are all trained by dry snow models, therefore the SPR algorithm is expected to be less accurate in areas affected by surface melting. Using reflectance data in seven SGLI channels (see Table 3), we developed the SPR algorithm to retrieve the snow grain size in two layers, the snow impurity concentration, the aerosol optical depth, and three broadband surface albedo values (visible (VIS), near-infrared (NIR), and shortwave (SW)) using SGLI 1 km ground resolution L1B data. To demonstrate its merits, we applied the SPR algorithm to SGLI data obtained each July from 2018 to 2021 over the GrIS. Table 3 provides a summary of the SGLI data employed by the SPR algorithm, the output standard data products from the algorithm, and the validation data.

We used field measurements for validation of the SPR algorithm. These validation data included snow grain size data measured at the East GReenland Ice core Project (EGRIP) (Vallelonga et al., 2014) camp (75.6N, 36.0W, 2,708 m a.s.l.) in the central GrIS in the summer season of 2018, and surface albedo measured at the SIGMA-A site (78.1N, 67.6E, 1,490 m a.s.l.) (Aoki et al., 2014a) in the northwest GrIS from 2012 to 2020. The SIGMA-A site observations were conducted by members of this study (Nishimura et al., 2023). Surface albedo has been measured since July 2012 by a net radiometer (CNR-4, Kipp and Zonen, Netherlands) installed to the AWS deployed at the SIGMA-A site on the GrIS in the northwestern part of Greenland. Although the SIGMA-A station is located at the snow accumulation area, it may still be subject to surface melt and runoff during summer. To investigate the impact of surface melting on the results produced by SPR algorithm, we included the shortwave surface albedo data from the PROMICE AWS product (Fausto et al., 2021) at the EGP (accumulation) and KAN-L (ablation) field station as a comparison.

Snow samplings were performed for two snow layers of 0–2 cm and 2–10 cm at a clean area in the EGRIP camp on 30 June, 5 July, 10 July, and 15 July 2018. The snow grain size was measured for snow samples collected from a topmost layer (0–1 cm) and a thick layer (0–5 cm) using the IceCube instrument (A2 Photonic Sensors, France, Gallet et al. (2009)). This measurement was performed mainly synchronized with the overpasses of GCOM-C/SGLI. The parameter measured by the IceCube instrument is the specific surface area (SSA) ] of snow grains, from which the corresponding sphere-equivalent snow grain radius (SGR) [m] is calculated by , where is the ice density (0.917 g ) (Gallet et al., 2009). The SGR values were used to validate snow grain size retrieved from SGLI data (see Figure 7). Snow grain shapes in the top 30 cm layer consisted of faceted crystals, depth hoar, rounded grains and ice layers during the observation period. The averaged snow densities were 0.128, 0.213, 0.298 and 0.352 g for the layers of 0–3, 3–6, 6–9 and 9–30 cm, respectively. Snow surface temperature varied in the range from C to C. The snow samples were stored in dust-free plastic bags on site and transported to Japan. Mass concentrations of snow impurities were measured from a quartz fiber filter and a nuclepore filter on which the snow impurities were collected, in a laboratory of Meteorological Research Institute in Tsukuba, Japan using the methods employed by (Aoki et al., 2011; Aoki et al., 2014b) and Kuchiki et al. (2015). The impurity components are elemental carbon (EC) and mineral dust (MD). For the quartz fiber filters EC mass concentrations were measured with the Lab OC-EC Aerosol Analyzer (Sunset Laboratory Inc., USA) using the thermal optical reflectance method (Chow et al., 1993). The MD mass concentrations were determined by gravimetric measurements of the nuclepore filter before and after the filtration. We assumed that the component of EC is equal to black carbon (BC) for comparison with satellite-derived snow impurity mass concentrations in the present study. As satellite-inferred snow impurity mass concentrations contain the effects of light absorption by BC and MD, we define the BC-equivalent snow impurity mass concentration by (Tanikawa et al., 2015).where and are the mass absorption coefficient (in units of ) and the impurity mass concentration (in units of g ), respectively. We used the hydrophilic values of and for the broad visible wavelength range (Aoki et al., 2011), in which the channels to retrieve snow impurity mass concentrations are contained.

FIGURE 7

3.1 Validation with field measurements obtained on the GrIS

We start our validation effort by comparing the satellite retrieved snow grain size to field measured SGR values in the EGRIP site. Recall that in our snow model, the upper snow layer thickness was set to 1.5 cm with snow density of , while the bottom layer was set to 98.5 cm with snow density of . These settings may not be consistent with the vertical structure of the snow during field measurements. The satellite-retrieved snow grain size depends on the light penetration depth for each channel. We calculated the penetration depth (Li et al., 2001) for a typical day (16 July 2018) from our forward RTM, using the retrieved snow and aerosol parameters as input, and found that the penetration depth is about 10 cm at the 3 shortest wavelengths, about 3.5 cm at 869 nm, about 1.5 cm at 1,240 nm, and cm at 1.6 and 2.2 m. Field IceCube measurements were obtained between 0–5 cm snow depths. Since longer wavelength ( nm) channels are more important for grain size retrieval, we consider our model setting (0–1.5 cm for the top layer, cm for the bottom layer) to be comparable with the IceCube measurements.

To compare with the measurements, we averaged retrieved upper layer and lower layer snow grain sizes and compared the average value with the measured 0–5 cm value. Overall during these 18 days, the retrieved upper layer snow grain sizes were between 30 and 100 m and the lower layer snow grain sizes were between 50 and 200 m.

The results are shown in Figure 7 with data selected based on matching satellite measurements within 1 h time difference and satellite viewing zenith angle smaller than 30°. In total, 16 averaged snow grain size samples from 0–5 cm were selected based on these criteria. Overall, the retrieved grain size shows a positive correlation with the measurements , but the limited sample size hampers statistical significance. Consequently, we shifted our validation efforts to in-situ snow albedo measurements, which are relatively abundant.

At the SIGMA-A site on the GrIS in the northwestern part of Greenland, albedo values were measured once a minute and 1-h averaged values were recorded. These hourly albedo values are compared with the nearest satellite-derived values ( hour) in Figure 8. The data points are color-coded by the associated averaged values from all satellite channels employed. The satellite-derived blue-sky albedo values generally agree with in-situ measurements for both SW and NIR albedo with relative errors smaller than 10%, especially for the NIR albedo which may indicate a higher accuracy in grain size retrieval. The relatively lower correlation coefficient of SW albedo may be caused by the difficulties in snow impurity retrieval, which could significantly affect the visible albedo of snow and create outliers that would significantly impact the results for a small sample size . To confirm this conjecture, we re-trained the algorithm with a dataset that has impurity concentration fixed to a level that is close to that determined by the field measurements (about ) and repeated the retrieval. The updated results displayed in Figure 9 show that fixing the impurity concentration indeed reduced the number of outliers in SW albedo retrieval whereas the overall correlation coefficients in both SW and NIR albedo are slightly decreased. Since the NIR albedo is mainly affected by the snow grain size this finding means that the accuracy of grain size retrievals are getting worse after fixing the impurity concentration. We also noticed that the data points with large values are often located at the bottom left corner of the albedo plots, which indicates larger grain size and tendency of melting. In summary, the results show that quantifying the presence of impurity and its interaction with other parameters (especially the lower layer snow grain size) is still a challenging problem that requires further attention. The albedo enhancement under cloudy (compared to clear) sky conditions due to the change of the spectral distribution of downward irradiance over snow surfaces (Liljequist, 1956; Aoki et al., 1999) can also increase the uncertainty as the field measurements and satellite overpasses are usually not exactly simultaneous.

FIGURE 8

FIGURE 9

To further investigate the relationship between and surface melting/bare ice, we compared our retrieval results with the surface shortwave albedo data from PROMICE automatic weather station data. In the snow accumulation area of the GrIS we chose the EGP field station as this site is the only one in the PROMICE dataset that is expected to be minimally affected by surface melting as it is located at 2,660 m above the sea level and close to the center of Greenland. On the other hand, we chose the KAN-L station located in the snow ablation area to examine the performance of the algorithm when applied to mostly melting snow/bare ice surfaces. We use the same time period as that adopted for our SIGMA-A site comparison and the results are shown in Figure 10. At the EGP station, the retrieved albedo data are tightly packed between 0.8–0.9 with relatively low averaged values (with all values below 20%). The lower correlation coefficient can be explained by the “restriction of range” effect, in which the measurements (or predictions) all fall within a narrow interval so that even small amounts of scatter/noise can overshadow the true relationship and reduce the correlation. There seems to be a slight bias of satellite-retrieved albedo compared to the field measurements of about 0.03–0.05 in Figure 10. This bias could be due to the platform obstruction effect as described by (Kokhanovsky et al., 2020) which we have not attempted to correct for. We also noticed a slight difference in shortwave albedo wavelength range (0.3–2.8 for the SIGMA-A site and 0.3–2.5 for the EGP site) which we consider to be minor but leaves room for further investigation. At the KAN-L station, however, the retrieved albedo show large dispersion with significantly increased RMSE and MAPE, which indicates downgraded retrieval quality from the algorithm. The average is the largest at this site (28.36%) compared to the other sites. Detailed statistics of shortwave albedo validation are listed in Table 4.

FIGURE 10

3.2 Application to SGLI data over Greenland

SGLI data have been available since early 2018. In Figures 11–15, we show applications of our MLCM and SPR algorithms to SGLI data obtained over the GrIS in July of 2018, 2019, 2020, and 2021. From these data we can infer not only how the snow properties evolved during the month of July for each year, but also how the snow surface of the GrIS changed from July 2018 to July 2021.

FIGURE 11

The retrieved snow grain size in the top and bottom layers of the two-layer snow configuration are shown in Figures 11, 12. The three columns from left show the 10-day average snow retrieval parameters for July 1–10, July 11–20, and July 21–30. The last column displays histograms of the 10-day averages and the overall July-averaged snow retrieval parameters, with their mode and mean values indicated in the legend. The retrieved values are smooth and stable and in basic agreement with ground truth data as demonstrated in Section 3.1. Over most of the GrIS, the upper layer snow grain size is smaller than about 100 m, and the lower layer snow grain size is a few hundreds of m. Although we already excluded pixels located at the edge of the GrIS identified as “non-snow’’ (by the MLCM algorithm), we still find pixels there with large grain size (m), which may be due to melting snow. At the beginning of July, the snow surface appears to be quite homogeneous across the GrIS. But later in July, a significant amount of pixels in the southern part of the GrIS has very large snow grain size (m), indicating that the snow has started to melt.

FIGURE 12

Among these 4 years, grain size values for July 2019 (especially for the bottom layer which is relatively stable under snow precipitation events) were found to be significantly higher with average bottom layer grain size being higher by about 54%, 25% and 27% than for the other 3 years (2018, 2020 and 2021) over the whole GrIS. This finding is consistent with Figure 15 showing that the surface temperature was also higher, implying that the GrIS was warmer in July 2019 than in July 2018, consistent with the results reported by Tedesco and Fettweis (2020) and Sasgen et al. (2020). A detailed description of the snow surface temperature product is provided elsewhere (Stamnes et al., 2007).

Consequently, both the blue sky and black sky surface albedo show lowest values for the year of 2019 compared to the other years (Figures 13, 14), with the mean blue sky and black sky albedo values for July of 2019 lowered by more than 0.02 compared to the other years. The albedo plots in Figures 13, 14, the temperature plots in Figure 15, and the plots of July 2019 in Figure 16 indicate a clear correlation between areas with lower albedo values, higher surface temperature and . These results suggest that can be a good indicator of surface melting.

FIGURE 13

FIGURE 14

FIGURE 15

FIGURE 16

4 Summary and conclusion

4.1 Uncertainty estimates and future work

References

• 1

  AmariS.ParkH.FukumizuK. (2000). Adaptive method of realizing natural gradient learning for multilayer perceptrons. Neural Comput.12, 1399–1409. 10.1162/089976600300015420

  • CrossRef
  • Google Scholar

• 2

  AndersonG. P.CloughS. A.KneizysF. X.ChetwyndJ. H.ShettleE. P. (1986). AFGL atmospheric constituent profiles (0-120km), AFGL-TR-86-0110 (OPI) Optical Physics Division, Air Force Geophysics Laboratory. MA 01736: Hanscom AFB.

  • Google Scholar

• 3

  AokiT.AokiT.FukaboriM.UchiyamaA. (1999). Numerical simulation of the atmospheric effects on snow albedo with a multiple scattering radiative transfer model for the atmosphere-snow system. J. Meteorological Soc. Jpn. Ser. II77, 595–614. 10.2151/jmsj1965.77.2_595

  • CrossRef
  • Google Scholar

• 4

  AokiT.KuchikiK.NiwanoM.KodamaY.HosakaM.TanakaT. (2011). Physically based snow albedo model for calculating broadband albedos and the solar heating profile in snowpack for general circulation models. J. Geophys. Res. Atmos.116, D11114. 10.1029/2010jd015507

  • CrossRef
  • Google Scholar

• 5

  AokiT.MatobaS.UetakeJ.TakeuchiN.MotoyamaH. (2014a). Field activities of the snow impurity and glacial microbe effects on abrupt warming in the arctic (SIGMA) project in Greenland in 2011-2013. Bull. Glaciol. Res.32, 3–20. 10.5331/bgr.32.3

  • CrossRef
  • Google Scholar

• 6

  AokiT.MatobaS.YamaguchiS.TanikawaT.NiwanoM.KuchikiK.et al (2014b). Light-absorbing snow impurity concentrations measured on Northwest Greenland ice sheet in 2011 and 2012. Bull. Glaciol. Res.32, 21–31. 10.5331/bgr.32.21

  • CrossRef
  • Google Scholar

• 7

  BairE. H.StillingerT.DozierJ. (2021). Snow property inversion from remote sensing (SPIReS): A generalized multispectral unmixing approach with examples from MODIS and Landsat 8 oli. IEEE Trans. Geoscience Remote Sens.59, 7270–7284. 10.1109/TGRS.2020.3040328

  • CrossRef
  • Google Scholar

• 8

  BakerN.AlexanderF.BremerT.HagbergA.KevrekidisY.NajmH.et al (2019). Workshop report on basic research needs for scientific machine learning: core technologies for artificial intelligence. Tech. rep., U.S. Dep. Energy Office Sci. Tech. Inf.10.2172/1478744

  • CrossRef
  • Google Scholar

• 9

  ChenN.LiW.GatebeC.TanikawaT.HoriM.ShimadaR.et al (2018). New neural network cloud mask algorithm based on radiative transfer simulations. Remote Sens. Environ.219, 62–71. 10.1016/j.rse.2018.09.029

  • CrossRef
  • Google Scholar

• 10

  ChenN.LiW.TanikawaT.HoriM.AokiT.StamnesK. (2014). Cloud mask over snow/ice covered areas for the GCOM-C1/SGLI cryosphere mission: validations over Greenland. J. Geophys. Res. Atmos.119. 10.1002/2014JD022017

  • CrossRef
  • Google Scholar

• 11

  ChenN.LiW.TanikawaT.HoriM.ShimadaR.AokiT.et al (2017). Fast yet accurate computation of radiances in shortwave infrared satellite remote sensing channels. Opt. express25, A649–A664. 10.1364/oe.25.00a649

  • CrossRef
  • Google Scholar

• 12

  ChenS.BillingsS. A.GrantP. M. (1990). Non-linear system identification using neural networks. Int. J. Control51, 1191–1214. 10.1080/00207179008934126

  • CrossRef
  • Google Scholar

• 13

  ChowJ. C.WatsonJ. G.PritchettL. C.PiersonW. R.FrazierC. A.PurcellR. G. (1993). The DRI thermal/optical reflectance carbon analysis system: description, evaluation and applications in US air quality studies. Atmos. Environ. Part A. General Top.27, 1185–1201. 10.1016/0960-1686(93)90245-t

  • CrossRef
  • Google Scholar

• 14

  D’AlimonteD.ZibordiG. (2003). Phytoplankton determination in an optically complex coastal region using a multilayer perceptron neural network. IEEE Trans. Geoscience Remote Sens.41, 2861–2868. 10.1109/TGRS.2003.817682

  • CrossRef
  • Google Scholar

• 15

  D’AlimonteD.ZibordiG.BerthonJ. F. (2004). Determination of cdom and nppm absorption coefficient spectra from coastal water remote sensing reflectance. IEEE Trans. Geoscience Remote Sens.42, 1770–1777. 10.1109/TGRS.2004.831444

  • CrossRef
  • Google Scholar

• 16

  D’AlimonteD.ZibordiG.BerthonJ.-F.CanutiE.KajiyamaT. (2012). Performance and applicability of bio-optical algorithms in different European seas. Remote Sens. Environ.124, 402–412. 10.1016/j.rse.2012.05.022

  • CrossRef
  • Google Scholar

• 17

  DeltzA. J.KuenzerC.GessnerU.DechS. (2012). Remote sensing of snow - a review of available methods. Int. J. Remote Sens.33 (13), 4094–4134. 10.1080/01431161.2011.640964

  • CrossRef
  • Google Scholar

• 18

  FanY.LiW.ChenN.AhnJ.-H.ParkY.-J.KratzerS.et al (2021). OC-SMART: a machine learning based data analysis platform for satellite ocean color sensors. Remote Sens. Environ.253, 112236. 10.1016/j.rse.2020.112236

  • CrossRef
  • Google Scholar

• 19

  FanY.LiW.GatebeC. K.JametC.ZibordiG.SchroederT.et al (2017). Atmospheric correction over coastal waters using multilayer neural networks. Remote Sens. Environ.199, 218–240. 10.1016/j.rse.2017.07.016

  • CrossRef
  • Google Scholar

• 20

  FaustoR. S.van AsD.MankoffK. D.VandecruxB.CitterioM.AhlstrømA. P.et al (2021). Programme for Monitoring of the Greenland Ice Sheet (PROMICE) automatic weather station data. Earth Syst. Sci. Data13, 3819–3845. 10.5194/essd-13-3819-2021

  • CrossRef
  • Google Scholar

• 21

  FreiA.TedescoM.LeeS.FosterJ.HallD.KellyR.et al (2012). A review of global satellite-derived snow products. Adv. Space Res.50, 1007–1029. 10.1016/j.asr.2011.12.021

  • CrossRef
  • Google Scholar

• 22

  GalletJ.-C.DomineF.ZenderC.PicardG. (2009). Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm. Cryosphere3, 167–182. 10.5194/tc-3-167-2009

  • CrossRef
  • Google Scholar

• 23

  Harris StuartR.FaberA.-K.WahlS.HörholdM.KipfstuhlS.VasskogK.et al (2023). Exploring the role of snow metamorphism on the isotopic composition of the surface snow at EastGRIP. Cryosphere17, 1185–1204. 10.5194/tc-17-1185-2023

  • CrossRef
  • Google Scholar

• 24

  HowP.LundM.AhlstrømA.AndersenS.BoxJ. E.CitterioM.et al (2022). PROMICE and GC-Net automated weather station data in Greenland. 10.22008/FK2/IW73UU

  • CrossRef
  • Google Scholar

• 25

  IshimotoH.AdachiS.YamaguchiS.TanikawaT.AokiT.MasudaK. (2018). Snow particles extracted from X-ray computed microtomography imagery and their single-scattering properties. J. Quantitative Spectrosc. Radiat. Transf.209, 113–128. 10.1016/j.jqsrt.2018.01.021

  • CrossRef
  • Google Scholar

• 26

  IshimotoH.MasudaK.ManoY.OrikasaN.UchiyamaA. (2012). Irregularly shaped ice aggregates in optical modeling of convectively generated ice clouds. J. Quantitative Spectrosc. Radiat. Transf.113, 632–643. 10.1016/j.jqsrt.2012.01.017

  • CrossRef
  • Google Scholar

• 27

  KingmaD. P.BaJ. (2017). Adam: a method for stochastic optimization. 10.48550/arXiv.1412.6980

  • CrossRef
  • Google Scholar

• 28

  KokhanovskyA.BoxJ. E.VandecruxB.MankoffK. D.LamareM.SmirnovA.et al (2020). The determination of snow albedo from satellite measurements using fast atmospheric correction technique. Remote Sens.12, 234. 10.3390/rs12020234

  • CrossRef
  • Google Scholar

• 29

  KokhanovskyA.LamareM.DanneO.BrockmannC.DumontM.PicardG.et al (2019). Retrieval of snow properties from the sentinel-3 ocean and land colour instrument. Remote Sens.11 (18), 2280–2323. 10.3390/rs11192280

  • CrossRef
  • Google Scholar

• 30

  KuchikiK.AokiT.NiwanoM.MatobaS.KodamaY.AdachiK. (2015). Elemental carbon, organic carbon, and dust concentrations in snow measured with thermal optical and gravimetric methods: variations during the 2007–2013 winters at sapporo, Japan. J. Geophys. Res. Atmos.120, 868–882. 10.1002/2014jd022144

  • CrossRef
  • Google Scholar

• 31

  LiW.StamnesK.ChenB.XiongX. (2001). Snow grain size retrieved from near-infrared radiances at multiple wavelengths. Geophys. Res. Lett.28 (1), 1699–1702. 10.1029/2000gl011641

  • CrossRef
  • Google Scholar

• 32

  LiljequistG. H. (1956). Energy exchange of an Antarctic snow-field: short-wave radiation. Oslo Nor. Polarinst2.

  • Google Scholar

• 33

  LyapustinA.TedescoM.WangY.AokiT.HoriM.KokhanovskyA. (2009). Retrieval of snow grain size over Greenland from MODIS. Remote Sens. Environ.113, 1976–1987. 10.1016/j.rse.2009.05.008

  • CrossRef
  • Google Scholar

• 34

  NeumaierA. (1998). Solving ill-conditioned and singular linear systems: a tutorial on regularization. SIAM Rev.40, 636–666. 10.1137/S0036144597321909

  • CrossRef
  • Google Scholar

• 35

  NishimuraM.AokiT.NiwanoM.MatobaS.TanikawaT.YamasakiT.et al (2023). Quality-controlled meteorological datasets from SIGMA automatic weather stations in northwest Greenland, 2012–2020. Earth Syst. Sci. Data15, 5207–5226. 10.5194/essd-15-5207-2023

  • CrossRef
  • Google Scholar

• 36

  PainterT.RittgerK.McKenzieP. S.DavisR.DozierJ. (2009). Retrieval of subpixel snow covered area, grain size, and albedo from MODIS. Remote Sens. Environ.113, 868–879. 10.1016/j.rse.2009.01.001

  • CrossRef
  • Google Scholar

• 37

  SasgenI.WoutersB.GardnerA. S.KingM. D.TedescoM.LandererF. W.et al (2020). Return to rapid ice loss in Greenland and record loss in 2019 detected by the GRACE-FO satellites. Commun. Earth and Environ.1, 8. 10.1038/s43247-020-0010-1

  • CrossRef
  • Google Scholar

• 38

  SirgueyP.MathieuR.ArnaudY. (2009). Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment. Remote Sens. Environ.113, 160–181. 10.1016/j.rse.2008.09.008

  • CrossRef
  • Google Scholar

• 39

  StamnesK.HamreB.StamnesS.ChenN.FanY.LiW.et al (2018). Progress in forward-inverse modeling based on radiative transfer tools for coupled atmosphere-snow/ice-ocean systems: a review and description of the AccuRT model. Appl. Sci.8, 2682. 10.3390/app8122682

  • CrossRef
  • Google Scholar

• 40

  StamnesK.StamnesJ. J. (2016). Radiative transfer in coupled environmental systems: an introduction to forward and inverse modeling. John Wiley and Sons.

  • Google Scholar

• 41

  StamnesK.WeiL.EideH.AokiT.HoriM.StorvoldR. (2007). ADEOS-II/GLI snow/ice products: Part I - Scientific basis. Remote Sens. Environ.111, 258–273. 10.1016/j.rse.2007.03.023

  • CrossRef
  • Google Scholar

• 42

  Steen-LarsenH. C.HörholdM.KipfstuhlS.FaberA.-K.FreitagJ.HughesA. G.et al (2022). 10 daily surface measurements over 90m transect, SSA, Density and Accumulation, from EastGRIP summer (May-August) of 2016-2019. 10.1594/PANGAEA.946763

  • CrossRef
  • Google Scholar

• 43

  StillingerT.RobertsD.CollarN.DozierJ. (2019). Cloud masking for Landsat 8 and MODIS terra over snow-covered terrain: error analysis and spectral similarity between snow and cloud. Water Resour. Res.55, 6169–6184. 10.1029/2019WR024932

  • CrossRef
  • Google Scholar

• 44

  TanikawaT.KuchikiK.AokiT.IshimotoH.HachikuboA.NiwanoM.et al (2020). Effects of snow grain shape and mixing state of snow impurity on retrieval of snow physical parameters from ground-based optical instrument. J. Geophys. Res. Atmos.125, e2019JD031858. 10.1029/2019JD031858

  • CrossRef
  • Google Scholar

• 45

  TanikawaT.LiW.KuchikiK.AokiT.HoriM.StamnesK. (2015). Retrieval of snow physical parameters by neural networks and optimal estimation: case study for ground-based spectral radiometer system. Opt. express23, A1442–A1462. 10.1364/oe.23.0a1442

  • CrossRef
  • Google Scholar

• 46

  TedescoM. (2015). Remote sensing of the cryosphere. New York, NY, USA: Wiley.

  • Google Scholar

• 47

  TedescoM.FettweisX. (2020). Unprecedented atmospheric conditions (1948–2019) drive the 2019 exceptional melting season over the Greenland ice sheet. Cryosphere14, 1209–1223. 10.5194/tc-14-1209-2020

  • CrossRef
  • Google Scholar

• 48

  VallelongaP.ChristiansonK.AlleyR. B.AnandakrishnanS.ChristianJ. E. M.Dahl-JensenD.et al (2014). Initial results from geophysical surveys and shallow coring of the Northeast Greenland Ice Stream (NEGIS). Cryosphere8, 1275–1287. 10.5194/tc-8-1275-2014

  • CrossRef
  • Google Scholar

• 49

  VandecruxB.BoxJ. E.AhlstrømA. P.AndersenS. B.BayouN.ColganW. T.et al (2023). The historical Greenland Climate Network (GC-Net) curated and augmented level-1 dataset. Earth Syst. Sci. Data15, 5467–5489. 10.5194/essd-15-5467-2023

  • CrossRef
  • Google Scholar

• 50

  VandecruxB.BoxJ. E.WehrléA.KokhanovskyA. A.PicardG.NiwanoM.et al (2022). The determination of the snow optical grain diameter and snowmelt area on the Greenland ice sheet using spaceborne optical observations. Remote Sens.14, 932. 10.3390/rs14040932

  • CrossRef
  • Google Scholar

• 51

  WarrenS. G. (2013). Can black carbon in snow be detected by remote sensing?J. Geophys. Res. Atmos.118, 779–786. 10.1029/2012jd018476

  • CrossRef
  • Google Scholar

• 52

  ZegeE.KatsevI.MalinkaA.PrikhachA.HeygsterG.WiebeH. (2011). Algorithm for retrieval of the effective snow grain size and pollution amount from satellite measurements. Remote Sens. Environ.115, 2674–2685. 10.1016/j.rse.2011.06.001

  • CrossRef
  • Google Scholar

• 53

  ZhouY.LiW.ChenN.FanY.StamnesK. (2023). A sensor-agnostic albedo retrieval method for realistic sea ice surfaces: model and validation. Cryosphere17, 1053–1087. 10.5194/tc-17-1053-2023

  • CrossRef
  • Google Scholar