3 D-VAR multilayer assimilation of X-band SAR data into a detailed snowpack model

Introduction Conclusions References


Introduction
Accurate knowledge of snowpack internal structure is essential for better understanding its evolution over time and provides greater benefit to snow forecasting, water resource monitoring and prediction of natural hazards, such as snow avalanches warning.For this purpose, snowpack evolution models are developed in order to simulate the evolution of snow based on meteorological variables.These models are currently limited due to the lack of snow in-situ stratigraphic observations.The new generation of spaceborne Synthetic Aperture Radar (SAR), with 1 m resolution and short revisit time provides information which can be assimilated into snow evolution models.Through reflectivity can be used to analyze the snowpack parameters calculated by a snow evolution model and therefore adjust these values according to modeling and observations error statistics.In this paper, the EBM initially developed by Longepe et al. (2009) based on Dense Media Radiative Transfer (DMRT) theory, which is compatible with dry snow measured at C-band (5 GHz), has been adapted for X (10 GHz) and Ku (14 GHz) frequency bands.The air-snow interface and snow-ground backscattering components are calculated using the Integral Equation Model (IEM) by Fung and Chen (2004).The dry snow permittivity is calculated based on Strong Fluctuation Theory (SFT) introduced by Stogryn (1984), which was tested and verified by Wang et al. (2000) for higher wave frequencies.The model calculates the total backscattering coefficient σ 0 pq for different polarization channels from the physical features of each snow layer (optical grain radius, density, thickness), the roughness of air-snow and snow-ground interface, and for given SAR acquisition conditions (frequency, incidence angle).
Dry snow is a dense and heterogeneous medium with strong variations of various physical properties such as grain optical diameter, density, and thickness.Many studies have been carried out on the retrieval of different snowpack properties from SAR data, such as snow water equivalent (Shi and Dozier, 2000), liquid water content (Shi et al., 1993), wet snow cover maps (Nagler and Rott, 2000), etc.In general, these studies concentrate on inversing the physical models, which enables the retrieval of such properties.For a multilayer snowpack, the number of observations, i.e. the SAR backscattering coefficients, is much smaller than the number of unknown variables, i.e. the snow properties of each layer.Classical estimation approaches based on the use of an inverse problem would not be viable.Instead, in our study, an adjoint operator of the direct EBM is developed to be used in a data assimilation scheme.
Data assimilation has been widely used in meteorology studies (Courtier et al., 1998;Uppala et al., 2005), land surface modeling (Slater and Clark, 2006;De Lannoy et al., 2010;Toure et al., 2011).Data assimilation using physically-based multi layer models has been initiated in recent studies, related to snow passive microwave radiance Introduction

Conclusions References
Tables Figures

Back Close
Full  (Toure et al., 2011) or albedo assimilation (Dumont et al., 2012).It has been proven to be an effective technique that allows the combination of observations into an a priori information on the state of a system, which then produce an a posteriori state.The a priori information is often referred to as guess parameters, whereas the a posteriori state is called the analysis.The guess parameters in this study are the physical properties of each snowpack layer simulated using the snow evolution model.The analysis is obtained by modifying the guess information based on the backscattering coefficient obtained from SAR acquisitions, according to the error statistics of both model and observations.The simulation of snowpack evolution is then continued with the analysis result.The intermittent assimilation algorithm is carried out each time a new SAR acquisition is available, therefore the assimilation is propagated over time, which allows to constrain the simulation using remote sensing observations.In our study, the three-dimensional variational analysis (3D-VAR) method (Courtier et al., 1998) is implemented.The case study is carried out using the snow evolution model SURFEX/ISBA-Crocus (Vionnet et al., 2012) and TerraSAR-X acquisitions performed over the French-Alps in winter 2009.The Argentière Glacier area has been chosen for this case study due to its large and rather uniformly snow-covered surface area.Introduction

Conclusions References
Tables Figures

Back Close
Full

