Snow and Ice Thickness Retrievals Using GNSS-R: Preliminary Results of the MOSAiC Experiment

The FSSCat mission was the 2017 ESA Sentinel Small Satellite (S^3) Challenge winner and the Copernicus Masters competition overall winner. It was successfully launched on 3 September 2020 onboard the VEGA SSMS PoC (VV16). FSSCat aims to provide coarse and downscaled soil moisture data and over polar regions, sea ice cover, and coarse resolution ice thickness using a combined L-band microwave radiometer and GNSS-Reflectometry payload. As part of the calibration and validation activities of FSSCat, a GNSS-R instrument was deployed as part of the MOSAiC polar expedition. The Multidisciplinary drifting Observatory for the Study of Arctic Climate expedition was an international one-year-long field experiment led by the Alfred Wegener Institute to study the climate system and the impact of climate change in the Arctic Ocean. This paper presents the first results of the PYCARO-2 instrument, focused on the GNSS-R techniques used to measure snow and ice thickness of an ice floe. The Interference Pattern produced by the combination of the GNSS direct and reflected signals over the sea-ice has been modeled using a four-layer model. The different thicknesses of the substrate layers (i.e., snow and ice) are linked to the position of the fringes of the interference pattern. Data collected by MOSAiC GNSS-R instrument between December 2019 and January 2020 for different GNSS constellations and frequencies are presented and analyzed, showing that under general conditions, sea ice and snow thickness can be retrieved using multiangular and multifrequency data.


Introduction
FSSCat was the winner of the "2017 ESA Sentinel Small Satellite (S^3) Challenge" and the overall winner of the 2017 "Copernicus Masters" competition [1]. The mission consists of two "federated" 6U (10 × 20 × 30 cm 3

) CubeSats ( 3 Cat5/A and 3 Cat-5/B) in support of the Copernicus Land and Marine
Environment services [2], which were successfully launched on 3 September 2020 onboard the Vega SSMS (Small Spacecraft Mission Service) maiden flight in a Sun Synchronous Orbit of 550 km height. The goal of the ESA S^3 Challenge was to develop, in a record time of one year, a novel small satellite mission to complement the Sentinel family missions. One of the instruments onboard FSSCat is the Flexible Microwave Payload-2 (FMPL-2) [3], a dual Global Navigation Satellite System reflectometer (GNSS-R) and microwave radiometer instrument implemented using a Software Defined Radio.
The GNSS signals produce a coherent reflection over the ice [4] that allows for the precise determination of the ice edge [5]. Moreover, other studies show the potential of GNSS-R measurements to infer sea-ice thickness, when the scattering is coherent and the complex reflected signals are downloaded to ground [6,7], based on the reflection between the sea-ice and the sea water underneath. However, this hypothesis has not been confirmed by in-situ measurements, as sea-ice thickness has been only retrieved at L-band by radiometers [8,9].
As part of the validation of the FMPL-2 GNSS-R instrument, the PYCARO-2 (P(Y) & C/A ReflectOmeter-2), an evolved version of PYCARO, [10] was designed to take part in the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition [11,12], the largest campaign ever aiming to study the Arctic Ocean.
The objective of MOSAiC was to understand the evolving Arctic climate system and the role it plays in changing the global climate. MOSAiC was a multidisciplinary campaign to study the Arctic Ocean from different points of view, and involving different teams: Ice, Ocean, Atmosphere, Ecology, Bio-geochemical systems, Satellite Remote Sensing, Modeling, and Aircraft teams. Most of them deployed a large suite of sensors to measure different parameters, such as ice thickness evolution, snow properties, the topography beneath the ice floe, etc. Most of the instruments run for over a year on top of a drifting ice floe. The diversity of simultaneous measurements and the extension of the campaign in time makes it an exceptional opportunity to cross-check and use the data retrieved as model validation for other sensors such as scatterometers, GNSS-R, or radiometer sensors on satellites.
Here, we will focus on the MOSAiC remote sensing site and its ancillary measurements, which include: • Three Remote Sensing Sites (alternating measurements every 2-3 weeks), with the following instruments: L-, C-, X-, Ka-, and Ku-band microwave scatterometers, P-to L-band, C-, X-, Ku-, K-, and W-band microwave radiometers, IR and hyperspectral cameras, and two multiconstellation and multiband GNSS-R instruments. The GNSS-R instruments ( Figure 1) are part of a joint effort from the Institut d'Estudis Espacials de Catalunya (IEEC): Institute of Space Sciences (ICE, CSIC), and Universitat Politècnica de Catalunya (UPC) IEEC sections.
• Regular transects (1 km length over different sea-ice types) to measure total sea ice thickness, snow depth and density, Ku-and Ka-band radar backscatter, L-band radiometry, and additionally surface albedo with the returning insolation from spring onwards.
• Other ice, snow, ocean, and atmospheric measurements.
This work describes the data collected by one of the GNSS-R instruments developed by IEEC/UPC. Section 2 describes the instrument, the experiment set up, and the ground truth data used, Section 3 describes the Interference Pattern captured by the instrument, produced by the sea-ice floe, modeled as a four-layer dielectric, and Section 4 analyzes the results of the data retrieved by the PYCARO-2 IEEC/UPC instrument and its comparison to the four-layer model is presented in Section 3. This work describes the data collected by one of the GNSS-R instruments developed by IEEC/UPC. Section 2 describes the instrument, the experiment set up, and the ground truth data used, Section 3 describes the Interference Pattern captured by the instrument, produced by the sea-ice floe, modeled as a four-layer dielectric, and Section 4 analyzes the results of the data retrieved by the PYCARO-2 IEEC/UPC instrument and its comparison to the four-layer model is presented in Section 3.

