1 D-Var multilayer assimilation of X-band SAR data into a detailed snowpack model

The structure and physical properties of a snowpack and their temporal evolution may be simulated using meteorological data and a snow metamorphism model. Such an approach may meet limitations related to potential divergences and accumulated errors, to a limited spatial resolution, to wind or topography-induced local modulations of the physical properties of a snow cover, etc. Exogenous data are then required in order to constrain the simulator and improve its performance over time. Synthetic-aperture radars (SARs) and, in particular, recent sensors provide reflectivity maps of snow-covered environments with high temporal and spatial resolutions. The radiometric properties of a snowpack measured at sufficiently high carrier frequencies are known to be tightly related to some of its main physical parameters, like its depth, snow grain size and density. SAR acquisitions may then be used, together with an electromagnetic backscattering model (EBM) able to simulate the reflectivity of a snowpack from a set of physical descriptors, in order to constrain a physical snowpack model. In this study, we introduce a variational data assimilation scheme coupling TerraSAR-X radiometric data into the snowpack evolution model Crocus. The physical properties of a snowpack, such as snow density and optical diameter of each layer, are simulated by Crocus, fed by the local reanalysis of meteorological data (SAFRAN) at a French Alpine location. These snowpack properties are used as inputs of an EBM based on dense media radiative transfer (DMRT) theory, which simulates the total backscattering coefficient of a dry snow medium at X and higher frequency bands. After evaluating the sensitivity of the EBM to snowpack parameters, a 1D-Var data assimilation scheme is implemented in order to minimize the discrepancies between EBM simulations and observations obtained from TerraSAR-X acquisitions by modifying the physical parameters of the Crocus-simulated snowpack. The algorithm then re-initializes Crocus with the modified snowpack physical parameters, allowing it to continue the simulation of snowpack evolution, with adjustments based on remote sensing information. This method is evaluated using multi-temporal TerraSAR-X images acquired over the specific site of the Argentière glacier (Mont-Blanc massif, French Alps) to constrain the evolution of Crocus. Results indicate that X-band SAR data can be taken into account to modify the evolution of snowpack simulated by Crocus.