Main components of the total backscattering coefficient
The Stokes vector, which contains the incoherent information related to the polarization of an electromagnetic wave, can be expressed as follows (Ulaby et al., 1981): where E h and E v represent the horizontal and vertical components of the Jones vector on the electric field, and .represents the expectation operator.
For given acquisition conditions, the Stokes vector scattered by a medium, g s , can be related to the incident one, g i , by a Stokes matrix M (Lee and Pottier, 2009), as g s = Mg i , with: reflection at the surface air-snow interface, volume scattering, volume-ground and ground-volume interactions, and reflection over the ground (Martini et al., 2003).Due to their small amplitude, the volume-ground and ground-volume contributions can be neglected (Floricioiu and Rott, 2001).The illustration of the three other mechanisms is shown in Fig. 1.The expression of the total polarimetric backscattered information can be written using the Mueller matrix corresponding to each mechanism: The air-snow interface and ground backscattering are modeled using the IEM introduced by Fung et al. (Fung and Chen, 2004), whereas the volume contribution is calculated using the Vector Radiative Transfer equation.

Surface backscattering
The matrix M as represents the second order polarimetric response backscattered by the air-snow interface.Its elements can be calculated from the air-snow interface roughness parameters, i.e. its correlation function w(x) and its root mean square (rms) height σ h , the incidence angle θ 0 and the emitted EM wave frequency f using the IEM (Fung and Chen, 2004).According to the IEM, the reflectivity may be expressed as: where p and q are equal to h or v , indicating a horizontal or vertical polarization, k 0 = 2πf c represents the wave number.The detailed mathematical expressions of the surface spectrum W (k) and the Fresnel reflection/transmission factor |I n pq | can be found in (Fung and Chen, 2004).Introduction

Conclusions References
Tables Figures

Back Close
Full

Volume backscattering
The volume backscattering M vol depends on various scattering mechanisms occurring during the propagation through a multilayer snowpack, which can be categorized into 4 types: related to transmission between two layers, absorption by the snow particles, scattering and coherent recombination.The amplitude of each mechanism depends largely on the dielectric properties of the snowpack medium.Therefore the permittivity of each layer, which characterizes its dielectric properties, needs to be calculated first.

Dry snow permittivity
Dry snow is considered as a dense and heterogeneous medium with strong variations of various physical properties such as grain optical diameter, density, and thickness.
Therefore the variance of permittivity across a snow layer is relatively high.Several snowpack characterization methods (Hallikainen et al., 1986) are largely based on the assumption that the scattering losses due to the correlation of EMW are negligible.However, at high frequency, the snowpack structure becomes bigger compared to the wavelength of X-or Ku-band EM waves (Huining et al., 1999) and the correlation between particles can no longer be ignored.The Strong Fluctuation Theory (SFT) introduced by Stogryn (Stogryn, 1984) can model the permittivity of such medium by using the effective permittivity ( eff ) that takes into account the scattering effect among ice particles at high frequencies.The expression of eff using the SFT is as follows (Huining et al., 1999): where g and δ g are the quasi-static permittivity and its variance, k 0 the wave number and L the correlation length, which is proportional to the average snow grain optical diameter and the snow density of the medium.Introduction

Conclusions References
Tables Figures

Back Close
Full

Transmission between two layers
The snowpack consists of layers with different physical properties.Therefore the model needs to take into account the energy loss due to transmission between two layers.
With the assumption of a smooth interface between two layers, the Fresnel transmission can be used.It is expressed through a matrix as follows (Ulaby et al., 1981): where 2 represents the Fresnel transmission coefficients of the pp channel, whereas g k(k−1) and h k(k−1) are the terms of Mueller matrix related to the co-polarized correlation (Longepe et al., 2009):

Attenuation
The particles in a snowpack are generally considered as spherical (Floricioiu and Rott, 2001;Koskinen et al., 2010).Due to the symmetry of the particle shape, the extinction of a wave propagating through the snowpack is independent of the polarization and may hence be represented by a scalar coefficient.The extinction is composed of an absorption and a scattering terms:

Conclusions References
Tables Figures