Instrument Description and Experiment Setup
The GNSS-R instruments belonging to the Remote Sensing site are divided in two. One of them is the ICE/CSIC instrument, which is a GNSS-R instrument working at linear polarization to collect the reflected signals off the sea-ice. Its aim is to provide sea-ice properties from scattered GNSS signals. The other one is the IEEC/UPC instrument (named PYCARO-2), which operates at Left hand Circular Polarization (LHCP) and aims to support the reflectivity modeling activities in support to the FMPL-2 data processing.
Both instruments share a common power and data interface, physical enclosure, and the Right Hand Circular Polarization (RHCP) antenna used for time and geo-reference the instrument. The primary objective of both GNSS-R instruments is to provide ground-truth data for reflection modeling and originally also to support the calibration/validation of the FMPL-2 [3] onboard the FSSCat mission [1,2]. This manuscript focuses only on the data retrieved by the IEEC/UPC PYCARO-2 instrument [13], which is equipped with an antenna working at the same polarization (LHCP) as the FMPL-2 instrument onboard FSSCat.
The down-looking antenna is a LHCP choke-ring antenna, tilted 45 • downwards to maximize

Instrument Description and Experiment Setup
The GNSS-R instruments belonging to the Remote Sensing site are divided in two. One of them is the ICE/CSIC instrument, which is a GNSS-R instrument working at linear polarization to collect the reflected signals off the sea-ice. Its aim is to provide sea-ice properties from scattered GNSS signals. The other one is the IEEC/UPC instrument (named PYCARO-2), which operates at Left hand Circular Polarization (LHCP) and aims to support the reflectivity modeling activities in support to the FMPL-2 data processing.
Both instruments share a common power and data interface, physical enclosure, and the Right Hand Circular Polarization (RHCP) antenna used for time and geo-reference the instrument. The primary objective of both GNSS-R instruments is to provide ground-truth data for reflection modeling and originally also to support the calibration/validation of the FMPL-2 [3] onboard the FSSCat mission [1,2]. This manuscript focuses only on the data retrieved by the IEEC/UPC PYCARO-2 instrument [13], which is equipped with an antenna working at the same polarization (LHCP) as the FMPL-2 instrument onboard FSSCat.
The down-looking antenna is a LHCP choke-ring antenna, tilted 45 • downwards to maximize the power of the reflected signals, as most GNSS satellites reach ∼50 • elevation angle in polar regions. The antennas have a very smooth radiation pattern, as shown in Figure 2a and their configuration (height and pointing) can be seen in Figure 2b.
The instruments were shipped to Tromsø, Norway, on July 2019, from where the Polastern icebreaker carrying the MOSAiC expedition departed on 20

Ground-Truth Data
In the vicinity of the GNSS instrument, different types of measurements were collected of both snow depth, and sea ice thickness. Measurements were repeated every few weeks by the combination of a Magnaprobe snow probe and broadband electromagnetic induction sensors (Geophex GEM-2). The Magnaprobe was used to retrieve snow thickness, with a precision better than 1 cm [16]. The GEM-2 was used to measure the combined thickness of the snow and the sea ice, as described by [17]. Finally, both measurements are combined to retrieve snow and sea-ice thickness separately. In this

Ground-Truth Data
In the vicinity of the GNSS instrument, different types of measurements were collected of both snow depth, and sea ice thickness. Measurements were repeated every few weeks by the combination of a Magnaprobe snow probe and broadband electromagnetic induction sensors (Geophex GEM-2). The Magnaprobe was used to retrieve snow thickness, with a precision better than 1 cm [16]. The GEM-2 was used to measure the combined thickness of the snow and the sea ice, as described by [17]. Finally, both measurements are combined to retrieve snow and sea-ice thickness separately. In this study, we use the quick-view GEM-2 thicknesses based on the 18 kHz in-phase channel and processed directly onboard and quality controlled against manual drill-hole observations. We chose this channel out of the available frequency set of the GEM-2 for optimal retrieval of sea ice thickness with an accuracy better than 10 cm.
The snow depth was also measured using an avalanche probe in the area surrounding the GNSS-R instrument mast. Figure 3a shows the evolution of the mean and standard deviation of the snow thickness in the area surrounding the GNSS instrument. Figure 3b shows a broad distribution of Remote Sens. 2020, 12, 4038 5 of 20 snow depths as a Probability Distribution Function (PDF) and points to a highly inhomogeneous snow surface, with developed snow dunes in two different time instants (blue, 12 January 2020, and green, 7 February 2020). Finally, Figure 3c shows the PDF of the sea-ice thickness measurements.
It can be observed that, for the period between December 2019 and January 2020, the PDF and the avalanche probe measurements are consistent, which is not the case for February 2020. The instrument performance is compared to these measurements, from where consistent ground-truth data is available. Furthermore, ice core drillings were performed during the first days of December, showing that the ice core during the week of the 21st of December 2019, was around 0.8 m and 1 m thick. The same ice core drillings were conducted in the time period shown in Figure 3c to validate the measurements taken by the Magnaprobe.
Remote Sens. 2020, 1, 0 5 of 21 study, we use the quick-view GEM-2 thicknesses based on the 18 kHz in-phase channel and processed directly onboard and quality controlled against manual drill-hole observations. We chose this channel out of the available frequency set of the GEM-2 for optimal retrieval of sea ice thickness with an accuracy better than 10 cm. The snow depth was also measured using an avalanche probe in the area surrounding the GNSS-R instrument mast. Figure 3a shows the evolution of the mean and standard deviation of the snow thickness in the area surrounding the GNSS instrument. Figure 3b shows a broad distribution of snow depths as a Probability Distribution Function (PDF) and points to a highly inhomogeneous snow surface, with developed snow dunes in two different time instants (blue, 12 January 2020, and green, 7 February 2020). Finally, Figure 3c shows the PDF of the sea-ice thickness measurements.

