Characterization of dry-snow sub-structure using GNSS reflected signals

Remote Sensing of Environment 124 (2012) 122–134

Contents lists available at SciVerse ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Characterization of dry-snow sub-structure using GNSS reflected signals Estel Cardellach a,⁎, Fran Fabra a, Antonio Rius a, Simone Pettinato b, Salvatore D'Addio c a b c

Institut de Ciències de l'Espai, ICE/IEEC-CSIC, Campus UAB, Ed. Ciències, Torre C5-parell, 08193 Bellaterra, Barcelona, Spain Nello Carrara Institute of Applied Physics (IFAC-CNR), Via Madonna del Piano, 10, 50019 Sesto Fiorentino (FI), Italy European Space Research and Technology Centre, ESTEC/ESA, Postbus 299, 2200 AG Noordwijk, The Netherlands

a r t i c l e

i n f o

Article history: Received 16 March 2012 Received in revised form 10 May 2012 Accepted 12 May 2012 Available online 8 June 2012 Keywords: GNSS reflectometry Dry snow Phase interferometry

a b s t r a c t Signals of Global Navigation Satellite Systems (GNSS) bi-statically reflected off pristine dry snow in DomeConcordia, Antarctica, were collected during an experimental campaign with the aim of investigating the potential capabilities of these signals to remotely sense dry snow areas. The data have been used to sense the sub-surface structure of the snow by means of a new holographic observable. The obtained reflected signals present interferometric patterns consistent with the coherent superposition of signals reflected in different layers of the snow sub-structure, down to a few hundred meters depth. A forward model has been developed to compare the amplitude and phase behavior of the interferences with the data. A new radio-holographic observable is also defined, to help identify the layer at which the internal reflections occur. This work represents the first step towards the potential use of a space-based GNSS-R mission for densely monitoring the depths of Antarctica's thick dry snow; and/or complementing L-band radiometric observations. © 2012 Elsevier Inc. All rights reserved.

1. Introduction The concept of GNSS-R (GNSS-Reflectometry) also known as PARIS (Passive Reflectometry and Interferometry System) was first suggested in 1993 in Martín-Neira (1993) as a multistatic radar technique to monitor ocean altimetry. The main advantages of PARIS are the larger number of simultaneous observations, giving wide swath capabilities, together with the possibility of achieving global coverage at a lower cost compared to active instrumentation. These advantages continue to increase as the updated GLONASS and the future GALILEO and COMPASS navigation constellations become available. On the other hand, the drawback of the concept arises from the fact that the signals were not initially conceived for that purpose: they are weaker than those obtained with dedicated active instruments, and the available code modulations (for civil uses) have relatively narrow bandwidth. Currently, some feasibility studies are being funded by the European Space Agency (ESA) to extend the PARIS concept for a future space-based mission (Martín-Neira et al., 2011). The altimetric capabilities of the GNSS-R have been tested in field campaigns (e.g. Cardellach et al., 2004; Lowe et al., 2002; Rius et al., 2002; Rius et al., 2010; Rius et al., 2011; Ruffini et al., 2004; Treuhaft et al., 2001). Besides the original altimetric purpose of the concept, a set of complementary purposes has also been proposed and experimentally checked, such as ocean roughness and surface

⁎ Corresponding author. Tel.: + 34 93 581 4358; fax: + 34 93 5814363. E-mail address: [email protected] (E. Cardellach). 0034-4257/$ – see front matter © 2012 Elsevier Inc. All rights reserved. doi:10.1016/j.rse.2012.05.012

wind estimates (e.g. Cardellach & Rius, 2008; Cardellach et al., 2003; Garrison et al., 2002; Zavorotny & Voronovich, 2000); soil moisture characterization (e.g. Masters et al., 2004); or sea-ice monitoring (Belmonte, 2007; Fabra et al., 2011; Semmling et al., 2011). Larson et al. (2009) determined that changes in snow depth surrounding geodetic GNSS stations can be clearly tracked in the corresponding multipath modulation of the amplitude of the GPS signal. In Larson's work, the change in the altitude of the reflecting surface is monitored by means of changes in the frequency of the interference between the direct radio-link, and the radio-link reflected off the snow–air interface. This technique could potentially be applied to a large network of geodetic GNSS ground stations to characterize the snow accumulation in their immediate surroundings (up to tens of meters distance from the site). Besides this work, the potential use of a GNSS-R approach for inferring other snow properties had not been checked. With this in mind, the GPS Sea-Ice and Dry-Snow (GPS-SIDS) project was an ESTEC/ESA contract to conduct two experimental campaigns at polar areas in which the Global Navigation Satellite System signals reflected off the ocean, sea-ice, and drysnow are used to infer geo-physical information. The goal of the project is to assess the applicability of the GNSS reflectometry technique to infer sea-ice and dry-snow properties. The first campaign was conducted in Greenland during winter 2008–2009, to focus on sea-ice products (Fabra et al., 2011; Semmling et al., 2011). This paper describes the work derived from the second campaign, intended to investigate the potential remote sensing capabilities of the GNSS-R over dry snow areas in Antarctica, a zone that would be densely covered by a future PARIS mission in a polar orbit.

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

After summarizing the main aspects of the experiment in Section 2, Section 3 describes some of the features detected in the reflected data. It includes the definition of a new observable, the lag-hologram, which uses spectral analysis to separate the different signal contributions collected by the receiver. The conclusion extracted from this preliminary analysis will be the hypothesis of the paper: the signals present patterns that result from interferences between reflections within the snow sub-structure layers. The hypothesis is evaluated under Section 4, where a forward model is presented, and then used to explain the interferometric features found in the data sets. Finally, the new technique is applied to the observations. Results and potential benefits of the new data products are presented in Section 5.

2. Experimental campaign and GNSS-R primary observables The GPS Open Loop Differential Real-Time Receiver, GOLD-RTR (Nogués-Correig et al., 2007) was used in this experiment. This is a dedicated open-loop receiver designed and manufactured at the Institute for Space Studies of Catalonia (IEEC). Two GNSS antennas were deployed: a zenith-looking Right-Hand Circular Polarized (RHCP) antenna for geodetic positioning, and a dual-polarized (RHCP + LHCP) horizon-looking antenna to collect the signals reflected off the snow. The dual-polarized antenna is a dual-frequency Antcom choke-ring antenna, hemispherical with b8 dBic peak gain at GPS L1 band, P/N: 123GM1215RL-A4A4-5NN-1. The GNSS antennas were mounted on top of the ∼ 45 m high American Tower at Dome-Concordia, Antarctica. This tower overlooks a wide area of pristine and well characterized dry snow (Fig. 1). Data were collected for twelve days in December 2009. The antenna points towards a flat horizon, and obstacles are not visible for hundred kilometers. A team from NelloCarrara Institute of Applied Physics (IFAC) installed and operated the equipments, and they also provided ancillary data, such as snow density profiles, snow temperature, and data obtained from the meteorologic station. Further information about the site and ancillary equipment can be found in Macelloni et al. (2006) and Macelloni et al. (2007). A low signal-to-noise (SNR) L-band electromagnetic field is transmitted by the GPS satellites orbiting at ∼ 20,000 km altitude. The signals are continuously transmitted at RHCP. Several codes modulate the signal by introducing 180° phase shifts. The Pseudo-Random

123