Introduction
Accurate knowledge of snowpack internal structure is critical for better understanding the snowpack evolution over time, and is essential for snow forecasting, water resource monitoring and prediction of natural hazards, such as avalanches.For this purpose, snow metamorphism models, such as Crocus (Brun et al., 1992;Vionnet et al., 2012), are developed in order to simulate the evolution of snowpack based on meteorological variables.These models are currently limited due to the lack of in situ snow stratigraphic measurements.For example, in the French Alps, the network of snow and meteorological observations contains about 150-180 stations, which is not enough to adjust a snow model to predict the state and the spatial variability of snowpack at small scale Published by Copernicus Publications on behalf of the European Geosciences Union.(20 m).This limitation results in potential divergences, accumulated errors and limited spatial resolution of the model.Therefore, exogenous data are crucial in order to constrain the simulator and improve its performance over time.
On the other hand, the radiometric properties of a snowpack measured at high frequencies depend strongly on its main physical parameters, like its depth, snow grain size and density.The electromagnetic backscattering model (EBM), initially developed by Longepe et al. (2009) based on dense media radiative transfer (DMRT) theory, allows for simulation of the backscattering coefficient σ 0 of dry snow from C band (5 GHz) to Ku band (14 GHz).The air-snow, σ as , and snow-ground, σ sg (or snow-ice, σ si ), interfaces backscattering components are calculated using the integral equation model (IEM) developed by Fung and Chen (2004).The snow permittivity is calculated using the strong fluctuation theory (SFT) (Stogryn, 1984).The SFT has been tested and verified in the literature (Wang et al., 2000;Tsang et al., 2007).It is also used in the DMRT model of multilayer snowpack developed by Longepe et al. (2009).This model is capable of simulating the interaction of electromagnetic waves with a layer of snow based on the physical parameters (thickness, optical diameter, snow density).The advantage of this model is the simple implementation and its moderate computation time, which is crucial in order to run the data assimilation process, where the electromagnetic model is repeatedly executed multiple times.With this model, we can calculate the total backscattering coefficient σ 0 pq for different polarization channels (p, q = H or V) from the physical features of each snow layer, the roughness of air-snow and snow-ice interfaces, and specific radar illumination (frequency, incidence angle).
The new generation of synthetic-aperture radar (SAR) satellite data provides images with metric resolution and short revisit time.The TerraSAR-X satellite, with 1.477 m ⇥ 2.44 m resolution and revisit time of 11 days, gives dense information both spatially and temporally on snowpack evolution.In this study, we propose a new process which uses these multi-temporal images of TerraSAR-X to constrain the Crocus model through data assimilation.
Data assimilation has been widely used in meteorological studies (Courtier et al., 1998;Uppala et al., 2005) and land surface modeling (Slater and Clark, 2006;De Lannoy et al., 2010;Toure et al., 2011).Data assimilation using physically based multilayer models has been initiated in recent studies, using passive microwave radiance (Toure et al., 2011) or albedo observations (Dumont et al., 2012).The advantages of the assimilation using SAR images are the quasiindependence with respect to atmospheric conditions, the high resolution of analysis, and the sensitivity of SAR responses to the presence and structure of volumetric mediums.The use of data assimilation on SAR data and meteorological models to predict certain physical properties of snowpack has been developed in the literature (Nagler et al., 2008;Takala et al., 2011).This study attempts to implement a data assimilation system which is capable of constraining a detailed snow metamorphism model at a layer scale (modification of the physical properties of each layer) using X-band SAR data.The assimilation techniques have proven effective in combining observations and a priori information to more realistically simulate snowpack conditions (i.e., 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 a 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 us to constrain the snowpack simulation using remote sensing observations.The adjustment made to the snowpack physical properties is based on error statistics of modeling (Crocus) and observation (SAR).This study reports, for the first time, on a new process based on the DMRT model and on the one-dimensional variational analysis (1D-Var) to assimilate TerraSAR-X data into the snow model Crocus.A global schematic of this process is presented in Fig. 1 the study of simulations and sensitivity of snowpack at X band.Section 6 presents the first results and discussion of data assimilation method in the particular case of the Argentière glacier, where the ground beneath the snow consists of ice.

Snowpack model Crocus
Crocus is a one-dimensional numerical model simulating the thermodynamic balance of energy and mass of a 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.It takes as inputs meteorological variables such as air temperature, relative air humidity, wind speed, solar radiation, long-wave radiation, and amount and phase of precipitation.In this study within the French Alps, these meteorological conditions are taken from the SAFRAN reanalysis, which combines ground-based, radiosondes 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, assumed to be homogeneous for a given mountain range and elevation in the French Alps region, 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, etc.) along with the internal physical properties for each layer (density, thickness, optical radius, etc.).This study uses the latest version of the detailed snowpack model Crocus, recently incorporated into 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.

Main components of the total backscattering coefficient
The Stokes vector, which contains the incoherent information related to the polarization of an electromagnetic wave (EMW), 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 h.i represents the expectation operator.For given acquisition conditions, the Stokes vector of radiation scattered by a medium, g s , can be related to the incident one, g i , by a Stokes matrix M (Lee and Pottier, 2009) (Lee and Pottier, 2009).
The first-order solution of the radiative transfer (RT) equation provides the total backscattered information from a snowpack that consists of a combination of five scattering mechanisms: reflection at the surface air-snow interface, volume scattering, volume-ice and ice-volume interactions, and reflection from the snow-ice interface (Martini et al., 2003).Due to their small amplitude, the volume-ice and ice-volume contributions can be neglected (Floricioiu and Rott, 2001).The illustration of the three other mechanisms is shown in Fig. 2. The expression of the total polarimetric backscattered information can be written using the Mueller matrix corresponding to each mechanism: (3) The air-snow interface (M as ) and snow-ice interface backscattering (M si ) are modeled using the IEM introduced by Fung and Chen (2004), whereas the volume contribution (M v ) is calculated using the vector radiative transfer equation.

Air-snow interface 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 rootmean-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, and 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).