Theoretical Background: IPT Applied to the Ice Floe
In navigation receivers, multipath is avoided as much as possible as it degrades positioning accuracy. For ground-based GNSS-R instruments, the signal collected by a GNSS antenna is a combination of the incident wave, i.e., the one coming directly from the GNSS satellite at RHCP and the signal reflected over a certain surface. The combination of the direct and reflected signals creates strong reflectivity fluctuations or fringes, called an "Interference Pattern". In the PYCARO-2 case, the reflection of the GNSS signal on the sea ice floe produced the same type of fluctuations.
The Interference Pattern Technique (IPT) was conceived by [18] to extract a number of geophysical parameters from the Earth by means of ground-based instruments. Most of its applications are soil moisture [19], vegetation height [20], water level [21] monitoring, sea state [22], or snow depth [23] estimation. However, most of those works are based on a three-layer model, with the first layer being air, a second thin layer of substrate (i.e., snow, crop field), and finally, a third layer which is a semi-infinite reflective surface (i.e., land). Interfering reflections from multiple layers of snow and ice in Antarctica were analyzed in [24], however, according to the authors' knowledge, the IPT has never been applied over an ice floe, which contains four dielectric layers: air, a snow layer on top of the sea ice, the sea ice itself, and the sea water underneath.
To study the Interference Pattern (IP) fringes created by the multiple reflections in the floe, the instrument retrieves the signal-to-noise ratio (or peak of the Delay-Doppler Map minus noise background measured a few lags before the correlation peak [25]) once every 10 s, collecting thousands of reflections from different bands and constellations. The total amount of reflections collected was 300.000 for Beidou B1D1, 125.000 for Beidou B2D1, 425.000 for Galileo E1C, 290.000 for Galileo E5b, 275.000 for GPS L1 C/A, and 235.000 for GPS L2 P(Y). The data used to perform this study was collected from 20 December 2019 to 27 January 2020. However, due to several operational procedures of the instrument, most of the data was obtained between 15 and 27 January 2020. This study is focused on the analysis of this period of data, from which results and methodology will then be applied to the rest of the campaign.

Four-Layer IPT Model: Theoretical Definition
The scattering geometry is modeled as shown in Figure 4, where the dielectric constant of each media as in Section 6.2 of [26], and in particular, the sea-ice dielectric constant model is based in [27], but for the direct model. This study is performed assuming the dielectric constant of the substrate is not varying. Furthermore, simulations have been performed to prove that the effect of varying the dielectric constant of the substrate produced a negligible effected as compared to the variation of either the snow or the ice thickness. In addition, the properties selected (i.e., snow density, ice type, and ice temperature) used to compute the dielectric constants shown in Figure 4 are the average values of what MOSAiC researchers' measured. Therefore, following the procedures described in [28], the power received by either the zenith-looking or the 45 • -looking antenna is proportional to (see Figure 2a).
where E i is the incident electric field, E r is the reflected electric field, F n (θ) is the antenna voltage pattern (amplitude and phase) for the up-looking and down-looking signals arriving to either the up-looking or down-looking antennas, with different patterns for each one (see Figure 2a), θ dir , and θ re f are the off-boresight arriving angles of the direct and reflected signals, θ i is the incidence angle of the signal arriving from a particular GNSS spacecraft, E 0 i is the direct signal, the one that would be received if there were no interference, R( r , θ) is the reflection coefficient as defined in (3), and ∆φ is the phase associated to the geometry. ∆φ is given by where λ is the electromagnetic wavelength, depending on the frequency band (L1 or L2), h is the antenna height, and θ is the satellite elevation angle. The reflectivity on top of the snow layer is computed iteratively as where r i,i+1 is the Fresnel coefficients between layers i and i+1, for the four-layer model being 1 = air, 2 = snow, 3 = ice, and 4 = saline water; S (see (4)) accounts for the reflectivity decrease due to the roughness at the interface [18], and Ψ is the phase produced by the reflection between the different layers, as detailed in (5) where t i+1 is the thickness of the i + 1 layer, r i is the permittivity of the i layer, and σ the surface roughness of the interface.
θ re f are the off-boresight arriving angles of the direct and reflected signals, θ i is the incidence angle of the signal arriving from a particular GNSS spacecraft, E 0 i is the direct signal, the one that would be received if there were no interference, R( r , θ) is the reflection coefficient as defined in (3), and ∆φ is the phase associated to the geometry. ∆φ is given by where λ is the electromagnetic wavelength, depending on the frequency band (L1 or L2), h is the antenna height, and θ is the satellite elevation angle. The reflectivity on top of the snow layer is computed iteratively as where r i,i+1 is the Fresnel coefficients between layers i and i+1, for the four-layer model being 1 = air, 2 = snow, 3 = ice, and 4 = saline water; S (see (4)) accounts for the reflectivity decrease due to the roughness at the interface [18], and Ψ is the phase produced by the reflection between the different layers, as detailed in (5) where t i+1 is the thickness of the i + 1 layer, r i is the permittivity of the i layer, and σ the surface roughness of the interface.

Four-Layer IP Model: PYCARO-2 Case
Each of the different surface permittivity values ( r i ) for this study have been modeled. The snow has been modeled as dry snow, as per in-situ measurements during the expedition, with a density of 296 kg/m 3 , and −25 • C [29]. The ice has been modeled as multiyear ice, with a mean temperature of −25 • C. Note that simulations have been carried out at different temperatures, and the impact of varying the mean temperature is negligible in front of the interference pattern produced by variations

Four-Layer IP Model: PYCARO-2 Case
Each of the different surface permittivity values ( r i ) for this study have been modeled. The snow has been modeled as dry snow, as per in-situ measurements during the expedition, with a density of 296 kg/m 3 , and −25 • C [29]. The ice has been modeled as multiyear ice, with a mean temperature of −25 • C. Note that simulations have been carried out at different temperatures, and the impact of varying the mean temperature is negligible in front of the interference pattern produced by variations on the layer thickness. Finally, the water layer has been modeled as saline water (32 psu) and at a temperature of −1.7 • C [29].

