## Abstract

A new retrieval algorithm for estimation of snow grain size and impurity concentration from spectral radiation data is developed for remote sensing applications. A radiative transfer (RT) model for the coupled atmosphere-snow system is used as a forward model. This model simulates spectral radiant quantities for visible and near-infrared channels. The forward RT calculation is, however, the most time-consuming part of the forward-inverse modeling. Therefore, we replaced it with a neural network (NN) function for fast computation of radiances and Jacobians. The retrieval scheme is based on an optimal estimation method with *a priori* constraints. The NN function was also employed to obtain an accurate first guess in the retrieval scheme. Validation with simulation data shows that a combination of NN techniques and optimal estimation method can provide more accurate retrievals than by using only NN techniques. In addition, validation with *in-situ* measurements conducted by using ground-based spectral radiometer system shows that comparison between retrieved snow parameters with *in-situ* measurements is acceptable with satisfactory accuracy. The algorithm provides simultaneous, accurate and fast retrieval of the snow properties. The algorithm presented here is useful for airborne/satellite remote sensing.

© 2015 Optical Society of America

## 1. Introduction

Snow/ice albedo is a fundamental component determining the radiation balance of the Earth-atmosphere system for climate studies, environmental physics and snow hydrology [1]. Snow generally has a high albedo that minimizes the absorption of shortwave radiation impinging on the snow surface, and thereby prevents the snow from melting. Snowpacks are an ensemble of aggregated ice microphysical particles, and the optical properties of these snow particles are governed by ice microphysical properties. Therefore, the most fundamental snow variables are ice microphysical properties that can influence the radiative transfer, snow grain size growth, and hydrological processes in snow [2]. Light absorbing snow impurity is also a fundamental variable which impacts the shortwave radiation impinging on the snow surface. The reason is that the impurity commonly consists of highly absorptive material which causes albedo reduction even when a small amount of impurity is embedded in a snowpack [3]. To fully understand these influences we must have reliable information on snow microphysics and impurity absorption properties, including effective snow grain size and mass concentration of snow impurity. Besides, we need to monitor these snow parameters on spatial and temporal scales.

Recent advances in multi-channel or hyper-spectral optical instruments have made quantitative monitoring possible. Retrievals of snow grain size using airborne optical instruments have been performed by several research groups [4–8]. High spatial resolution images such as the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data [9] and Airborne Multi-Spectral Scanner (AMSS) data [10] provide accurate snow information related to the water resources in basin as well as mountain areas. Such imagery can be also used to evaluate the accuracy of the retrieved snow grain size by comparing retrieved values with *in-situ* field measurement data.

Satellite remote sensing adaptation from multispectral to hyperspectral sensors has been tested through various studies. The grain size before and after snow fall over mountain areas was mapped using two sets of the Moderate Resolution Imaging Spectroradiometer (MODIS) data [11]. One set was obtained from the morning pass of the Terra satellite and the other from the afternoon pass of the Aqua satellite. They detected new snow from two images, with the retrieval error of 5% to 40%. For global scale monitoring, the snow grain size over Antarctica was monitored using a normalized difference approach between MODIS band 1 and band 2 radiances [12]. The snow grain size over Antarctica, obtained from MODIS data, showed similar but large seasonal variation in all years between 2000 and 2005 [13]. A snow/ice retrieval algorithm for application to Global Imager (GLI) data was developed [14], and used to infer snow grain size distributions and impurity concentrations over the northern hemisphere from GLI data [6]. The latter study also used snow grain size and temperature distributions retrieved from GLI data to construct a map of snow metamorphism potential, which is a useful indicator of environmental change. Monitoring of snow grain size was used as an indicator of snowmelt, and to identify the occurrence of snowmelt over large areas in Greenland in 2004 [15]. Another application of the snow grain size retrievals demonstrates the potential for snow fall detection. Time series of monitoring the snow grain size made it possible to detect snow fall events after sudden grain size decreases [16].

The remote sensing of light absorbing contaminants in snow is generally a challenging task because it is difficult to discriminate between absorption in the atmosphere and in the snow. So this area of research is a very challenging one in Earth observation by the remote sensing. However, the Landsat Thematic Mapper (TM) can be used to identify areas where the snow impurity concentration is high [17]. The mass fraction of snow impurities were retrieved from AMSS data, and a single experiment comparing mass fractions inferred from AMSS data and from field measurements was conducted [5]. The comparison, assuming mineral dust as the snow impurities, revealed good agreement between AMSS-retrieved and *in-situ* measured values. The monitoring of snow impurity concentrations from GLI data indicated reasonable results by comparison with *in-situ* data obtained in Hokkaido, Japan and Alaska [18].

It is well recognized that both snow physical and radiative processes vary significantly with the vertical position in the snow [19]. This recognition is the motivation for the new algorithms to be developed for the Global Change Observation Mission-Climate/Second-generation Global Imager (GCOM-C/SGLI) which is scheduled to launch in the 2016 fiscal year by the Japan Aerospace Exploration Agency (JAXA) [20]. Near-infrared channels at wavelengths such as *λ* = 0.86,1.24 *μ*m and 1.64 *μ*m can be used to retrieve vertical information about snow grain size [21]. Because longer wavelengths carry information about shallower layers of the snow, the results derived from several near-infrared wavelengths give information about the vertical distribution of the snow grain size. This finding indicated that several near-infrared channels available on the AMSS or GLI sensor can be used to retrieve information about the grain size distribution [5,18]. However, most previous retrieval studies employed a one-layer snow model which is unrealistic because real snow has significant vertical variation from the top to the bottom of the snow layer. Thus, there could be a difficulty in the interpretation because the different channels provide a different grain size from the one-layer snow model. The radiance measured in visible and near infrared channels originates from different snow depths, and the snow grain size varies with snow depth. In order to retrieve snow parameters accurately, we need a multi-channel and a multi-angle optical instrument in conjunction with a multi-layer snow model to take into account the vertical variation of the snow pack.

In general, a multi-channel and multi-angle (ground/satellite/airborne) observation system provides a significant amount of information about the target material [22]. Hence, use of such a system in conjunction with a multi-layer snow model is expected to improve the accuracy of retrieved snow parameters. Additionally we also expect to obtain detailed information about the target material, including parameters that have not been retrieved before, leading to new retrieval products. However, in order to achieve this goal, the techniques used to retrieve system parameters (for example, snow grain size and impurity concentration embedded in snow or oceanic/atmospheric parameters, and so on) from the ground/airborne/satellite-measured data may be need to be modified, adapted, and enhanced. Well-known techniques to retrieve information about environmental parameters employ a radiative transfer (RT) model as a forward model and look-up-tables (LUTs) as an inversion model. This approach has the advantage of being simple and computationally fast. But increased information provided by the instruments makes this approach difficult to use because it cannot easily be generalized to solve forward/inverse problems which in general involve complicated, non-linear relationships between the observed radiances and the retrieval parameters. Instead this situation suggests that we should explore employing more powerful statistical/analytical optimizations in the data/image processing. The neural network (NN) technique is a very promising mathematical tool to solve the forward/inverse problems in remote sensing accurately [23]. There are many practical advantages. If the forward/inverse model can each be replaced by a simple neural network function, the NN technique may provide increased computational speed, improved accuracy, and robustness. There are many examples of the use of NN techniques to solve the inverse problem encountered in the remote sensing applications [24–27]. One example is the retrievals of total ozone column amounts and cloud effects from ultraviolet irradiance instruments by using the NN method [27]. Comparison of methodology for the retrievals shows that the NN method retrieved more valid results than LUT method. The NN method can take various environmental factors into account, implying that the NN method has a high flexibility and robustness in the retrievals compared to the LUT method.

The NN technique can also provide the entire Jacobian matrix with practically no additional computational effort. Access to the Jacobian matrix is a big advantage in forward/inverse modeling because it is needed in non-linear optimal estimation methods used to retrieve physical parameters. We expect a combination of NN techniques and optimal estimation to give more accurate retrievals than by using only NN techniques or the traditional LUT approach.