Noise (PRN) codes are designed to isolate the signal transmitted by one particular GPS space vehicle from the others received simultaneously. To achieve it, the codes are sequences of arbitrary phase jumps that result in a quasi-null cross-correlation function when correlated with any other transmitter PRN (orthogonal-like behavior). Therefore, the GPS receiver needs to cross correlate the signal captured by the antenna with replicas of the PRN codes, to discern, identify and separate them. These cross-correlation functions are also called waveforms. A further detailed explanation about GNSS/GPS can be found in e.g. Misra and Enge (2006). Unlike standard GNSS receivers, which give the observables at the estimated peak of the waveform, the GOLD-RTR system provides the complete waveform for each of the front-ends (zenith, LHCP-limb, and RHCP-limb antenna ports). The reflected signals are crosscorrelated with a model of the signal that uses phase and group-delay information obtained in real-time by the zenith-looking antenna (direct signals). This can be done because of the low altitude of the receiver above the snow surface, and the static conditions of the observations. Under such conditions, the delays and Doppler frequencies of the reflected signals are similar enough to the direct signals to operate the open-loop tracking using direct signal's information. The waveforms are not sampled in the continuous delay-space, but in time or delay intervals. We call these intervals lags. Each waveform is sampled in 64 complex delay-lags (in-phase/quadrature, I/Q), at 15-meter inter-lag resolution. Up to 10 waveforms can be captured simultaneously (different GNSS satellites from different front-ends), at 1 KHz rate. The system also provides standard navigation observables from the zenith antenna. Further details about the system and the data can be found in Cardellach et al. (2011). The campaign continuously ran for twelve days, from December 10 to 21, 2009, taking up to 129 GB of data every day. The observation window was set from 5 ∘ to 65 ∘ elevation and 255 ∘ to 30 ∘ azimuth, to constrain the measurements to an area of pristine snow. The ground tracks in polar coordinates can be seen in Fig. 2.

3. Properties of the reflected signals The reflected signals present coherence times longer than 1 s. This is a very long coherence time compared to GNSS signals reflected

Fig. 1. (Left) Photography of the GNSS antennas which participated in this experiment, mounted on top of the American Tower at Dome-Concordia, Antarctica, overlooking the horizontal pristine snow area (antenna pointing Nord-West). (Right) Map of the Concordia Station around the American Tower: the station buildings are at ∼ 800 m in the South-East direction from the American Tower. Another small building (Helene shelter) is at ∼ 100 m behind the antenna's pointing direction. Note that South points upward in this projection.

124

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

20

28 24

16

PRN

20

Waveform I/Q (a.u.)

32

10

0

−10

8 4 25m 50m 100m 300m

Fig. 2. Satellite ground tracks in polar coordinates, for one particular day of campaign (December 16, 2009). The color scale states for the GNSS transmitter identifier, the Pseudo-Random-Noise code number (PRN). The tracks repeat every day (except for a daily ∼4 min time shift; possible changes in the constellation–maneuvers, malfunctions, and receiver scheduling issues).

Waveform amplitude (a.u.)

12

20

15

10

5

0 10

20

30

40

50

60

Delay (15−m lag) Fig. 3. Amplitudes of a sequence of 1-s integrated complex waveforms collected with the GOLD-RTR receiver (PRN 13, December 16, 2009), between 44.5∘ and 45.5∘ elevation (only 1 out of every 10 waveforms are here shown, to avoid overloading the plot). (Top) In-phase and Quadrature components (in gray and black respectively). (Bottom) Total amplitude.

∼44.65 ∘, ∼44.9 ∘, or ∼ 45.28 ∘. We could think that those are effects of the noise (lower SNR levels at the end of the trailing edge, lag 37), however, some of these new peaks have significant SNR levels, and

Amplitude (a.u.)

20

10

0

20

Amplitude (a.u.)

off other types of surfaces. For instance, reflections off the rough ocean present several‐millisecond coherence only. Part of this long coherence time can be understood by the fact that the receiver crosscorrelates the reflected signals with a signal model that includes real-time information obtained by the direct radio-link. In other words, the receiver tends to stop the reflected signal using the dynamics of the direct one. The other factor to understand the long coherence is the fact that the snow surface is very smooth and the scattering is essentially specular. The snow surface topography in the area is essentially flat, with global slope b0.2°, and roughness characterized by very long auto-correlation length (compared to GPS electromagnetic wavelengths) and ∼1 cm RMS vertical dispersion (e.g. Petit et al., 1982; Six et al., 2004). That is, the surface roughness is not big enough to induce diffuse scattering, which would have introduced fluctuation in the phase. Specular reflections tend to generate waveforms with the shape of the signal modulation's auto-correlation function, with no further deformations (in opposition to diffuse scattering). The signal modulation used in this experiment, the civil available GPS C/A, has an autocorrelation function shaped as a triangle, of total time-width given by 2 ⋅ LC/A. The modulation chip length is LC/A ∼ 0.98 μs, which corresponds to ∼ 293 m in the space domain. The received waveforms, nevertheless, do not show the expected triangle shape, but a series of distorted triangles, with added tails and secondary peaks. Moreover, these shapes change gradually in time, as displayed in Fig. 3. From Fig. 3 we can see that the time series of the amplitude of the reflected field is not constant, but oscillates, and that these oscillations differ between distant lags. We have also found that the fluctuations build a pattern that repeats every day. An example of the amplitude oscillation pattern found is given in Fig. 4: a short time series is chosen to perceive the high rate components of the interference and their repeatability. It shows a sequence of 1-second integrated amplitudes for two different lags and four different days: lag 22, which approximately corresponds to the peak of the direct waveform (nominal zero delay); and lag 37, delayed by ∼ 225 m (see Fig. 3 to locate both lags within the waveform). The oscillation patterns at lag 22 and lag 37 present similarities, but do not perfectly match with each other. Further delayed lags tend to slightly increase the rate: double peaks appear sometimes where only one peak was detected in lag 22. For example, it happens around elevation

10

0 44.5

45.0

45.5

Elevation angle (deg) Fig. 4. Amplitude of lag-22 (top) and lag-37 (bottom) of the 1-s integrated waveforms, between 44.5∘ and 45.5∘ elevation, for PRN 13. Days 16 to 19 December 2009 have been plotted in different hues of gray. Note that some of these points are explicitly contained in Fig. 3-bottom (December 16).

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

a hint of them seemed to emerge in lag 22. Moreover, all days present the same patterns, thus suggesting that they cannot be just noise, but some signal. In the GNSS geodesy community the term multipath designates a particular type of reflections, near the receiving system (e.g. Byun et al., 2002; Elósegui et al., 1995). Multipath is usually an undesired effect which might mask or deteriorate the GNSS observables, because it interferes with the main ray generating oscillating patterns in the amplitude and the phase. The frequency of these patterns is given by fm ¼

−1 dρ λ dt

ð1Þ

where ρ is the delay between the reflected and the direct signals, expressed in units of length; and λ is the electromagnetic wavelength of the signal. Assuming planar horizontal reflectors, the phase with which the multipath signal reaches the receiver, measured in cycles with respect to the phase of the main signal, can be modeled as

holographic observable which uses each of the lags of the received waveform. The generation of the lag-hologram follows the steps below: • time series of N (typically we used 128) complex (I/Q) waveforms at 1 s sampling obtained from the front-end of the receiver connected to the horizon-looking antenna are taken: wf(τi, tj); • the phases of each τi lag of these waveforms are then counter-rotated by the phase of the direct signal. The direct signal is here defined as the peak of the waveform (lag-22, τ22) obtained by the front-end connected to the zenith-looking antenna: wf(τi, tj)e− iϕdir(τ22, tj); • a Fourier analysis (by FFT) is conducted independently on the time series of each lag: LH(τi, f) = FFT{wf(τi, tj)e − iϕ(τ22, tj)}; • since the geometric parameter that changes with time (and thus forces the potential interference to change) is the elevation angle, it is more practical to express the frequency in terms of elevation rate (oscillation cycles/degree-elevation) rather than Hz (oscillation cycles/second). The conversion between them is given by
f

ρ 2Hm Φm ¼ ¼ sinðeÞ λ λ