Interference Pattern in the RHCP Zenith-Looking Antenna
This study is focused on two frequencies: 1575.42 MHz, used by GPS L1 C/A code and Galileo E1C; and 1207.14 MHz, used by Galileo E5b and Beidou B2D1 signals. As a first example, the four-layer IP model simulation results are shown in Figure 5a for the up-looking RHCP antenna at 1575.42 MHz for different values of snow, ice thickness, and surface roughness. Note that the antenna heights used to Remote Sens. 2020, 12, 4038 8 of 20 simulate the model are the ones indicated in Figure 2b. Moreover, the contribution to the interference pattern from the reflected path to the up-looking antenna starts to vanish at elevation angles larger than 25 • due to the antenna pattern (see Figure 2a). Therefore, the model is only presented up to 30 • . Note that the presented IPT is the normalized version of P(θ) from Equation (1) This study is focused on two frequencies: 1575.42 MHz, used by GPS L1 C/A code and Galileo E1C; and 1207.14 MHz, used by Galileo E5b and Beidou B2D1 signals. As a first example, the four-layer IP model simulation results are shown in Figure 5a for the up-looking RHCP antenna at 1575.42 MHz for different values of snow, ice thickness, and surface roughness. Note that the antenna heights used to simulate the model are the ones indicated in Figure 2b. Moreover, the contribution to the interference pattern from the reflected path to the up-looking antenna starts to vanish at elevation angles larger than 25 • due to the antenna pattern (see Figure 2a). Therefore, the model is only presented up to 30 • . Note that the presented IPT is the normalized version of P(θ) from Equation (1)   As it can be seen, the position of the different peaks as well as their shape depends on the thickness of the snow layer. The position of both peaks and notches depends mainly on the snow thickness: from 10 to 15 cm thickness, the position of the notches shifts towards the left (lower elevation angles), but for a thicker snow layer (20 cm) the notch position moves towards higher elevation angles (e.g., the notch at 15 • of elevation in Figure 5a). In addition, there is also a small dependence on the ice thickness. Comparing the top (ice thickness = 1.2 m) and bottom (ice thickness = 1.6 m) plots in Figure 5a, the notch position is shifted slightly towards lower elevation angles, where thickness variation has less impact in the notch position.
A similar phenomenon happens at 1207.14 MHz (i.e., Galileo E5b and Beidou B2D1) where, depending on both the snow and the ice thickness in the four-layer model, the position where the minima occur fluctuates. As shown in Figure 5b, the four-layer model is also affected by the snow thickness in both the peaks' shape and the notch position. However, the ice thickness has a larger impact on the notch position as compared to the 1575.42 MHz case as, due to the longer wavelength, the penetration depth into the ice is larger.
Finally, Figure 6 presents a set of curves representing different values of both ice and snow thickness at 1575.42 MHz and 1207.14 MHz. The 1207.14 MHz signal is more sensitive to variations in the ice thickness than to snow thickness (10 to 15 cm), in terms of the positions of peaks and notches. In the 1575.42 MHz case, the nulls' position in the IP move ∼0.5 • in elevation, while in the 1207.14 MHz case they move ∼1.5 • .
Remote Sens. 2020, 1, 0 9 of 21 As it can be seen, the position of the different peaks as well as their shape depends on the thickness of the snow layer. The position of both peaks and notches depends mainly on the snow thickness: from 10 to 15 cm thickness, the position of the notches shifts towards the left (lower elevation angles), but for a thicker snow layer (20 cm) the notch position moves towards higher elevation angles (e.g., the notch at 15 • of elevation in Figure 5a). In addition, there is also a small dependence on the ice thickness. Comparing the top (ice thickness = 1.2 m) and bottom (ice thickness = 1.6 m) plots in Figure 5a, the notch position is shifted slightly towards lower elevation angles, where thickness variation has less impact in the notch position.
A similar phenomenon happens at 1207.14 MHz (i.e., Galileo E5b and Beidou B2D1) where, depending on both the snow and the ice thickness in the four-layer model, the position where the minima occur fluctuates. As shown in Figure 5b, the four-layer model is also affected by the snow thickness in both the peaks' shape and the notch position. However, the ice thickness has a larger impact on the notch position as compared to the 1575.42 MHz case as, due to the longer wavelength, the penetration depth into the ice is larger.
Finally, Figure 6 presents a set of curves representing different values of both ice and snow thickness at 1575.42 MHz and 1207.14 MHz. The 1207.14 MHz signal is more sensitive to variations in the ice thickness than to snow thickness (10 to 15 cm), in terms of the positions of peaks and notches. In the 1575.42 MHz case, the nulls' position in the IP move ∼0.5 • in elevation, while in the 1207.14 MHz case they move ∼1.5 • .

Interference Pattern in the LHCP Down-Looking Antenna
As the LHCP signal is received by the 45 • -tilted down-looking antenna, its reception is not limited by the lower back lobes of the radiation pattern of the RHCP antenna. For this reason, model predictions are presented for elevations from 5 • to 60 • . Figure 7a presents the IP at LHCP and 1575.42 MHz. Using the same parameters as for the RHCP signal case, the sensitivity to snow thicknesses is even larger than for the RHCP case (Figure 5a). Note that around 30 • and 45 • of elevation, the ripples produced by the IP have the same shape for different snow thickness, i.e., the curve for 10 cm has the same shape as that for 20 cm, both for an ice thickness of 1.2 m and 1.6 m.