Back Close
Full It can also be computed through the effective permitivity eff (Huining et al., 1999): The attenuation matrix represents the gradual loss in EMW intensity while penetrating through a multilayer snowpack, composed of layers with different physical properties.It takes into account the energy loss by absorption and scattering mechanisms based on the extinction coefficient κ e and thickness d of the layer, as well as the loss by transmission effect while an EM propagates through different layers: The Att down is the intensity loss when propagating from the surface to layer k , whereas Att up represents the intensity loss from layer k to the surface.The exponential factor, which takes into account the gradual loss of energy throughout the layer, is deduced from the basic radiative transfer equation dI = Iκ e dr where r = d/ cos θ.

Scattering by the particles
The phase matrix P k , under the hypothesis of spherical particles has the form shown in Eq. ( 2) where the cross-polarization terms P 12 and P 21 , which correspond to σ hv and σ vh , are null.In the backscattering case, with the assumption of spherical particles, the SFT phase matrix can be simplified to P k = Introduction

Conclusions References
Tables Figures

Back Close
Full

Calculation of the volume backscattering
Considering a snowpack made of n distinct layers, where θ k is the incidence angle and d k is the thickness of layer k , the total contribution of the volume backscattering mechanism M vol can be written as follows:

Ground backscattering
The backscattering M g of the snow-ground interface may be computed as: where R(θ n ) represents the contribution of the underlying ground backscattering and can be determined using the IEM.

The detailed snowpack model Crocus
Crocus is a one-dimensional numerical model simulating the thermodynamic balance of energy and mass of the snowpack.Its main objective is to describe in detail the evolution of internal snowpack properties based on the description of the evolution of morphological features of snow grains during their metamorphism.the general scheme of Crocus.It takes as inputs the meteorological variables air temperature, relative air humidity, wind speed, solar radiation, long wave radiation, amount and phase of precipitation.When used in the French mountain ranges (Alps, Pyrenees and Corsica), these quantities are provided by the SAFRAN system, which combines ground-based, radioprobes and remote sensing (cloudiness) observations with an a priori estimate of meteorological conditions from a numerical weather prediction (NWP) model (Durand et al., 1993;Durand, 2009).SAFRAN meteorological fields are assumed to be homogeneous within a given mountain range and provide a description of the altitude dependency of meteorological variables.The output of Crocus includes scalar physical properties of the snowpack (snow depth, Snow Water Equivalent -SWE, surface temperature, albedo, . . . ) along with the internal physical properties for each layer (density, thickness, optical grain radius, . . .).This study uses the latest version of the detailed snowpack model Crocus, recently incorporated in the land surface scheme ISBA within the SURFEX interface (Vionnet et al., 2012).Among other advantages over previous versions of Crocus, this allows seamless coupling of the snowpack to the state of the underlying ground.

Introduction to the assimilation method
Variational assimilation aims to integrate observation data into guess variables through the use of an observation operator.It is widely used in meteorological studies in order to relate observations and modeling (Courtier et al., 1998;Dumont et al., 2012).
The method concentrates on searching a solution that minimizes simultaneously the distance between observations and simulation results and the distance between initial guess variables and the analysed variables.A scheme of this process is presented in Fig. 3.The outputs of the EBM described in the previous section, such as backscattering coefficient at HH and VV polarizations, are used as elements of the observation operator H ebm (x):

Conclusions References
Tables Figures

Back Close
Full where x represents the set of variables describing the snowpack properties (here, density and optical grain diameter for each snow layer).The 3D-VAR (Courtier et al., 1998) algorithm is based on the minimization of a cost function J(x), defined as: where x is called the state vector, and can be modified after each iteration of the minimization, x g is the initial guess of the state vector and remains constant during the whole process.Therefore x − x g 2 serves as a distance between the modified profile and the starting point.The observed polarimetric response, y obs , contains calibrated values of the backscattering coefficients σ 0 for different polarimetric channels.
Therefore, y−H ebm (x) 2 represents the distance between simulated and observed radiometric quantities.The process also requires the estimation of the error covariance matrices of observations/simulations R and of the model B, the guess error covariance.