125

 cycle f ½Hz f ½Hz 2π ¼ de ¼ de deg−el ½ deg−el=s  ½rad=s 360 dt dt

ð3Þ

ð2Þ

being e the elevation angle of observation, and Hm the vertical distance at which the reflecting surface is located with respect to the receiving antenna. The multipath field sums coherently with the main field, with a rotated phase with respect to the main signal as given in Eq. (2). As the conditions change (either Hm or the elevation angle of observation e), the multipath phase in Eq. (2) changes too, introducing a rotation of the added field with respect to the main one. This phenomenon introduces oscillations in both amplitude and phase. In our experimental set-up, the only dynamic parameter is the elevation angle. As will be shown later, the rate of change of the elevation angle of observation is too low for a near-by reflector (or even from the shelter at the base of the tower) to introduce the scintillating patterns observed in the data. A near-by multipath phenomenon, thus, cannot be the source of interference. In the GNSS radio-occultation community, the term multipath, or tropospheric multipath, is applied to the phenomena that occur under certain atmospheric conditions, for which the GNSS signals split in several rays. A technique called radio-holography is then used in GNSS radio-occultations to identify and separate atmospheric multipath (e.g. Igarashi et al., 2000). Similarly, a new observable is next defined to shed some light on the source of the patterns. It is based on the holographic approach, but taking into account that different lags within the waveform might contain complementary information. 3.1. New observable: lag-hologram Radio holography uses coherent properties of the signals propagating through a medium. These properties arise due to the high stability of the GNSS signal and its high sensitivity to layered structures. This approach seeks to obtain the maximum spatial compression of the main ray separate from that of the other rays trajectories. This makes it possible to evaluate the intensity of radio waves at each ray trajectory and to determine the corresponding frequency displacement from a reference ray. A reference wave field (reference ray) is used to reveal the spectra from the total received field. We choose the direct signal (with no reflection) as a reference field, aiming to see the rest of possible contributing rays present in the data (atmospheric multipath, surface, sub-surface reflections … as sketched in Fig. 5). Once the reference field has been used to counter-rotate the phase of the reflected signal, a spectral analysis is performed. In radiooccultation applications, the holography is applied at the peak of the received waveform solely (because this is the only data provided by standard and radio-occultation receivers). We here present a new

• finally, each lag (τi) of the lag-hologram is normalized (independently to the other lags):

LH

norm

  τi ; f j ¼

  1 LH τi ; f j : ∑k ∥LH ðτi ; f k Þ∥

ð4Þ

• Another possible normalization would have been a single factor for the entire lag-hologram. We have chosen the lag normalization to give more relative power to the features at the end of the waveform, on those lags of the waveform where the overall power is weak, unmasking frequency contributions that otherwise would be too low compared to the values around the lag-hologram's peak power. The disadvantage of the lag normalization is that the weaker frequencies in powerful lags might also become masked. The resulting lag-holograms look like the example given in Fig. 6, where we can see the set of frequency components (rays) contributing to each lag. These contributions are asymmetric (only negative with respect to the direct signal frequency), grouped in a certain region of the frequency domain, and gradually taking more negative values as the lags further delay. The asymmetry of the spectrum is an important aspect of the spectral analysis, since it indicates that the frequency contributions are not induced by modulations in the signal power (which introduce symmetric features), but rotating phasors (consistent oscillations of phase and amplitude). Moreover, the lag-holograms present very high repeatability, day after day. This is displayed in Fig. 7, where the lag-holograms of observations occurring at the same geometric conditions as the example given in Fig. 6 are presented for the next four subsequent days (17 to 20 December 2009, PRN 13). The features are almost identical every day. Similar interferences and their repeatability are found in all observations during the campaign (every PRN, at every azimuth and elevation). 3.2. Preliminary conclusion: work hypothesis A reflection off the horizontal air–snow interface reaches the receiver with a phase (with respect to the direct signal) given by Eq. (2), where Hm is now the altitude of the receiver above the snow surface, Hr. The time derivative of Eq. (2) gives the interferometric frequency (in Hz) at which the oscillations occur, and the factor in Eq. (3) can be used to express it in units of cycle per degree-elevation. Then, f surf ½Hz ¼

−1 dρ −2H r de ¼ cosðeÞ λ dt dt λ

ð5Þ

126

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

Fig. 5. The lag-hologram observable is meant to separate and identify the different signals being captured by the horizon-looking antenna. This sketch contains some of them: snow external surface reflections; sub-structure reflections; multipath off external reflecting elements; and atmospheric multipath. The external reflecting elements could be the shelter at the base of the tower, or the tower itself. The atmospheric multipath is generated by strong gradients of the atmosphere, and usually inslant geometries (unlike this sketch, in low elevation angles of observation, here added for completeness).

f surf ½cycle= deg−el ¼

−2H r 2π : cosðeÞ 360 λ

ð6Þ

Fig. 6 corresponds to observations around 45 ∘ elevation. The frequency that corresponds to the snow surface is thus ∼−5.8 cycle/ deg-el. Any reflection off planar-elements above the snow surface 30

0.11

20

0.09 0.08

10

0.07 0.06

0 0.05

Power (a.u.)

Interf.Freq.(cycle/deg−elev.)

0.10

0.04

−10

0.03 0.02

−20

0.01 −30

0.00 8

12 16 20 24 28 32 36 40 44 48 52 56 60

would correspond to frequencies computed with Hm smaller than Hr, thus frequencies slower (less negative) than fsurf. Fig. 8 compiles the first negative peak of the lag-holograms observed during December 16 2009, which clearly follow Eq. (6). Therefore, it seems that the main bulge of frequency contributions generally captured by the lag-holograms, which are faster (more negative) than fsurf, might come either from depths below the snow external surface, or from reflectors which do not follow Eq. (2). Such reflectors could be tilted reflecting surfaces located above the snow, but at some horizontal distances. Fig. 1-right displays the environment of the observation tower: the Concordia Station buildings are at ∼ 800 m distance from the tower, in the direction opposite to the main beam of the antenna. These distances are too long to be captured within our waveform (42 lags after the direct signal, which is 630 m delay with respect to the reception of the direct radio-link). A small shelter is closer to our antennas, at ∼100 m in the back-lobe direction. The fact that it is in the blind area of the antenna, together with its small size (which hinders the generation of continuous multipath, from any incidence angle), makes it difficult to believe that it might generate the strong and continuous patterns observed in the data. Thus, our hypothesis is that the presence of multiple reflections within the dry snow sub-structural layers is captured in the data as interferometric patterns. A model is next presented to evaluate this hypothesis.

Delay (15−m lags) 4. Forward model: multiple-ray single-reflection (MRSR) Fig. 6. Example of lag-hologram for a time series of N = 128 1-s complex waveforms. The waveforms are the same ones used in Fig. 3. Here the frequency is given as interferometric cycles per degree-elevation. The zero-frequency corresponds to the reference ray, the direct one. According to Eq. (6), frequencies more negatives than − 5.8 cycle/deg-el correspond to scattering off reflecting-elements below the snow– air interface.

Rays may be reflected off both the external snow surface and internal snow boundaries. Given the clear multiple interference patterns observed in the data, we have neglected volumetric scattering (Wiehl et al., 2003), which would not produce interference patterns,

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

0

Interf.Freq.(cycle/deg−elev.)

Interf.Freq.(cycle/deg−elev.)

0

−10

−20

−30

−10

−20

−30 8

12 16 20 24 28 32 38 40 44 48 52 56 60

8

12 16 20 24 28 32 38 40 44 48 52 56 60

Delay (15−m lags)

Delay (15−m lags)

0

Interf.Freq.(cycle/deg−elev.)

0

Interf.Freq.(cycle/deg−elev.)