Snow volume backscattering
The volume backscattering M v depends on various scattering mechanisms occurring during the propagation through a multilayer snowpack, which can be categorized into four types: (1) transmission between two layers, (2) attenuation by the snow particles, (3) scattering and (4) 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 to be a dense and heterogeneous medium with strongly variable physical properties.Therefore the variance of permittivity across a snow layer is relatively high.The SFT, introduced by Stogryn (1984), can model the permittivity of such a medium by using the effective permittivity (✏ eff ) that takes into account the scattering effects among ice particles at high frequencies.The expression of ✏ eff using the SFT is as follows (Huining et al., 1999): where j is the imaginary unit; ✏ g and δ ✏ g are the quasistatic permittivity and its variance; k 0 is the wave number; and L = 0.85 D/3 is the correlation length, with D being the snow optical diameter.

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 k is the layer number and k − 1 is the layer above it, t 2 represents the Fresnel transmission coefficients of the pp channel, and 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 to be 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 term: It can also be computed through the effective permittivity ✏ 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: Att down is the intensity loss (attenuation) 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 ✓ and I is the EMW intensity.

Scattering by the particles
The phase matrix P k , under the hypothesis of spherical particles, has the form shown in Eq. ( 2), where the crosspolarization terms P 12 and P 21 , which correspond to σ hv and σ vh , are equal to 0. In the backscattering case, with the assumption of spherical particles, the SFT phase matrix can be simplified to P k = 3 s 8⇡ I 4 , where I 4 is the (4 ⇥ 4) identity matrix (Tsang et al., 2007).The assumption of spherical particles can simplify the modeling problem; however, it prevents the simulations of the backscattering coefficient over crosspolarization channels (HV and VH).

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 (Fig. 2), the total contribution of the volume backscattering mechanism M v can be written as follows:

Snow-ice interface backscattering
The backscattering M si of the snow-ice interface is computed as where R(✓ n ) represents the contribution of the snow-ice interface backscattering and can be determined using the IEM., optical diameter 1 mm, and snow depth 30-400 cm.The glacier roughness is fixed at σ si = 0.9 cm and l si = 8.6 cm.

Sensitivity of the EBM to snowpack parameters
In order to assess the sensitivity of the EBM outputs with respect to the different properties of a snowpack, a set of simulations were run for various snowpack structures.A random data set was generated corresponding to a snow height varying from 30 to 400 cm (SWE from 75 to 1000 mm with snow density set at 250 kg m −3 ).Measurements of the roughness parameters of air-snow interface and snow-ice interface are not available; therefore, empirical values for the correlation length l and the rms height σ from Oh et al. (1992) have been used.The values of σ as = 0.4 cm and l as = 8.4 cm, equivalent to a slightly rough surface, are used for the air-snow interface; however σ si = 0.9 cm and l si = 8.6 cm, corresponding to a rough surface, are chosen for the snow-ice interface due to the characteristics of ice beneath the snowpack over the study area.
The results of EBM simulations are plotted vs. SWE in Figs. 3 and 4. In Fig. 3, snow density is fixed at 250 kg m −3 , while the optical diameter is varied from 0.2 to 1 mm.The backscattering contribution at the air-snow interface, being inferior to −40 dB, is not represented here.As the SWE increases, the volume backscattering coefficient becomes more important until it reaches a value comparable to the snow-ice interface backscattering.The vertical dispersion of the volume backscattering represents the sensitivity of the EBM to optical diameter.Lowest values correspond to an optical diameter of 0.2 mm, whereas the highest ones correspond to an optical diameter of 1 mm.
In Fig. 4, where the optical diameter is fixed at 1 mm and snow density varies from 200 to 600 kg m −3 , the vertical dispersion of the volume backscattering represents the sensitivity of the EBM to snow density.By comparing Figs. 3 and 4, we can observe that the EBM is strongly sensitive to the optical diameter and moderately sensitive to the snow density.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), and wet snow mapping (Nagler and Rott, 2000).In general, these studies concentrate on inverting the EBM, which enables the retrieval of such snowpack properties.For a multilayer snowpack, the number of observations, i.e., the number of 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.

Introduction to data assimilation
The aim of variational assimilation is to integrate observational data with guess variables through the use of an observation operator.The method concentrates on searching for a solution that minimizes simultaneously the distance between observations and simulation results and the distance between initial guess variables and the analyzed variables.A schematic of this process is presented in Fig. 1.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): where x represents the set of variables describing the snowpack properties (here, density and optical diameter for each snow layer).The 1D-Var algorithm is based on the minimization of a cost function J(x), defined as where x is called the state vector, which can be modified after each iteration of the minimization, and x g is the initial guess of the state vector and remains constant during the whole process.Therefore, kx − x g k 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, ky − H ebm (x)k 2 represents the distance between simulated and observed radiometric quantities.The process also requires the estimation of the error covariance matrix of observations/simulations (R) and of the guess error covariance matrix (B).