The purpose of this paper is to discuss a new forward/inverse modeling approach that relies on the use of NN techniques combined with the non-linear optimal estimation method. In order to achieve this goal, we show a validation of the NN function used as a forward/inverse model. Then we confirm the non-linear optimal estimation method by using simulated data as a whole retrieval algorithm. Finally, we construct a simple retrieval procedure based on a 2-layer snow model and apply it to a ground-based spectral data set to retrieve snow parameters of topmost-and subsurface-layer snow grain sizes and light absorbing snow impurity concentrations which are again the fundamental parameters controlling the spectral albedo of snow.

## 2. Retrieval algorithm

Our proposed algorithm to retrieve the physical parameters from the observed data consists of three key parts or tasks, *i.e.* (i) compilation of a data set of corresponding pairs of (snow) physical parameters with geometry and spectral radiances (or albedo) obtained by the RT model for the coupled atmosphere-snow system, (ii) the NN training of forward and inverse models, respectively, *i.e.* determination of the coefficients of the NN function, (iii) installation of forward/inverse NN functions into the retrieval scheme based on a non-linear, iterative optimal estimation method with *a priori* constraints (Fig. 1). We briefly explain each step in the following section. It should be noted that in general, the retrieval of snow physical parameters from remote sensing data is an ill-posed problem. So, we need the regularization using some additional *a priori* information to solve the ill-posed problem. Besides, there are additional problems to be addressed before the retrievals such as the influence of the atmosphere and the restriction of the spectral information. We will also discuss these problems in the following section.

#### 2.1. Radiative transfer model for the coupled atmosphere-snow system

The first task is to make the data set
${\{\mathbf{P},\mathbf{R}\}}_{s}{}_{=1,\dots ,{N}_{\mathrm{set}}}$ corresponding to the radiant fluxes **R** and the retrieval parameters **P** for training of the NN function. For this purpose, we used the discrete-ordinate radiative transfer (DISORT) model for the atmosphere-snow system as the forward model [28,29]. A multi-layer snow model was assumed to obtain accurate retrievals of snow parameters. Since the snow grain size at the snow surface changes significantly on short time scales [30, 31], we set the snow layer structure as follows: a topmost layer of snow depth *d* = 0− 1 cm and a subsurface layer of *d* = 1 − 100 cm. In the forward model of snow, input parameters are snow physical parameters (the topmost- and the subsurface-layer snow grain size, *r*_{top}, *r*_{sub} and the light absorbing snow impurity concentration, *c*_{imp}) and geometry parameters (solar zenith angle, *θ*_{0}), and the output parameters are the radiant fluxes. Besides, the light absorbing snow impurity concentration *c*_{imp} was assumed to be homogenious in the whole layer. Scattering properties can be computed from the classical light-scattering theory, the Mie theory. A smooth phase function such as the Henyey-Greenstein function is applied. For the definition of snow grain size, we used the effective snow grain radius, or “area-weighted mean radius,” instead of geometric radius since spherical droplets scatter the amount of light in proportion to their cross-sectional area [32]. A log-normal size distribution of snow grain sizes was used to calculate scattering properties [33]. The wavelength dependent refractive indices of ice were used in the calculation [34]. The snow impurity was assumed to be black carbon model [35], externally mixed with the snow particles. For the atmospheric model, since analysis is confined to the Hokkaido in Japan, we adopted a standard mid-latitude winter model and divided it into 24 layers to determine the optical properties and adequately take into account gaseous absorption by H_{2}O, O_{3}, SO_{2} and the uniformly mixed gases. We used continental average aerosol and a cumulus (continental, clean) cloud model [35] and set the aerosol vertical distribution from 0 to 2 km, and the cloud from 2 to 5 km. According to Hori et al. [6], a sensitivity analysis of the atmospheric aerosol model on the snow parameter retrievals showed that the snow grain sizes were less sensitive to the selection of aerosol model except for the absorptive aerosol model cases (e.g., urban aerosol with 23 km visibility) or unrealistically heavy aerosol loading (e.g., rural aerosolwith 5 km visibility). The result implys that the snow parameters retrievals are not considered to be seriously in error for our test site without any unrealistically heavy/absorptive aerosol loading. For the cloudness, the sensitivity analysis showed that the retrieval algorithm does not support the case for the thick optical thickness due to the less sensitivity of spectral information on the snow. Hence, the retrieval algorithm performs under the almost clear sky.

#### 2.2. Neural network training

In order to retrieve multiple-layer snow parameters simultaneously, we need to use an inversion method. The traditional LUT interpolation method is one option but this method is not suitable for simultaneous and accurate retrieval of multiple parameters [27]. Instead, we used an iterative optimal estimation method with loosely constrained *a priori* information [36, 37]. But this approach is very time consuming because we need to call the forward RT model repeatedly to compute radiant fluxes and Jacobians required in the non-linear iterative optimal estimation method described in the next section. In order to solve this problem, we replaced the forward RT model with a forward radial basis function neural network (RBF-NN) to establish a relationship between the retrieval parameters and the spectral radiant fluxes. The RBF-NN functions consist of two layers with a radial basis layer of hidden neurons, and a linear layer of output neurons.

### 2.2.1. Parameters to spectral radiant fluxes training (forward NN training)

We start by training a RBF-NN to provide a relationship between the retrieval parameters **P** and the radiant fluxes **R** and denote this forward NN function **R** = **m**_{NN}(**P**). We use this analytic relationship to derive an analytic expression for the Jacobians **J**. Both the radiant fluxes and the Jacobians are required in the non-linear optimal estimation described in Section 2.3. The model **m** is used to compute a large set of
${\{\mathbf{P},\mathbf{R}\}}_{s}{}_{=1,\dots ,{N}_{\mathrm{set}}}$. This set is used for training/testing of the forward NN. The input data **P** are snow and atmospheric parameters (*i.e.* snow grain size, snow impurity, cloud and aerosol optical depth) and solar zenith angle. The output data **R** are the radiant fluxes. The spectral radiant flux *R _{i}* in channel

*i*obtained from the RBF-NN is given as:

*N*and

*N*are the total number of neurons and input retrieval parameters, respectively. The

_{p}*P*’s are the input parameters, and

_{k}*a*,

_{i j}*b*,

*c*and

_{jk}*d*are the optimized coefficients.

_{i}The Jacobians (partial derivatives) **J** are also required for retrieval by the non-linear optimal estimation. These Jacobians can be calculated by analytical differentiation as follows

Note that we use the same optimized coefficients of *a _{i j}*,

*b*,

*c*obtained in Eq. (1) to calculate the Jacobians

_{jk}**J**.

### 2.2.2. Spectral radiant fluxes to parameters training (inverse NN training)

We train the RBF-NN as the inverse NN function
$\mathbf{P}={\mathbf{m}}_{\mathrm{NN}}^{-1}(\mathbf{R})$ in which the spectral radiant fluxes, computed for a large set of aerosol and cloud optical depths, and solar zenith angles, are used in the accurate forward RT model to generate input data **R**, whereas the snow parameters are the output data **P**. The RBF-NN function is given by

*N*is the total number of neurons, and

*N*is the number of input parameters. The

_{r}*R*’s are the input parameters, and

_{i}*e*,

_{k j}*f*,

*g*and

_{ji}*h*are the optimized coefficients. The snow parameters

_{k}**P**can be obtained directly from the Eq. (3). This inverse training is very useful, because it provides approximate retrieval parameters that can be used as starting values in the optimal estimation. Since optimal estimation methods rely on Gauss-Newton type inversion, they need starting values that are reasonably close to the solution in order to converge. For routine observations, previous measurement data or retrieved parameters are often used as initial parameters in the iterative inversion process. However, this approach may not be suitable because the previous measurement data could not be used to start a retrieval after a snowfall. Thus, the inverse NN function is progressive idea.

#### 2.3. Non-linear optimal estimation

We use cost-function minimization techniques appropriate to non-linear iterative spectral fitting. The retrieval parameters are assembled into a state vector **x** = [*x*_{1}, *x*_{2},*…*, *x _{n}*]

*, and the measured radiant fluxes into another vector*

^{T}**Y**

_{meas}= [

*Y*

_{1},

*Y*

_{2},

*…*,

*Y*]

_{m}*. In optimal estimation, the update of the state vector*