Adjoint operator and minimization algorithm
In order to minimize the cost function J, one needs to calculate its gradient: Once the adjoint operator is developed, the minimization of J can be achieved using a gradient descent algorithm.Each iteration consists of modifying the vector x according to the Newton method:

Conclusions References
Tables Figures

Back Close
Full where ∇ 2 J(x n ) is the gradient of second order (Hessian) of J:

General comments on the chosen analysis method
In general, modeling techniques are used to establish the relationship between the physical properties of natural environment and observations measured by a specific equipment (such as SAR or optical sensors).An inverse approach may then be developed to characterize the environment using the observations.However, such problems often require solving an underdetermined system, with a number of unknown quantities higher than the number of equations.
In our case, the length of the input state vector x can reach 100 (in the case of snowpack with 50 layers, frequently generated by Crocus), whereas the output of the model only consists of backscattering coefficients corresponding to the polarimetric channels of SAR data.Therefore the realization of an inverse model is theoretically highly ill-conditionned and practically pointless.Data analysis methods, on the other hand, require a vector of guess variables relatively close to the actual values, allowing to add an a priori information.The snowpack variables calculated by Crocus are used as guess in our assimilation scheme.The fundamental goal is to modify the initial guess variables, while balancing the errors of guess, modeling and measurements.It should be noted that, as the problem remains underdetermined, the analysis scheme only serves as a method to improve the initial guess variables using the new observations from SAR data.The quality of improvement is based on the estimation of the initial guess vector x g and on the precision of the EBM.

TCD Introduction Conclusions References
Tables Figures

Back Close
Full  In order to carry out the comparison between the backscattering coefficients σ sim (obtained from the EBM using Crocus snowpack profile as an input) and σ TSX (obtained from TerraSAR-X reflectivity), the images were multi-looked to (20 m × 20 m) wide pixels and Frost filter was applied with window size 5 × 5 pixels.
Figure 5 shows the evolution in time of TerraSAR-X backscattering coefficients on the same area at altitude 2700 m over the Argentière Glacier.It can be observed that for the period from 6 January to 17 January (blue triangle) and from 8 February to 19 February (red circle), the sets of comparison values are well below the equality line, which means the backscattering coefficients decreased over these periods.The opposite effect can be noted for the period from 17 January to 8 February (green cross).The ground medium beneath the snowpack consists of glacier ice and its roughness can be considered as constant as it is not affected by snowfall or wind induced transportation, therefore these changes in backscattering may be due to the evolution of snowpack structure.Introduction

Conclusions References
Tables Figures

Back Close
Full

Crocus snowpack data
The intrinsic parameters of a snowpack needed for EBM simulations are simulated by Crocus, which consists of the number of snow layers, and their density, optical grain diameter and thickness.These quantities are used as inputs for the simulation of volume backscattering mechanism.The relation between open-loop (i.e., without assimilation) Crocus data and the TerraSAR-X reflectivity for different altitudes over the Argentière Glacier is shown in Fig. 6.The optical thickness τ is the product of snow depth by the relative permittivity (Tsang et al., 2007).In the case of multilayer snowpack, it is defined as: According to the SFT, the effective permitivity eff is largely depends on the correlation length L = 0.85D/3 where D is the optical grain diameter.It can be observed that the vertical profile description provided by Crocus may be used to separate TerraSAR-X reflectivity at different altitudes.Figure 6 also indicates the sensivity of the reflectivity to volume-related snow parameters.

Sensitivity of the EBM to snowpack parameters
In order to assess the sensitivy of the EBM outputs with respect to the different properties of a snowpack, a set of simulations were run for varying snowpack structures.The random dataset correspond to a thin to fairly thick snowpack (from 30 cm to 200 cm), with low to high density (250 to 500 kg m −3 ) and small to large values of the optical diameter (0.5 to 1 mm).Measurements of the roughness parameters of air-snow interface and ground are not available, therefore empirical values for the correlation length and the rms height from (Oh et al., 1992)  for the air-snow interface; whereas σ h = 1.12 cm and l = 8.8 cm, corresponding to very rough surface, are chosen for snow-ground interface due to the characteristic of ice beneath the snowpack over the study area.The results of EBM simulations are plotted vs SWE in Fig. 7.As it can be observed, as the SWE exceeds 500 mm, which is the case for Argentière Glacier in winter, the contribution of the ground backscattering mechanism is considerably lower than the volume one.It shows that in our case, the EBM output are sensitive to the evolution of the internal snowpack structure.