Interference Pattern in the LHCP Down-Looking Antenna
As the LHCP signal is received by the 45 • -tilted down-looking antenna, its reception is not limited by the lower back lobes of the radiation pattern of the RHCP antenna. For this reason, model predictions are presented for elevations from 5 • to 60 • . Figure 7a presents the IP at LHCP and 1575.42 MHz. Using the same parameters as for the RHCP signal case, the sensitivity to snow thicknesses is even larger than for the RHCP case (Figure 5a). Note that around 30 • and 45 • of elevation, the ripples produced by the IP have the same shape for different snow thickness, i.e., the curve for 10 cm has the same shape as that for 20 cm, both for an ice thickness of 1.2 m and 1.6 m. Figure 7b presents the IP at LHCP and 1207.14 MHz. As compared to Figure 5b, the simulated IPs are now sensitive to both ice and snow thickness variations. Comparing both IPs (Figure 7b top and Figure 7b bottom) around 20 • of elevation, it can be appreciated that the valleys are modulated by the snow thickness, and their depth (amplitude) and shape depend on the ice thickness. Figure 7b presents the IP at LHCP and 1207.14 MHz. As compared to Figure 5b, the simulated IPs are now sensitive to both ice and snow thickness variations. Comparing both IPs (Figure 7b top and 7b bottom) around 20 • of elevation, it can be appreciated that the valleys are modulated by the snow thickness, and their depth (amplitude) and shape depend on the ice thickness. Comparing the overlaid curves from Figure 7a,b and for different ice thicknesses, in Figure 8 it can be clearly seen that the ice thickness variation does not affect the interference pattern produced between 30 • and 42.5 • (R2 in Figure 8) of elevation, although at low elevation angles, between 10 • and 30 • , the IP is sensitive to both snow and ice thickness. Note that, for the sake of simplicity, two ice Comparing the overlaid curves from Figure 7a,b and for different ice thicknesses, in Figure 8 it can be clearly seen that the ice thickness variation does not affect the interference pattern produced between 30 • and 42.5 • (R2 in Figure 8) of elevation, although at low elevation angles, between 10 • and 30 • , the IP is sensitive to both snow and ice thickness. Note that, for the sake of simplicity, two ice thickness are shown, but simulations have been carried out for the complete range between 0.6 m. and  In order to summarize, the IPT with the four-layer model is sensitive to changes in the snow and ice thickness in both RHCP and LHCP signals. However, at 1575.42 MHz, the RHCP has a lower sensitivity to ice thickness. In the 1207.14 MHz case, as the wavelength is longer and thus able to penetrate more into the ice, it shows more sensitivity to ice thickness variations. The LHCP case shows a mixture of both phenomena depending on the elevation angle: for low elevation angles, the IPT is affected by both snow and ice thickness, but for angles between 30 • and 42.5 • , the IPT mainly depends on the snow thickness. Note that this model does not take into account any antenna orientation, as antenna pattern of the received signal is compensated for.

Data Analysis
The PYCARO-2 instrument receives the direct and reflected GNSS signals. In this section, the waveforms' signal-to-noise ratio (SNR) are presented as radar plots for the different frequencies and constellations. The fringes produced by the reflection geometry can be easily detected by the angle of arrival of the signal (i.e., azimuth and elevation angles). Figure 9a shows the IPT corresponding to Galileo E1C code, and Figure 9b to the GPS L1 C/A code. It can be clearly identified, in both RHCP signals, where the reflections are coming from, located in the azimuth range from ∼180 • to ∼300 • . However, in the LHCP signal collected by the down-looking antenna, other interference patterns can be identified in both GPS and Galileo, most of them grouped around 220 • and 250 • . Moreover, at 1207.14 MHz (Figure 9c,d) Galileo E5B and Beidou B2D1 signals are also exhibiting an IP. It can be seen that the pattern can be appreciated around 210 • and 260 • of azimuth, whereas for other angles the fringes are slightly changed or even lost. In order to summarize, the IPT with the four-layer model is sensitive to changes in the snow and ice thickness in both RHCP and LHCP signals. However, at 1575.42 MHz, the RHCP has a lower sensitivity to ice thickness. In the 1207.14 MHz case, as the wavelength is longer and thus able to penetrate more into the ice, it shows more sensitivity to ice thickness variations. The LHCP case shows a mixture of both phenomena depending on the elevation angle: for low elevation angles, the IPT is affected by both snow and ice thickness, but for angles between 30 • and 42.5 • , the IPT mainly depends on the snow thickness. Note that this model does not take into account any antenna orientation, as antenna pattern of the received signal is compensated for.

Data Analysis
The PYCARO-2 instrument receives the direct and reflected GNSS signals. In this section, the waveforms' signal-to-noise ratio (SNR) are presented as radar plots for the different frequencies and constellations. The fringes produced by the reflection geometry can be easily detected by the angle of arrival of the signal (i.e., azimuth and elevation angles). Figure 9a shows the IPT corresponding to Galileo E1C code, and Figure 9b to the GPS L1 C/A code. It can be clearly identified, in both RHCP signals, where the reflections are coming from, located in the azimuth range from ∼180 • to ∼300 • . However, in the LHCP signal collected by the down-looking antenna, other interference patterns can be identified in both GPS and Galileo, most of them grouped around 220 • and 250 • . Moreover, at 1207.14 MHz (Figure 9c,d) Galileo E5B and Beidou B2D1 signals are also exhibiting an IP. It can be seen that the pattern can be appreciated around 210 • and 260 • of azimuth, whereas for other angles the fringes are slightly changed or even lost. Because of the antenna pattern of the 45 • down-looking antenna, most of the IPs are detected around 220 • and 250 • , where the fringes are stable at both frequencies. In order to compare the measurements from PYCARO-2 and the four-layer model described in Section 3, the signals have been filtered, and only those received in the 220 • -250 • of azimuth are used.
In order to ease the visualization of the signal in a 2D form, each SNR measurement retrieved by the instrument has been binned depending on the elevation angle of the transmitting satellite. Each measurement has been divided depending on the receiving antenna (i.e., RHCP or LHCP) and band (1575.42 MHz or 1207.14 MHz). Figure 10 shows, after calibration and azimuth filtering, a 2D histogram of the PYCARO-2 measurements overlaid with the RHCP (a) and the LHCP (b) four-layer model for different snow and ice thicknesses. In order to ease the representation and comparison of the interference pattern, the model curve has been rescaled to have its range in the same span as the PYCARO-2 measurement.
Note that, as seen in the previous section, both snow and ice thickness have an impact on the 1575.42 MHz band. However, it is not large enough to clearly distinguish which model parametrization is closer to the observed pattern. For the 1207.14 MHz case, the number of reflections received for low elevation angles is very small, and therefore, the ripples caused by the interference pattern cannot be clearly identified.
When moving to the LHCP case, as seen in Figure 10b, the retrieved signal is noisier at 1575.42 MHz, where there is only a small range of elevation angles (20 • to 35 • ) with ripples present. However, at 1207.14 MHz, the signal presents notches that can be easily identified as the ones described by the red curve in Figure 10b. In this case, the signal at 1207.14 MHz is able to penetrate deeper into the substrate (i.e., the snow), and therefore it is slightly less sensitive to snow variations. As shown in Figure 3b,c, the snow depth has a larger dispersion than the sea-ice thickness. In order to ease the visualization of the signal in a 2D form, each SNR measurement retrieved by the instrument has been binned depending on the elevation angle of the transmitting satellite. Each measurement has been divided depending on the receiving antenna (i.e., RHCP or LHCP) and band (1575.42 MHz or 1207.14 MHz). Figure 10 shows, after calibration and azimuth filtering, a 2D histogram of the PYCARO-2 measurements overlaid with the RHCP (a) and the LHCP (b) four-layer model for different snow and ice thicknesses. In order to ease the representation and comparison of the interference pattern, the model curve has been rescaled to have its range in the same span as the PYCARO-2 measurement.
Note that, as seen in the previous section, both snow and ice thickness have an impact on the 1575.42 MHz band. However, it is not large enough to clearly distinguish which model parametrization is closer to the observed pattern. For the 1207.14 MHz case, the number of reflections received for low elevation angles is very small, and therefore, the ripples caused by the interference pattern cannot be clearly identified.
When moving to the LHCP case, as seen in Figure 10b, the retrieved signal is noisier at 1575.42 MHz, where there is only a small range of elevation angles (20 • to 35 • ) with ripples present. However, at 1207.14 MHz, the signal presents notches that can be easily identified as the ones described by the red curve in Figure 10b. In this case, the signal at 1207.14 MHz is able to penetrate deeper into the substrate (i.e., the snow), and therefore it is slightly less sensitive to snow variations. As shown in Figure 3b,c, the snow depth has a larger dispersion than the sea-ice thickness. Remote Sens. 2020, 1,