127

−10

−20

−30

−10

−20

−30 8

12 16 20 24 28 32 38 40 44 48 52 56 60

8

12 16 20 24 28 32 38 40 44 48 52 56 60

Delay (15−m lags)

Delay (15−m lags)

Fig. 7. Repeatability of the lag-holograms for different days: lag-holograms obtained with the same time series of observations (PRN 13, and elevation angles) as in Fig. 6, but for the next four subsequent days: 17 to 20 December 2009 (sequentially from left to right, and top to bottom). Note that the plots zoom into the range of frequencies of interest [−30,5], whereas Fig. 6 displayed the range [−30,30]. Color scale in Fig. 6.

Interf.Freq.(cycle/deg−elev.)

and focused on scattering off internal layers. We have taken a geometrical optics approach, where the different contributions are modeled as rays bouncing in different layers. Different aspects of the model we implemented are sketched in Figs. 9, and 10. We assume locally −2 −3 −4 −5 −6 −7 −8 −9 −10 10

20

30

40

50

60

Elevation (deg) Fig. 8. Averaged location of the first negative peak in the lag-holograms collected during December 16, 2009. The solid line is the theoretical interferometric frequency corresponding to a reflection of the snow–air interface, fsurf in Eq. (6). The dispersion (error bars) is consistent with the frequency resolution of the lag-hologram, given by the length of the time series (128 s) and Eq. (3): with de/dt in the range of values between 0.003 and 0.008 deg/s, the frequency resolution lies between 1.0 and 2.6 cycle/ deg − el.

horizontal layers, parallel incidence, and propagation/reflection through the snow layers following the Snell's law. The model considers a set of rays contributing to the total received signal, where each ray might suffer single-reflections solely, so we call it MultipleRays Single-Reflection (MRSR) model. The model was later tested with multiple reflections per ray, but the outcome was not satisfactory: the contributions generated with 3-way reflections within each layer resulted too weak and delayed, and just by themselves they did not match the data behavior as nicely as the single reflection model. The snow sub-structure is given by the vertical profile of the snow density, which can be transformed into vertical profile of the snow permittivity. We have followed the formulation in Ulaby et al. (1990). The real part of the relative permittivity, ′, is obtained by solving Eq. (7); then, its imaginary part, ″, can be calculated as in Eq. (8): ′ −1 ′ −1 ¼ νi ′ i ′ 3 i  i þ 2′ ″ ″  ¼ 3νi  i 

ð7Þ ′ 2

′

ð Þ ð2 þ 1Þ     i þ 2′ Þ ′ i þ 2 ′ Þ2 ′

ð8Þ

where ′i and ″i are the real and imaginary parts of the permittivity of the ice. These values have been set to 2.93 and 0.001 respectively in

128

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

zero-delay (lag 22 in our real data), and it represents the contamination by the direct signal being captured by the horizon-looking antenna. It should not be a big contribution, since this antenna is cross-polar; LC/A is the length pffiffiffiffiffiffiffiof ffi the GPS C/A triangle function, and j is the complex unity, j ¼ −1. Eq. (11) is thus the MRSR model, where Ui and ρi will be given by Eqs. (27) and (39) respectively, as explained below. 4.1. Amplitude of the i-layer contribution, Ui The next step in the MRSR model is to estimate the amplitude with which the ray reflected off of the bottom of the i-layer would reach the receiver. All the amplitudes will be normalized to the incident field at the surface level. We will take into account: • the transmission and reflection (Fresnel coefficients, T and R); and • the propagation attenuation within each layer (attenuation coefficient α). The equations are extracted from Ulaby et al. (1990) and modified to depend on the angle in the incident medium. This modification uses trigonometry, Snell's law, and assumes non-magnetic mediums: k1 sin θ1 ¼ k2 sin θ2

ð12Þ

pffiffiffiffiffi pffiffiffiffiffi 1 sin θ1 ¼ 2 sin θ2

ð13Þ

(where kj is the wave-number in medium j and θj the angle from the normal direction in the j-medium). The modified expressions read as follows: Fig. 9. Sketch of the procedure to build the total received complex waveform, once the amplitude Ui and the delay ρi of the contribution from each i-layer have been calculated (Sections 4.1 and 4.2 respectively). The phases ϕi depend on the delays ρi as indicated in Eq. (10). The signal modulation is accounted, by means of the triangle functions.

T⊥ ¼

our implementation (Ulaby et al., 1990). And νi is the volume fraction of ice in the snow, which depends on the snow density, ρ, as:

T∥ ¼

νi ¼

ρ : 0:916

ð9Þ

The following sub-sections will describe the equations of the amplitude Ui (Section 4.1) and the delay ρi (Section 4.2) with which a field that incises into the snow, propagates down to the i-layer, rebounds, and propagates back to the snow–air interface, finally reaches the receiver. The received complex waveform is not simply the complex sum of all these i-fields, but the modulation of the GNSS signals must also be taken into account. The i-layer contribution can be obtained as the GPS C/A auto-correlation function (triangle function) shifted by the delay at which the i-layer reflected signal arrives (with respect to the direct signal), ρi (Eq. (39)), and multiplied by the amplitude corresponding to the i-layer, Ui (Eq. (27)). Moreover, the i-layer contribution to the electromagnetic field at the receiver coordinates has a different phase than the phase of the direct ray contribution, that is: there is a rotation in complex space between the two contributions. Fig. 9 sketches this approach of building the total complex waveform. The phase is given by ρi : λ The total received waveform will thus be:

ϕi ¼ 2π

wf ðτ Þ ¼ wf

dir

X

N layers

ðτ Þ þ

i¼0

"

!

1:0−

jρi −τ j jϕ Ui e i LC=A

ð10Þ

2 cos θ1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : cos θ1 þ 2 − sin2 θ1 1

cos θ1 þ

ð11Þ

where τ is the delay-variable along the waveform; wfdir is the GPS C/A auto-correlation function (triangle function) centered at the nominal

2 cos θ1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi: 1 1 2 − 2

sin θ1

2

ð15Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 cos θ1 − 1 2 −ð1 sin θ1 Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi R⊥ ¼  1 cos θ1 þ 1 2 − 1 sin θ1 Þ2 :

ð16Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 cos θ1 − 1 2 −ð1 sin θ1 Þ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi R∥ ¼  2 cos θ1 þ 1 2 − 1 sin θ1 Þ2 :

ð17Þ

The transmitted GPS signals are essentially Right-Hand Circular Polarized (RHCP). In order to obtain the reflection and transmission coefficients for circularly polarized signals, the following combinations are required: T co ¼

1 ðT þ T ⊥ Þ: 2 ∥

T cross ¼

Rco ¼

1 ðT −T ⊥ Þ: 2 ∥

1 ðR þ R⊥ Þ: 2 ∥

Rcross ¼

#

ð14Þ

1 ðR −R⊥ Þ: 2 ∥

ð18Þ ð19Þ

ð20Þ ð21Þ

where sub-index co and cross relate to the circular co-polar and crosspolar transitions respectively. Both down-ward and up-ward propagating signal will need to be modeled in the internal reflection scenario. Our convention states that the incident layer/medium is the first sub-index, whereas the

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

129

Fig. 10. Sketch of the single-reflection approach implemented to model the snow internal reflections. The amplitude factors affecting each segment are plotted in blue, red for the rest of parameters used in the computations. R is the location of the receiving antenna, S is the specular point of a reflection off the snow–air interface. T is not included in the figure, it represents the GNSS transmitter position, very far away, so it illuminates the area with parallel incidence. The incident RHCP power at the surface level is set to 1.