^{T}**x**

*at iteration step*

_{i}*i*is given by [22]:

The measurement vector **Y**_{meas} has covariance matrix **S*** _{ε}* corresponding to noise

*ε*,

**Y**

*=*

_{i}**F**(

**x**

*,*

_{i}**b**) are simulated radiant fluxes generated by the forward model

**F**(

**x**

*,*

_{i}**b**) which is a (non-linear) function of

**x**

*and*

_{i}**b**(vector of model parameters not retrieved but sources of error).

**K**

*is the Jacobian matrix of simulated radiant flux partial derivatives with respect to*

_{i}**x**

*. The*

_{i}*a priori*state vector is

**x**

*, with covariance*

_{a}**S**

*. We used the inverse NN function to set ${\mathbf{x}}_{a};{\mathbf{x}}_{a}={\mathbf{m}}_{\mathrm{NN}}^{-1}({\mathbf{Y}}_{\mathrm{meas}})$.*

_{a}**G**

*is the gain matrix of contribution functions. The inverse process starts from an initial guess*

_{i}**x**

_{0}; set to

**x**

*. A convergence criterion checks progress towards the solution $\widehat{\mathbf{x}}$ that minimizes the cost function. Besides measurement noise and*

_{a}*a priori*covariance (smoothing error), other sources of error are uncertainties in the elements of

**b**, and forward-model uncertainties due to physical and/or mathematical simplifications.

#### 2.4. Snow layer structure and the data set ${\{\mathbf{P},\mathbf{R}\}}_{s}{}_{=1,\dots ,{N}_{\mathrm{set}}}$

In order to train the forward and inverse NN functions, we generate a large data set of ${\{\mathbf{P},\mathbf{R}\}}_{s}{}_{=1,\dots ,{N}_{\mathrm{set}}}$ in the ground-based spectral radiometer system for albedo and flux (GSAF) [38] by using the RT model for the atmosphere-snow system described in the section 2.1. The GSAF has multiple spectral channels in the visible and near-infrared spectral range, and measures three spectral radiant flux components: downward diffuse flux, downward global flux, and upward flux. We apply our proposed algorithm to the GSAF to retrieve the snow parameters as its validation and application. Details of the GSAF will be described in Section 4.

We used *N*_{set} = 10,000 forward computations with randomly selected geometry and atmosphere/snow parameters: *θ*_{0} in the range [40*°* 70*°*]; aerosol optical thickness, *τ _{a}* in the range [0.01 − 1.3]; cloud optical thickness,

*τ*in the−range [0 − 1];

_{c}*r*

_{top}in the range [10 − 2000

*μ*m];

*r*

_{sub}in the range [10 − 2000

*μ*m] and

*c*

_{imp}in the range [0−.01 − 1.5 ppmw]. The spectral albedo, which is the ratio of the upward flux and the downward global flux, and direct and diffuse fraction of the downward global flux at center wavelengths

*λ*= 0.443,0.87 and 1.225

*μ*m were simulated for a total of

*N*

_{set}= 10,000 combinations of the components.

## 3. Validation of the retrieval algorithm by using the accurate RT model to create a synthetic data set

In this section, we use synthetic data to confirm the quality of the forward and the inverse NN functions as well as the retrieval algorithm based on the non-linear optimal estimation.

#### 3.1. The forward NN function: **R** = **m**_{NN}(**P**)

We compared spectral albedos calculated by the forward NN function **R** = **m**_{NN}(**P**) with those calculated by the spectrally detailed DISORT model (accurate reference model) {**R}** as a function of the geometry and atmosphere/snow parameters (Fig. 2). The data set
${\{\mathbf{P},\mathbf{R}\}}_{s}{}_{=1,\dots ,{N}_{\mathrm{set}}}$ was used as the accurate reference model since the data set was generated by the DISORT model. The plotted albedos were calculated for a two-layer homogeneous snow model. As the results, the two sets of albedos agreed almost completely for all wavelengths with determination coefficients *R*^{2} close to 1. The results indicate that the forward NN function perfectly performs very well and produces accurate spectral albedo values. When we use non-linear optimal estimation method to do the retrievals, we will call the forward RT model many times in the iterative retrieval process. Use of the the forward NN function will speed up the forward calculation as well as the snow parameter retrievals by a factor of 100 or more compared to the LUT approach.

#### 3.2. The inverse NN function: $\mathbf{P}={\mathbf{m}}_{\mathrm{NN}}^{-1}(\mathbf{R})$

Next, the inverse NN function was validated by using the same data set
${\{\mathbf{P},\mathbf{R}\}}_{s}{}_{=1,\dots ,{N}_{\mathrm{set}}}$. Figures 3(a)–(c) show the comparision of snow parameters as the first guess
${\mathbf{x}}_{0}={\mathbf{m}}_{\mathrm{NN}}^{-1}(\mathbf{R})$ calculated by the inverse NN function with those {**P**} as an input data. The results give a high correlation with a slope of almost 1. The first guess is highly correlated with the input (reference) data in all three parameters. The root mean square errors (RMSE) are also very small. It follows from these results that the inverse NN function can provide a good first guess to the retrieval algorithm. The accuracy of the retrieval algorithm is highly dependent on the first guess, implying that having a good initial value could greatly improve the retrieval results. In fact, the optimal estimation may diverge if the first guess is poor, as alluded to previously.

#### 3.3. The retrieval algorithm as a whole system

Finally, we validate our retrieval algorithm containing the forward/inverse NN functions and the optimal estimation as a whole system. The same data set
${\{\mathbf{P},\mathbf{R}\}}_{s}{}_{=1,\dots ,{N}_{\mathrm{set}}}$ was used. Figures 3(d)–(f) show a comparison of input snow parameters {**P**} and the corresponding retrieved parameters
$\widehat{\mathbf{x}}$ by the non-linear optimal estimation. The determination coefficients *R*^{2} in all three parameters were very high (> 0.98) and were higher than those obtained from only the inverse NN function used to get first guesses [Figs. 3(a)–(c)]. The RMSE values are also improved after the iterations by the non-linear optimal estimation, and most of the data are within a relative error of 5%. The reason is that the first guess obtained by the inverse NN function provides a very accurate initial value for the non-linear optimal estimation. Additionally, the Jacobians derived from the forward NN function provide also very accurate values. Therefore, our algorithm used to the non-linear optimal estimation with the RBN-NN method gives more accurate retrievals than by using only inverse NN techique.

In summary, our algorithm proposed here has several unique features:

- The information available in all channels is used to solve the inverse problem by the optimal estimation method to produce simultaneous retrieval of all desired parameters.
- The forward NN function is used to compute the simulated spectral radiant fluxes and the Jacobians, which speeds up the retrieval and saves computing time.
- The inverse NN function is used to obtain a good first guess of our retrieval, implying that this approach will effectively decrease the number of iterations required for the optimal estimation algorithm to converge, which save computing time in the retrieval process as well, and thereby makes the algorithm more efficient.

## 4. Validation and application of the retrieval algorithm by using in-situ measurements at Sapporo and Memuro in Japan

In order to further validate the retrieval algorithm, we apply it to *in-situ* measurement data taken by the GSAF [38]. The channels installed in the GSAF are close to those of the MODIS and SGLI sensors. So, this is the ground-based system for remote sensing of snow parameters. We retrieve snow parameters from the GSAF data, and compare them with data obtained from snow pit work to validate the retrieval algorithm. In this section, we provide a introduction of the instrumentation, the observation site, and the snow conditions. Then we discuss the validation results and provide comparisons with data obtained from the snow pit work.

#### 4.1. Instrumentation and observation site

The GSAF has multiple spectral channels in the visible and near-infrared spectral range and measures three spectral radiant flux components: downward diffuse flux, downward global flux, and upward flux (Fig. 4). So, the GSAF provides the spectral albedo and the direct/diffuse fraction of downward flux each minute during only the daytime.