Snow and Ice Thickness Retrievals: Results and Discussion
In order to compare the model and the 2D histogram, the median curve shown in Figure 11 has been computed over a running window with a width of 0.5 • in elevation angle. In order to compare those curves and the model, a nonlinear least square (NLS) minimization is performed in different steps. Note that the error function used for the process is the mean-squared error (MSE) of the SNR (dB).

Snow and Ice Thickness Retrievals: Results and Discussion
In order to compare the model and the 2D histogram, the median curve shown in Figure 11 has been computed over a running window with a width of 0.5 • in elevation angle. In order to compare those curves and the model, a nonlinear least square (NLS) minimization is performed in different steps. Note that the error function used for the process is the mean-squared error (MSE) of the SNR (dB). In order to properly retrieve snow and ice thickness, a first NLS minimization is performed for the LHCP signal at elevation angles between 30 • and 42.5 • , where the minimization parameter is the snow thickness. As seen in the previous section, the LHCP signal around 30 • and 42.5 • elevation angle is insensitive to ice thickness variations but not to snow thickness variations. In this case, the notch position of the different ripples can be used to estimate the snow thickness.
When the snow thickness is retrieved, a second NLS minimization is performed for the RHCP up-looking signal, between 5 • and 25 • of elevation. In this second case, the variable used for the second minimization is the ice thickness, using the snow thickness retrieved by the first NLS minimization. Note that in the RHCP case, both snow and ice thickness have an impact on the notch position, and therefore it is very advisable to first retrieve the snow thickness, as in the LHCP case the signal is mainly sensitive to the snow. Moreover, as the RHCP signal is coming from the up-looking antenna, the IPT can only be applied to very low elevation angles, where the antenna pattern has similar values for both the direct and the reflected signal. In this case the reflected signal is also RHCP, as the incidence angle is smaller than the Brewster angle.

Error Function Analysis and Ambiguity Removal
Both the snow and ice thickness IP produce periodic notches. Therefore, almost equal interference patterns are generated by different combinations of snow and ice thickness. To resolve this uncertainty, data from the two frequency bands (1575.42 MHz and 1207.14 MHz) can be used. The error function with respect to snow and ice thickness is shown in Figure 12. It can be clearly seen that the periodicity of the error function at 1575.42 MHz and 1207.14 MHz are different. In our case, there is a point where both error functions exhibit a minimum at the same time (blue dashed line). This point is the actual snow-ice thickness. In order to properly retrieve snow and ice thickness, a first NLS minimization is performed for the LHCP signal at elevation angles between 30 • and 42.5 • , where the minimization parameter is the snow thickness. As seen in the previous section, the LHCP signal around 30 • and 42.5 • elevation angle is insensitive to ice thickness variations but not to snow thickness variations. In this case, the notch position of the different ripples can be used to estimate the snow thickness.
When the snow thickness is retrieved, a second NLS minimization is performed for the RHCP up-looking signal, between 5 • and 25 • of elevation. In this second case, the variable used for the second minimization is the ice thickness, using the snow thickness retrieved by the first NLS minimization. Note that in the RHCP case, both snow and ice thickness have an impact on the notch position, and therefore it is very advisable to first retrieve the snow thickness, as in the LHCP case the signal is mainly sensitive to the snow. Moreover, as the RHCP signal is coming from the up-looking antenna, the IPT can only be applied to very low elevation angles, where the antenna pattern has similar values for both the direct and the reflected signal. In this case the reflected signal is also RHCP, as the incidence angle is smaller than the Brewster angle.

