## 1 Introduction  
Monitoring carbon and water cycles is key to understanding changes in the Earth's ecosystem and climate. Multiple geophysical parameters affect these cycles, including soil moisture, soil freeze/thaw state, snow water equivalent (SWE), snow depth, aboveground biomass, vegetation water content, and permafrost active layer properties. These geophysical parameters also affect the reflections of global navigation satellite system (GNSS) signals and other signals of opportunity (SoOp). Unlike active microwave remote sensing modalities, GNSS-reflectometry (GNSS-R) provides a low-cost solution to observe the geophysical variables. Multiple spaceborne GNSS-R systems have been launched, including the UK Disaster Monitoring Constellation (UK-DMC) [1] satellite, the TechDemoSat-1 (TDS-1) [2] satellite, the Cyclone GNSS (CYGNSS) [3] constellation, the Flexible Microwave Payload (FMPL)-2 instrument on a CubeSat [4], the Spire GNSS-R CubeSats [5], and the BuFeng-1 A/B twin satellites [6]. The use of GNSS-R for remote sensing has proven to be less expensive than active monostatic radar. Thus, with a fraction of the cost of a single active microwave remote sensing satellite, multiple GNSS-R receiver satellites can be launched [7]. Using multiple transmit-receive satellite pairs can also improve the revisit time. Moreover, the same GNSS-R constellation can be used to retrieve multiple geophysical parameters. For instance, CYGNSS was originally designed for ocean wind speed retrieval, but it is now also used for soil moisture retrievals and flood detection, to name a few uses of the data [8; 9; 10]. However, GNSS-R signal modeling and geophysical parameter retrieval can be more challenging compared to monostatic radar systems. Estimating geophysical parameters from GNSS-R with high accuracy requires accurate modeling of the reflected signal. Multiple models have been developed for both ocean and land surfaces [11; 12; 13; 14; 15; 16; 17; 18; 19]. For example, the Soil And VEgetation Reflection Simulator (SAVERS) [12] models scattering from both the Earth land surface and its overlying vegetation. This model is based on the bistatic radar equation integrated over a digital elevation model (DEM). The scattering process is divided into two components: near-specular and incoherent. The near-specular component is derived from a physical optics model that includes wavefront sphericity and a parameter for the effective beamwidth of the DEM surface facets. The incoherent component is based on the advanced integral equation method (AIEM) and is minimal except for very rough surfaces. Another example is the analytical Kirchhoff solutions (AKS) model, which calculates the scattered signal from the ground surface only [17]. This model is the successor to the fine scale partially coherent patch model (FPCP) [14]. The surface in AKS is modeled as physical patches defined by a DEM with superimposed random roughness. Furthermore, the analytical solutions for coherent and incoherent scattering are expressed in terms of the spectrum of the random roughness. Thus, AKS does not require division of the random roughness into discrete scales. Finally, the geometric optics with topography (GOT) model [11], which is the basis for the model in this paper, also models the ground surface only. This model is based on the geometrical optics limit of the Kirchhoff integral. All three models (SAVERS, AKS, and GOT) are based on the Kirchhoff integral equation and use a DEM to parameterize their formulations. Furthermore, they are suitable for spaceborne GNSS-R signals as well as airborne and ground GNSS-R signals. A comparison of these models, presented in [20], was performed on a single CYGNSS track over a Soil moisture Sensing Controller and oPitmal Estimator (SoilSCAPE) site located in the San Luis Valley (SLV), CO, USA [21]. Moreover, the input parameters were unified for all of the models in this intercomparison. For example, all models used the one-arc-second Shuttle Radar Topography Mission (SRTM) DEM.