Simulations of direct and reflected wave trajectories for ground-based GNSS-R experiments

The detection of Global Navigation Satellite System (GNSS) signals that are reflected off the surface, along with the reception of direct GNSS signals, offers a unique opportunity to monitor water level variations over land and ocean. The time delay between the reception of the direct and reflected signals gives access to the altitude of the receiver over the reflecting surface. The field of view of the receiver is highly dependent on both the orbits of the GNSS satellites and the configuration of the study site geometries. A simulator has been developed to determine the location of the reflection points on the surface accurately by modeling the trajectories of GNSS electromagnetic waves that are reflected by the surface of the Earth. Only the geometric problem was considered using a specular reflection assumption. The orbit of the GNSS constellation satellites (mainly GPS, GLONASS and Galileo), and the position of a fixed receiver, are used as inputs. Four different simulation modes are proposed, depending on the choice of the Earth surface model (local plane, osculating sphere or ellipsoid) and the consideration of topography likely to cause masking effects. Angular refraction effects derived from adaptive mapping functions are also taken into account. This simulator was developed to determine where the GNSS-R receivers should be located to monitor a given study area efficiently. In this study, two test sites were considered: the first one at the top of the 65 m Cordouan lighthouse in the Gironde estuary, France, and the second one on the shore of Lake Geneva (50 m above the reflecting surface), at the border between France and Switzerland. This site is hidden by mountains in the south (orthometric altitude up to 2000 m), and overlooking the lake in the north (orthometric altitude of 370 m). For this second test site configuration, reflections occur until 560 m from the receiver. The planimetric (arc length) differences (or altimetric difference as WGS84 ellipsoid height) between the positions of the specular reflection points obtained considering the Earth’s surface as an osculating sphere or as an ellipsoid were found to be on average 9 cm (or less than 1 mm) for satellite elevation angles greater than 10 , and 13.9 cm (or less than 1 mm) for satellite elevation angles between 5 and 10 . The altimetric and planimetric differences between the plane and sphere approximations are on average below 1.4 cm (or less than 1 mm) for satellite elevation angles greater than 10 ◦ and below 6.2 cm (or 2.4 mm) for satellite elevation angles between 5 and 10 . These results are the means of the differences obtained during a 24 h simulation with a complete GPS and GLONASS constellation, and thus depend on how the satellite elevation angle is sampled over the day of simulation. The simulations highlight the importance of the digital elevation model (DEM) integration: average planimetric differences (or altimetric) with and without integrating the DEM (with respect to the ellipsoid approximation) were found to be about 6.3 m (or 1.74 m), with the minimum elevation angle equal to 5 . The correction of the angular refraction due to troposphere on the signal leads to planimetric (or altimetric) differences of an approximately 18 m (or 6 cm) maximum for a 50 m receiver height above the reflecting surface, whereas the maximum is 2.9 m (or 7 mm) for a 5 m receiver height above the reflecting surface. These errors Published by Copernicus Publications on behalf of the European Geosciences Union. 2262 N. Roussel et al.: GNSS-R simulations increase deeply with the receiver height above the reflecting surface. By setting it to 300 m, the planimetric errors reach 116 m, and the altimetric errors reach 32 cm for satellite elevation angles lower than 10 . The tests performed with the simulator presented in this paper highlight the importance of the choice of the Earth’s representation and also the nonnegligible effect of angular refraction due to the troposphere on the specular reflection point positions. Various outputs (time-varying reflection point coordinates, satellite positions and ground paths, wave trajectories, first Fresnel zones, etc.) are provided either as text or KML files for visualization with Google Earth.

GLONASS and Galileo), and the position of a fixed receiver are used as input.Four different simulation modes are proposed depending on the choice of the Earth surface model (local plane, osculating sphere or ellipsoid) and the consideration of topography likely to cause masking effects.Angular refraction effects derived from adaptive mapping functions are also taken into account.This simulator was developed to determine where the GNSS-R receivers should be located to monitor efficiently a given study area.In this study, two test sites were considered.The first one at the top of the 65 m Cordouan lighthouse in the Gironde estuary, France, and the second one in the shore of the Geneva lake (50 m above the reflecting surface), at the border between France and Switzerland.This site is hidden by mountains in the south (orthometric altitude up to 2000 m), and overlooking the lake in the north (orthometric altitude of 370 m).For this second test site configuration, reflections occur until 560 m from the receiver.The planimetric (arc length) differences (resp.altimetric difference as WGS84 ellipsoid height) between the positions of the specular reflection points obtained consider-35 ing the Earth's surface as an osculating sphere or as an ellipsoid were found to be on average 9 cm (resp.< 1 mm) for satellite elevation angles greater than 10°and 13.9 cm (resp.< 1 mm) for satellite elevation angle between 5°and 10°.The altimetric and planimetric differences between the 40 plane and sphere approximations are on average below 1.4 cm (resp < 1 mm) for satellites elevation angle greater than 10°and below 6.2 cm (resp.2.4 mm) for satellite elevation angle between 5°and 10°.These results are means of the differences obtained during a 24h simulation with a com-45 plete GPS and GLONASS constellation, and thus depends on how the satellite elevation angle is sampled over the day of simulation.The simulations highlight the importance of the Digital Elevation Model (DEM) integration: average planimetric differences (resp.altimetric) with and without inte-50 grating the DEM (with respect to the ellipsoid approximation) were found to be about 6.3 m (resp.1.74 m) with the minimum elevation angle equal to 5°.The correction of the angular refraction due to troposphere on the signal leads to planimetric differences (resp.altimetric) about 18 m (resp.6