Error Function Analysis and Ambiguity Removal
Both the snow and ice thickness IP produce periodic notches. Therefore, almost equal interference patterns are generated by different combinations of snow and ice thickness. To resolve this uncertainty, data from the two frequency bands (1575.42 MHz and 1207.14 MHz) can be used. The error function with respect to snow and ice thickness is shown in Figure 12. It can be clearly seen that the periodicity of the error function at 1575.42 MHz and 1207.14 MHz are different. In our case, there is a point where both error functions exhibit a minimum at the same time (blue dashed line). This point is the actual snow-ice thickness.
In the case of snow thickness (Figure 12, top), the signal presents a periodicity every 10 cm, and therefore both signals will not present a minimum at the same time in the observation range (i.e., 5 to 35 cm). Therefore, the actual snow thickness can be estimated by overlapping the two normalized error functions. By looking to the sum of both normalized curves, the minimum value is the actual thickness (14.4 cm), which is very close to the measurements by the MOSAiC Remote Sensing team, which were 14.5 ± 0.2 cm.
As for the ice thickness case, as shown in Figure 12 (bottom), a minimum in the MSE function is present every ∼20 cm. Similar to the snow thickness case, the periods are not the same for each band. However, in this case, as the observation range (i.e., 0.5 to 2.5 m) is larger, there are different combinations of thickness values producing a minimum at the same time: 0.58 m. and 1.24 m. Since some a priori knowledge of the ice is required to resolve the ambiguity, we rely on the measured value of 1.21 m on 15 January 2020. Therefore, from both possible solutions, the correct value is likely 1.24 m.
Remote Sens. 2020, 1, 0 15 of 21 Figure 12. NLS minimization output for the data set between the 15th and the 27th of January, 2020. MSE is represented for both snow and ice thickness, for different signal frequencies and polarization, and it has been normalized (z-score normalization) to ease its visualization. Note that the exact place where both signals present a minimum at the same thickness is represented by the blue dashed line thickness and the ice thickness have an impact on the notch positions; therefore, it is very advisable to first retrieve the snow thickness. As in the LHCP case, the signal is mainly sensitive to the snow.
In the case of snow thickness (Figure 12, top), the signal presents a periodicity every 10 cm, and therefore both signals will not present a minimum at the same time in the observation range (i.e., 5 to 35 cm). Therefore, the actual snow thickness can be estimated by overlapping the two normalized error functions. By looking to the sum of both normalized curves, the minimum value is the actual thickness (14.4 cm), which is very close to the measurements by the MOSAiC Remote Sensing team, which were 14.5 ± 0.2 cm.
As for the ice thickness case, as shown in Figure 12 (bottom), a minimum in the MSE function is present every ∼20 cm. Similar to the snow thickness case, the periods are not the same for each Figure 12. NLS minimization output for the data set between the 15th and the 27th of January 2020. MSE is represented for both snow and ice thickness, for different signal frequencies and polarization, and it has been normalized (z-score normalization) to ease its visualization. Note that the exact place where both signals present a minimum at the same thickness is represented by the blue dashed line thickness and the ice thickness have an impact on the notch positions; therefore, it is very advisable to first retrieve the snow thickness. As in the LHCP case, the signal is mainly sensitive to the snow.