medium into which the signal propagates is the second one (e.g. Tij is the transmission coefficient from i into j). Similarly, Rij refers to the coefficient of a reflection signal incident from medium/layer i off the interface with the j. Note that the opposite notation is also found in the literature. The coefficients given above relate the amplitude of the transmitted/ reflected field with respect to the incident field. In addition to these effects, the medium might also attenuate the signal. The attenuation constant α relates the total amplitude of the field after traveling certain length z across the medium, with an exponential decay given by e− αz (Ulaby et al., 1990). The attenuation depends on the imaginary part of the permittivity: α¼

pffiffi 2π Ima g λ

ð22Þ

(being  the relative permittivity). Given the little amount of signal transmitted in the cross-polar component (less than 0.1% of the power transmitted in the co-polar component), and the fact that only 1-reflection per ray is being considered in this model, we will not take into account: • the cross-polar component of the transmitted signals; and • the co-polar component of the reflected signals. • We neither consider geometrical loss due to the small differences between the distances traveled through the open-air by different rays.

Fig. 10 sketches the multi-layered snow with the amplitude factors affecting each segment of the propagation. As mentioned above, we will consider only Rcross − polar and Tco − polar. Therefore, the R in the figure and in the following analysis must be understood as Rcross − polar and the Ts as Tco − polar. This model is intended for the signal received with the LHCP antenna. To obtain the model for the signal received with the RCHP antenna we only need to replace Rcross by Rco in the analysis below. The LHCP ray reflected at the external surface will have an amplitude relative to the incident RHCP field as U 0 ¼ R01 :

ð23Þ

The one coming from a single reflection off of the bottom of the first internal layer will have an amplitude (w.r.t. incident RHCP field at the surface) −αd1

U 1 ¼ T 01 e

−αd1

R12 e

−2αd1

T 10 ¼ R12 T 01 T 10 e

:

ð24Þ

where d1 is the distance traveled inside the first layer: d1 = H1/cos(θ1). Note that T01 ≠ T10, but they both must be computed interchanging the corresponding angles and permittivities. The amplitude of the LHCP signal that has been single-reflected in the second layer would be: −αd1

U 2 ¼ T 01 e

−αd2

T 12 e

−αd2

R23 e

−αd1

T 21 e

T 10

ð25Þ

130

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

where ρint − 2 is the delay-contribution from the internal propagation through the layer-2

which can be expressed as: i¼2

−2αdi

U 2 ¼ R23 Π i¼1 T i−1i T ii−1 e

:

ð26Þ

It is straightforward to derive the amplitude of the contribution coming from the kth layer: −2αdi U k ¼ Rk;kþ1 Π i¼k i¼1 T i−1i T ii−1 e Hi : di ¼ cosðθi Þ

ð29Þ

ð30Þ

and A1 S ¼ D1 sinðθ0 Þ:

ð31Þ

D1 being the horizontal extent of the propagation inside the 1layer of snow (see Fig. 10): D1 ¼ 2H 1 tanðθ1 Þ:

ð32Þ

A compact way to express it is: ρ1 ¼ ρ0 þ ρ int−1 −D1 sinðθ0 Þ:

ð33Þ

Similarly, the ray that manages to propagate down to layer-2, get reflected off its bottom, and propagate upward to reach the receiver is: h i ρ2 ¼ ρ int−2 þ ρ int−1 þ ρSR þ ρTS −A2 S −ρTR

ð36Þ

ρ3 ¼ ρ int−3 þ ρ int−2 þ ρ int−1 þ ρSR þ ½ρTS −ðD1 þ D2 þ D3 Þ sinðθ0 Þ−ρTR

ð34Þ

ð37Þ

with ρ int−3 ¼ 2n3

H3 cosðθ3 Þ

ð38Þ

and D3 = 2H3 tan(θ3). Therefore, the general expression for the i-layer is:

ð28Þ

where H1 cosðθ1 Þ

and the distance between the specular point S and the point A2, A2 S, is:

being D2 = 2H2 tan(θ2). To complete these examples, the 3rd layer would read:

where subscripts TS, SR, and TR mean Transmitter–Specular reflection point, Specular reflection point–Receiver, and Transmitter–Receiver (direct radio-link), respectively. The ray which propagates into the first layer of snow, and gets reflected off of its bottom, is delayed with respect to the direct-link by ρ1. This delay has several contributions: (1) the one given by the internal propagation through the 1-layer, ρint − 1; (2) the distance from the specular point to the receiver, ρSR; (3) the distance from the transmitter to the point in which this ray enters the snow. This distance is equal to the distance between the transmitter and the specular reflection point, except for the segment between A1 and the specular point S; (4) finally, in order to reference the delay to the direct radio-link reception, we need to subtract the direct distance between the transmitter and the receiver, ρTR. This four contributions are summarized in the following equation:

ρ int−1 ¼ 2n1

ð35Þ

A2 S ¼ ðD1 þ D2 Þ sinðθ0 Þ

The second step of the model is to compute, for each i-layer, the delay of the ray such that manages to propagate down into the i-layer, is reflected off the bottom of that layer, and propagates upward towards the receiver. These delays, ρi, are given with respect to the direct reception of the signal (radio-link from the GNSS transmitter to the receiver with no reflection, ρTR in Fig. 10). As described below, most of the contributions to the i-delay are also common to the rays that have been reflected from the layers above it. Therefore, an iterative approach can easily solve the problem. In our notation, 0-layer is the external air, so the 0-delay corresponds to the reflection off the snow's external surface (point S in Fig. 10):

h i ρ1 ¼ ρ int−1 þ ρSR þ ρTS −A1 S −ρTR

H2 cosðθ2 Þ

ð27Þ

4.2. Delay of the i-layer contribution, ρi

ρ0 ¼ ρTS þ ρSR −ρTR

ρ int−2 ¼ 2n2

ρi ¼ ρ0 þ

k¼i X k¼1

k¼i X Hk − 2nk Dk cosðθk Þ k¼1

! sinðθ0 Þ

ð39Þ

being Dk ¼ 2Hk tanðθk Þ:

ð40Þ

4.3. Validation The model given above has been applied to a vertical profile of snow densities, provided by the IFAC team. The profile was measured in-situ at the experimental site, from either snow pits (first 10 m) or ice cores (EPICA community members, 2004; Macelloni et al., 2006) (down to 700 m depth). The vertical sampling of the profile is discrete in a non-homogenous set of depths, from sub-meter intervals in the upper part of the snow medium, to 1–2 m in the deeper parts. We have taken each of the depths of the discrete sampling as an interface of the layered snow structure upon which the model is applied. As shown in Fig. 11, the densities are very low at the surface layer, of the order of ∼ 0.3 g/cm 3, present several sharp transitions of ±0.2 g/cm 3 during the first 10 m depth, to gradually increase up to a saturation value of 0.92 g/cm 3 starting at ∼ 250 m depth. We have used this density profile for running the MRSR model and compare to our data, regardless of the location of the reflection and the time of the day, that is, we assume that the snow density is homogeneous across the observation area and has no time variations within the twelve days of experiment. An example of several waveforms synthesized with this model are compiled in Fig. 12. The series represents sequential complex waveforms (their real and imaginary parts) generated under identical geometric conditions as those real waveforms in Fig. 3. The two figures present similar features, folding and unfolding of the waveform defining oscillating peaks and troughs, consistent with interferometric patterns given by the layered-model. If anything, the model results are more conservative than the extreme features found in the real data. Varying the contamination level of the direct signal (wf dir in Eq. (11)) might improve this performance. The lag-hologram corresponding to 128-s observations around the geometric conditions of the synthetic data in Fig. 12 has been also produced. To do so, time series of 128 synthetic complex waveforms, which include the ones presented in the example above, generated

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

30

131

0.22 0.20

20

0.18 0.16

10

0.14 0.12

0.6

0 0.10 0.08

Power (a.u.)