Adjoint operator and minimization algorithm
In order to minimize the cost function J , one needs to calculate its gradient: If the model is denoted H ebm : B ! R, where B and R are the domain of definition of x and y, then the function r H t ebm satisfying 8 x, y, hr H t ebm y, xi B =hy, r H ebm xi 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.
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 until J is converged to its minimum: where r 2 J(x n ) is the gradient of second order (Hessian) of J :

Estimation of error covariance matrices
With preset air-snow interface and snow-ice interface parameters, the original model input vector x = [x Crocus x air-snow x snow-ice ] 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 optical diameter and the density of the ith layer of the snowpack.This means that the analysis process does not modify directly the thickness of each layer; however this parameter can be changed indirectly in the subsequent simulations by Crocus.At the first iteration of the algorithm, x is equal to 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 simulation, is a (2 n ⇥ 2 n) definite positive matrix.Each element of B is computed as where σ i and σ j represent the standard deviation of the errors on x i and x j , which have been experimentally estimated to 0.3 mm and 65 kg m −3 , respectively, for optical diameter and snow density.
The coefficient γ ij represents the correlation between errors on x i and x j and is modeled as These values are issued from an ensemble of slightly perturbed Crocus runs, obtained by varying their meteorological inputs over one winter season.The deviations between these runs, considered to be elementary perturbations, were then statistically studied and fitted with the model of Eq. ( 21) for each pair of variables.
In this case study, SAR data are only available for the HH channel; therefore the error covariance matrix R reduces to a scalar, deduced from the radiometric uncertainty of TerraSAR-X (0.5 dB) and the error of the EMB (inferred from the sensitivity of the EBM).The calculations at several altitudes over the Argentière glacier gives the average value of R = 0.03.

General comments on the chosen analysis process
In general, modeling techniques are used to establish the relationship between the physical properties of a natural environment and observations measured by 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 a 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 highly impractical.
Data analysis methods, on the other hand, require a vector of guess variables relatively close to the actual values.The snowpack variables calculated by Crocus are used as guess variables in our assimilation scheme.The fundamental goal is to modify the initial guess variables, while balancing the errors of the guess variables, 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.

Sensitivity of TerraSAR-X data
For this study, TerraSAR-X descending acquisitions over 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.Table 1 provides the main features of TerraSAR-X data sets.Figure 5 shows the location and a TerraSAR-X image of Argentière glacier captured on 6 January 2009.Meteorological forcing data provided by SAFRAN from 2400 to 3000 m altitude in steps of 100 m elevation on horizontal terrain were used to drive the detailed snowpack model Crocus throughout the entire season 2008-2009(starting on 1 August 2008)).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 a Frost filter (Frost, 1981) was applied using window size of 5 ⇥ 5 pixels.
In order to study the sensitivity of TerraSAR-X data to the changes in snow properties, Fig. 6 shows the comparison of TerraSAR-X backscattering coefficients (σ TSX ) on different dates at the altitudes of 2400, 2700 and 3000 m on Argentière glacier.For the period from 6 to 17 January (blue triangles) and from 8 to 19 February (red circles), it can be observed that the sets of comparison values are well below the equality line, which means the backscattering coefficients decreased between successive observations.The opposite effect  can be noted for the period from 17 January to 8 February (green crosses).The medium beneath the snowpack consists of glacier ice, and its roughness can be considered to be constant; therefore these increases and decreases in backscattering suggest that the σ TSX can be related to the modification of the snow condition.As can be observed in the snow precipitation chart on the bottom right, the green period has significantly more snowfall than the other two periods.

Simulation of Crocus snowpack data
The intrinsic parameters of a snowpack needed for EBM simulations are simulated by Crocus, which consist of a number of snow layers, their density, optical diameter, and thickness.These quantities are used as inputs for the simulation of the 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. 7.The optical thickness (⌧ ) is the product of snow depth and the extinction coefficient (Tsang et al., 2007).In the case of multilayer snowpack, it is defined as The Cryosphere, where the extinction coefficient  e is calculated using Eq. ( 9) and d the thickness of the snow layer.It can be observed that the snowpack stratigraphy provided by Crocus may be used to separate TerraSAR-X reflectivity at different altitudes.Figure 7 also indicates the sensitivity of the reflectivity to volume-related snow parameters.