Model and Measured Signal Overlay
Once the thickness parameters have been retrieved, they are used as input for the four-layer IPT model. The modelled IPT is then overlaid on top of the PYCARO-2 data, as shown in Figure 13. Note that the IPT model has been normalized in range (i.e., added a bias and multiplied by a scale factor) in order to better compare it with the SNR values of PYCARO-2.
It can be noticed that the IPT curves are not 100% coincident with PYCARO-2 data in all elevation angles. Despite that, the estimated thickness values are extracted from an MSE minimization process, therefore giving a minimum error for all elevation angles with respect to the median curve used. It can be noticed that the IPT curves are not 100% coincident with PYCARO-2 data in all elevation angles. Despite that, the estimated thickness values are extracted from an MSE minimization process, therefore giving a minimum error for all elevation angles with respect to the median curve used.
Note that in the 1575.42 MHz RHCP signal the notches are also very clear, and therefore we can check that the IPT model and the median curve present the same notches at 6 • , 8 • , 12 • , 14.5 • , and 17.5 • of elevation angle. In the same way, the median curve and the IPT model for the 1207.14 MHz RHCP signal is presenting two notches at the same position, around 12 • and 15.5 • of elevation angle.
In order to validate the model for other data sets, another period of data has been used, between the 23rd and 27th of December, 2019. In this case, the snow thickness measured from MOSAiC researchers was 11-13 cm, and the estimated ice thickness by means of ice core drilling performed some days before, on December 20, 2019, was around 90 cm. Figure 14 shows the error functions for both ice and snow thickness. In this case, for the minimization process at LHCP, both curves exhibit a minimum at 12 cm of snow thickness, being consistent with the measurements carried out during the campaign. Furthermore, the error curves for RHCP present once more a combination of two values with a minimum of 0.79 m and 1.68 m. In this case, thanks to the knowledge of the actual ice floe thickness, the selected solution is 0.79 m. Note that in the 1575.42 MHz RHCP signal the notches are also very clear, and therefore we can check that the IPT model and the median curve present the same notches at 6 • , 8 • , 12 • , 14.5 • , and 17.5 • of elevation angle. In the same way, the median curve and the IPT model for the 1207.14 MHz RHCP signal is presenting two notches at the same position, around 12 • and 15.5 • of elevation angle.
In order to validate the model for other data sets, another period of data has been used, between the 23rd and 27th of December 2019. In this case, the snow thickness measured from MOSAiC researchers was 11-13 cm, and the estimated ice thickness by means of ice core drilling performed some days before, on 20 December 2019, was around 90 cm. Figure 14 shows the error functions for both ice and snow thickness. In this case, for the minimization process at LHCP, both curves exhibit a minimum at 12 cm of snow thickness, being consistent with the measurements carried out during the campaign. Furthermore, the error curves for RHCP present once more a combination of two values with a minimum of 0.79 m and 1.68 m. In this case, thanks to the knowledge of the actual ice floe thickness, the selected solution is 0.79 m.
The overlaid version is presented in Figure 15. In this case, the notches at 1207.14 MHz, both RHCP and LHCP, are well-matched, as well as the minima and maxima at 1575.42 MHz, LHCP. Despite that, the modeled IPT for the RHCP signal at 1575.42 MHz, and the actual measured data by PYCARO-2 has some mismatch. Even then, the snow and ice thicknesses selected are the ones presenting a minimum MSE between the median curve and the IPT curve.
Even though the estimated ice thickness in this second case differs by 10 cm from the actual measured value, the accuracy of the GEM-2 instrument is 10 cm as well. In addition, as pointed out by the MOSAiC researchers, the ice floe is not evenly growing. To prove that, underwater pictures were taken by a remotely operated vehicle (ROV) during MOSAiC expedition. Furthermore, it is known that small (e.g., 10 cm) and local (meter-scale) undulations are present during the sea-ice freeze up [30], and therefore a single GEM-2 measurement cross-validated using an ice drilling may not be representative of the overall sea-ice thickness but an approximation. It is important to note that this second example of measurements was conducted during December, which is the time period when the ice floe is patched together to form new ice. In this case, some regions may present thinner or thicker spots of ice, which can affect the GNSS-R data measured by PYCARO-2. This phenomena may not be present in the data presented for January, as the different ice sheets have been already merged together.
Remote Sens. 2020, 1, 0 17 of 21 Figure 14. NLS minimization output for the data set between the 23rd and the 27th of December, 2019. MSE is represented for both snow and ice thickness, for different frequencies and signal polarization, and it has been normalized (z-score normalization) to ease its visualization. Note that the blue-dashed line shows the exact thickness at which both signals present a minimum.
The overlaid version is presented in Figure 15. In this case, the notches at 1207.14 MHz, both RHCP and LHCP, are well-matched, as well as the minima and maxima at 1575.42 MHz, LHCP. Despite that, the modeled IPT for the RHCP signal at 1575.42 MHz, and the actual measured data by PYCARO-2 has some mismatch. Even then, the snow and ice thicknesses selected are the ones presenting a minimum MSE between the median curve and the IPT curve.
Even though the estimated ice thickness in this second case differs by 10 cm from the actual measured value, the accuracy of the GEM-2 instrument is 10 cm as well. In addition, as pointed out by the MOSAiC researchers, the ice floe is not evenly growing. To prove that, underwater pictures were taken by a remotely operated vehicle (ROV) during MOSAiC expedition. Furthermore, it is known that small (e.g., 10 cm) and local (meter-scale) undulations are present during the sea-ice freeze up [30], and therefore a single GEM-2 measurement cross-validated using an ice drilling may not be representative of the overall sea-ice thickness but an approximation. It is important to note that this second example of measurements was conducted during December, which is the time period when the ice floe is patched together to form new ice. In this case, some regions may present thinner or thicker spots of ice, which can affect the GNSS-R data measured by PYCARO-2. This phenomena may not be present in the data presented for January, as the different ice sheets have been already merged together. Figure 14. NLS minimization output for the data set between the 23rd and the 27th of December 2019. MSE is represented for both snow and ice thickness, for different frequencies and signal polarization, and it has been normalized (z-score normalization) to ease its visualization. Note that the blue-dashed line shows the exact thickness at which both signals present a minimum.

Conclusions
This work has presented the GNSS-R theory and techniques required to estimate the snow and the ice thickness of an ice floe. Multifrequency and dual-polarization GNSS-R can be used to determine both snow and ice thickness by means of the IPT. A peculiarity is that when nonlinear least-square

Conclusions
This work has presented the GNSS-R theory and techniques required to estimate the snow and the ice thickness of an ice floe. Multifrequency and dual-polarization GNSS-R can be used to determine both snow and ice thickness by means of the IPT. A peculiarity is that when nonlinear least-square minimization is applied to both RHCP and LHCP signals, it shows a periodicity in the error function, leading to multiple solutions. These uncertainties can be solved by comparing the error functions from different bands. Snow thickness can be easily estimated from this process, as the periodicity of the error with the snow thickness has the same order than the observation range (i.e., a period every 10 cm, observation range ∼40 cm). However, for the ice thickness case, the periodicity of the error with the ice thickness is one order of magnitude smaller than the observation range (i.e., a period every 20 cm, observation range ∼2 m), which leads to multiple solutions for ice thickness that, from the data available for PYCARO-2, cannot be resolved without a priori information. It has been seen that the snow thickness can be more easily retrieved, as it is the one inducing larger changes in the LHCP signal, which led to centimeter errors as compared to the in-situ measurements. However, this is not the case for the ice thickness, where ambiguities may appear, and a priori information is required. Furthermore, it has been seen that the error in the ice thickness retrieval for the second example differs 10 cm from the in-situ measurements.
The interference pattern produced by the four-layer model exhibits some sensitivity to ice thickness variations, especially at the lower frequency band. This means that signals do reflect off the bottom of the ice and reach the air media with still enough power to induce an interference with the direct signal. This behavior supports the hypothesis in [6], from which measurements of the sea ice draft could be possible with GNSS-R. However, there is a wide range of elevation angles where the snow thickness is dominant in the reflection (i.e., between 30 • and 42.5 • of elevation angle, the ice thickness cannot be measured). In this range of angles, the reflection takes place mostly in the snow-ice interface, while for other elevation angles, the interference pattern is produced in the air-snow, snow-ice, and ice-water interfaces.
From the results presented, it is clearly seen that the snow layer on top of the sea-ice is dominant, while the ice thickness seems to have a smaller impact to GNSS reflections. This phenomena can be also translated to other sensors, where the snow layer shall be properly modeled and taken into account to, for instance, provide an ice thickness measure.
The MOSAiC campaign ended on 30 September 2020, and data from different sensors will be released in the future. For this reason, the acquired data processed through the described algorithm can be used for validation of other MOSAiC instruments such as microwave radiometers. Furthermore, the combination of GNSS-R measurements and measurements from other instruments may increase the accuracy of the GNSS-R retrieved snow and ice thickness, thus helping to remove the ambiguities for different ice thicknesses.