Introduction
The Global Navigation Satellite System (GNSS), which includes the American GPS, the Russian GLONASS, and the European Galileo (which is getting denser) uses L-band microwave signals to provide accurate 3-D positioning on any 75 point of the Earth surface or close vicinity.Along with the space segment development, the processing techniques have also improved considerably, with a better consideration of the various sources of error in the processing.Among them, multipaths still remain a major problem, and the mitigation of their influence has been widely investigated (Bilich A.L. , 2004).ESA (European Space Agency) first proposed the idea of taking advantage of the multipaths phenomenon in order to assess different parameters of the reflecting surface (Martin-Neira M. , 1993).This opportunistic remote sensing technique, known as GNSS-Reflectometry (GNSS-R), is based on the analysis of the electromagnetic signals emitted continuously by the GNSS satellites and detected by a receiver after reflection on the Earth's surface.Several parameters of the Earth surface can be retrieved either using 90 the time-delay between the signals received by the upper (direct signal) and the lower (reflected signal) antennas, or by analyzing the waveforms (temporal evolution of the signal power) corresponding to the reflected signal.This technique offers a wide-range of applications in Earth sciences.The time-delay can be interpreted in terms of altimetry as the difference of height between the receiver and the surface.Temporal variations of sea (Lowe S.T. et al. , 2002;Ruffini G. et al. , 2004;Löfgren J.S. et al. , 2011;Semmling A.M. et al. , 2011;Rius A. et al. , 2012) and lakes level (Treuhaft P. 100 et al. , 2004;Helm A. , 2008) were recorded with an accuracy of a few cm using in situ and airborne antennas.Surface roughness can be estimated from the analysis of the Delay-Doppler Maps (DDM) derived from the waveforms of the reflected signals.They can be related to parameters such as soil 105 moisture (Katzberg S. et al. , 2006;Rodriguez-Alvarez N. et al. , 2009, 2011) over land, or wave heights and wind speed (Komjathy A. et al. , 2000;Zavorotny A.U. et al. , 2000;Rius A. et al. , 2002;Soulat F. et al. , 2004) over the ocean, or ice properties (Gleason S. , 2006;Cardellach E. et al. , 110 2012).GNSS-R technique presents two main advantages: (1) a dense spatial and temporal coverage, not only limited to a single measurement point or a non repetitive transect as using classical GNSS buoys, (2) a guarantee of service for the next decades (because of the strategic role played by these 115 systems).GNSS-R altimetric accuracy is today at the level of few cm.But this technique will benefit, in the future, from improved processing technique and from the densification of the GNSS constellation.The commonly-used GNSS-R system consist of two antennas (figure 1): the first one is right-120 hand circular polarized (RHCP) and zenith-facing to receive the direct waves.The second one, left-hand circular polarized (LHCP) and nadir-facing to receive the reflected waves.These reflected waves will predominantly change their polarization from RCHP to LHCP by reflecting at near-normal 125 incidence.The reflected signals have an additional path delay with respect to the direct ones.The analysis of the path difference between these direct and reflected signals is used to estimate the relative height difference between the two antennas.In order to anticipate the impact of the geometric con-130 figuration of the experiment, a simulator has been developed to estimate the positions of reflection points using a specular reflection point assumption.Four different methods were implemented: approximating the Earth's surface as a local plane, as an osculating sphere, as an ellipsoid or integrat-135 ing a Digital Elevation Model (DEM).In addition, the signal bending due to the neutral part of atmosphere is taken into account using the Adaptive Mapping Functions (AMF) from Gegout et al. (2011) and made available by GRGS (Groupe de Recherche en Géodésie Spatiale).Simulations were per-140 formed for different configurations: variations in the reflectometer height, mask effects due to terrain, satellite network geometry.
This article is composed by three main parts following the logical structure of the figure 2. The first part presents the 145 datasets used for initiating simulations, the second one concerns the methodologies for the determination of the reflection points while the last one deals with the simulator performances and simulation results.

Design of the simulator 150
The simulator has been developed in the GNU R language, generally used for data processing and statistical analysis.A user manual and a description of the R language can be found on the website http://www.r-project.org/.The main interest of such a language remains in that it is distributed under GNU 155 GPL license which does R routines an open source program, available on various platforms (i.e.GNU/Linux, FreeBSD, NetBSD, OpenBSD, Mac OS and Windows).The simulator is composed by three main blocks (figure 2): an input block which contains the different elements manda-160 tory for the processing; a processing block where the user can choose which algorithm to be used, and an output block containing the different results of the simulation.
As inputs, this simulator requires the receiver coordinates, the satellite ephemeris and a set of optional environmental 165 parameters such as a DEM in order to take the possible masking of the terrestrial topography into account, as well as adap-tive mapping functions to integrate atmospheric delays and bending effects.
As outputs, the simulator provides the time-varying reflection point coordinates, but also various KML files (Keyhole Markup Language -standard format used by Google Earth) such as satellites positions and ground paths, waves trajectories and Fresnel first surfaces which can be opened using the Google Earth visualization tool.

GNSS orbit parameters
The simulations are based on the determination of the positions of the specular reflection points, once the receiver and the satellites positions are known Simulations are performed for a given receiver position in the WGS84 coordinates system and height above the ground.It is possible to apply an elevation or azimuthal angles mask to the simulations to avoid satellites with low elevation angle for instance.The elevation angle mask commonly used is set to 10°min and 90°max and no mask is set in azimuth.

SRTM Digital Elevation Model
The most realistic simulation needs the integration of a Digital Elevation Model (DEM) in order not to only take the possible masking of satellites into account, but to get more ac-200 curate and exact positions of the specular reflection points as well.The hole-filled version 4 of the Shuttle Radar Topography Mission (SRTM) DEM, with a spatial resolution of 90 m at the equator is used (Jarvis J. et al. , 2008).The altitudes are given with reference to the EGM96 geoid model.Uncertainty 205 on altitude is around 16 m over mountainous areas (Rodriguez E. et al. , 2005).It is made available by files of 5°x 5°f or land areas between 60°N and 60°S by the Consortium for Spatial Information (CGIAR-CSI): http://srtm.csi.cgiar.org/.

Earth Gravitational Model EGM96
In order to be able to convert between ellipsoidal heights (with respect to the WGS84 ellipsoid) and altitudes (with respect to the EGM96 geoid model) when producing KML files or when integrating a DEM, the knowledge of the geoid undulation is mandatory.In this study, we interpolate a 15 x 215 15-Minute Geoid Undulation Grid file derived from EGM96 model in a tide-free system released by the U.S. National Geospatial-Intelligence Agency (NGA) EGM Development Team: http://earth-info.nga.mil/GandG/wgs84/gravitymod/.The 220 error on the interpolation is lower than 2 cm (NASA and NIMA , 1998).