The GSAF is a unique radiometer system: (1) downward and upward fluxes are measured with the same sensor (one sensor) by turning the sensor unit implying that radiometric calibration of only one sensor is required to measure albedo without any calibration and comparison of output signals of two separate radiometers to measure albedo, and (2) the GSAF can measure the diffuse fraction of the downward flux which is used for the correction of the cosine response of a sensor required for accurate measurements of spectral albedo and the diffuse component of the downward flux. The deviation from an ideal cosine response shows that the GSAF underestimates the radiant flux by about 20% at incident angle *θ* = 60*°* and 40% at *θ* = 80*°*. The spectral dependence is however very small, less than a few percent, for all values of *θ*. Thus, it is easy to systematically correct for this deviation.

The spectral channels of the GSAF are at center wavelengths *λ* = 0.40, 0.443, 0.675, 0.87, 1.02, 1.225 and 1.60 *μ*m and are close to the GCOM-C/SGLI or Terra(Aqua)/MODIS channels. For the measurements taken in 2007–2008, five spectral channels at *λ* =0.40, 0.443, 0.87, 1.225 and 1.60 *μ*m were available, while for 2008–2009, seven channels at *λ* = 0.40, 0.443, 0.675, 0.87, 1.02, 1.225 and 1.60 *μ*m were available. From these channels, we used 0.443, 0.87 and 1.225 *μ*m to retrieve snow parameters in the validation and the application of the algorithm because these wavelengths are sensitive to the snow parameters defined by the GCOM-C/SGLI cryospheric products [20].

We used the GSAF data measured at two sites in Hokkaido, Japan. These two sites are the meteorological observation field at the Institute of Low Temperature Science, Hokkaido University (Sapporo site; 43*°*04′56″N, 141*°*20′30″E; 15m a.s.l) and the National Agricultural Research Center for the Hokkaido Region at the Memuro research station (Memuro site; 42*°*53′20″N, 143*°*04′26″E; 80 m a.s.l). One site is located at the west side of Hokkaido close to the Sea of Japan and the other at the east side close to the Pacific ocean. The difference in location leads to a contrast in snow accumulation due to meteorological differences (more details are described in the section 4.2). The GSAF instrument was mounted at the height of 1.8 m and deployed on flat grasslands without any surface slope and any shadow by forests. The snow surface within the field of view of the radiometer is homogeneous for both sites. Note that in the winter season of 2008–2009 there are no data available during a short period at the Sapporo site, and after February 8, 2009 at the Memuro site due to a malfunction of the GSAF operational system.

Figure 5 shows daily variations in the spectral albedo for three spectral channels at *λ* = 0.443,0.87 and 1.225 *μ*m during two winters. Hereafter, we refer to the daily 30 min-averaged data, collected from 1131 to 1200 h local time, around local solar noon. The Sapporo site has many cloudy or snowy days because Sapporo is located on the Sea of Japan side of Hokkaido where cold northwesterly winds off the Asian Continent bring frequent snowfalls, whereas at the Memuro site there are many clear days although heavy snowfalls related to extratropical cyclones sometimes occur. So, the contrast between the two sites can be seen in the spectral albedo.

When snow has accumulated on the grass land, the spectral albedo is drastically changed. In early December the visible albedo jumps up from approximately 0.1 to 0.9. In the accumulation season, there are frequent snowfalls at the Sapporo site so the albedo oscillates approximately between 0.9 and 1.0 in the visible channels, and between 0.5 and 0.8 at the *λ* = 1.225 *μ*m channel. In contrast, the spectral albedo at the Memuro site increased steeply and then gradually decreased with elapsed time after the snowfall. With the start of the melting season in March, the albedo decreases at both sites.

#### 4.2. Snow condition

*In-situ* measurements by snow pit work were conducted twice a week for the Sapporo site and twice a month for the Memuro site through the winter season. Snow depth, snow type, snow grain size, snow density, snow temperature, and mass concentration of the snow impurities were measured. In general the measured snow depth increases until the end of February or the beginning of March, and decreases in March at those regions (Fig. 5). We refer to the former period as the “accumulation period” and the latter as the “melt season”. Unfortunately no data from snow-pit-work are available in 2008–2009 for the Memuro site.

For the Sapporo-site, the observed snow types in the upper layers in the accumulation season consist of precipitation particles, decomposing and fragmented precipitation particles and rounded grains while the lower layers mainly consist of melt forms. The formation of lower melt forms is due to the upward soil heat flux from the underlying ground. In the melt season, all layers consist of melt forms with ice layer. For the Memuro-site, internal snow types in the accumulation season were quite different from those observed at the Sapporo-site though results of snow-pit works were limited at the Memuro-site. Solid-type depth hoar (faceted crystals) and depth hoar were dominant in the middle and lower layers of the snowpack. The snow shapes employed in this paper are in accordance with ‘The International Classification for Seasonal Snow on the Ground’ published by the International Hydrological Programme of the UNESCO [39].

The snow grain size is estimated by using a handheld lens with a scale of 10 *μ*m resolution. We define grain size here as half the branch width of dendrites, as half of the width of narrower part of broken crystals or the radius of each spherical particle (*r*_{2} of Aoki et al. [40]). We estimated the minimum, median and maximum values of *r*_{2} in each snow layer because its use in radiative transfer modeling led to a good correlation between spectral hemispherical reflectance and *r*_{2} [18, 40–42]. The averaged values of measured snow grain mode sizes in the snow surface
${r}_{2}^{\mathrm{sfc}.-1}$ (*d* = 0−1 cm depth) and those in the sub-surface
${r}_{2}^{\mathrm{sfc}:1-3}$ (*d* = 1−3 cm depth) during the two winters are shown in Fig. 6. During the accumulation season, grain sizes were small while during the melting season, grain sizes were large with an air temperature rise.

The snow from two snow layers of *d* = 0 − 2 cm depth and of *d* = 0 − 10 cm depth was sampled to estimate the mass concentration of snow impurities. Since the technical details of the method used to estimate it were described elsewhere [38, 43, 44], we briefly introduce the methodology in this paragraph. The snow samples were collected in stainless bags in winter 2007–2008 and in dust-free plastic bags in winter 2008–2009. The mass concentration of snow impurities was estimated by filtering the snow samples with a Silver Membrane filter of a pore size of 0.45 *μ*m. Total mass concentration of snow impurities (
${c}_{\mathrm{all}}^{\mathrm{sfc}-2}$ for *d* = 0−2 cm depth, and
${c}_{\mathrm{all}}^{\mathrm{sfc}-10}$ for *d* = 0−10 cm depth) was estimated from the difference in the filter weight before and after filtering [19]. Afterward, the mass concentration of elemental carbon (EC) (
${c}_{\mathrm{EC}}^{\mathrm{sfc}-2}$ for *d* = 0 − 2 cm depth and
${c}_{\mathrm{EC}}^{\mathrm{sfc}-10}$ for *d* = 0 10 cm depth) and organic carbon (OC) (
${c}_{\mathrm{OC}}^{\mathrm{sfc}-2}$ for *d* = 0−2 cm depth and
${c}_{\mathrm{OC}}^{\mathrm{sfc}-10}$ for *d* = 0 −10 cm depth) was obtained for each sampled snow layer by a Carbon Aerosol Laboratory Instrument (Sunset Laboratory Inc.) based on the thermal optical reflectance method [45]. For this measurement we adopt the Interagency Monitoring of Protected Visual Environments (IMPROVE) thermal evolution protocol [46]. The instrument was calibrated by a standard solution of sucrose. The reproductive experiment of this instrument by the standard showed good agreement within 3%. Finally we estimated the dust concentration (
${c}_{\mathrm{dust}}^{\mathrm{sfc}-2}$ for *d* = 0 − 2 cm depth and
${c}_{\mathrm{dust}}^{\mathrm{sfc}-10}$ for *d* = 0 10 cm depth) by subtracting the total carbon concentration (*c*_{EC} + *c*_{OC}) from the total mass concentration of snow impurity (*c*_{all}).

#### 4.3. Comparison of snow parameters between measurements and retrievals

