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.

 

Fig. 1 Conceptual illustration of snow retrieval algorithm. The RBF-NN is trained as the forward and inverse NN function, respectively, based on a large data set of {P,R} computed by the DISORT model. Then, the non-linear optimization loop did the retrieval. Y_meas is set to the measurement data and the solar zenith angle as input data. The vector $\widehat{\mathbf{x}}$ contains the retrieved snow parameters.

Download Full Size | PDF

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, θ₀), 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₂O, O₃, SO₂ 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:

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_jk obtained in Eq. (1) to calculate the Jacobians 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

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₁, x₂,…, x_n]^T, and the measured radiant fluxes into another vector Y_meas = [Y₁, Y₂, …, Y_m]^T. In optimal estimation, the update of the state vector x_i at iteration step 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_i and b (vector of model parameters not retrieved but sources of error). K_i is the Jacobian matrix of simulated radiant flux partial derivatives with respect to x_i. The a priori state vector is x_a, with covariance S_a. We used the inverse NN function to set ${\mathbf{x}}_a;{\mathbf{x}}_a={\mathbf{m}}_{\mathrm{NN}}^{-1}({\mathbf{Y}}_{\mathrm{meas}})$. G_i is the gain matrix of contribution functions. The inverse process starts from an initial guess x₀; set to x_a. A convergence criterion checks progress towards the solution $\widehat{\mathbf{x}}$ that minimizes the cost function. Besides measurement noise and 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: θ₀ in the range [40° 70°]; aerosol optical thickness, τ_a in the range [0.01 − 1.3]; cloud optical thickness, τ_c in the−range [0 − 1]; 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² 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.

 

Fig. 2 Comparison of spectral albedo derived by the forward NN function R = m_NN(P) with those of the reference data set {R}. R² and RMSE are a determination coefficient and a root mean square error, respectively.

Download Full Size | PDF

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.

 

Fig. 3 (top) Comparison of the first guess x₀ of snow parameters derived by the inverse NN function $\mathbf{P}={\mathbf{m}}_{\mathrm{NN}}^{-1}(\mathbf{R})$ with those by the reference data set {P}, and (bottom) comparison of retrieved snow parameters $\widehat{\mathbf{x}}$ with those by the reference data set {P}. Color bar is the relative error for each point in the scatter plot, which is [(P – x₀)/P] ×100 for the first guess (top), and $[(P-\widehat{x})/P]\times 100$ for the retrieved parameter (bottom). Unit of x- and y-axis is μm for r_top and r_sub, and ppmw for c_imp. R² and RMSE are the determination coefficient and the root mean square error, respectively.

Download Full Size | PDF

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² 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:

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.

 

Fig. 4 Photographs show GSAF instrument and three measurement modes [38]. (a) The GSAF was mounted at a height of 1.8 m above the ground surface. (b) The diffuse component of the downward flux is measured with the shadow blade blocking the direct component. The downward global (upward) radiant flux is measured by turning the sensor unit upward (downward); (c) global mode and (d) upward mode. (e) Schematic illustration of the GSAF.

Download Full Size | PDF

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.

 

Fig. 5 Snow albedo measurements at Sapporo [(a) and (b)] and Memuro [(c) and (d)] during two winters of 2007–2008 [(a) and (c)] and 2008–2009 [(b) and (d)]. These measured values are 30 min averages from 1131 to 1200 h local time, around local solar noon. There are no data available during the period marked by arrows. The snow depth data are also plotted in each figure.

Download Full Size | PDF

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₂ of Aoki et al. [40]). We estimated the minimum, median and maximum values of r₂ in each snow layer because its use in radiative transfer modeling led to a good correlation between spectral hemispherical reflectance and r₂ [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.

 

Fig. 6 Daily variations of snow grain size and snow impurity concentration measured in the two winter season 2007–2008 and 2008–2009 and those retrieved by the non-linear optimal estimation method; Sapporo site [(a) and (b)], and Memuro site [(c) and (d)]. ${\widehat{r}}_{\mathrm{top}}$ and ${\widehat{r}}_{\mathrm{sub}}$ are retrieved topmost- and subsurface-grain size, respectively. $r_2^{\mathrm{sfc}.-1}$ is the in-situ measured snow grain radius at a snow surface (d = 0 − 1 cm depth) and $r_2^{\mathrm{sub}:1-3}$ is the same but averaged in 1 to 3 cm layer. The upper and lower extents of r₂-error bar represent the maximum and minimum of the in-situ measured snow grain radius, respectively. ĉ_imp is retrieved snow impurity. $c_{\mathrm{eq}.\mathrm{BC}}^{\mathrm{sfc}.-2}$ are in-situ measured the BC-equivalent impurity concentration for d = 0 − 2 cm depth. The upper extent of $c_{\mathrm{eq}.\mathrm{BC}}^{\mathrm{sfc}.-2}$-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. There are no in-situ data available at Memuro in 2008/09 winter seasion.

Download Full Size | PDF

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

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₂. The determination coefficient R² 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² = 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.

 

Fig. 7 Relationship between retrieved snow grain size and in situ measured one: (a) topmost-layer grain size ${\widehat{r}}_{\mathrm{top}}$ versus in-situ measured depth-averaged grain size $r_2^{\mathrm{sfc}.-1}$ and (b) subsurface-layer snow grain size ${\widehat{r}}_{\mathrm{sub}}$ versus $r_2^{\mathrm{sub}:1-3}$. The solid line represents a 1:1 relationship. R² is the determination coefficient. In Fig. (a), two open circles are for $r_2^{\mathrm{sfc}.-1}$ ≥ 1000 μm.

Download Full Size | PDF

Table 1. Summary of statistical analysis between retrieved snow physical parameters and in situ measurements.

View Table

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² 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² 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.

 

Fig. 8 Relationship between retrieved mass concentration ĉ_imp and in situ measured BC-equivalent mass concentration of snow impurity c_eq.BC. (a) The snow layers of the in situ measurements are 0 to 2 cm (bullet) and 0 to 10 cm (open circle), and (b) 0 to 2 cm (bullet) when ${\widehat{r}}_{\mathrm{top}}$ ≤ 80 μm and 0 to 10 cm (open circle) when ${\widehat{r}}_{\mathrm{top}}$ > 80 μm. The solid line represents a 1:1 relationship. R² is the determination coefficient.

Download Full Size | PDF

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:

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² 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² ≥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² ≥ 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² 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² = 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.