Evaluation of the process and discussions
Crocus snow stratigraphic profiles were computed for seven different altitudes over the Argentière glacier, from 2400 to 3000 m.The level of liquid water content per volume (LWC v ) at the times and locations of analysis is 0 %; therefore the condition of dry snow is satisfied.Figure 5 shows the approximate locations of each study area on the glacier.
Figure 8 presents the implementation of the SAR data assimilation process into Crocus.The top part of the figure shows the Crocus simulation of snowpack without assimilation of SAR data.At instant t, Crocus simulates the snow stratigraphic profile from the previous state of snowpack (instant t − 1) and the meteorological data hourly provided from SAFRAN.The time lag between instant t and instant t − 1 is therefore one hour.We call this simulation "open loop".The bottom part of the figure shows the implementation of data assimilation into the execution of Crocus.Every 11 days, a TerraSAR-X acquisition is used to modify the snowpack stratigraphic profile of Crocus through an assimilation process.The snow profile before assimilation is called "guess", and the analyzed snow profile after assimilation is called "assimilated".Consequently, at the date of the first TerraSAR-X acquisition (6 January 2009), open-loop and guess profiles are identical.Once this first SAR acquisition is assimilated into Crocus, guess and assimilated profiles differ.This modification permits the constrainment of a physical snowpack simulation using external information acquired at different dates.Figure 9 shows the results of simulation and analysis using the TerraSAR-X time series from 6 January to 24 March 2009.The reflectivity on the glacier crevasse area (2600 m elevation) has a very high standard deviation due to the cracks and has therefore been masked.The red line corresponds to the TerraSAR-X reflectivities along the glacier, whereas the cyan diamond shape indicates the EBM simulations for the Crocus open-loop profiles.The blue triangles indicate the EBM simulation of the guess profiles.These guess profiles are in turn modified by the assimilation process to become the assimilated profiles.The EBM simulations of the assimilated profiles are shown in green circles.The assimilated profile is used to reinitialize Crocus for the next iteration, which then produces the guess profile for the next assimilation when a new SAR acquisition is available.
The agreement between TerraSAR-X reflectivity and the output of the EBM using Crocus simulated profiles can be observed in Fig. 9, where EBM simulations of assimilated profiles converge gradually over time toward the TerraSAR-X backscattering coefficient.The graph corresponding to 2 March 2009 shows that the convergence has been reached at all altitudes, as EBM simulations of guess and assimilated profiles are much closer to the TerraSAR-X measurements than the open-loop profiles.
Table 2 shows a comparison of root-mean-square error (RMSE) between simulated and measured reflectivities for different types of profile: open loop, guess and assimilated.It can be observed that the σ snow values converge gradually toward the σ TSX for the guess and assimilated profiles.At the last date of acquisition (24 March), the RMSEs for guess and assimilated profiles are below 1 dB, while the open-loop profile still gives an RMSE higher than 3 dB.
Figure 10 shows a detailed analysis of the modifications of the optical diameter and density of each layer due to data assimilation on 6 January, 8 February and 13 March 2009 at the altitude of 2400 m.It can be observed that the assimilation algorithm tends to modify the optical diameter and density in the deep layers which have a strong influence on the backscatter intensity and whose slight modification significantly reduce the discrepancy between TerraSAR-X  of optical diameter and snow density made by data assimilation also indirectly modify others physical properties of the Crocus-simulated snowpack.
These results show that we have combined three models (Crocus, EBM, adjoint model) and the TerraSAR-X data to constrain spatially and temporally the snowpack evolution.The use of data assimilation on SAR data to predict certain physical properties of snowpack has been developed in Nagler et al. (2008) and Takala et al. (2011).However, it is the first time that active X-band radar data have not been used directly to perform an assessment of snowpack properties but instead used to estimate physical parameters of each snow layer through a data assimilation algorithm.This algorithm needs to be further validated in the future using in situ measurements and advanced 3-D imaging techniques (Ferro-Famil et al., 2012).