Results and discussions
With preset air-snow interface and ground parameters, the original model input vector x = [x Crocus x surface x ground ] t may be reduced to the Crocus variables consisting in density and optical diameter for each snow layer: where d i and ρ i are respectively the grain optical diameter and the density of the ith layer of the snowpack.This means that the analysis process does not modify the thickness of each layer, therefore the total snow depth is unchanged but not the total SWE.
At first iteration of the algorithm, x = x g , given by the Crocus snow profile.The covariance matrix B, which represents the error of the guess profile, i.e. of the Crocus calculation, is a (2n × 2n) definite positive matrix.Each element of B is computed as: where σ i represents the standard deviation of the error on x i , which has been experimentally estimated to 0.3 mm and 65 kg m −3 respectively for grain optical diameter and snow density.Introduction

Conclusions References
Tables Figures

Back Close
Full The coefficient γ ij represents the correlation between errors on x i and x j and are modelled as: where ∆h ij is the distance in cm between layer i and layer j.The values of α and β depend on different types of correlations and can be splitted into 3 cases 5: - These values are issued from an ensemble of slightly perturbed Crocus runs, i.e. obtained by differences in their meteorological inputs, over one winter season.The deviation between these runs, considered as elementary perturbations, have been then statistically studied and fitted with the model of Eq. ( 22) for each pair of variables only.
In this case study, SAR data is available for the HH channel, therefore the error covariance matrix R reduces to a scalar, equal to the variance of the SAR image intensity over the studied area.The calculations of the variance at several altitudes over the Argentière Glacier gives the average value of R = 0.03.
Crocus snow stratigraphic profiles have been computed for 7 different altitudes over the Argentière Glacier, from 2400 m to 3000 m.The level of Liquid Water Content per volume (LWC v ) at the time and location of experiment is 0 percent, therefore the condition of dry snow is satisfied.Figure 4 shows the approximate locations of each study area on the glacier.
Figure 8 shows the results of simulation and analysis using two TerraSAR-X acquisitions.The reflectivity on the glacier crevasse area (around 2600 m) has very high standard variation due to the cracks, and therefore has been masked.The red line corresponds to the TerraSAR-X reflectivities along the glacier, whereas the cyan diamond- assimilation, i.e. independent of TerraSAR-X observations.The blue triangles indicate the EBM simulation of the "guess states" of snowpack, which is the Crocus output that took into account all the modifications of previous dates.These "guess states" are in turn modified by the assimilation process to become the "analysed states".The EBM simulations of the analysed states are shown in green circles.The analysed snow parameters are used to reinitialize Crocus for the next iteration, which then produce the guessed state for the next assimilation when new SAR acquisition is available.The agreement between TerraSAR-X reflectivity and the output of EBM using Crocus simulated profiles can be observed in Fig. 8. Overall the EBM is able to follow important changes like those between the altitudes 2500 m and 2900 m.There is a considerable increase of reflectivity at the altitude of 3000 m, probably due to firn, which accounts for very high range of density.The gap between EBM simulations and TerraSAR-X reflectivities in this zone may be due to the exceptional structure of snowpack with potentially extremely rough internal layers, for which the EBM is not adapted.
Figure 9 shows a detailed analysis of the modifications of the snow stratigraphic profiles done by the data assimilation process.The input parameters contain the snow grain optical diameter from Crocus and the density of each snow layer.It can be observed that the assimilation algorithm tends to modify the grain optical diameter and density in the deep layers in order to affect the backscattering intensity and reduce the discrepancy between TerraSAR-X observations and Crocus simulations.The speed of densification process is therefore modified in the Crocus simulations with assimilation.