Adaptive Mapping Functions
The neutral atmosphere bends the propagation path of the GNSS signal and retards the speed of propagation.The range 225 between the satellite and the tracking site is neither the geometric distance nor the length of the propagation path, but the radio range of the propagation path (Marini J.W. , 1972).
For GNSS-R measurements, the tropospheric effects induced by the neutral part of the atmosphere are an impor-230 tant source of error.Indeed, GNSS-R measurements are often done at low elevation angle where the bending effects are maximal.Accurate models have to be used to mitigate signal speed decrease and path bending.It is commonly accepted to model tropospheric delays by calculating the zenith 235 tropospheric delay and obtaining the slant tropospheric delays with a mapping function.New mapping functions have been developed in the 2000's (Boehm J. et al , 2006;Niell A. , 2001) and significantly improve the geodetic positioning.Although modern mapping functions like VMF1 (Boehm J. 240 et al , 2006b) and GPT2/VMF1 (Lagler K. et al. , 2013) are derived from Numerical Weather Models (NWM), most of these mapping functions ignore the azimuth dependency which is usually introduced by two horizontal gradient parameters -in north-south and east-west directions -estimated 245 directly from observations (Chen G. et al. , 1997).More recently, the use of ray-traced delays through NWM directly at observation level has shown an improvement on geodetic results (Hobiger T. et al , 2008;Nafisi V. et al , 2012;Zus F. et al , 2012).The Adaptive Mapping Functions (AMF) are 250 designed to fit the most information available in NWM -especially the azimuth dependency -preserving the classical mapping function strategy.AMF are thus used to approximate thousands of atmospheric ray-traced delays using a few tens of coefficients with millimetre accuracy at low elevation 255 (Gegout P. et al. , 2011).AMF have a classical form with terms which are function of the elevation.But, they also include coefficients which depend on the azimuth to represent the azimuthal dependency of ray-traced delays.In addition, AMF are suitable to adapt to complex weather by chang- In order to assess the ocean tide influence on the positions of the reflection points estimated at an offshore experimental site located at the top of the Cordouan lighthouse (45°35'11"N ; 1°10'24"W), we use 24 hours of REFMAR (Réseau de Référence des Observations Marégraphiques) tide gauge observations, with a sampling frequency of 5 minutes.The tide gauge records of the station of Royan (45°37'14.07"N;1°01'40.12",located 12 km from the lighthouse) are the property of MEDDE (Ministère de l'Ecologie, du Développement Durable et de l'Energie), and they are available on the REFMAR website (http://refmar.shom.fr)).
3 Methodology : determination of the positions of reflection points The difference of phase between the two antennas (A-RHCP and B-LHCP on figure 1) at an epoch t for the i th GNSS satellite can be seen as a classical single difference between two receivers used for relative positioning as follows : where λ is the wavelength of the GNSS wavelength carrier, ∆φ i AB the measured carrier phase difference between the direct and received signals expressed in cycles, ∆δ i AB the difference in distance between the direct and received signals, ∆N i AB is the difference of phase ambiguity between the direct and received signals, c the speed of light in vacuum, ∆t AB the receivers clock bias difference.As the baseline between the two receivers is short (a few cm to a few tenth of cm), and in the case of low altitude of the receivers, both tropospheric and ionospheric effects are neglected due to the spatial resolution of the current atmospheric and ionospheric models.Besides, when both antennas are connected to same receiver, the receiver clock bias difference is also cancelled out.In this study, we only consider the difference in distance between direct and reflected signals as illustrated in figure 1.
The processing block contains four algorithms for determining the positions of the specular reflection points: the first considering the Earth as a local plane in the vicinity of the reflection point, the second as an osculating sphere, the third as an ellipsoid which corresponds to the WGS84 ellipsoid which has been expanded until the ellipsoid height of the receiver equals the height of the receiver above the reflecting surface (see subsection 3.3), and the last one uses the ellipsoid approximation but takes the Earth's topography into account: see figure 3. Comparisons between the different approximations of the Earth shape will be performed in subsection 4.1.All of them are based on iterative approaches to solve the Snell-Descartes law for reflection: the unique assumption is that the angle of incidence is equal to the angle of reflection on a plane interface separating two half-space media (a locally planar approximation is adopted when the surface is not everywhere planar).In the plane, sphere and ellipsoid approximations, the specular reflection point of a given satellite is contained within the plane defined by the satellite, the receiver and the center of the Earth.With regards to the DEM 320 integration, reflection can occur everywhere.In order to be able to compare the specular reflection points positions obtained by integrating a DEM, and to simplify the problem, we will only consider the reflections occurring within the plane, even while integrating a DEM.

Local plane reflection approximation
Refering to figure 3, let us consider the projection of the receiver R0 on the osculating sphere approximation (see subsection 3.2).We define the local plane P as the plane tangent to the sphere at R0.Let T 0 be the projection of the satellite 330 on P and R the symmetry of R0 relative to P .We look for the positions of the specular reflection points on P .Considering the Thales theorem in the triangles R SR0 and ST T 0, we have (see figure 3): And so:

Local sphere reflection approximation
The model we consider is an osculating sphere.Its radial direction coincides with the ellipsoidal normal and its center is 340 set at an ellipsoidal height equal to the negative value of the Gaussian radius of curvature defined as: with ϕ the latitude of the receiver, a and b the semi-major and semi-minor axis of the modified ellipsoid (see subsection 345 3.3).Please refer to (Nievinski F.G. and Santos M.C. , 2010) for further information on the different approximations of the Earth, particularly on the osculating sphere.
J. Kostelecky and C. Wagner already suggested an algorithm to retrieve the specular reflection point positions by 350 approximating the Earth as a sphere in (Kostelecky J. et al. , 2005;Wagner C., Klokocnik J. , 2003).Their algorithm is based on an optimized iterative scheme which is equivalent to make the position of a fictive specular point vary until verifying the first law of Snell-Descartes.A similar approach will 355 be used in this paper in the subsection 3.3 with the ellipsoid approximation.Here we chose to adopt a more analytical algorithm, first proposed by (Helm A. , 2008).In order to validate this algorithm, comparisons between it and the iterative one developed for the ellipsoid approach will be performed, by setting the minor and major axis of the ellipsoid equal to the sphere radius (see part 4.2.1).
Let us consider the vertical plane formed by the transmitter (GNSS) satellite (T), the receiver (R) and O', the centre of the Earth (figure 4).We assume that the specular reflection point (S) will be included in that plane.Let us consider the following orthonormal reference systems of coordinates: 1997), with O the centre of the Earth.WGS84 has Z polar and X,Y equatorial.The receiver and transmitter coordinates are known in this system.
-(O , x, y) R2 : a local two-dimensional system, obtained by the rotation of the (O,X,Y,Z) system around the Z axis, in such a way that x r = 0, and a translation 00 with 0 the center of the osculating sphere.
-(S, x , y ) R3 : a local two-dimensional system, obtained by a rotation around the z axis and a r E translation of the (O',x,y) system in such a way that x' and the local vertical are colinear, and that the system origin coincides with the specular reflection point S.
If H is the height of the receiver above the ground, the position of the receiver is: with r E the Gaussian radius of curvature at the latitude of the receiver ϕ r .
The position of the GNSS satellite transmitter considering ε the elevation angle of the satellite (considering zenith angle reckoned from the ellipsoidal normal direction) and τ the angle RT O is given by: Using the trigonometric sine formula in the R-T-0' triangle: We finally obtain: The Snell-Descartes law for reflection can be expressed as the ratios of the coordinates of the receiver and the transmitter in (S, x', y'): The coordinates in R3 can be derived from the coordinates 400 in R2 from: where γ is the rotation angle between the two systems (figure 4) .So (9) becomes: Following (Helm A. , 2008), we proceed to the substitution t = tan( γ 2 ), and (11) becomes: And finally becomes: with: Equation ( 13) is solved to determine the roots of this polynom using an iterative scheme based on the Newton method (Nocedal J. et al. , 2006).

Ellipsoid reflection approximation 420
We consider an ellipsoid corresponding to the WGS84 one extended such as the ellipsoid height of the receiver is equal to the receiver height above the reflecting surface.In other words, the WGS84 ellipsoid is expanded until its surface coincides with the reflecting surface, at nadir of the receiver 425 (surface base point).The cartesian coordinates of this surface base point must remain unchanged, when computed either from the original geodetic coordinates (λ, φ, h) W GS84 and ellipsoid constant (a, b) W GS84 on the one hand, or their N. Roussel et al.: GNSS-R simulations modified values (λ , φ , h , a , b ) on the other hand, where λ = λ is the longitude, φ = φ is the latitude, and h = 0.The ellipsoid thus remains geocentric and its axes are scaled as follows: where c = a 2 cos(φ This ellipsoidal extension is only done once as long as the receiver position remains unchanged with respect to the reflecting surface; it is redone if the reflecting surface changes (e.g., tidal waters) but is not with changes in the satellite direction.
We define the two normalized anti-incident r t and scattering r s vectors.When the Snell-Descartes law is verified, the sum of these two vectors (bisecting vector dr) coincides with the local vertical.The determination of the location of the reflection point is based on iterative process proposed earlier by (Gleason S. et al. , 2009), and enhanced with a dichotomy process.Let us consider three points on the ellipsoid: -S1 the projection of the receiver on the ellipsoid -S3 the projection of the transmitter on the ellipsoid -S2 the projection of the middle of [S1S3] on the ellipsoid We calculate dr, the correction in direction, considering the location of each the three points: We consider then the direction of the correction dr.If the correction is in the satellite direction, the sign is considered as positive, and negative if the correction is in the receiver direction.If the signs of dr S1 and dr S2 are different, it means that the specular reflection point is located between S1 and S2.We thus consider a new iteration with S1 = S1, S3 = S2 and S2 the projection on the ellipsoid of the middle of the new S1 and S3 points.We thus eliminate the part between the initial S2 and S3 points.Else if the signs of dr S2 and dr S3 are different, we consider a new iteration with S1 = S2 and S3 = S3 (and S2 the projection on the ellipsoid of the middle of the new S1 and S3 points).The iterative process stops when the difference between incident and reflected angle (with respect to the local vertical) is close to zero with a fixed tolerance of 10 −7 °.

Ellipsoid reflection approximation combined with a DEM
The two first approaches presented above are well adapted in the case of an isolated receiver, located on the top of a lighthouse, for instance.In most of the cases, the receiver is located on a cliff, a sand dune, or a building overhanging the sea surface or a lake.It can however be really appropriate and necessary to incorporate a Digital Elevation Model (DEM) into the simulations, in order not to only take the 480 mask effects (e.g., a mountain occulting a GNSS satellite) into account, but also to get more accurate and realistic positions of specular reflection points.The method we propose here consists of three steps later detailed in subsections 3.4.1,3.4.2and 3.4.3.485 1.A "visibility" determination approach to determine if the receiver is in sight of each GNSS satellite.
2. A determination of the specular reflection point position.
3. A "visibility" determination approach to determine if 490 the determined specular point is in sight from both receiver and satellite.
It is important to keep in mind that a DEM gives altitudes above a reference geoid.For consistency purpose, the positions of the receiver and the transmitter, and the DEM grid 495 points have all to be in the same reference system.So it is absolutely mandatory to convert the EGM96 altitudes from the SRTM DEM into WGS84 ellipsoidal heights by adding the geoid undulation interpolated from EGM96.

Visibility of the GNSS satellite from the receiver 500
This algorithm aims to determine the presence of mask between the receiver and the satellite.The visibility of the satellite and of the receiver, both from the specular point will be checked once the potential specular point position will be found.

505
Let R, S, and T be the locations of the receiver, the specular point and the satellite/transmitter on the ellipsoid.We interpolate the ellipsoidal heights along the path [T SR] with a step equal to the DEM resolution, with a bivariate cubic or bilinear interpolation.Cubic interpolation is used when 510 the gradient is big, linear interpolation otherwise.Tests show millimetric differences between cubic and linear interpolation for flat zones but can reach 1 m for mountainous areas.We thus obtain a topographic profile from R to T .For each segment of this topographic profile, we check if it intersects 515 the path [T R].If it does, it means that the satellite is not visible from the receiver.If not, we check the next topographic segment, until reaching the end of the path (i.e.T ).

Position of the specular point
Once the satellite visibility from the receiver is confirmed, 520 the next step consists in determining the location of the specular reflection point S along the broken line defined as in subsection 3.4.1.In order to simplify the process, we only consider the specular points located into the plane formed algorithm is similar to the one used for the ellipsoid approximation and is based on a dichotomous iterative process.
The segments formed by the points of the 2D DEM (see figure 5) are all considered susceptible to contain a specular reflection point.For each of this segment the sign of the correction to apply at each of the two extremities of the segment is checked following the same principle that for the ellipsoid approximation (see subsection 3.3), but with a local vertical component defined as the normal of the considered segment.If the signs are equal, no reflection is possible on this segment.Otherwise, we apply the dichotomous iterative method presented in subsection 3.3 until convergence with respect to the tolerance parameter (fixed to 10 −7 °).

Visibility of the determined specular reflection
point from the satellite and the receiver 540 Once the position of the specular reflection point is determined, we check if it is visible from the satellite and the receiver thanks to the algorithm presented in subsection 3.4.1.

Corrections of the angular refraction due to the troposphere 545
Our goal is to determine the location of the reflection point.
Only the angular refraction will be considered.The reflected minus direct range is left as future work.In order to correct the anisotropy of propagation of radio waves used by the GNSS satellites, we use AMF calculated from the 3-hourly delayed cut-off in model levels computed by the ECMWF (European Centre for Medium-Range Weather Forecasts).AMF tropospheric corrections were computed following (Gegout P. et al. , 2011) and provided by GRGS for this study.Given the geometric specificities of the specular reflection point, two paths have to be checked for propagation error: the first one from the satellite to the surface, and the second from the surface to the receiver.The main steps of the process are the following: 1 We consider the position of the specular reflection point 560 without any correction of the angular refraction; 2 We calculate the corrections to apply to this specular point knowing the incident and reflecting angle corresponding to the considered reflection point.We thus obtain a corrected incident angle.Figure 6 shows the correction to apply as a function of the elevation angle; 3 From the corrected incident angle, a corrected position of the specular point is calculated, making the reflecting angle being equal to the corrected incident angle; 4 With the new position of the specular point and to reach 570 a better accuracy of the point position, a second iteration is performed computing the corrections to apply to this new incident angle.

Correction of the satellite-surface path
First and foremost, the parallax problem for the wave emit-575 ted by a known GNSS satellite is solved.At first sight, the position of the specular reflection point calculated without any correction of the angular refraction is considered, given by the algorithm approximating the Earth's shape as a sphere given in paragraph 3.2.We use here AMF calculated from the 580 projection of the receiver on the surface, considering that the AMF planimetric variations are negligible for ground-based observations (i.e.we consider that we can use the same AMF for every specular reflection points, which is valid only if the specular reflection points are less than few tens of km 585 from the receiver and that the specular points lie on an equalheight surface).We thus obtain the corrected incident angle of the incident wave.Considering the law of Snell-Descartes, the reflecting angle must be equal to the corrected incident angle, for the specular reflection point position.590

Correction of the surface-receiver path
The aim here is to adjust the surface-receiver path to accommodate for the consequences of angular refraction.With the corrected reflection angle, we can deduce the corrected geometric distance between the reflection point and the receiver, 595 using this time AMF calculated from the receiver, assuming that the AMF altimetric variations are non-negligible (i.e. the part of the troposphere corresponding to the receiver height will have a non-negligible impact on the AMF).Considering the corrected geometric distance between the reflection 600 point and the receiver, the corrected position of the reflection point is obviously determined.It is indeed obtained as the intersection of a circle whose radius is equal to the correct geometric distance, with the surface of the Earth assimilated as a sphere, an ellipsoid, or with a DEM, depending on which 605 approximation of the Earth is taken into account.The whole process is iterated a second time to reach a better accuracy of the reflection point location.In fact, the first corrections were not perfectly exact since computed from an initially false reflection point location, and the second iter-610 ation brings the point closer to the true location.More iterations are useless (corrections to apply are not significant).Figure 6 shows an example of elevation corrections to apply as a function of the satellite elevations.This figure has been computed from simulations done on a receiver placed 615 on the Geneva Lake shore (46°24'30N" ; 6°43'6"E ; 471m): see subsection 4.1 page 8.

Footprint size of the reflected signal
The power of the received signal is mostly due to coherent reflection and most of scattering is coming from the first Fres-620 nel zone (Beckmann P. and Spizzichino A. , 1987).The first Fresnel zone can be described as an ellipse of semi-minor axis (r a ) and semi-major axis (r b ) equal to (Larson K.M. and Nievinski F.G. , 2013): With λ the wave length (m), h the receiver height (m) and the satellite elevation seen from the specular reflection point (rad) (i.e.corresponds to the reflection angle).
the shore of the Geneva lake (46°24'30N";6°43'6"E).This site is hidden by mountains in the South (orthometric altitude up to 2000 m), and overlooks the lake in the North (orthometric altitude of 370 m).
For both sites, precise GPS and GLONASS ephemeris at a 15-minute time-sampling come from IGS standard products (known as "SP3 orbit").

Validation of the surface models
Simulations were performed in the case of the Geneva Lake 645 shore, for a 24-hour experiment, on the 4 th october 2012.

Cross-validation between sphere and ellipsoid approximations
Local sphere and ellipsoid approximation algorithms have been compared by putting the ellipsoid semi-major and minor axis equal to the sphere radius.Planimetric and altimetric differences between both are below 6.10 −5 m for a receiver height above reflecting surface between 5 and 300 m and are then negligible.The two algorithms we compare are completely different: the first is analytical and the second is based on a iterative scheme and both results are very similar, which confirms their validity.

Cross-validation between ellipsoid approximation and DEM integration
The algorithm integrating a DEM has been compared to the 660 ellipsoid approximation algorithm by putting a flat DEM as input (i.e. a DEM with orthometric altitude equal to the geoid undulation).Results for satellite elevation angles above 5°a re presented in table 1.
As we can see in table 1, planimetric and altimetric mean 665 differences are subcentimetric for a 5 and 50 m receiver height and centimetric for a 300 m receiver height.However, some punctual planimetric differences reach 20 cm in the worst conditions (reflection occurring at 3449 m from the receiver corresponding to a satellite with a low elevation 670 angle), which can be explained with the chosen tolerance parameters but mainly because due to the DEM resolution, the algorithm taking a DEM into account approximating the ellipsoid as a broken straight line, causing inaccuracies.For a 50 m receiver height, planimetric differences are below 4 cm 675 (reflections occurring until 557 m from the receiver).With regards to the altimetric differences, even for reflections occurring far from the receiver, the differences are negligible (submillimetric).The map of the reflected points obtained for a big receiver height above the reflecting surface will in fact be the same 695 as the one obtained for a smaller receiver height, but more stretched.Henceforth, the higher the receiver height, the bigger the "measurable" area, but the less dense the ground coverage of the data (less reflection points per surface unit).

Study case: the influence of tides 700
As an illustration of a possible application of the simulator, tide influence on the position of the specular reflection points was assessed.Simulations at the Cordouan lighthouse were achieved integrating ocean tide from the tide gauge in Royan, by time-varying the receiver height above the sea surface in 705 order to simulate the tide.The vertical visibility mask was set to 10-90°, in order to avoid the weaker accuracy of determination of the specular reflection points positions for satellites with low elevation angle, as highlighted in subsection 4.3.2.By comparing the results with simulations made with a fixed-710 receiver height of 60 m above the sea surface, it appears that the 3D offsets reach values higher than 12 m for the maximum tide values (< 3 m) (figure 10).We can expect even higher discrepancies by taking into account satellites whose elevation angle would be lower than 10°.

Geneva Lake
Three sets of simulation have been performed in the case of the Geneva Lake shore, for a 24-hour experiment, on the 4 th october 2012: first configuration considering a receiver height of 5 m above lake level second configuration considering a receiver height of 50 m above lake level third configuration considering a receiver height of 300 m above lake level as for an airborne experiment (e.g.hovering helicopter).
Each series has been computed using the four algorithms of determination of the reflection points (local planimetric approximation, local osculating sphere approximation, ellipsoid approximation and the algorithm taking a DEM into ac-730 count).Results are presented on figures 11 to 14 and in table 2. They show the distances between the specular points and the receiver (arc lengths), and the differences between the positions given by each algorithm.
Influence of the receiver height above the reflecting sur-735 face It appears that both planimetric and altimetric differences between the method used increase with the receiver height above the reflecting surface.This is explainable by the fact that the higher the receiver is, the farther the reflection points will be from the receiver, and the bigger the impact of the Earth approximation will be.For a 5-m receiver height, reflection occurs up to approximately 60 m from the receiver, whereas for a 300-m receiver height, it occurs up to 3400 m (6700 m when integrating the DEM).It means that, in the second case, reflections occur in the mountains in the south of the receiver hence big differences between the sphere algorithm and the algorithm taking the DEM into account.For a 5 m receiver height above the reflecting surface and considering satellites with elevation angles above 5°, mean planimetric differences are below 1.3 cm between the osculating sphere and ellipsoid approximation and below 1.3 mm between the sphere and plane approximations.With regards to the comparison between the plane and ellipsoid approximations, the mean planimetric differences are about 1.4 cm.Altimetric differences are negligible for all of them.
With a 50 m receiver height above the reflecting surface, mean planimetric (resp.altimetric) differences reach 14 cm (resp.< 1 mm) between the sphere and ellipsoid approximations, 6.2 cm (resp. 2 mm) between the sphere and plane approximations, and 15 cm (resp. 2 mm) between the plane and ellipsoid approximations.
With a 300 m receiver height above the reflecting surface, mean planimetric (resp.altimetric) differences reach 83 cm (resp.< 1 mm) between the sphere and ellipsoid approxi-765 mations, 2.19 m (resp.8 cm) between the sphere and plane approximations, and 2.35 m (resp.8 cm) between the plane and ellipsoid approximations.
It is worth noticing that the sphere approximation is closer to the plane than the ellipsoid approximation when reflec-770 tions occur not too far from the receiver (below 560 m), and conversely if reflection occurs far from the receiver.

Influence of the satellite elevation angle
Secondly, by plotting the differences as functions of the satellite elevation angles, we can observe that the lapses between 775 the different algorithms vary in an inversely proportional way than the satellite elevation angle (and so, proportionally to the point distance from the receiver).The lower the satellite elevation angle is, the farther the specular reflection points from the receiver and the larger the impact of the Earth approxima-780 tion is.The choice of the algorithm used to perform the simulations becomes thus really important for the farthest reflection points (i.e., for low satellite elevation angles, and high receiver height above the reflecting surface).For example, mean planimetric differences between the local sphere and 785 ellipsoid approximation with a 50 m receiver height are about 14 cm considering satellites with elevation angles above 5°a nd are about 9 cm considering only satellites with elevation angles above 10°.With a 300 m receiver height above the reflecting surface, mean planimetric difference between 790 sphere and ellipsoid approximations is about 83 cm for satellites with elevation angle above 5°and 54 cm for a minimum elevation angle set to 10°.

Influence of the DEM integration
For continental surfaces the full integration of the DEM in 795 the simulation play a crucial role for a good calculation of the reflection points.The integration of a DEM leads to the suppression of 245 specular reflection points out of the 905 points determined during 24 hours the 4 th of October 2012 with the sphere approximation algorithm (figure 15).These 800 245 points came from a wave emitted by a satellite hidden by a mountain located in the south part of the area.In the north part, any reflection point is valid when taking a DEM into account, because in that direction, the topography is flat over the Geneva Lake, and so, satellites are all visible and reflections are possible.Moreover, the points positions have been rectified while taking a DEM into account, since the others algorithms consider that reflections occur (in first approximation) in a plane around the projection of the receiver and without integrating the problem of the presence of topogra-810 phy.

Comparison between the different models of the Earth surface
For a 5-m receiver height, and for satellite elevations greater than 10°, the mean planimetric difference between the ellipsoid and the sphere algorithm is equal to 1.4 cm whereas for a 300-m receiver height it is equal to 83 cm.The approximation done by considering the Earth as a sphere, an ellipsoid or a plane does not really affect the precision of the specular reflection point determination when reflections do not occur too far from the receiver i.e. for low receiver height and high satellite elevation.For example, if we consider that we need an uncertainty on the determination of the specular reflection position below 20 cm, the choice of the approximation of the Earth shape will have no influence if reflections occur until Concerning the algorithm taking the DEM into account, the differences obtained with respect to the sphere or ellipsoid algorithms are quite big even if the specular reflection point is close enough from the receiver.For instance, the mean planimetric (resp.altimetric) difference between the ellipsoid algorithm and the one integrating the DEM is equal to 70 cm (resp.18 cm) for a 5-m receiver height, and is equal to 78 m (resp.25 m) for a 300-m receiver height, and with satellite elevation angle above 5°.It is worth noticing that these differences will highly depend on the flatness of the considered area.

Angular refraction due to the troposphere
Given the geometric configuration of the satellite, the reflection point and the receiver, the same elevation angle correction will have a different effect according to the receiver height above the reflecting surface.It turns out that considering a same satellite at a given time, the corresponding reflection point will be farther for a big receiver height above the reflecting surface than for a smaller one.Consequently, for the same elevation angle correction, the resulting correction of the reflection point position will be higher in the first case than in the second one.Figure 16 shows the differences, in terms of geometric distances, between the reflection points positions obtained with and without correcting the angular refraction and for different receiver heights.It appears that for low satellite elevation angle and high receiver height, the angular refraction has a non-negligible influence on the specular point positions (116 m (resp.32 cm) for a 300-m receiver height and satellites elevation angle lower than 10°).

5 Conclusions
In this paper, we presented a simulator based on real GNSS satellites ephemeris, as a user-friendly tool, for modelling the trajectories of GNSS electromagnetic waves that are reflected on the surface of the Earth and therefore preparing GNSS-R 865 campaigns more efficiently.The originality of this simulator remains mainly in the integration of a DEM and of the correction of the angular refraction due to the troposphere.The results of simulations led us to a better understanding of the influence of some parameters on the reflection geom-870 etry, namely by quantifying the impact of the receiver height but also the influence of the satellite elevations, the natural topography (DEM), and the troposphere perturbation.
The different simulations realized near to quite rugged topography lead us to the following conclusions:

875
the DEM integration is really important for mountainous areas: planimetric differences as arc length (resp.altimetric differences as ellipsoid height) can reach 5.4 km (resp.1.0 km) for a 300-m receiver height, considering satellite with elevation angle greater than 5°.

880
differences between sphere and ellipsoid approximations are negligible for specular reflection points close from the receiver (closer than 50-60 m) i.e. small receiver height and/or high satellites elevations.For instance, planimetric differences are smaller than 11 cm 885 for a 5-m receiver-height, considering satellites with elevation angle greater than 10°.Altimetric differences are negligible.
sphere and plane approximations show really small differences in the vicinity of the receiver (smaller differ-890 ences than between the sphere and ellipsoid approximations): maximum differences are about 1.5 cm (resp.3 mm) with a 5 m receiver height (i.e.reflections occurring until 56 m from the receiver).
with regards to the plane and ellipsoid approxima-895 tions, differences are bigger than between the plane and sphere approximations when reflections occur closer than 550 m from the receiver.For farther reflections, differences between plane and ellipsoid become smaller than between plane and sphere.

900
the angular refraction due to troposphere can be negligible with regards to the position of the specular reflection point when the receiver height is below 5 m, but is absolutely mandatory otherwise, particularly for satellites with low elevation angle where the correction to apply 905 is exponential.
As a final remark, it is worth reminding that the farther the specular reflection point from the receiver is, the more important the influence of the different error sources will be: Earth approximation, DEM integration, angular refraction.
The farthest specular reflection points will be obtained for high receiver height and low satellite elevation.This simulator is likely to be of great help for the preparation of in situ experiments involving the GNSS-R technique.Further developments of the simulator will be soon implemented, such as receiver installed on a moving platform in order to map the area covered by airborne GNSS-R measurements campaigns and on-board a LEO satellite.The impact of the tide on the size of the reflecting area is non-negligible (decametric 3D-differences), and it is worth noticing that the gaps would have been even bigger integrating satellites with low elevation angle.Note also that the periodic variations of the 3D differences are only linked to the tide, since the mean of the satellite elevation angles does not show periodic variations during the day of simulation (43.3 ± 3.5°over the period).(Topography amplified by a factor 3) Yellow lines: direct waves, sphere approximation algorithm ; Green lines: direct waves, taking a DEM into account ; Blue lines: reflected waves, sphere approximation algorithm ; Red lines: reflected waves, taking a DEM into account.It can be noticed that some yellow and blue lines (direct and reflected waves, sphere approximation algorithm) go through the moutain (reflection points having been calculated inside the moutain), whereas any red or green line (direct and reflected waves, intergrating a DEM) go through it.

4N.
Roussel et al.: GNSS-R simulations 2.6 Data used for a simulator usage illustration 325 study casesSimulations and tests of parameters have been performed on two main sites: for simulations in the case of the Cordouan lighthouse are presented in figure7and figure 8 .These simulations were performed considering the sphere 685 approximation algorithm and a 15 minute time-step.Figure9(a) shows the variation of the distance between reflected points and receiver, as a function of the satellite elevation angle, and for several receiver heights above the reflecting surface and figure 9 (b) shows the variation of the 690 area of the first Fresnel surface.Such figures have been produced by performing simulations on the Cordouan lighthouse and varying the receiver height above the reflecting surface.
(figure 14 (b)).In order to get reflections below 125 m from the receiver, considering satellites with elevation angle above 5°, the receiver height above the reflecting surface should not exceed 25-30 m (figure9 (a)), which would corresponds to a first Fresnel zone area between 830 300 and 400 m 2 .

Fig. 1 .
Fig. 1.Principle of GNSS-Reflectometry.T : satellite/transmitter, S: specular reflection point, : satellite elevation, δAB(t) : additional path covered by the reflected wave, d : interdistance between the LHCP and RHCP antennas and h: height of the receiver above the reflecting surface.

Fig. 2 .
Fig. 2. Data flow chart of the simulator.Three main blocks: an input block which contains the different elements mandatory for the processing; a processing block where the user can choose which algorithm to be used, and an output block containing the different results of the simulation, namely KML files to be opened with Google Earth.

Fig. 3 .
Fig. 3. Determination of the specular reflection point in a local plane approximation and local difference with the sphere and ellipsoid approximations and DEM integration.S: specular reflection point position.R: receiver position.T: transmitter/satellite position.h: height of the receiver above the ground surface.

Fig. 4 .
Fig. 4. Local osculating sphere approximation : the three different reference systems of coordinates.S: specular reflection point position.R: receiver position.T: transmitter/satellite position.(0, X, Y, Z)R1: WGS84 Cartesian system.(0 , x, y)R2: local two-dimensional system, centred on the center of the osculating sphere, obtained by the rotation of the R1 system around the Z axis, in such a way that xr = 0.(S, x , y )R3: a local two-dimensional system, obtained by a rotation around the z axis and a rR translation of the R2 system in such a way that x' and the local vertical are colinear and that the system origin coincides with the specular reflection point S.

Fig. 5 .
Fig. 5. Determination of the specular reflection point integrating a DEM S: specular reflection point position.R: receiver position.T: transmitter/satellite position.A dichotomous process is applied for each topographic segment of the DEM to find if there is a point where the bisecting angle (equal to the sum of the anti-incident and scattering vectors) is colinear with the local normal vector.

Fig. 6 .
Fig. 6.Effect of the neutral atmosphere on the elevation angle.An exponential correction must be made for satellites with low elevation angle.

Fig. 7 .
Fig. 7. Positions of the specular reflection points and first Fresnel zones for one week of simulation on the Cordouan lighthouse with a 15-minute sampling rate (i.e.satellites positions actualized each 15 minutes).Only GPS satellites with elevation angle greater than 5°have been considered.Note the gap in the North direction.

Fig. 8 .
Fig. 8. First Fresnel zones and some direct and reflected waves display: 24h Cordouan lighthouse simulation with GPS constellation.

Fig. 10 .
Fig.10.Assessment of the tide influence.The impact of the tide on the size of the reflecting area is non-negligible (decametric 3D-differences), and it is worth noticing that the gaps would have been even bigger integrating satellites with low elevation angle.Note also that the periodic variations of the 3D differences are only linked to the tide, since the mean of the satellite elevation angles does not show periodic variations during the day of simulation (43.3 ± 3.5°over the period).

Fig. 11 .
Fig. 11.Planimetric and altimetric differences between the specular reflection points obtained with the different algorithm.Receiver height above the reflecting surface: 5 m.a) Planimetric differences as arc length (m).b) Altimetric differences as ellipsoid height (m).Note the dispersion within results for a fixed elevation angle which is a consequence of the azimuth variability in the ellipsoidal radius of curvature.

Fig. 12 .
Fig. 12. Planimetric and altimetric differences between the specular reflection points obtained with the different algorithm.Receiver height above the reflecting surface: 50 m.a) Planimetric differences as arc length (m).b) Altimetric differences as ellipsoid height (m).Note the dispersion within results for a fixed elevation angle which is a consequence of the azimuth variability in the ellipsoidal radius of curvature.

Fig. 13 .
Fig. 13.Planimetric and altimetric differences between the specular reflection points obtained with the different algorithm.Receiver height above the reflecting surface: 300 m. a) Planimetric differences as arc length (m).b) Altimetric differences as ellipsoid height (m).Note the dispersion within results for a fixed elevation angle which is a consequence of the azimuth variability in the ellipsoidal radius of curvature.