Conclusions
This study presents a new system using data assimilation and a multilayer snowpack backscattering model based on the radiative transfer theory to constrain the evolution of a snowpack simulated by the snow model Crocus.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 1D-Var data assimilation based on the linear tangent and adjoint operator of the EBM, we are able to modify, in a physically consistent way, the snowpack profiles calculated using the snowpack evolution model Crocus.This process has been applied to a time series of TerraSAR-X images and Crocus simulations during the winter of 2008-2009 over the Argentière glacier.Results show that SAR data can be taken into account to efficiently modify the evolution of snowpack simulated by Crocus.This process can be further developed and used in real applications such as large-scale snow cover monitoring or snowpack evolution over a long period of time.
This system, however, does have some limitations, like the inability to simulate and assimilate under wet snow conditions due to the hypothesis used in the EBM.Another important hypothesis made in this study concerns the spherical shape of snow grains.On the one hand, this assumption highly simplifies the modeling problem but, on the other hand, prevents the simulations over cross-polarization channels (HV and VH).The discussion on how to resolve these limitations should be addressed in another study on the modeling of electromagnetic waves interactions with a snowpack.
Future studies will concentrate on calibrating the assimilation process using in situ measurements.Direct field measurements of the optical diameter using recently developed methods (Gallet et al., 2009;Arnaud et al., 2011) allow for a direct comparison to Crocus output (Morin et al., 2013).Future developments will also benefit from the recently finalized prognostic representation of optical diameter in Crocus (Carmagnola et al., 2014).

Figure 1 .
Figure 1.Global schematic of the data analysis used in this study.The inputs 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.
. Section 2 introduces the Crocus snowpack evolution model.Section 3 describes the DMRT electromagnetic backscattering model.The 1D-Var data assimilation method is presented in Sect. 4. Section 5 contains

Figure 2 .
Figure 2. Main backscattering mechanisms occurring within a multilayer snowpack obtained from the radiative transfer equation at first order (Longepe et al., 2009): air-snow reflection (M as ), volume scattering (M v ) and reflection over the snow-ice interface (M si ).

Figure 4 .
Figure 4. Test of EBM simulations on X-band, HH polarization for varying snow depth and density: snow density 200-600 kg m −3 , optical diameter 1 mm, and snow depth 30-400 cm.The glacier roughness is fixed at σ si = 0.9 cm and l si = 8.6 cm.
is the distance in centimeters between layer i and layer j .The values of ↵ and β depend on different types of correlations and can be split into three cases: correlation D − D: ↵ = 0.11 and β = 1; correlation ⇢ − ⇢: ↵ = 0.13 and β = 1; correlation D − ⇢: ↵ = 0.15 and β = 0.66.

Figure 5 .
Figure 5. Top panels: location of the TerraSAR-X acquisition in the French Alps.Bottom panel: a cropped image on the Argentière glacier area.The approximate positions of different altitudes on the Argentière glacier: 2400, 2700 and 3000 m on the TerraSAR-X images are indicated.The red line represents the continuous trail on the glacier where the SAR data will be used in the case study; the marks on this line delineate each 100 m of altitude.

Figure 6 .
Figure 6.Comparison of TerraSAR-X reflectivities between two different dates of winter season 2008-2009 at the altitudes of 2400, 2700 and 3000 m on Argentière glacier.The small graph on the bottom right shows the snow precipitation level for each period of comparison.

Figure 7 .
Figure 7. TerraSAR-X reflectivity plotted as function of optical thickness derived from Crocus output.Each point corresponds to a date of acquisition TerraSAR-X.

Figure 9 .Figure 10 .
Figure9.Results of simulation and analysis using eight TerraSAR-X acquisitions performed in winter 2009.σ TSX (red) are mean values obtained from the SAR images over the Argentière glacier (corresponding to the red line of Fig.5).σ sim (blue) represents the output of simulations using Crocus snowpack variables as inputs.Simulations obtained after data analysis are shown in green.Error bars show the standard deviation of the measured reflectivities.

Sensitivity and simulation of snowpack at X band 5.1 Study site: Argentière glacier
The area of interest covers the Argentière glacier (altitude: 2771 m; 45.94628 • N, 7.00456 • E).The size of the domain is approximately 5 km ⇥ 6 km.Over the glacier, altitude varies from 2400 to 3200 m, and the snowpack is essentially composed of dry snow.