Figure 6 shows the snow grain sizes and snow impurity concentrations measured at Sapporo and Memuro for the two winter seasons of 2007–2008 and 2008–2009, and those retrieved by the non-linear optimal estimation method. Since we used a two-layer snow model that takes into account the vertical variation of the snow grain size, we retrieve a topmost-layer grain size
${\widehat{r}}_{\mathrm{top}}$ as well as subsurface-layer grain size
${\widehat{r}}_{\mathrm{sub}}$. For the snow impurity retrieval, since we used one-layer snow model, a depth-averaged snow impurity is retrieved. First, we discuss the retrieved snow parameters, and then we show the validation results by comparing the retrieved ones with the *in-situ* measured ones.

### 4.3.1. Variations of snow parameters during the two winters

The results clearly show that at the Memuro site the retrieved snow grain size decreased when snowfalls occurred and then gradually increased afterwards whereas at the Sapporo site, short-term variations (several-day cycle) were measured due to the frequent snowfalls. This different pattern is due to the meteorological differences between the two sites. The subsurface-layer grain size ${\widehat{r}}_{\mathrm{sub}}$ was generally larger than that of the topmost-layer ${\widehat{r}}_{\mathrm{top}}$. Furthermore, the variation of ${\widehat{r}}_{\mathrm{sub}}$ was larger than that of ${\widehat{r}}_{\mathrm{top}}$ throughout the season. The snow grain sizes in the melting season increase gradually from February to April and become the largest at the end of March.

On the snow impurity concentration *ĉ*_{imp}, a low snow impurity concentration was retrieved at both sites due to the new snow after snow is accumulating. The value of *ĉ*_{imp} at the Sapporo site was generally higher than that at the Memuro site during the two winter seasons. It could be because Sapporo is an urban area while Memuro is a rural area. High concentrations of snow impurities were retrieved in early December. Since the measured snow depth is less than 30 cm in early December, the optical depth of the snow is not large enough compared to the snow optical depth used in the retrieval algorithm. The algorithm does not support such a shallow snow case as visible light could reach the ground resulting in the retrieval of a lower visible albedo. We did not include these values in the validation of retrieved snow physical parameters.

We compare the retrieved snow grain sizes and the *in-situ* measurement data. In Fig. 6, the error bars show the measured snow-grain-radius range based on the measured minimum and maximum snow grain sizes. Comparison between the retrieved snow grain size and the measured one shows that at the Sapporo site, both
${\widehat{r}}_{\mathrm{top}}$ and
${\widehat{r}}_{\mathrm{sub}}$ in the accumulation season were roughly consistent with measured snow grain sizes and most retrieved snow grain sizes fall within the range of the *in-situ* snow grain size. In the melting season, the values of
${\widehat{r}}_{\mathrm{sub}}$ were qualitatively consistent with the *in-situ* measurements, whereas the values of
${\widehat{r}}_{\mathrm{top}}$ were somewhat underestimated. This underestimation particularly for large snow grain was already reported [18, 38, 47–49]. More detailed discussion of this issue is provided in the next section. At the Memuro site, the retrieved results were roughly consistent with measured snow grain sizes but revealed a tendency for both
${\widehat{r}}_{\mathrm{top}}$ and
${\widehat{r}}_{\mathrm{sub}}$ to be underestimated.

The retrieved concentration of snow impurities and the *in-situ* measurement data are also shown in Fig. 6. We assumed the snow impurity to be the black carbon (BC) in the retrieval algorithm described in the Section 2.1. Thus, the concentration of optically equivalent BC *c*_{eq.BC} was calculated by using *c*_{EC} and *c*_{dust} obtained in the filtering and the carbon analyses, and was used for the comparison between the retrieved snow impurity and the *in-situ* measurement data. The BC-equivalent impurity concentration *c*_{eq.BC} is defiend by

*d*= 0 − 2 cm depth with error bars. The upper extent of the error bar represents the concentration of dust ${c}_{\mathrm{dust}}^{\mathrm{sfc}.-2}$, whereas the lower extent is the concentration of EC ${c}_{\mathrm{EC}}^{\mathrm{sfc}.-2}$ for

*d*= 0 − 2 cm depth obtained in the filtering and the carbon analyses. Most of

*ĉ*

_{imp}fell between ${c}_{\mathrm{EC}}^{\mathrm{sfc}.-2}$ and ${c}_{\mathrm{dust}}^{\mathrm{sfc}.-2}$, and were generally closer to ${c}_{\mathrm{eq}.\mathrm{BC}}^{\mathrm{sfc}.-2}$ The results indicate that EC would be the main component of light-absorbing snow impurities at Sapporo and Memuro. More detailed discussion is provided in the section 4.3.3.

### 4.3.2. Comparison of snow grain size

We retrieve the topmost-layer grain size
${\widehat{r}}_{\mathrm{top}}$ and the subsurface-layer grain size
${\widehat{r}}_{\mathrm{sub}}$ from the two-layer snow model considering the vertical variation of the snow grain size. Though we used three channels for simultaneous retrieval of snow parameters,
${\widehat{r}}_{\mathrm{top}}$ is mainly retrieved from the albedo at *λ* = 1.225 *μ*m. Because the light absorption by ice is strong at *λ* = 1.225 *μ*m, only the topmost-layer snow contributes to the reflected light (upward flux) at 1.225 *μ*m. Furthermore, the snow impurity has almost no effect at this wavelength [3]. Therefore, the snow albedo at *λ* = 1.225 *μ*m has information about the topmost-layer grain size. Meanwhile
${\widehat{r}}_{\mathrm{sub}}$ is retrieved from three wavelengths at *λ* = 0.443, 0.87 and 1.225 *μ*m due to the contribution of *ĉ _{imp}* and
${\widehat{r}}_{\mathrm{top}}$. The light absorption by ice at

*λ*= 0.87

*μ*m is several hundred times weaker than that at

*λ*= 1.225

*μ*m. So, we expect the channel at

*λ*= 0.87

*μ*m to give information about the grain size in deeper layers. Figure 7 shows the relationship between the retrieved snow grain sizes and the

*in situ*measured depth-averaged value

*r*

_{2}. The determination coefficient

*R*

^{2}for ${\widehat{r}}_{\mathrm{top}}$ is 0.653, implying that the retrieved snow grain size and the

*in situ*measured one are well correlated. The statistics was roughly consistent with or better than that reported in the previous stidies by the GSAF snow parameter retrievals and the satellite products [18, 38]. However, ${\widehat{r}}_{\mathrm{top}}$ is underestimated for the large grain sizes occurring in the melting season (two open circles in Fig. 7(a)). ${\widehat{r}}_{\mathrm{top}}$ and ${r}_{2}^{\mathrm{sfc}-1}$ were 478

*μ*m and 1000

*μ*m, respectively, on March 11, 2008 and 499

*μ*m and 1500

*μ*m, respectively, on March 21, 2008 at the Sapporo site. So, the value of RMSE for ${\widehat{r}}_{\mathrm{top}}$ is 256

*μ*m (Table 1) which is not acceptable for small grains such as a precipitation particles and decomposing and fragmented precipitation particles. Similar results were reported in previous studies [18, 38, 47–49]. One possible reason of this underestimation is the “sun crust” observed on wet snow surfaces under a clear sky. It can be frequently seen in the melting season. Tanikawa et al. [49] confirmed the high reflectance in the forward scattering direction of the sun crust over the granular snow. This increased reflectance due to sun crust would cause an underestimate of the retrieved snow grain size. Another probable reason is less sensitivity to the larger grain size at

*λ*= 1.225

*μ*m. The wavelength at

*λ*= 1.225

*μ*m has a strong light absorption by ice, implying that the sensitivity to the larger grain size is too small [50]. When selecting the data for ${r}_{2}^{\mathrm{sfc}.-1}$ < 1000

*μ*m, the values of

*R*

^{2}= 0.941 with RMSE = 47

*μ*m are obtained. This result indicates that the wavelength at

*λ*= 1.225

*μ*m would be suitable for the small and middle grain sizes retrievals.