Density (g/cm^3)

0.8

0.06

0.4

0.04 0.02 0

100

200

300

Fig. 11. Density profile to which the MRSR model has been applied. The model runs down to 700 m (density saturates at 0.92 g/cm3 after ∼250 m depth).

Waveform I/Q (a.u.)

under identical geometric conditions as those in Fig. 6, are analyzed by the FFT algorithm, applied to each lag-sequence. The resulting laghologram is shown in Fig. 13. The lag-hologram clearly shows a discrete set of interference frequency stripes, rather than a continuous or broad spectrum. The structure of stripes is not symmetric, indicating that the frequency components are phasor-rotations rather than amplitude modulations (which would generate symmetric positive and negative stripes). The zero-frequency corresponds to the direct signal, while ∼−5.8 cycle/deg corresponds to the theoretical interferometric frequency of a reflection off the snow external surface (from Eq. (6)). This figure can be compared to the lag-hologram in Fig. 6, computed with real data under identical geometry. The model relies on a discrete set of layers, given by the discrete sampling of the density profile. This may induce fake interfaces, just

0.10

0.05

0.00

Fig. 13. Lag-hologram of a series of 128 synthesized complex waveform around ∼ 45∘ elevation. The geometry of each of the 128 time series of synthetic observations is the same as each of the 128 real observations in Fig. 6.

because the permittivity is not continuously sampled. In spite of the similarity between our real data and synthetic data lag-holograms, we have run an exercise to identify the effects of the discretization of the density profile into the final holographic image. A smooth analytical expression of the density profile is used, resulting from an exponential fit to the ground truth discrete density profile. The exercise consists of generating two lag-holograms, one in which the smooth analytical profile is sampled at 1 cm-depth, and another run at 1‐m-depth intervals. Both synthetic lag-holograms presented a main contribution stripe coming from the reflection off the snow surface. The 1-cm discrete profile did not present any further holographic stripes, whereas the 1-m resolution one shown some weak stripes in those lags delayed more than 375 m with respect to the direct signal. These stripes, byproducts of the discretization, are very weak compared to the stripes we tend to find in real data. Therefore, we considered that the discretization effects were negligible. In this section we have shown that the MRSR model, applied to realistic snow density profiles, produces synthetic data with similar features as the real data. These similarities are observed both in the waveforms and in lag-holograms, and thus validate the hypothesis of the paper. Regarding the model, the effects of the discretization have also been checked. We conclude that they are much weaker than the signals, thus negligible. 5. Application

−0.05

Waveform amplitude (a.u.)

0.00 8 12 16 20 24 28 32 36 40 44 48 52 56 60

Depth (m)

0.10

0.05

0.00 10

20

30

40

50

60

Delay (15−m lag) Fig. 12. (Top) Real and imaginary parts of the complex waveforms synthesized using the model in Eq. (11). An amplitude of 0.05 has been assigned to the LHCP leakage of the direct signal. The geometric conditions reproduce those in Fig. 3. (bottom). Total amplitude of the waveforms on top.

Once we have evidence that the interferometric features found in the data can be explained by reflections of the signal inside the snow layered structure, we can propose the use of such signal features to sense the snow interiors. The concept is to associate the frequencies at which the lag-hologram present intense stripes to the different snow-layer depths that are mostly reflecting the L-band signals. The MRSR model can be used for converting interferometric frequencies into snow depths. To that end, we use the delays of the rays bouncing off each i-layer as provided by the model. These are shown in Fig. 14. This figure also indicates the limit of our system, which in the current set-up cannot capture delays greater than 630 m from the direct radio-link (42 15-m lags, number of lags after the direct signal, currently allocated at lag-22). This means that signals reflected off layers deeper than ∼350 m would be out of our detection window. Moreover, as seen in Fig. 3, power levels around lag 60 drop to nearly noise level. These are 38 15-m lags after the direct signal, or equivalently 570 m delay. This delay corresponds to depths between 200 and 300 m, deepest depths from which we are getting signal.

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

Because the incidence angle of the observation constantly evolves in time, the delay between the direct and reflected signals also changes in time. Similarly to the explanation given for the multipath, a receiver tracking the direct signal but also receiving a contribution from a reflection off the i-snow layer, does not lock on this other component because the latter arrives with a different frequency. This interferometric frequency (with respect to the direct one) can be obtained from Eqs. (1) and (3), where ρ is now the delay of the layer-reflected signal with respect to the direct one. The interferometric frequencies corresponding to each snow layer, deduced from values as in Fig. 14 and Eqs. (1) and (3), are displayed in Fig. 15. With this conversion method, we transform the frequencyaxis of the lag-hologram into a depth-axis, and thus the spectral stripes relate to the snow layers that reflect signal towards the receiver. The sensitivity of this relationship to uncertainties and errors in the vertical profile of snow density used in the model is analyzed in Appendix A. Appendix A shows that the vertical location of a particular frequency stripe is not very sensitive to non-modeled perturbations in the snow density profile, especially when the errors in the a-priori profile correspond to deep snow layers. A perturbation of 0.1 gr/cm 3 at 50 m depth, or at deeper layers, introduces less than 1% error in the vertical location of the stripes. As an intermediate step to identify the reflecting layers, we first integrate the lag-holograms along the lag-axis, to obtain the total spectral power as a function of the snow depth. An example is given in Fig. 16. The profile shows four clear echoes located at ∼5, 90, 130, and 240 m depth. The figure also displays the integrated spectra obtained with the MRSR model and the given density profile. Some–but not all–of the reflecting layers agree with the data. The discrepancies could be due to inaccuracies in the density profile assumed by the model, or by locally tilted interfaces (not considered in this initial model). Some snow layers appear consistently in many of the lag-holograms as reflecting elements. Those are identified in Table 1. However, in spite of the high temporal repeatability found in each GPS satellite data, the layers identified by different satellites are not always coincident, that is, the results present high temporal repeatability, but poor geographic consistency. This could be due to a diversity of causes, among them possible inhomogeneities in the snow across the scanned area, ∼500 m long; and the limited capability of the simple MRSR model to explain all the features in the data. To account for tilted layers in the sub-surface structure could be a possible way to improve the forward model.

50

Frequency (cycle/deg−el)

132

40 30

300 m 200 m

20

10 deg 100 m

10

80 deg 0 0

100

200

300

400

500 0

Snow depth (m)

30 60 90

Elevation (deg)

Fig. 15. A ray propagating into the snow, down to the x-meter depth layer, reflected at that layer, and propagated upward to reach the receiver, arrives with a y-frequency with respect to the direct signal (Eqs. (1) and (3)). It depends on the elevation angle of observation: gray hues, from 10∘ to 80∘ elevation, in steps of 10∘, lighter to darker respectively. The dependency on the elevation angle is also shown on the right frame, where the interferometric frequency in the signals reflected by layers 100, 200, and 300 m deep are plotted as function of the elevation angle. All values from the MRSR model.

The vertical resolution of the identified layers is mainly given by the length of the time series used to generate the spectral analysis. Other secondary factors are the geometry (see different slopes in Fig. 15 as function of the elevation angle). In this study, where 128 samples of 1-s integration each were used, the approximate vertical resolution is between 5 and 10 m. This could be improved by increasing the length of the time series, as shown in Fig. 16, which shows an example of enhanced vertical resolution obtained with lagholograms made out of 256 1-s samples. The integrated spectrum computed with 256 modeled waveforms resolves more layers than the model computed with 128. Unfortunately, our analysis of lagholograms made out of 256-samples was not fully satisfactory, presenting poorer agreement with the data. This is probably due to the significant geometric changes suffered by the observation along the event. The horizontal resolution of the measurement in the direction across the line-of-sight is of the order of 5–10 m, and it is given by the first Fresnel zone. The resolution along the line-of-sight direction is given by the displacements of the specular points in deep layers (related to Dk in Fig. 10 and Eq. (40)). Assuming that reflections might occur down to ∼ 300 m depth, the resolution along the line-of-sight is of the order of ∼350 m. Note that line-of-sight resolution worsens with depth.