Fig. 14 .
Fig. 14. 3D differences between the specular reflection points obtained with the different algorithm.Receiver height above the reflecting surface of 5 m (a), 50 m (b) and 300 m (c).Note the dispersion within results for a fixed elevation angle which is a consequence of the azimuth variability in the ellipsoidal radius of curvature.

Fig. 15 .
Fig. 15.Influence of the topography -Direct and reflected waves display.(Topographyamplified by a factor 3) Yellow lines: direct waves, sphere approximation algorithm ; Green lines: direct waves, taking a DEM into account ; Blue lines: reflected waves, sphere approximation algorithm ; Red lines: reflected waves, taking a DEM into account.It can be noticed that some yellow and blue lines (direct and reflected waves, sphere approximation algorithm) go through the moutain (reflection points having been calculated inside the moutain), whereas any red or green line (direct and reflected waves, intergrating a DEM) go through it.

Fig. 16 .
Fig. 16.Angular refraction correction as a function of the satellite elevation angle and for different receiver height above the reflecting surface.a) Planimetric differences as arc length (m).b) Altimetric differences as ellipsoid height (m).

Table 1 .
Cross-validation between ellipsoid approximation and DEM integration

Table 2 .
Maximum differences between the positions of the specular reflection points obtained with the different algorithms and for different receiver heights above the reflecting surface.For each cell of this table, the first number is result obtained with minimum satellite elevation angle set to 5°, and the second number is the result obtained with minimum satellite elevation angle set to 10°.