For the grain size in the subsurface layer [Fig. 7(b)],
${\widehat{r}}_{\mathrm{sub}}$ is relatively well correlated with
${r}_{2}^{\mathrm{sub}:1-3}$ and *R*^{2} is 0.771. The value of RMSE = 173 *μ*m with a small bias is acceptable for melt forms with large grains (Table 1). These statistics were roughly consistent in the previous studie [18, 38]. However, it is insufficient for the small grains such as precipitation particles and decomposing and fragmented precipitation particles. There are both underestimates and overestimates of
${\widehat{r}}_{\mathrm{sub}}$ around
${r}_{2}^{\mathrm{sub}:1-3}$ *∼* 100 *μ*m. The penetration depth for the mean effective radius of a spherical snow particle *r*_{eff} = 100 *μ*m is approximately 2 cm and the critical snow depth is up to 5 cm [2]. Since snow grain size often changes sharply with depth near the surface layer [30, 31], the assumption of a homogeneous snow layer in the subsurface layer might be cause an insufficient RMSE. We need to add more snow layers in the retrieval algorithm in the future.

### 4.3.3. Comparison of snow impurity

In this section, we show the validation results for the concentration of snow impurities. From the results in the Fig. 6, the main component of light-absorbing snow impurities is the EC at both Sapporo and Memuro sites. The mineral dust is another main component of snow impuritues which is less absorptive than EC [3, 40]. Hereafter, we validate the retrieved snow impurity concentration *ĉ*_{imp} with *c*_{eq.BC} values.

Figure 8(a) depicts scatter plots between the retrieved *ĉ*_{imp} and *in-situ c*_{eq.BC} for *d* = 0 − 2 cm and *d* = 0 − 10 cm, respectively. The results of the statistical analysis are shown in Table 1. The both statistics were basically insufficient and the values were different between
${c}_{\mathrm{eq}.\mathrm{BC}}^{\mathrm{sfc}.-2}$ and
${c}_{\mathrm{eq}.\mathrm{BC}}^{\mathrm{sfc}.-10}$. According to Kuchiki et al. [38], the penetration depths depend on the values of snow grain size and snow impurities which are retrieved together. The penetration depths at *λ* = 0.443 *μ*m and 0.87 *μ*m are 2.9 and 1.1 cm for the case of the small grain size *r*_{eff} = 20 *μ*m (with a mass concentration of a black carbon *c*_{BC} = 0.1 ppmw), respectively, while those for large snow grains *r*_{eff} = 2000 *μ*m (with *c*_{BC} = 0.1 ppmw) were 32 and 9.7 cm, respectively. This situation indicates that for the validation we need to vary the snow layer of the *in-situ* measurements according to the values of the topmost snow grain radius
${\widehat{r}}_{\mathrm{top}}$. Thus, we employed
${c}_{\mathrm{eq}.\mathrm{BC}}^{\mathrm{sfc}.-2}$ values when
${\widehat{r}}_{\mathrm{top}}$ *≤* 80 *μ*m and
${c}_{\mathrm{eq}.\mathrm{BC}}^{\mathrm{sfc}.-10}$ when
${\widehat{r}}_{\mathrm{top}}$ > 80 *μ*m. The threshold snow grain radius (
${\widehat{r}}_{\mathrm{top}}$ = 80 *μ*m) was determined so as to obtain the highest correlation coefficient between the retrieved values of *ĉ*_{imp} and *c*_{eq.BC} based on the method of Kuchiki et al. [38]. Figure 8(b) shows the relationship between the retrieved snow impurity concentration *ĉ*_{imp} and the *in-situ* measured BC equivalent snow impurity concentration *c*_{eq.BC} based on the rule described above. Retrieved mass concentrations of impurity *ĉ*_{imp} are relatively correlated with *c*_{eq.BC} with *R*^{2} of 0.529 and RMSE of 0.197 ppmw. These statistics were totally improved and were better than those reported in previous studies [18,38] which relied on a one-layer homogeneous snow model to retrieve snow physical parameters. A plausible reason is that we used a two-layer snow model considering vertical variation in the snow grain size. The value of *ĉ*_{imp} was retrieved from a combination of three channels at *λ* = 0.46,0.86 and 1.225 *μ*m. The channel at *λ* = 0.46 *μ*m is the most sensitive of the three channels to the concentration of snow impurity, but it is also sensitive to snow grain size [3, 40]. In this study, we employed a two-layer snow model to retrieve accurate snow grain size, resulting in improved accuracy of the retrieved snow impurity concentrations.

In spite of the improved accuracy of the snow impurity retrievals, the retrieved *ĉ*_{imp} values are generally higher than the *c*_{eq.BC} values with a bias of +0.050 ppmw and RMSE of 0.197 ppmw (Table 1). This discrepancy is probably due to our assumption that the snow impurity employed in the retrieval algorithm consisted entirely of black carbon; that is, the retrieved *ĉ*_{imp} is defined as the mass concentration of impurities optically equivalent to that of black carbon. Scanning electron microscope photographs of snow impurities taken in Sapporo site show irregularly-shaped particles (dust) as well as coagulated particles (EC), and the averages of EC, OC, and dust mass fraction were 5.4%, 8.5%, and 86.1% from the filtering and the carbon analyses for the six winter seasons of 2007–2013 [44]. The mass concentrations of dust were considerably larger than that of EC and OC. This finding implies that besides the black carbon assumed in the retrieval algorithm, absorption by mineral dust and organic carbon should also be taken into account. Mineral dust has optical properties substantially different from those of black carbon, and can absorb and scatter solar radiation, resulting in a reduction of the visible albedo of snow when the dust particles are deposited in the snow surface. For the organic carbon, carbon molecules contain in organic matter in the snow samples. Since a wide variety of organic carbon forms are present [51], there are uncertainties in the inherent optical properties of OC: we do not know with certainty the snow impurity type, the particle size/shape distributions and the chemical composition. So, for the validation we did not consider the effect of OC on the BC-equivalent impurity concentration defined by Eq. (6). Determining the inherent optical properties of the mineral dust and organic carbon is important in order to improve the accuracy of snow impurity retrievals. Another possible reason for the discrepancy discussed above is the snow model employed in the retrieval algorithm. We ignored vertical inhomogeneity in the retrieval of the snow impurity concentration. The field measurements showed that the mass concentration of snow impurities in the snow layer *d* = 0 − 2 cm depth is higher than that in *d* = 0−10 cm depth. Since the retrieval result is a depth-averaged snow impurity concentration, it may be desirable to use a multi-layer model for snow impurity retrieval.

## 5. Summary

This paper describes a new snow retrieval algorithm based on optimal estimation. We used a radial basis function neural network (RBF-NN) for fast computations of radiant fluxes, Jacobians and retrieval parameters. Both a forward and an inverse RBF-NN were employed. A non-linear optimal estimation method was employed as the retrieval scheme. Our algorithm approach has several unique features:

- The forward NN function is used to provide accurate spectral radiant quantities and the Jacobians, which reduces the computation time of the retrieval by a factor of 100 or more compared to the LUT approach.
- The inverse NN function is used to obtain a good first guess of our retrieval. This approach will prevent the optimal estimation algorithm from diverging, and decrease the number of iterations required for it to converge.
- The information available in all channels is used to solve the inverse problem by the optimal estimation method providing a simultaneous retrieval of all desired parameters.

We confirm the quality of the forward and the inverse RBF-NN function as well as the retrieval scheme based on the non-linear optimal estimation method by using the dataset derived by DISORT model as an accurate reference model. The forward RBF-NN function yields spectral albedos in close agreement with the reference model for all wavelength with determination coefficients *R*^{2} close to 1. The inverse RBF-NN function, which can be used for the first guess in the retrieval algorithm, produces accurate physical parameters with *R*^{2} *≥*0.97 from the result of the comparison with the reference model. The comparison between retrieved physical parameters derived from our proposed algorithm containing the forward and inverse RBF-NN functions, and the corresponding reference physical parameters shows a high correlation in all three parameters with *R*^{2} *≥* 0.98 and RMSE < 10%. These values were improved after the iteration by the non-linear optimal estimation, compared to only the inverse NN function. The reason is that both the Jacobians and the first guess derived from the forward and inverse RBF-NN functions, respectively, provide very accurate velues in the non-linear optimal estimation.