1400

80 deg

Delay (m)

1200 1000

300 m

800

10 deg

600

200 m

400 200

100 m

0 0

100

200

300

Snow depth (m)

400

500 0

30 60 90

Elevation (deg)

Fig. 14. A ray propagating into the snow, down to the x-meter depth layer, reflected at that layer, and propagated upward to reach the receiver, arrives with y-delay with respect to the direct signal (according to the model in Eq. (39)). The delay depends on the elevation angle of observation: gray hues, from 10∘ to 80∘ elevation, in steps of 10∘, lighter to darker respectively. The dependency on the elevation angle is also shown on the right frame, where the delays in the signals reflected by layers 100, 200, and 300 m deep are plotted as function of the elevation angle. All values from the MRSR model.

Accumulated spectrum (a.u.)

1600 3.5 3.0

Data−128 Model−128

2.5

Model−256

2.0 1.5 1.0 0.5 0.0 0

50

100

150

200

250

300

Snow depth (meters) Fig. 16. (Circles) The frequencies in the lag-hologram shown in Fig. 6 have been converted into snow depth using Eqs. (1) and (3). This figure shows the sum of all its lag (integration of the lag-hologram along the lag-axis). Observation corresponding to PRN13, 16 December 2009, at 45∘ elevation. (Triangles) Integrated spectral power obtained with the MRSR model run in 128-steps, the geometry of which is identical to the geometry in the 128 samples used to generate the data in circles. (Inverted triangles) Same as the triangles, but using a series of 256 synthetic observations to improve the frequency resolution.

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134 Table 1 List of snow sub-structural layers that reflect signal towards the receiver producing interference patterns, as consistently appear in most of the lag-holograms at ∼ 45∘ elevation angle of observation. Depth (m)

10

70

130

133

MIUR, and ESA C.N. 22046/08/NL/EL and 20066/06/NL/EL. The authors are grateful to the Italian–French logistic team at CONCORDIA Station for their kind assistance during the experimental campaign. Estel Cardellach is under the Spanish Ramón y Cajal program.

240

6. Summary and conclusions

Appendix A. Sensitivity of the depth to inaccuracies in the snow density profile The technique presented in the paper identifies the sub-surface reflecting layers as frequency stripes in a holographic analysis. The spectral analysis is independent of any forward model. However, each frequency stripe is later assigned to a given depth of the reflecting layer. This step does require to run a model. In particular, an a-priori estimate of the snow density profile is needed. This appendix tries to understand how sensitive the vertical location of the layer is to inaccuracies of the snow density profile. To assess this question, a set of synthetic runs have been performed. Based on a realistic but smooth analytic expression of the density profile, a set of perturbations to this smooth profile has been added. A perturbation is here a Gaussian function added to the smooth profile, of relatively large intensity (0.1 gr/cm 3, representing more than 10% of the highest density) and vertical size (∼ 20 m half-width), and variable depth-location. It is illustrated in Fig. A.17top. For each of these perturbations, the model has been run to find the relationship between interferometric frequency and snow depth. This has been done for a geometry around 45 ∘ elevation angle. Different geometries would yield different results, but of similar order of magnitude. The error introduced by these uncertainties is given in Fig. A.17bottom. It is clear that even large uncertainties such as the Gaussian perturbations added in the smooth profile introduce errors below 2.5%. In 5 out of the 6 cases the error is lower than 1%.

Acknowledgments

Fig. A.17. (Top) A smooth analytical expression of the snow density is used to estimate a generic relationship between interferometric frequency and snow depth. The smooth profile is then perturbed by Gaussian bulks, of large intensity (> 10% of the highest density), ∼20 m vertical half-size, and located at different depths. (Bottom) Error in the snow depth location introduced by each Gaussian perturbation, with respect to the depths obtained with the non-perturbed reference profile.

This work has been funded by ESA C.N. 21793/08/NL/ST and AYA2008-05906-C02-02/ESP Spanish grants. The work conducted in Dome-C is partially covered by the Project on Glaciology of the PNRA-

a

1.1

Snow density [gr/cm^3]

Interference fringes found in GNSS-R observables collected during an experimental field campaign at Dome Concordia, Antarctica, December 2009, cannot be explained by near-by multipath or external reflecting elements. The patterns show strong temporal repeatability. All indications point to reflections off internal layers of the snow as sources of these interference fringes. A simple model, MRSR, of L-band scattering and propagation through sub-surface layers of snow has been implemented, and its qualitative agreement with the data is shown. The main conclusion is that the L-band GNSS signals penetrate into the snow, they are reflected by some of its layers, and the reflected components interfere with each other and with the direct radio-link to generate the observed beating patterns. The depth of penetration is estimated to be down to 200–300 m, and a few of the layers are more reflecting than the rest. These foremost reflecting layers are identified by means of a radioholographic technique, applied independently on each lag of the GNSS waveform, thus building what we have called a lag-hologram observable. Further delayed lags of the waveform contain information about signal being more delayed with respect to the reference–directsignal. Therefore, the complete lag-space must be inspected rather than just the peak of the signal, in order to obtain information from reflections occurring at deeper layers. While the lag-hologram identifies bright frequency stripes, the MRSR can also be used to link these frequencies into the depths of the reflecting layers that induced them. It has been shown that this translation is rather independent of inaccuracies in the snow density profile used by the model (perturbations in the density of the layers greater than 10% of the highest density, typically introduce 1% error in the vertical location of the layer). The sole knowledge of the most reflecting layers might be of interest to complement L-band radiometric experiments, such as the spaceborne SMOS and Aquarius missions. Despite the high temporal repeatability of the interferometric features found in the real data, the geographic consistency is rather low. Perhaps we are sensing inhomogeneities in the snow, but this has not been proven. Several other reasons that might explain this performance are contained in the validity of the model's assumptions: (1) locally horizontal snow layers, (2) rays internally reflected once solely, and (3) one single density profile for the entire area, ∼500 m wide. All these assumptions should be carefully revised, and improvements in the model could be made to attempt more complete retrievals of the snow sub-surface contents. For instance, estimates of the snow density or permittivity could be refined with the appropriate forward model; or tomographic approaches implemented to solve for the 3-D structures, including tilted layers and spatial inhomogeneities. Having proved that the GNSS signals are sensitive to snow substructures down to a few hundreds of meters depth is the first step toward enhanced applications that might be relevant to a possible future PARIS space-based mission (Martín-Neira et al., 2011). Such a mission, currently at the feasibility assessment level, will possibly follow a near-polar orbit, which would densely sample the vast and thick dry-snow cover of the Antarctica with 200–300 m deep observations.

1.0 0.9

Reference

0.8 0.7 0.6 0.5 0.4 0.3 0

50

100

150

200

250

300

350

Snow depth [meter]

Relative depth error [%]

b

2.5 2.0 1.5

Depth = 20m

Depth = 150m

Depth = 50m

Depth = 200m

Depth = 100m

Depth = 250m

1.0 0.5 0.0 −30

−25

−20

−15

−10

−5

Interf.freq.[cycles/degree−elev.]

134

E. Cardellach et al. / Remote Sensing of Environment 124 (2012) 122–134