Conclusions
The results of this study show the potential of using data analysis method and a multilayer snowpack backscattering model based on the radiative transfer theory in order to improve snowpack detailed simulation.The proposed new backscattering model adapted to X-band and higher frequencies enables a fairly accurate calculation of EMW losses in each layer of the snowpack.Through the use of 3D-VAR data assimilation based on the linear tangent and adjoint operator of the EBM, we have the possibility to modify the snowpack profiles calculated by the snowpack evolution model Crocus.
The output of this process shows interesting results, and could be further developed and used in real application such as snow cover area monitoring on massif scale or snowpack evolution through a period of time.Future studies will be concentrated on calibrating the assimilation process using in-situ measurements.Direct field measurements of the optical grain diameter using recently developed methods (Gallet et al., 2009;Arnaud et al., 2011) allow a direct comparison to Crocus output (Morin et al., 2013).Improvements will also stem from the fact that the optical grain diameter will become a prognostic variable of Crocus (currently it is a quantity diagnosed from several semi-empirical microstructure properties, see Vionnet et al., 2012).This work is currently under development.Introduction

Conclusions References
Tables Figures

Adjoint operator and
In order to minimize the cos late its gradient:  According to the SFT, the effective permitivity ef f is largely    According to the SFT, the effective permitivity ef f is largely depends on the correlation length L = 0.85D/3 where D is an Electromagnetic Backscattering Model (EBM) and a data assimilation method, SAR Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2) where M 11 = |σ 0 vv | 2 and M 22 = |σ 0 hh | 2 represent the co-polarized backscattering coefficients, M 12 = |σ 0 vh | 2 and M 21 = |σ 0 hv | 2 the cross-polarized backscattering terms, whereas M 33 = e(σ 0 vv σ 0 hh +|σ 0 hv | 2 ), M 44 = e(σ 0 vv σ 0 hh −|σ 0 hv | 2 ), M 34 = − m(σ correlation terms.Due to the reflection symmetry, the other terms of M are equal to zero (Lee and Pottier, 2009).The solution of the RT equation at order 1 provides a total backscattered information from a snowpack that consists of a combination of five scattering mechanisms: Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure
Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

)
If the model is denoted H ebm : B → R, with B and R are the domain of definition of x and y, then the function ∇H t ebm satisfying: ∀x, y, ∇H t ebm y, x B = y, ∇H ebm x R is the adjoint operator of H ebm .In our case, due to the complexity of the EBM, an analytical solution of the gradient is time consuming and unreliable.Therefore, numerical differentiation has been used to calculate the adjoint model.
Discussion Paper | Discussion Paper | Discussion Paper | For this study, TerraSAR-X descending acquisitions on the region of Chamonix Mont-Blanc, France from 6 January 2009 to 24 March 2009 are available for continuous assimilation, with a revisit time of 11 days.
Discussion Paper | Discussion Paper | Discussion Paper | have been taken into consideration.The values of σ h = 0.4 cm and l = 8.4 cm, equivalent to a slightly rough surface, are used Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 3 .
Fig. 3. Global scheme of the data analysis used in this study.The input of the process are the SAR reflectivities, σ 0 , (observation) and the snowpack stratigraphic profile calculated by Crocus (guess).The output is the analyzed snowpack profile x that minimizes the cost function.

Fig. 4 .Fig. 5 .
Fig. 4. (Top) Location of the TerraSAR-X acquisition in the French Alps.(Bottom) The crop image on the Argentière Glacier area.The three rectangles show the approximate positions of different altitudes on the Argentière Glacier: 2400 m, 2700 m and 3000 m on the TerraSAR-X images.The red line indicates the continuous trail on the glacier where the SAR data will be used in the case study.
Table shows the main features of TerraSAR- steps of 100 m elevation on horizontal terrain were used to drive the detailed snowpack model Crocus throughout the whole season 2008-2009 (starting on 1 August 2008).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | assimilation of X-band SAR data into a detailed snowpack model