A simple retrieval procedure for two-snow-layer model was constructed and applied to the GSAF data taken in Sapporo and Memuro, Japan during two winter seasons 2007–2008 and 2008–2009. Snow grain sizes and snow impurity concentrations were retrieved from the spectral albedo at *λ* = 0.443, 0.87 and 1.225 *μ*m. We found that retrieved snow parameters were well correlated with snowfall and snow aging. Validation with *in-situ* measurements showed that the determination coefficients *R*^{2} were approximately 0.529 for the snow impurity concentration, and 0.653 and 0.771 for the topmost- and the subsurface-layer snow grain size, respectively. In particular, the topmost-layer snow grain size derived mainly from the wavelength at *λ* = 1.225 *μ*m, was well correlated with the small and middle grain size (*R*^{2} = 0.971). These statistics were roughly consistent with or better than those reported in the previous studies [18, 38]. Therefore, our algorithm proposed here is useful for monitoring local variations in snow physical parameters and as validation data for satellite snow products with high precision and accuracy.

In order to retrieve more detailed information about snow vertical structure, we will try adding more layers in the snow model, and using more spectral channels because radiation measured in different channels originates from different depths in the snow [21, 38]. This is a common issue for both snow grain size and impurity retrievals. Since there are seven GSAF channels available, we will explore this issue in the future. In addition to this, determining the inherent optical properties of the mineral dust and organic carbon is important for improving the snow impurity retrievals.

Finally, JAXA starts the GCOM-C mission which aims to conduct global, long-term observation of the carbon cycle and radiation budget [20]. The SGLI aboard GCOM-C is an optical sensor, having two major new features: 250 m spatial resolution and polarization/multidirectional observation capabilities. We define the retrievals of the snow grain size and the concentration of snow impurity as well as snow cover extent as a cryospheric product, and prepare the algorithms for the cryospheric product of the GCOM-C/SGLI [52]. The retrieval algorithm proposed here is now implemented for the snow parameter retrievals. The algorithm presented in this study is useful for the airborne/satellite remote sensing.

## Acknowledgments

We thank T. J. Yasunari, M. Takahashi, T. Kuno, Y. Sawada, K. Shimoyama, T. Nakai, T. Sueyoshi, and J. Mori of the Institute of Low Temperature Science, Hokkaido University for performing the snow pit work at Sapporo throughout two winters, and E. Tanaka of the Japan Meteorological Agency and M. Niwano of the Meteorological Research Institute for analysis of snow impurities. We also thank Y. Iwata and T. Hirota of the National Agricultural Research Center for Hokkaido Region for providing the observational site and meteorological data in Memuro. This work was supported in part by (1) the GCOM-C/SGLI Mission, JAXA, (2) the Ministry of Education, Culture, Sports, Science and Technology, Grant-in-Aid for Scientific Research (S), 23221004, 2011, and (3) the Grant for Joint Research Program, the Institute of Low Temperature Science, Hokkaido University.

## References and links

**1. **J. Hansen and L. Nazarenko, “Soot climate forcing via snow and ice albedos,” Proc. Natl. Acad. Sci. USA **101**, 423–428 (2004). [CrossRef]

**2. **W. J. Wiscombe and S. G. Warren, “A model for the spectral albedo of snow, I, Pure snow,” J. Atmos. Sci. **37**, 2712–2733 (1980). [CrossRef]

**3. **S. G. Warren and W. J. Wiscombe, “A model for the spectral albedo of snow, II, Snow containing atmospheric aerosols,” J. Atmos. Sci. **37**, 2735–2745 (1980). [CrossRef]

**4. **A. W. Nolin and S. Liang, “Progress in bidirectional reflectance modeling and applications for surface particulate media: snow and soils,” Remote Sens. Rev. **18**, 307–342 (2000). [CrossRef]

**5. **T. Tanikawa, T. Aoki, and F. Nishio, “Remote sensing of snow grain size and snow impurities from airborne multispectral scanner data using a snow bidirectional reflectance distribution function model,” Ann. Glaciol. **34**, 74–80 (2002). [CrossRef]

**6. **M. Hori, T. Aoki, K. Stamnes, and W. Li, “ADEOS-II/GLI snow/ice products: Part III - Retrieved results,” Remote Sens. Environ. **111**, 291–336 (2007). [CrossRef]

**7. **A. A. Kokhanovsky, V. V. Rozanov, T. Aoki, D. Odermatt, C. Brockmann, O. Krüger, M. Bouvet, M. Drusch, and M. Hori, “Sizing snow grains using backscattered solar light,” Int. J. Remote Sens. **22**, 6975–7008 (2011). [CrossRef]

**8. **A. Mary, M. Dumont, J. P. Dedieu, Y. Durand, P. Sirguey, H. Milhem, O. Mestre, H. S. Negi, A. A. Kokhanovsky, M. Lafaysse, and S. Morin, “Intercomparison of retrieval algorithms for the specific surface area of snow from near-infrared satellite data in mountainous terrain, and comparison with the output of a semi-distributed snow-pack model,” The Cryosphere **7**, 741–761 (2013). [CrossRef]

**9. ** AVIRIS homepage, http://aviris.jpl.nasa.gov/links/index.html

**10. ** AMSS homepage, http://sharaku.eorc.jaxa.jp/ADEOS2/COMMON/AMSS/index.html

**11. **M. Tedesco and A. A. Kokhanovsky, “The semi-analytical snow retrieval algorithm and its application to MODIS data,” Remote Sens. Environ. **111**, 228–241 (2007). [CrossRef]

**12. **T. Scambos, T. Haran, M. Fahnestock, T. H. Painter, and J. Bohlander, “MODIS-based mosaic of Antarctica (MOA) data sets: Continent-wide surface morphology and snow grain size,” Remote Sens. Environ. **111**, 242–257 (2007). [CrossRef]

**13. **Z. Jin, T. P. Charlock, P. Yang, Y. Xie, and W. Miller, “Snow optical properties for different particle shapes with application to snow grain size retrieval and MODIS/CERES radiance comparison over Antarctica,” Remote Sens. Environ. **112**, 3563–3581 (2008). [CrossRef]

**14. **K. Stamnes, W. Li, H. Eide, T. Aoki, M. Hori, and R. Storvold, “ADEOS-II/GLI snow/ice products: Part I -Scientific basis,” Remote Sens. Environ. **111**, 258–273 (2007). [CrossRef]

**15. **A. Lyapustin, M. Tedesco, Y. Wang, T. Aoki, M. Hori, and A. A. Kokhanovsky, “Retrieval of snow grain size over Greenland from MODIS,” Remote Sens. Environ. **113**, 976–987 (2009). [CrossRef]

**16. **H. Wiebe, G. Heygster, E. Zege, T. Aoki, and M. Hori, “Snow grain size retrieval SGSP from optical satellite data: Validation with ground measurements and detection of snow fall events,” Remote Sens. Environ. **128**, 11–20 (2013). [CrossRef]

**17. **J. G. Winther, “Landsat Thematic Mapper (TM) derived reflectance from a mountainous watershed during the snow melt season,” Nordic Hydrol. **23**, 273–290 (1992).

**18. **T. Aoki, M. Hori, H. Motoyohi, T. Tanikawa, A. Hachikubo, K. Sugiura, T. J. Yasunari, R. Storvold, H. Eide, K. Stamnes, W. Li, J. Nieke, Y. Nakajima, and F. Takahashi, “ADEOS-II/GLI snow/ice products: Part II - Validation results using GLI and MODIS data,” Remote Sens. Environ. **111**, 274–290 (2007). [CrossRef]

**19. **T. Aoki, A. Hachikubo, and M. Hori, “Effects of snow physical parameters on shortwave broadband albedos,” J. Geophys. Res. **108**, 4616 (2003). [CrossRef]

**20. **K. Imaoka, M. Kachi, H. Fujii, H. Murakami, M. Hori, A. Ono, T. Igarashi, K. Nakagawa, T. Oki, Y. Honda, and H. Shimoda, “Global Change Observation Mission (GCOM) for monitoring carbon, water cycles, and climate change,” Proc. IEEE **98**(5), 717–734 (2010). [CrossRef]

**21. **W. Li, K. Stamnes, and B. Chen, “Snow grain size retrieved from near-infrared radiances at multiple wavelengths,” Geophys. Res. Lett. **28**, 1699–1702 (2001). [CrossRef]