References Belmonte, M. (2007), Bistatic scattering of global positioning system signals from arctic sea ice, Ph.D. thesis, University of Colorado. Byun, S. H., Hajj, G. A., & Young, L. E. (2002). Development and application of GPS signal multipath simulator. Radio Science, 37(6), 1098. http://dx.doi.org/10.1029/ 2001RS002549. Cardellach, E., Ao, C., de la Torre-Juárez, M., & Hajj, G. (2004). Carrier phase delay altimetry with GPS-reflection/occultation interferometry from Low Earth Orbiters. Geophysical Research Letters, 31(10), L10402. Cardellach, E., Fabra, F., Nogués-Correig, O., Oliveras, S., Ribó, S., & Rius, A. (Oct. 2011). GNSS-R ground-based and airborne campaigns for ocean, land, ice and snow techniques: Application to the GOLD-RTR datasets. Radio Science, 46(RS0C04). http: //dx.doi.org/10.1029/2011RS004683. Cardellach, E., & Rius, A. (2008). A new technique to sense non-Gaussian features of the sea surface from L-band bi-static GNSS reflections. Remote Sensing of Environment, 112, 2927–2937. http://dx.doi.org/10.1016/j.rse.2008.02.003. Cardellach, E., Ruffini, G., Pino, D., Rius, A., Komjathy, A., & Garrison, J. L. (2003). Mediterranean balloon experiment: Ocean wind speed sensing from the stratosphere using GPS reflections. Remote Sensing of Environment, 88(3), 351–362. Elósegui, P., Davis, J. L., Jaldehag, T. K., Johansson, J. M., Niell, A. E., & Shapiro, I. I. (1995). Geodesy using the Global Positioning System: The effects of signal scattering on estimates of site position. Journal of Geophysical Research, 100(B7), 9921–9934. EPICA community members (2004). Eight glacial cycles from an Antarctic ice core. Nature, 429, 623–628. Fabra, F., Cardellach, E., Rius, A., Ribó, S., Oliveras, S., Nogués-Correig, O., et al. (Nov 2011). Phase altimetry with dual polarization GNSS-R over sea ice. IEEE Transactions on Geoscience and Remote Sensing, 50(6). http://dx.doi.org/10.1109/TGRS.2011.2172797. Garrison, J. L., Komjathy, A., Zavorotny, V. U., & Katzberg, S. J. (2002). Wind speed measurement using forward scattered GPS signals. IEEE Transactions on Geoscience and Remote Sensing, 40(1), 50–65. Igarashi, K., Pavelyev, A., Hocke, K., Pavelyev, D., Kucherjavenkov, I. A., Matyugov, S., et al. (2000). Radio holographic principle for observing natural processes in the atmosphere and retrieving meteorological parameters from radio occultation data. Earth Planets and Space, 52(11), 893–899. Larson, K. M., Gutmann, E. D., Zavorotny, V. U., Braun, J. J., Williams, M. W., & Nievinski, F. G. (2009). Can we measure snow depth with GPS receivers? Geophysical Research Letters, 36(L17502). http://dx.doi.org/10.1029/2009GL039430. Lowe, S. T., Zuffada, C., Chao, Y., Kroger, P., Young, L. E., & LaBrecque, J. L. (2002). 5-cm Precision aircraft ocean altimetry using GPS reflections. Geophysical Research Letters, 29(10), 1375. Macelloni, G., Brogioni, M., Pampaloni, P., & Cagnati, A. (2007). Multifrequency microwave emission from the Dome-C area on the East Antarctic plateau: Temporal and spatial variability. IEEE Transactions on Geoscience and Remote Sensing, 45(7), 2029–2039. http://dx.doi.org/10.1109/TGRS.2007.890805. Macelloni, G., Brogioni, M., Pampaloni, P., Cagnati, A., & Drinkwater, M. R. (2006). DOMEX 2004: An experimental campaign at Dome-C Antarctica for the calibration of

spaceborne low-frequency microwave radiometers. IEEE Transactions on Geoscience and Remote Sensing, 44(10), 2642–2653. http://dx.doi.org/10.1109/TGRS.2006.882801. Martín-Neira, M. (1993). A Passive Reflectometry and Interferometry System (PARIS): Application to ocean altimetry. ESA Journal, 17, 331–355. Martín-Neira, M., D'Addio, S., Buck, C., Floury, N., & Prieto-Cerdeira, R. (2011). The PARIS ocean altimeter in-orbit demonstrator. IEEE Transactions on Geoscience and Remote Sensing, 49(6), 2209–2237. http://dx.doi.org/10.1109/TGRS.2010.2092431. Masters, D., Axelrad, P., & Katzberg, S. (2004). Initial results of land-reflected GPS bistatic radar measurements in SMEX02. Remote Sensing of Environment, 92, 507–520. http://dx.doi.org/10.1016/j.rse.2004.05.016. Misra, P., & Enge, P. (2006). Global Positioning System: Signals, measurements, and performance. Lincoln, Massachusetts: Ganga-Jamuna Press 01773, ISBN: 0-9709544-1-7. Nogués-Correig, O., Cardellach Galí, E., Sanz Campderròs, J., & Rius, A. (2007). A GPS-reflections receiver that computes Doppler/delay maps in real time. IEEE Transactions on Geoscience and Remote Sensing, 45(1), 156–174. Petit, J. R., Jouzel, J., Pourchet, M., & Merlivat, L. (1982). A detailed study of snow accumulation and stable isotope content in Dome C (Antarctica). Journal of Geophysical Research, 87(C6), 4301–4308. http://dx.doi.org/10.1029/JC087iC06p04301. Rius, A., Aparicio, J. M., Cardellach, E., Martín-Neira, M., & Chapron, B. (2002). Sea surface state measured using GPS reflected signals. Geophysical Research Letters, 29(23), 2122. Rius, A., Cardellach, E., & Martín-Neira, M. (2010). Altimetric analysis of the sea surface GPS reflected signals. IEEE Transactions on Geoscience and Remote Sensing, 48(4), 2119–2127. http://dx.doi.org/10.1109/TGRS.2009.2036721. Rius, A., Cardellach, E., Oliveras, S., Valencia, E., Park, H., Camps, A., et al. (June 2011). Altimetry with GNSS-R interferometry: First proof of concept experiment. GPS Solutions. http://dx.doi.org/10.1007/s10291-011-0225-9. Ruffini, G., Soulat, F., Caparrini, M., Germain, O., & Martín-Neira, M. (2004). The Eddy Experiment: Accurate GNSS-R ocean altimetry from low altitude aircraft. Geophysical Research Letters, 31(12), L12306. Semmling, A. M., Beyerle, G., Stosius, R., Dick, G., Wickert, J., Fabra, F., et al. (2011). Detection of Arctic Ocean tides using interferometric GNSS-R signals. Geophysical Research Letters, 38, L04103. http://dx.doi.org/10.1029/2010GL046005 4 pp. Six, D., Fily, M., Alvain, S., Henry, P., & Benoist, J. P. (2004). Surface characterisation of the Dome Concordia area (Antarctica) as a potential satellite calibration site, using Spot 4/vegetation instrument. Remote Sensing of Environment, 89(1), 83–94. http://dx.doi.org/10.1016/j.rse.2003.10.006. Treuhaft, R. N., Lowe, S. T., Zuffada, C., & Chao, Y. (2001). Two-cm GPS altimetry over Crater Lake. Geophysical Research Letters, 28(23), 4343–4346. Ulaby, F., Moore, R., & Fung, A. (1990). Microwave remote sensing active and passive, vol. I-III, : Artech House, Inc. Wiehl, M., Légrésy, R., & Dietrich, R. (2003). Potential of reflected GNSS signals for ice sheet remote sensing. Progress in Electromagnetics Research, 40, 177–205. Zavorotny, V. U., & Voronovich, A. G. (2000). Scattering of GPS Signals from the ocean with wind remote sensing application. IEEE Transactions on Geoscience and Remote Sensing, 38(2), 951–964.