**22. **C. D. Rodgers, *Inverse Methods for Atmospheric Sounding - Theory and Practice* (World Scientific Publishing Co. Ple. Ltd., 2000).

**23. **P. M. Atkinson and A. R. L. Tatnall, “Neural networks in remote sensing - introduction,” Int. J. Remote Sens. **18**, 699–709 (1997). [CrossRef]

**24. **L. Tsang, Z. Chen, S. Oh, and R. J. Marks, “Inversion of snow parameters from passive microwave remote sensing measurements by a neural network trained with a multiple scattering model,” IEEE Trans. Geosci. Remote Sens. **GE-30**, 1015–1024 (1992). [CrossRef]

**25. **A. Abdelgadir, “Forward and inverse modeling of canopy directional reflectance using a neural network,” Int. J. Remote Sens. **19**, 453–471 (1998). [CrossRef]

**26. **L. Meng, Y. He, J. Chen, and Y. Wu, “Neural network retrieval of ocean surface parameters from SSM/I data,” Mon. Weather Rev. **135**, 586–597 (2007). [CrossRef]

**27. **L. Fan, W. Li, A. Dahlback, J. J. Stamnes, S. Stamnes, and K. Stamnes, “New neural-network-based method to infer total ozone column amounts and cloud effects from multi-channel, moderate bandwidth filter instruments,” Opt. Express **22**, 19595–19609 (2014). [CrossRef] [PubMed]

**28. **K. Stamnes, S. C. Tsay, W. J. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. **27**, 2502–2509 (1988). [CrossRef] [PubMed]

**29. **K. Stamnes, B. Hamre, J. J. Stamnes, G. Ryzhikov, M. Biryulina, R. Mahoney, B. Hauss, and A. Sei, “Modeling of radiation transport in coupled atmosphere-snow-ice-ocean systems,” J. Quant. Spectrosc. Radiat. Transfer **112**, 714–726 (2011). [CrossRef]

**30. **M. Schneebeli and S. A. Sokratov, “Tomography of temperature gradient metamorphism of snow and associated changes in heat conductivity,” Hydrol. Process. **18**, 3655–3665 (2004). [CrossRef]

**31. **M. Matzl and M. Schneebeli, “Measuring specific surface area of snow by near-infrared photography,” J. Glaciol. **52**, 558–564 (2006). [CrossRef]

**32. **J. E. Hansen and L. D. Travis, “Light scattering in planetary atmosphere,” Space Sci. Rev. **16**, 527–610 (1974). [CrossRef]

**33. **T. C. Grenfell and S. G. Warren, “Relationship of a nonspherical ice particle by a collection of independent spheres for scattering and absorption of radiation,” J. Geophys. Res. **104**, 31697–31709 (1999). [CrossRef]

**34. **S. G. Warren and R. E. Brandt, “Optical constants of ice from the ultraviolet to the microwave: a revised compilation,” J. Geophys. Res. **113**, D14220 (2008). [CrossRef]

**35. **M. Hess, P. Koepke, and I. Schult, “Optical properties of aerosols and clouds: the software package OPAC,” Bull. Am. Meteorol. Soc. **79**, 831–844 (1998). [CrossRef]

**36. **R. Spurr, K. Stamnes, H. Eide, W. Li, K. Zhang, and J. J. Stamnes, “Simultaneous retrieval of aerosols and ocean properties: a classic inverse modeling approach. I. Analytic Jacobians from the linearized CAO-DISORT model,” J. Quantum Spectrosc. Radiat. Transfer **104**, 428–449 (2007). [CrossRef]

**37. **W. Li, K. Stamnes, R. Spurr, and J. J. Stamnes, “Simultaneous retrieval of aerosol and ocean properties by optimal estimation: SeaWiFS case studies for the Santa Barbara Channel,” Int. J. Remote Sens. **29**, 5689–5698 (2008). [CrossRef]

**38. **K. Kuchiki, T. Aoki, T. Tanikawa, and Y. Kodama, “Retrieval of snow physical parameters using a ground-based spectral radiometer,” Appl. Opt. **48**, 5567–5582 (2009). [CrossRef] [PubMed]

**39. **C. Fierz, R. L. Armstrong, Y. Durand, P. Etchevers, E. Greene, D. M. McClung, K. Nishimura, P. K. Satyawali, and S. A. Sokratov, *The International Classification for Seasonal Snow on the Ground* (International Association of Cryosphere Sciences, UNESCO-International Hydrological Programme, 2009).

**40. **Te. Aoki, Ta. Aoki, M. Fukabori, A. Hachikubo, Y. Tachibana, and F. Nishio, “Effects of snow physical parameters on spectral albedo and bi-directional reflectance of snow surface,” J. Geophys. Res. **105**, 10219–10236 (2000). [CrossRef]

**41. **A. A. Kokhanovsky, T. Aoki, A. Hachikubo, M. Hori, and E. P. Zege, “Reflective properties of natural snow: approximate asymptotic theory versus in situ measurements,” IEEE Trans. Geosci. Remote Sens. **43**, 1529–1535 (2005). [CrossRef]

**42. **T. Tanikawa, T. Aoki, M. Hori, A. Hachikubo, O. Abe, and M. Aniya, “Monte Carlo simulations of spectral albedo for artificial snowpacks composed of spherical and nonspherical particles,” Appl. Opt. **45**, 5310–5319 (2006). [CrossRef] [PubMed]

**43. **T. Aoki, K. Kuchiki, M. Niwano, Y. Kodama, M. Hosaka, and T. Tanaka, “Physically based snow albedo model for calculating broadband albedos and the solar heating profile in snowpack for general circulation models,” J. Geophys. Res. **116**, D11114 (2011). [CrossRef]

**44. **K. Kuchiki, T. Aoki, M. Niwano, S. Matoba, Y. Kodama, and K. Adachi, “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. **120**, 868–882 (2015).

**45. **J. C. Chow, J. G. Watson, L. C. Pritchett, W. R. Pierson, C. A. Frazier, and R. G. Purcell, “The DRI thermal/optical reflectance carbon analysis system: description, evaluation and applications in U.S. air quality studies,” Atmos. Environ. **27**, 1185–1201 (1993). [CrossRef]

**46. **J. C. Chow, J. G. Watson, D. Crow, D. H. Lowenthal, and T. Merrifield, “Comparison of IMPROVE and NIOSH carbone measurements,” Aerosol Sci. Technol. **34**, 23–34 (2001). [CrossRef]

**47. **M. Fily, B. Bourdelles, J. P. Dedieu, and C. Sergent, “Comparison of in situ and Landsat Thematic mapper derived snow grain characteristics in the Alps,” Remote Sens. Environ. **59**, 452–460 (1997). [CrossRef]

**48. **J. E. Kay, A. R. Gillespie, G. B. Hansen, and E. C. Pettit, “Spatial relationship between snow contaminant content, grain size, and surface temperature from multispectral images of Mt. Rainier, Washington (USA),” Remote Sens. Environ. **86**, 216–231 (2003). [CrossRef]

**49. **T. Tanikawa, T. Aoki, M. Hori, A. Hachikubo, and M. Aniya, “Snow bidirectional reflectance model using non-spherical snow particles and its validation with field measurements,” EARSeL eProceedings **5**, 137–145 (2006).

**50. **T. Tanikawa, M. Hori, T. Aoki, A. Hachikubo, K. Kuchiki, M. Niwano, S. Matoba, S. Yamaguchi, and K. Stamnes, “In situ measurements of polarization properties of snow surface under the Brewster geometry in Hokkaido, Japan and northwest Greenland ice sheet,” J. Geophys. Res. **119**, 13946–13964 (2014).

**51. **M. O. Andreae and A. Gelencsér, “Black carbon or brown carbon? The nature of light-absorbing carbonaceous aerosols,” Atmos. Chem. Phys. **6**, 3131–3148 (2006). [CrossRef]

**52. **N. Chen, W. Li, T. Tanikawa, M. Hori, T. Aoki, and K. Stamnes, “Cloud mask over snow-/ice-covered areas for the GCOM-C1/SGLI cryosphere mission: Validations over Greenland,” J. Geophys. Res. **119**, 12287–12300 (2014).