Microwave Imaging Radiometers by Aperture Synthesis—performance Simulator (part 1): Radiative Transfer Module

The Synthetic Aperture Interferometric Radiometer Performance Simulator (SAIRPS) is a three-year project sponsored by the European Space Agency (ESA) to develop a completely generic end-to-end performance simulator of arbitrary synthetic aperture interferometric radiometers. This means, on one side, a generic radiative transfer module from 1 to 100 GHz, including land and ocean covers, as well as a fully 3D atmosphere and Faraday ionospheric rotation based on variable TEC. On the other hand, the instrument can have an arbitrary array topology (number of antenna elements, and their time-dependent position and orientation). Receivers' topology can also be modified, starting from a very generic one to connecting and disconnecting subsystems, whose parameters can be individually configured. These parameters can be defined either by mathematical functions or by input data files, including the frequency and temperature dependence. Generic calibration and image reconstruction algorithms that are suitable for arbitrary array topologies have also been implemented, as well as tools to compute the instrument performance metrics, i.e., radiometric accuracy, sensitivity, angular resolution, etc. This manuscript presents the generic architecture of the SAIRPS, the algorithms implemented in the Radiative Transfer Module, and simulation results showing its performance. A companion manuscript (Part II) describes the instrument and calibration modelling, the image reconstruction algorithms, and the validation tests that were performed.


Introduction
Synthetic aperture radiometry was developed in the 1950s to obtain high-resolution radio images of the sky.In 1983, LeVine and Good [1] first proposed its use for Earth observation as a way to increase the angular resolution of microwave radiometers.Until the advent of synthetic aperture microwave radiometry, brightness temperature maps were obtained by either a mechanical or electrical scan of a large antenna.In synthetic aperture interferometric radiometers (SAIRs), brightness temperature images are formed through a Fourier synthesis process in a snapshot basis, after cross-correlating all the signal pairs collected by the array elements.The mean value of the image is also measured using at least one real aperture radiometer connected to one of the antenna elements.For simplicity of operation, array elements are located in a plane, and the Z-axis is orthogonal to it, but this is not a restriction of the Fourier synthesis process itself.
So far, only the MIRAS (Microwave Imaging Radiometer by Aperture Synthesis) instrument on board the ESA (European Space Agency) SMOS (Soil Moisture and Ocean Salinity) mission launched in 2009 is using successfully this new technique to generate global and frequent soil moisture and ocean surface salinity maps (e.g., [2]).The interested reader is referred to [3] for a summary of the principles of operation of this type of instrument.The SEPS (SMOS End-to-end Performance Simulator) was developed to study the MIRAS instrument performance and to develop the geophysical parameters retrieval algorithms needed before SMOS was launched [4].
After SMOS' success, today, a number of instruments are planned or under study: the GeoStar instrument is the baseline payload for the Precipitation and All-weather Temperature and Humidity (PATH) mission from NASA (USA) [5], the Geostationary Atmospheric Sounder (GAS) instrument is under study for post-MSG operational satellites observations (Europe) [6], and the Geostationary Interferometric Microwave Sounder (GIMS) instrument from NSSC-CAS (China) [7].The study of the instrument performance in terms of angular resolution and radiometric performance (

Main Building Blocks of a Synthetic Aperture Interferometric Radiometer Performance Simulator
The main purpose of a generic SAIRPS is to simulate and compute figures of merit of the performance for arbitrary SAIRs, with arbitrary receiver and array topologies, to assist in the definition of future instruments and missions.As a generic instrument simulator, it is required to have a flexible architecture and functions to simulate various instrument conditions.
Figure 1 presents the overall architecture of a SAIRPS, which consists of five modules: -The Geometry Module, which is capable of simulating arbitrary antenna positions and orientations (i.e., polarization or rotation) in the array, that are also defined individually as time-dependent variables (even within the integration time, which is subdivided in a number of subintervals to allow the simulation of image blurring, for example, due to instrument motion).The most common array configurations are provided as pre-defined standards and include: uniform linear array, Minimum (or Low) Redundancy Linear Array (MRLA or LRLA), star-like arrays with arbitrary number of arms, hexagonal arrays, circular arrays, U-and T-arrays . . .and the antenna spacing can be arbitrarily selected, or even be non-uniform in which case it will follow a geometric progression.It is also in this module that the orbital parameters are set, with a wide range of choices given to the user.-The RTM and BT Maps Module that either computes the BT maps (4 Stokes parameters) at a given frequency and polarization, as a function of the incident and azimuth angles, and geophysical parameters, or ingests user defined BT maps, computing the FOV of the instrument and the surface area captured in each snapshot, according to the inputs given by the Geometry Module, and provide snapshot-based, interpolated, BT maps to the Instrument module.The RTM and BT maps module is detailed as shown in Figure 2, and it is the object of this manuscript.-The Instrument Module, which can simulate a generic receiver configuration.It includes the basic SAIR instrument components and their corresponding errors; the noise impact induced by instrument components; the physical temperature change impact; hardware non-idealities and uncertainties in the characterization of the instrument.-The Calibration Module is able to calculate the calibrated complex visibility samples from the raw observables produced by the Instrument Module.This module is able to simulate both internal and external calibrations.-The Image Reconstruction Module is able to produce a reconstructed brightness temperature map from the calibrated visibility samples, from the most basic image reconstruction method (a non-uniform FFT) [9], to the more sophisticated G-matrix method [10] (solved using either the Truncated Singular Value Decomposition, or the Conjugate Gradient), and the CLEAN method.Finally, -The Performance Module is used to evaluate the instrument's performance in angular and radiometric resolutions.
The Instrument, Calibration, and Image Reconstruction Modules are the object of the second part of this manuscript.
geophysical parameters, or ingests user defined BT maps, computing the FOV of the instrument and the surface area captured in each snapshot, according to the inputs given by the Geometry Module, and provide snapshot-based, interpolated, BT maps to the Instrument module.The RTM and BT maps module is detailed as shown in Figure 2, and it is the object of this manuscript.

-
The Instrument Module, which can simulate a generic receiver configuration.It includes the basic SAIR instrument components and their corresponding errors; the noise impact induced by instrument components; the physical temperature change impact; hardware non-idealities and uncertainties in the characterization of the instrument.

-
The Calibration Module is able to calculate the calibrated complex visibility samples from the raw observables produced by the Instrument Module.This module is able to simulate both internal and external calibrations.
The Image Reconstruction Module is able to produce a reconstructed brightness temperature map from the calibrated visibility samples, from the most basic image reconstruction method (a non-uniform FFT) [9], to the more sophisticated G-matrix method [10] (solved using either the Truncated Singular Value Decomposition, or the Conjugate Gradient), and the CLEAN method. Finally,

-
The Performance Module is used to evaluate the instrument's performance in angular and radiometric resolutions.
The Instrument, Calibration, and Image Reconstruction Modules are the object of the second part of this manuscript.A SAIRPS must also support different simulation modes to support different studies' objectives: -Snapshot Mode to process a single BT snapshot, corresponding to the instrument elementary integration period, in order to compare against the reference BT maps, -Monte-Carlo Mode to evaluate the instrument performance by analyzing a sequence of consecutive snapshots for a given constant brightness temperature map.A set of system/instrument random errors can be included in the simulation and the user has the possibility to select which of the components parameters are fixed and which ones are allowed to vary according to its configuration.-Time Evolution Mode to analyze instrument's drifts or to generate synthetic multi-look observations to simulate the actual instrument performance and develop geophysical parameters retrieval algorithms.

RTM block Design
The main blocks of the RTM are: -Land surface contribution computation, including bare soil emission, vegetation-covered soils emission, ice and snow, and other effects, such as rocky soil, urban areas, inland water bodies, and topography effects, -Ocean surface contribution computation, including flat sea emission, and wind effects and the induced azimuthal signature in the four Stokes parameters, -Atmosphere attenuation computation (scattering effects are not included), and -Sky contribution computation, including direct and reflection contributions over land and over sea.
In Figure 3, the schematic block diagram of the RTM algorithms and how they interact with each other is shown.The most important inputs for each algorithm block are also shown: A SAIRPS must also support different simulation modes to support different studies' objectives: -Snapshot Mode to process a single BT snapshot, corresponding to the instrument elementary integration period, in order to compare against the reference BT maps, -Monte-Carlo Mode to evaluate the instrument performance by analyzing a sequence of consecutive snapshots for a given constant brightness temperature map.A set of system/instrument random errors can be included in the simulation and the user has the possibility to select which of the components parameters are fixed and which ones are allowed to vary according to its configuration.-Time Evolution Mode to analyze instrument's drifts or to generate synthetic multi-look observations to simulate the actual instrument performance and develop geophysical parameters retrieval algorithms.

RTM block Design
The main blocks of the RTM are: -Land surface contribution computation, including bare soil emission, vegetation-covered soils emission, ice and snow, and other effects, such as rocky soil, urban areas, inland water bodies, and topography effects, -Ocean surface contribution computation, including flat sea emission, and wind effects and the induced azimuthal signature in the four Stokes parameters, -Atmosphere attenuation computation (scattering effects are not included), and -Sky contribution computation, including direct and reflection contributions over land and over sea.
In Figure 3, the schematic block diagram of the RTM algorithms and how they interact with each other is shown.The most important inputs for each algorithm block are also shown: The generic Radiative Transfer Model (RTM) follows the approach described in Figure 2 of [11], including the Earth's surface emission and the atmospheric emission.In that example, the Earth's surface is the ocean, the atmosphere is a non-scattering one, the input variables are the Earth Incidence Angle (EIA), the Sea Surface Temperature (SST), the SSS (Sea Surface Salinity), the output variables are the up-welling and down-welling atmospheric brightness temperatures (TUP and TDN) at Top Of the Atmosphere (TOA), and τ the atmospheric opacity.Finally, Ω is an intermediate variable that accounts for the atmospheric path correction in the scattered down-welling radiation.
In the general case, the RTM has to solve the radiative transfer equation (RTE) with the appropriate boundary conditions in the surfaces of the Earth, and in the top of the atmosphere ( [12]): where B T is the Stokes vector. �  in Equation ( 1) is the extinction matrix given by: where  0 is the number of particles per unit of volume, k is the wave number, fpq are the forward scattering amplitudes, the operator < > stands for the average over the ensemble of scatterers in the volume (raindrops, each one with a canting angle), and ℜ and ℑ are the real and imaginary parts operators, respectively.The emission vector  � is given by: The generic Radiative Transfer Model (RTM) follows the approach described in Figure 2 of [11], including the Earth's surface emission and the atmospheric emission.In that example, the Earth's surface is the ocean, the atmosphere is a non-scattering one, the input variables are the Earth Incidence Angle (EIA), the Sea Surface Temperature (SST), the SSS (Sea Surface Salinity), the output variables are the up-welling and down-welling atmospheric brightness temperatures (T UP and T DN ) at Top Of the Atmosphere (TOA), and τ the atmospheric opacity.Finally, Ω is an intermediate variable that accounts for the atmospheric path correction in the scattered down-welling radiation.
In the general case, the RTM has to solve the radiative transfer equation (RTE) with the appropriate boundary conditions in the surfaces of the Earth, and in the top of the atmosphere ( [12]): where T B is the Stokes vector.k e in Equation ( 1) is the extinction matrix given by: where n 0 is the number of particles per unit of volume, k is the wave number, f pq are the forward scattering amplitudes, the operator < > stands for the average over the ensemble of scatterers in the volume (raindrops, each one with a canting angle), and and are the real and imaginary parts operators, respectively.The emission vector F is given by: Fpθ, ϕq " " k a1 pπ ´θ, π `φq k a2 pπ ´θ, π `φq ´ka3 pπ ´θ, π `φq ´ka4 pπ ´θ, π `φq ı , k a1 pθ, φq " k e11 pθ, φq ´2π and the phase matrix P is given by: In Equations ( 2) and ( 4) the f pq coefficients are the vector-scattering amplitude functions that provide the amplitude, phase and polarization information of the scattered field at q-polarization Ñ E q " êq f pq `θ, θ 1 , ϕ, ϕ 1 ˘e´jkr {r, when a plane wave at p-polarization r is incident on each scatterer, and as in Equation (2) the operator < > stands for the average over the ensemble of scatterers in the volume.
Assuming a non-scattering atmosphere a simplified RTM expression for the top of the atmosphere (TOA) brightness temperature (T B,p ) at p = v, h polarizations can be written as: T B,p " T UP `Tatm,UP{DN ¨ep ¨Ts `Tatm,UP{DN ¨Γp " T DN `Tatm,DN ¨Tsky where T S is the surface's temperature, e p is the surface's emissivity, T atm,UP{DN " exp ´´τ atm,UP{DN ¯, and τ atm,UP{DN are the up-welling and down-welling atmospheric opacities, Г p is the Earth's surface reflection coefficient, T UP{DN are the up-and down-welling atmospheric brightness temperatures, and T sky is the contribution from outer space radiation sources (e.g., the cosmic background T cos ~2.7 K, the galaxy plane and the Sun, which are very much frequency-dependent . . .etc.).Actually, Γ p is not a reflection coefficient since it integrates all the bistatic contributions over the upper hemisphere scattered in the direction of the sensor.These effects will be discussed at the end of this document.
In the next sections, the models used to compute e p and Γ p will be discussed, leaving for the last part the discussion of the modeling of the atmospheric contributors.

Earth's Surface Contribution
Land/Cryosphere Contributions As commented before, different cases have been considered: bare soil, soil covered by vegetation, and soil covered by ice or snow, inland water bodies (fresh or salt water), urban areas, and mountainous regions (large topography effects).Since most of current models for vegetation and snow covered soils have a limited frequency range of validity, and require very detailed information which is usually not available, simplified models have been adopted and some parameters, e.g., the albedos, have been adjusted to fit modelled brightness temperatures with reported experimental data.Land emission closely follows the models used in the direct model of the SMOS retrieval algorithm [13], and have been extended in their range of validity whenever general models exist.

A. Bare Soil Emission
Dependence on the soil moisture (SM) is introduced through the dielectric permitivity of the soil used in the Fresnel reflection coefficient.In addition to the soil moisture data, mixing models require data about the sand (S) and clay (C) fractions, the soil bulk and solid soil material densities and the soil dielectric constant.A number of models exist for the soil dielectric constant [14][15][16][17]) even though the most recent measurements with SMOS seem to indicate that [15] is more accurate.
For bare soils the emission is given by: e bare p pθq " 1 ´Γbare p pθq where T S is the soil effective temperature, and the emissivity e bare p pθq depends on the soil's characteristics (moisture, texture, roughness and -eventually-salinity, not included here).

‚
Flat soil emission: for smooth bare soils, the reflection coefficients are given by the (power) Fresnel reflection coefficients as given by: where ε r is the dielectric constant and the soil is assumed to be a non-magnetic medium (µ r = 1).In this work, the widely used Dobson model (Equations ( 10)-( 12)) is selected for the soil dielectric constant [16], because of its wide frequency range of validity (1-18 GHz), despite the higher accuracy at L-band of other dielectric constant models such as the Mironov et al. one [15], which is in better agreement with measurements for a larger range of soil texture types: In the above equation, ρ b is the soil bulk density, which is a function of the soil texture (default value = 1.3 g/cm 3 ), ρ s is the soil particle density (assumed to be constant and equal to 2.664 g/cm 3 ), ε pa is the dielectric constant of solid particles (ε pa « 4.7), α = 0.65, β " β 1 ´j¨β 2 is an empirically-calculated complex function of soil texture parameter usually calculated as in [16,17], SM is the soil moisture value in m 3 /m 3 , and ε s f w " ε 1 s f w ´j¨ε 2 s f w is the dielectric constant of free water included in the soil with: where ε w0 is the static dielectric constant of water, ε w8 is the high frequency limit of the dielectric constant of water, f is the frequency [Hz], rτ w is the relaxation time of water, ε 0 " 8.854¨10 ´12 F{m, where the coefficients SGEF i , BERE i , and BEIM i are provided in [18], and S and C are the sand and clay fractional contents.
In the case of dry sand, since it has almost no bound water and hence has specific dielectric constant behavior, it has specific water capacities and can be very dry, leading to large penetration depths.Hence, it is often considered that the dielectric constant of sand can be expressed as ε r dry sand « 2.53 -j¨0.05, which is fairly constant over a wide range of frequencies, including our frequencies of interest.

‚
Surface's roughness effects are accounted for by using a "rough" surface reflectivity Γ bare, rough p that is computed using the following empirical relationship: where Q is the so-called polarization mixing factor Q " 0.35 ı , HR is an effective surface's roughness (no units), and NR p is a polarization-dependent integer between 0 and 2 used to parameterize the roughness effects of the incidence angle [19].
The expression for HR is given by: where the transition moisture point (XMVT): is a function of the wilting point (WP): WP pC, Sq " 0.06774 ´0.00064¨S `0.00478¨C and the field capacity (FC) is given by: where the S and C are expressed in %.
Finally HR max is a function of the land cover type, k is the wave number, and σ is the RMS height.

‚
The effective soil temperature (T s ) depends on the soil properties and moisture content profile within the soil volume.The effective temperature is usually computed using a surface temperature and the temperature at a depth where it is almost constant.The actual profile and depth are dependent upon the soil type actual profile and the level at which the deep soil temperature is obtained.The effective soil temperature is written as a function of the soil temperature at depth (T soil depth , approximately at 0.5 m depth), and surface soil temperature (T soil surf , approximately at 5 cm) as follows: T s " T soil depth `Ct ¨´T soil sur f ´Tsoil depth ¯ (18) where: is a parameter depending mainly on the frequency and the soil moisture.If the soil is very dry, soil layers at depth (deeper than one meter for dry sand) contribute significantly to the soil emission, and the value of C t is lower than 0.5.Conversely, if the soil is very wet, the soil emission originates mainly from layers at the soil surface and C t « 1.In Equation ( 19) w 0 and bw 0 are parameters that depend mainly on the soil texture and structure.Often it can be assumed that T s " T soil sur f .

B. Vegetation-Covered Soils Emission
Is modelled by means of the three layer approach described in [20].Its impact on brightness temperature is two-fold: a) it absorbs and scatters the direct bare soil radiation, and attenuates and reflects above surface radiation directly and indirectly, through bare soil reflectivity, and b) it introduces its own upward and downward radiation; the latter leading to an indirect contribution through soil reflectivity and self-attenuation.
Whenever the air-vegetation reflection coefficient is negligible (Γ air-veg ~0) [20], and neglecting scattering, the vegetation emission model simplifies to the widely known as the τ-ω model [19] (τ: optical depth, ω: single scattering albedo) generally used at low microwave frequency bands.At higher frequency bands, even though multiple scattering starts playing a role, the same model can be used provided the parameters are properly tuned.Several types of vegetation have to be considered: low vegetation (grassland, crops, etc.), and forest vegetation (coniferous, evergreen and deciduous forests).

‚
For low vegetation, the τ-ω model is given by: where the first term corresponds to the soil emission attenuated by the vegetation canopy, the second one corresponds to the vegetation upwelling emission including the first order scattering through the single scattering albedo (ω canopy p ), and the third one to the vegetation down-welling emission reflected on the soil surface, and then attenuated and scattered in the vegetation in the upwelling path.The vegetation layer attenuation is L canopy p " e τ canopy p {cospθq , where τ canopy p is the nadir optical depth of the vegetation layer.
From the full-polarimetric RTE, τ canopy p can be computed as: where k e11 and k e22 are the extinction coefficients at v and h polarizations, elements 11 and 22 of the extinction matrix k e (Equation ( 2)), but the above equations are very difficult to apply in practice.
Actually, the computation of the optical depth must account not just for green vegetation, but for litter and intercepted water as well: The effect of vegetation structure on the optical depth can be modelled by: where τ N AD is the attenuation at nadir, which is estimated as a function of the Leaf Area Index (LAI), and C pol is the polarization correction factor (C pol ą 1 for vertical structures).However, for "large" footprints including a variety of vegetation types, polarization and incidence angle effects are averaged.
Recent results have shown that the effects of water interception by the vegetation canopy due to rain or dew may be very significant, i.e., optical depth may increase by a factor of two or three during and after rainfalls over a fallow for instance.However, estimating the fraction of intercepted water, which depends on the intensity of the rainfall events, vegetation type and evaporation fluxes, would be very difficult, and an accurate modelling of these effects is not known to date to authors' knowledge.
Even though it is not well known, it is likely that the effect of litter may be significant, and may explain why very high values of b (b ~0.4 in Equation ( 27)) must be used over natural vegetation covers such as prairies.The opacity due to the litter can be modelled as [21,22]: where c L " 0.24 and LWC is the Liquid Water Content in the litter in [kg/m 2 ], computed as: where M gL " a L ¨SM `bL (0 ď M gL ď 0.8) , with default parameters B sL = 0.3 kg/m 2 , a L " 2.33, and b L = 0.At 1.4 GHz many studies have related the vegetation opacity to the Vegetation Water Content (VWC) as: Good correlation has also been obtained for green vegetation between τ p and the Leaf Area Index (LAI), although this parameterization is less accurate during the senescence phase, during which the opacity might be underestimated from low LAI values over some vegetation types.A possible parameterization is given by: where the tt p are parameters that allow accounting for the dependence of τ p with the incidence angle.Even though all vegetation parameters b 1 s , b 2 s , tt h and tt v are function of the canopy type, the dependence of b 1 s and b 2 s on the canopy hydric status and the change of the vegetation structure in relation with phenology is neglected (b 1 s " 0.06, b 2 s " 0), as well as the dependence of τ p with the incidence angle and polarization, which can also be neglected (tt h " tt v " 1).
In Equation ( 20), the single-scattering albedo ω canopy p , as well as the opacities, should be given as inputs, and are not computed by integration of the phase matrix.At L-band, by default, it can be safely assumed that ω canopy p « 0.

‚
Over forests, a simple τ-ω empirical approach is not entirely appropriate because the emission/scattering processes are more complex.Since trunks and branches are comparable in size to the electromagnetic wavelength, multiple scattering effects cannot be neglected, and a simple first-order approach is not reliable [23][24][25].Despite the above considerations, from an operational point of view, a simple τ-ω approach is kept here, as it is in the SMOS L2 Soil Moisture processor.Three forest categories are aggregated: (1) needle leaf and broadleaf (including Tropical forests and woodland); (2) mixed forests; and (3) woodlands.The same general procedure is applied for the 3 categories as in the low vegetation case, although the parameters are specific of each category: where F V is the fraction of the LAI corresponding to the understory contribution.As a result of the variability in orientation of branches and leaves, a simple τ N AD constant independent on the polarization and incidence angle may be used: b 1 F " 0.290 for deciduous broadleaf, evergreen broadleaf, and woodlands, b 1 F " 0.360 for needle leaf, b 1 F " 0.325 for mixed forests, and b 2 F " 0.06¨FV (FV forest-type dependent).
The single scattering albedo may also be considered constant (i.e., independent on angle, polarization and time), but it is not negligible, since it actually is ω canopy p 0.08 at L-band.At other frequencies, suitable values were used as input parameters.
Alternatively, the vegetation attenuation could be computed from the vegetation height and the extinction constant, to be computed from the vegetation dielectric constant, but these parameters -to authors' knowledge-are not available at global scale.

‚
Frozen soils cover large areas at high latitudes and sometimes altitudes.At mid latitudes, frozen soil can also be expected in winter, especially for morning-pass orbits (e.g., 6 AM).Experience shows that the dielectric properties of frozen soil are very close to those of dry soil, while vegetation is almost fully transparent [26].It is often considered that for frozen soils the dielectric constant can be written [27].ε r " 5.0 ´j¨0.5 ‚ Ice covered soils emission is computed in a similar way [28], but the dielectric constant of the ice layer should be used instead.The input parameters are the temperature (T in [K]) and the frequency (f in [GHz]): α " ˆ0.00504 `0.0062¨ˆ3 00 T ´1˙˙¨e ´22.1¨p 300 T ´1q (32a) It can be noted that for pure ice, ε i is very small and, therefore, it is rather transparent, with penetration depths up to several hundreds of meters.

‚
Snow covered soils account for about 40% of the Northern hemisphere land mass, and it varies seasonally.Dielectric properties depend on its history: while fresh, dry snow is transparent to microwave radiation; as snow melts its dielectric constant increases, snow grain size and liquid water content increase and may be totally opaque (at 0 ˝C) when wet.
The emissivity of the soil surface covered by snow is modelled in a similar way as the vegetation covered soil.Dry and wet snow models are used.For dry snow, the input parameters are the temperature (T in [K]) between ´40 ˝C and 0 ˝C, and the dry snow density (P s in [g/cm 3 ]) or the ice volume fraction (v i ), for frequencies between 0.8 and 37 GHz [28].
For wet snow, at a temperature near freezing, the input parameters are the volumetric water content (m v ) between 1% and 12%, and the dry snow density (P s in [g/cm 3 ]) between 0.09 and 0.38 g/cm 3 , and for frequencies between 3 and 37 GHz [28].
This behavior is shown in Figure 4 below as a function of frequency for both dry and wet snow.

D) Other Effects
• Rocks and rocky outcrops areas are not well modelled at present.They are assumed to behave as very dry soils.Field measurements do not show a significant effect from rocks, which are usually encountered in barren areas or in mountains regions, etc.In [20] permittivity values are given for rocks at 400 MHz and 35 GHz.They range from 2.4 to 9.6.Approximate expressions do exist for rocks, but a default constant dielectric constant value is usually assumed: • Urban Areas are usually largely covered by concrete and asphalt, therefore a bare soil model is used, but using the dielectric constant of dry asphalt.

‚
Rocks and rocky outcrops areas are not well modelled at present.They are assumed to behave as very dry soils.Field measurements do not show a significant effect from rocks, which are usually encountered in barren areas or in mountains regions, etc.In [20] permittivity values are given for rocks at 400 MHz and 35 GHz.They range from 2.4 to 9.6.Approximate expressions do exist for rocks, but a default constant dielectric constant value is usually assumed: ‚ Urban Areas are usually largely covered by concrete and asphalt, therefore a bare soil model is used, but using the dielectric constant of dry asphalt.
‚ Inland Water Bodies and Rivers are modelled using a "flat" "bare soil model" with (Q = 0) and using the fresh water dielectric constant model, which depends on the frequency and temperature (salinity is assumed to be equal to zero).In [29] a model is presented that has been validated from 19 to 86 GHz data, and in its derivation it included data up to 500 GHz.The implementation used in this project follows that of [28], which is applicable for frequencies (f ) up to 1 THz, physical temperature (T) from 0 ˝C to 30 ˝C, and salinities (S) from 0 to 40 psu.Please note that the formulas below will be used in the sea surface emissivity calculation with the appropriate value of salinity S (‰0).

‚
Topography effects produce a change of the local incidence angle, polarization mixing, and eventually a shadowing and multiple scattering, which are not accounted for in this model.The equations to compute the scattered field components from the illuminating wave components are summarized below [30]: where the effective (electric field amplitude) reflection coefficients are given by: and R v `θl,k ˘and R h `θl,k ˘are the Fresnel reflection coefficients for the electric field (or the modified ones to account for surface roughness, vegetation, etc.) at v and h polarizations.The power reflection coefficients were derived by taking the square of the module of the (amplitude) electric field reflection coefficients.The angle between the facet and the observation direction θ l,k is given by: where nk is the normal to the facet, which exhibits a tilt angle α, oriented by an azimuth angle β with respect to the global plane of incidence.The other unitary vectors are defined as follows: ĥi,k " ki,k ˆẑ vs " ĥs ˆk s ĥs " ks ˆẑ ˇˇks ˆẑ ˇˇ( 44) qk " ks ˆn k ˇˇks ˆn k ˇˇ" ps,k " qk ˆk s (47) ki,k " ´I ´2¨n k ¨n k ¯k s (48) With the above quantities, the local reflectivities Г v,l and Г h,l can be represented in the global reference system, taking into account the polarization rotation [30]: Γ h pθq " Γ v,l pθ l q ¨sin 2 pϕq `Γh,l pθ l q ¨cos 2 pϕq (50) where: sin pϕq " sin pβq ¨sin pθ l q (51) The above implementation is extremely cumbersome for practical purposes.To overcome these difficulties, using a Geometric Optics model, and a 30 m resolution Digital Elevation Model (DEM), Utku and Le Vine [31] computed the net brightness temperature change at vertical and horizontal polarizations, as a function of the incidence angle, for different scenarios with different rms slopes.
Note that these changes in the surface's emission will also be affected by the atmospheric attenuation (Equation ( 5)).
Figure 5 top shows these variations for low and high topography scenes (increasing at h-polarization, and decreasing at v-polarization), and the 4th order polynomial fit. Figure 5, bottom, shows the residual fits, which are negligible for all practical purposes. where: The above implementation is extremely cumbersome for practical purposes.To overcome these difficulties, using a Geometric Optics model, and a 30 m resolution Digital Elevation Model (DEM), Utku and Le Vine [31] computed the net brightness temperature change at vertical and horizontal polarizations, as a function of the incidence angle, for different scenarios with different rms slopes.Note that these changes in the surface's emission will also be affected by the atmospheric attenuation (Equation ( 5)).
Figure 5 top shows these variations for low and high topography scenes (increasing at hpolarization, and decreasing at v-polarization), and the 4th order polynomial fit. Figure 5, bottom, shows the residual fits, which are negligible for all practical purposes.In SAIRPS a ~450 m resolution global DEM has been used instead.This lower resolution translates into a reduced rms slopes value (mean square slopes), but since the locations of the scenarios used in [31] are known, an adjustment has been made between the rms slopes computed at 30 m and the ones at 450 m resolution.Equations ( 52) and ( 53), together with the computed coefficients listed in Table 1, provide the polynomial fits for the topography effects assuming Ts = 300 K [31]: In SAIRPS a ~450 m resolution global DEM has been used instead.This lower resolution translates into a reduced rms slopes value (mean square slopes), but since the locations of the scenarios used in [31] are known, an adjustment has been made between the rms slopes computed at 30 m and the ones at 450 m resolution.Equations ( 52) and ( 53), together with the computed coefficients listed in Table 1, provide the polynomial fits for the topography effects assuming T s = 300 K [31]:  The change in the reflection coefficient can then be estimated as: Other secondary effects induced by topography are a variation of the atmospheric up-welling temperature (T atm, up ) term due to the path variations through the atmosphere, and a variation of the atmospheric down-welling temperature scattered over the Earth's surface from directions away from the main specular direction in the global reference frame (T atm,sc ).These two effects will be treated later on, when dealing with the atmospheric effects.

Ocean Surface Contributions
The polarimetric emissivity of an irregular surface in the XY plane at h and v polarizations, observed from an observation angle θ (relative to 0 ˝at nadir) and an azimuth angle φ (relative to 0 ˝in the direction of the wind), are related to its scattering properties by [32,33].
where the polarimetric bistatic scattering coefficients γ mnpq (θ, φ, θ i , φ i ), describe the scattering in the direction (θ, φ) of a wave incident from the direction (θ i , φ i ).The coefficient γ mnpq describes the resulting cross-correlation of the scattered waves at m and p polarizations due to incident waves at n and q polarizations respectively.
Since the above ocean surface's emission model is very computationally intensive, because it requires a two-fold integral over 2π stereo-radians for each direction, the Meissner and Went's model is used instead [11].It has been validated from 6 to 90 GHz using SSM/I and WindSat data for a wide range of wind speeds, and includes the variations with both the incidence and azimuth angles.In this formulation, the ocean's emission is modelled as: e p " e 0,p `∆e w,p `∆e ϕ,p ( The first term corresponds to the flat surface's emissivity, the second term corresponds to the isotropic wind-induced emissivity, which depends on wind speed w, and the third term corresponds to the wind direction signal, which contains the dependence on wind direction ϕ relative to the azimuthal look.

A) Flat Sea Surface Emissivity
The flat sea surface emissivity is governed by the Fresnel (power) reflection coefficient at the interface air-sea, that depends on the polarization, incidence angle, and dielectric constant, which is a function of the frequency, salinity and temperature (given in Equation ( 37)).e 0,p pθq " 1 ´Γ0,p pε r p f , T s , Sq , θq (57)

B) Isotropic Wind-induced Emissivity Signal
The equation that fits the isotropic emission due to the wind speed is listed below (Equation ( 60)).Typical θ re f is 55.2 ˝, and T s,ref = 20 ˝C.For θ i ě θ re f the behaviour is extrapolated linearly [11].

C) Wind-directional Signal: Stokes Parameters
The wind-directional signal is given by Equation (64) for the four Stokes parameters, although -to our knowledge-, no reliable, published measurements for S 3 and S 4 exist below 10.7 or above 37.0 GHz.The coefficients A p, f i ´θre f , w ¯are given in Equations ( 62)-(69).In Equation (67) they are expressed in terms of the true Stokes parameters (as opposed to the modified Stokes parameters v and h) [11].
Above these frequencies, the authors are not aware of models as accurate and validated as the ones described in [11].Below these frequencies (e.g., 1.4 GHz) the model in [34], obtained from a fit of SMOS data has been implemented (Equations ( 70) and ( 71)).This model is, in principle, more precise, although it does not predict any azimuthal or polarimetric (S3 and S4) dependence of the Stokes emission vector, although these signatures are almost negligible at these lower frequencies.The coefficients w 1j , w 2j, b j , W j , B, a and b in Equations ( 70) and (71) can be found in [34].r T B,p " e p ¨Ts "

D) Sea ice
Sea ice in the Arctic is the top layer of the 'halocline', a 200-300 m thick layer of low salinity, but also very low temperatures, close to ´2 ˝C [35].It is modelled as a layer of salty ice (20-30 psu), typically 2 to 3 m thick, and in some regions up to 4 to 5 m thick, at ´2 ˝C on top of the sea water.Since the information on the ice thickness is not available, and there is only sensitivity to the ice thickness at low microwave frequencies, i.e., L-band, and up to a maximum of 0.6-1 m thickness [36], for practical reasons it can be considered as an infinite layer.
Antarctic sea ice is typically first-year 1 to 2 meters thick, and it is modelled as an infinite layer as well, because due to the higher salt content, the penetration depth is even smaller.Continental ice cover over Antartica is more challenging to model at low microwave frequencies because, despite is thickness (several kilometers) the ice is almost pure ice, with very large penetration depths, and the microwave signatures depend on the presence of underground lakes of fresh water, as it has been observed with ESA SMOS and NASA Aquarius missions [37].At higher frequencies, the penetration depth is smaller and the ice layer appears as an infinite layer of pure ice.

Earth's Atmosphere Contribution
Up to 1 THz, the specific attenuation due to dry air and water vapor, can be most accurately evaluated at any pressure, temperature and humidity by summation of the individual resonance lines from oxygen and water vapor, together with small additional factors for the non-resonant Debye spectrum of oxygen below 10 GHz, pressure-induced nitrogen attenuation above 100 GHz, and a wet continuum to account for the excess water vapor-absorption found experimentally.Below 100 GHz, the absorption due to gas constituents is mainly contributed by water vapor and oxygen, around 22 and 55-60 GHz respectively (Figure 6).Near 60 GHz, at sea-level pressures, many oxygen absorption lines merge together to form a single, broad absorption band.At higher altitudes, the individual lines of the oxygen attenuation become resolved at lower pressures.Some additional molecular species (e.g., oxygen isotopic species, oxygen vibrationally excited species, ozone, ozone isotopic species, and ozone vibrationally excited species, and other minor species) are not included in the line-by-line prediction method.These additional lines are insignificant for typical atmospheres, but may be important for a dry atmosphere.model of the vertical structure of the atmosphere with measured temperature, pressure and water vapor density profiles [20,28]: (77) In the RTM module of SAIRPS the atmospheric density vertical profile is actually computed from the atmospheric temperature, pressure, and relative humidity profiles as described in [38].The integration of Equations ( 72)-(74) suffer from serious stability problems, especially around the oxygen absorption bands, where the predicted Tup and Tdn vary significantly with the grid step and frequency.The problem was solved using an efficient 96 Gaussian quadrature integration up to 72 km height.The values of the geophysical parameters in the grid points are interpolated from the values given by the pressure levels in the auxiliary data files.
When computing the up-and down-welling atmospheric brightness temperatures at low frequencies, the attenuation is very small and, paradoxically, it is prone to large relative (as opposed to absolute) errors.At L-band (1.4 GHz), the following alternative formulations for the atmospheric attenuation and up-welling atmospheric contribution derived from NASA/Aquarius observations [39] are used instead, which accurately predicts the up-welling atmospheric temperature as well as the atmospheric losses up to 70° incidence angle.
with: In the scope of an atmospheric model to be included in the SAIRPS framework, a model that includes gas attenuation (water vapor and oxygen), rain and water or ice clouds is implemented, while scattering effects have been neglected.
The atmospheric upwelling and down-welling emitted radiation can be computed as [20]: T dn pθ, Hq " where k a pzq is the absorption coefficient at a height z 1 (equal to the emissivity in thermodynamic equilibrium), and can be decomposed in three terms: the contribution of the different gas constituents, k gas and eventually the extinction coefficients of clouds and rain precipitation, k clouds and k rain : Except for water vapor variations, the relative composition of the atmosphere is quite stable up to 90 km above sea level.For mid-latitudes, the 1962 US Standard Atmosphere gives a generalized model of the vertical structure of the atmosphere with measured temperature, pressure and water vapor density profiles [20,28]: P pzq " P 0 ¨e´z rkms 7.7 (77) In the RTM module of SAIRPS the atmospheric density vertical profile is actually computed from the atmospheric temperature, pressure, and relative humidity profiles as described in [38].The integration of Equations ( 72)-( 74) suffer from serious stability problems, especially around the oxygen absorption bands, where the predicted T up and T dn vary significantly with the grid step and frequency.The problem was solved using an efficient 96 Gaussian quadrature integration up to 72 km height.The values of the geophysical parameters in the grid points are interpolated from the values given by the pressure levels in the auxiliary data files.
When computing the up-and down-welling atmospheric brightness temperatures at low frequencies, the attenuation is very small and, paradoxically, it is prone to large relative (as opposed to absolute) errors.At L-band (1.4 GHz), the following alternative formulations for the atmospheric attenuation and up-welling atmospheric contribution derived from NASA/Aquarius observations [39] are used instead, which accurately predicts the up-welling atmospheric temperature as well as the atmospheric losses up to 70 ˝incidence angle.
T up " t2.3058 ´3.2699¨10 ´3¨pT ph DEM q ´273.15q`4.2328¨10 ´3 ¨pP ph DEM q ´900q `1.4417¨10Since the above model does not provide an estimate for T dn , it can be assumed that T dn = T up , which is very accurate at L-band.Finally, to account for the Earth's curvature effects, the effective incidence angle is computed as follows: with R Earth the radius of the Earth.Figure 7a,b show the evolution of the effective incidence angle and the difference between the true and the effective angles as a function of the geometrical incidence angle.As it can be appreciated, the difference is smaller than 1 ˝up to ~83 ˝, and smaller than 2 ˝up to ~87 ˝.Therefore, except for a small corona at near grazing angles, this correction is negligible.
with  ℎ the radius of the Earth.Figure 7a,b show the evolution of the effective incidence angle and the difference between the true and the effective angles as a function of the geometrical incidence angle.As it can be appreciated, the difference is smaller than 1° up to ~83°, and smaller than 2° up to ~87°.Therefore, except for a small corona at near grazing angles, this correction is negligible.

A) Water Vapor Attenuation
In the range 1 < f < 1000 GHz, temperature ´100 ˝C < T < 50 ˝C, pressure 10 ´5 mb < P < 1013 mbar, and water vapor density 0 < ρ v < 20 g/m 3 , the expressions providing the water vapor and oxygen attenuation are ( [28], Section 8-2.3 and Equation (8.19c)): e " ρ 0 0.7223¨ϑ (84) k s " 0.444¨ϑ 7.5  (86) f " 0.0145¨ϑ 4.5  (87) where P d is the dry air pressure, e is the water vapor partial pressure, S i is the strength of the i-th line, F i is the line shape factor, d is the width parameter for the Debye spectrum, N 2 d is the dry continuum due to pressure-induced nitrogen absorption and the Debye spectrum, and the sum extends over all the lines.The above equations are evaluated with the parameters given in Table A1.

B) Oxygen Attenuation
In the range 1 < f < 1000 GHz, temperature ´100 ˝C < T < 50 ˝C, pressure 10 ´5 mb < P < 1013 mbar, and water vapor density 0 < ρ v < 20 g/m 3 , the expressions providing the oxygen attenuation are [28] δ i " pe i `fi ¨ϑq ¨P¨ϑ 0.8 (103) These equations are very similar to those for the water vapor absorption, but the δ factor in the line shape factor (F i ) is included to correct for the interference effects in the oxygen lines.The above equations are evaluated with the parameters given in Table A2.

C) Rain Attenuation
The vertical distribution of rain roughly extends up to 4 km and a logarithmic model is used to compute the extinction coefficient for h/v polarizations vs. height at each frequency.The extinction coefficient (scattering plus absorption) due to rain at any frequency 0 < f < 1000 GHz and rain rate 0 < R r < 150 mm/h, and T > 0 ˝C can be computed using the following formulas [28], Sections 8-8.2.
However, these formulas do not account for the differential attenuation at horizontal and vertical polarizations.Instead the Rec.ITU-R P.838-2 "Specific attenuation model for rain for use in prediction methods" could be used, but the formulas provided are only deemed to be valid up to 55 GHz.Instead, we interpolate the coefficients given in the Table A3 for any particular frequency.

D) Water Clouds Attenuation
Water clouds extend roughly up to 10 km, but their structure varies.The attenuation coefficient of water clouds is computed using Equation (110), using the absorption coefficient of the pure water (as a function of the frequency f in [GHz], and the temperature T in [ ˝C]) and the cloud water content m v in [g/m 3 ] [28], Section 8-7.3.

E) Ice Clouds Attenuation
Similarly, the attenuation coefficient of ice clouds is computed using Equation (108), using the absorption coefficient of the pure water (as a function of the frequency f in [GHz], and the temperature T in [ ˝C]) and the cloud water content m v in [g/m 3 ] [28], Section 8-7.3.
If the cloud and rain structure are known (vertical profiles of m v , T, and R r , e.g., Figure 8) then, Equations ( 108), ( 110), (112) can be directly applied, if they are not known or not available, a default atmosphere model with a "homogeneous" cloud from 1 to 4 km height and a concentration of m v = 5 g/m 3 is used, and for the rain cell, a "homogeneous" rain cell from 0 to 1 km.absorption coefficient of the pure water (as a function of the frequency f in [GHz], and the temperature T in [°C]) and the cloud water content mv in [g/m 3 ] [28], Section 8-7.3.
If the cloud and rain structure are known (vertical profiles of mv, T, and Rr, e.g., Figure 8) then, Equations ( 108), ( 110), (112) can be directly applied, if they are not known or not available, a default atmosphere model with a "homogeneous" cloud from 1 to 4 km height and a concentration of mv = 5 g/m 3 is used, and for the rain cell, a "homogeneous" rain cell from 0 to 1 km.

Cosmic, Isotropic, and Galactic Noises
The sky noise includes three different contributions [40,41]): • The cosmic background is fairly constant and equal to   = 2.275 K.

•
The isotropic component (extra-galactic term) decreases with frequency roughly as: • The galactic component also decreases with frequency roughly as: where β lies between 2.5 and 2.6 between 45 and 408 MHz, and it is equal to 2.81 ± 0.16 from 1.375 to 7.5 GHz.A Global Sky Model derived from all publicly available unpolarized large-area radio surveys can be found at [42].In the framework of this simulation tool the  0 = 1420 MHz map used in the SMOS processors, scaled according to Equation (118) with β = 2.81, is used (Figure 9).Once the central frequency is defined, and a Sky radio map is known at a given frequency f0, the three contributions can be immediately computed.

‚
The cosmic background is fairly constant and equal to T cos " 2.275 K.

‚
The isotropic component (extra-galactic term) decreases with frequency roughly as:

‚
The galactic component also decreases with frequency roughly as: where β lies between 2.5 and 2.6 between 45 and 408 MHz, and it is equal to 2.81 ˘0.16 from 1.375 to 7.5 GHz.A Global Sky Model derived from all publicly available unpolarized large-area radio surveys can be found at [42].In the framework of this simulation tool the f 0 " 1420 MHz map used in the SMOS processors, scaled according to Equation (118) with β = 2.81, is used (Figure 9).Once the central frequency is defined, and a Sky radio map is known at a given frequency f 0 , the three contributions can be immediately computed.In the case of a non-flat land surface, radiation emitted by some elevated parts of the relief incident on another part of the surface may scatter towards the space-borne radiometer.This radiation is shadowing the radiation from the "hidden" sky.The local angle of incidence θH of the horizon for a given surface facet determines the range of angles from which the sky radiation arrives to the facet (θl ≤ θH).For larger values, the incident radiation comes from the land at the brightness temperature Tsr.
An accurate computation of the effective reflectivity can be very difficult and cumbersome for facets with general bidirectional scattering.However, assuming a Lambertian terrain (with reflectivity Гd) then the emission is not polarized and direction independent.The total upwelling brightness temperature Tup,relief is the sum of the radiation from a flat horizon, Tup,flat plus an increase ΔTup ( [30], Chapter 4).
An upper limit Δ , can be estimated by assuming that the elevated surface is a black body at constant temperature (Th = T0) whereas the flat surface facet is a non-black Lambert scatterer at the same physical temperature.With this simplification where c is the azimuthal average of cos 2   .The effective reflectivity Γ eff can now be evaluated as: Finally,

B) Over the ocean
The approach over the ocean is different, because this effect has already been numerically evaluated in [29]: it is called the atmospheric path length correction.It is typically parameterized as where Ω  ( ,  = 0 / ) = 0 and Ω  (  = 0,  ) = 0 .This automatically guarantees that Δ , vanishes for a smooth surface (w = 0 m/s) and for a completely opaque (τ = 0) and a completely In the case of a non-flat land surface, radiation emitted by some elevated parts of the relief incident on another part of the surface may scatter towards the space-borne radiometer.This radiation is shadowing the radiation from the "hidden" sky.The local angle of incidence θ H of the horizon for a given surface facet determines the range of angles from which the sky radiation arrives to the facet (θ l ď θ H ). For larger values, the incident radiation comes from the land at the brightness temperature T sr .
An accurate computation of the effective reflectivity can be very difficult and cumbersome for facets with general bidirectional scattering.However, assuming a Lambertian terrain (with reflectivity Г d ) then the emission is not polarized and direction independent.The total upwelling brightness temperature T up,relief is the sum of the radiation from a flat horizon, T up,flat plus an increase ∆T up ( [30], Chapter 4).
T up, relie f " T up, f lat `∆T up (119) An upper limit ∆T up,max can be estimated by assuming that the elevated surface is a black body at constant temperature (T h = T 0 ) whereas the flat surface facet is a non-black Lambert scatterer at the same physical temperature.With this simplification where c is the azimuthal average of cos 2 θ H .The effective reflectivity Γ eff can now be evaluated as:

B) Over the Ocean
The approach over the ocean is different, because this effect has already been numerically evaluated in [29]: it is called the atmospheric path length correction.It is typically parameterized as ∆T SC,p " Ω p pτ, wq where Ω p pτ, w " 0 m{sq " 0 and Ω p pτ " 0, wq " 0. This automatically guarantees that ∆T SC,p vanishes for a smooth surface (w = 0 m/s) and for a completely opaque (τ = 0) and a completely transparent (τ = 1) atmosphere.Opaque and transparent atmospheres are isotropic, and therefore no atmospheric path length correction exists.The values of Ω p for different Earth incidence angles, frequencies, polarization and atmospheric opacities, for a reference atmosphere at 281 K can be found in [11].

RTM Model Outputs
Finally, once the brightness temperatures for the different covers (land, ocean, mixed), the atmospheric contributions, and the scattered down-welling radiation over the land and/or the ocean have been computed for each pixel in the brightness temperature map, they are weighted (summed) with a weight equal to the fractional area.It is important that the pixel size be much smaller that the antenna footprint (in a real aperture radiometer) or the spatial resolution (in a synthetic aperture radiometer) in order to properly integrate their effects within the beam (either real or synthetic).
Finally, the total signal collected by the radiometer antenna at a given polarization is the beam-weighted sum over the radiation from all visible facets within the antenna footprint: with dΩ " A¨cos pθ l q R 2 " A h ¨cos pθ l q R 2 ¨cos pαq (125) being A the true area of the facet, A h the area of the facet projected in the horizontal plane, and R the distance from the radiometer antenna to the facet.Actually the complete antenna pattern, including the co-and cross-polar antenna patterns are included [43].At this moment, the apparent brightness temperatures at TOA in the Earth's surface (local) reference frame have been calculated, including both the Earth's surface and the atmospheric contributions.Now, these brightness temperatures have to be interfaced with the BT Maps module to produce observables as similar as possible to those that will be produced by the sensors.
Due to the relative orientation (α) between the local reference frame (over the Earth's surface: vor hpolarizations) and the antenna reference frame (X, Y), and (eventually) the Faraday rotation (φ Faraday ) the electric fields incident in the antennas are rotated by an angle: where A " cos ´α `ϕFaraday ¯, B " sin ´α `ϕFaraday ¯.Therefore, since the polarimetric brightness temperatures T pq are proportional to , the brightness temperatures in the antenna reference frame at a given polarization (X or Y) can be expressed as a linear combination of the brightness temperatures in the Earth's pixel reference frame (v or h).
where T vh " T hv " pT 3 `j¨T 4 q {2 and T 3 and T 4 are the brightness temperatures associated to the third and fourth Stokes parameters in the Earth's reference frame.As compared to T hh and T vv , T 3 and T 4 have negligible values (except for a few Kelvin in the ocean).Therefore, if (T xy )= (T yx )=0 and (T xy )= (T yx ), and T vh " T hv " T 3 " 0. Equation (124) reduces to: and a non-zero third Stokes parameter appears at the TOA even if it does not exist at the BOA.

Input Data Collection and Strategy for Time-Space-Interpolation
Regarding the generation of the Earth's surface emission and scattering models described in the previous sections, a subset of the ECMWF Auxiliary Data Files used for SMOS is used.Among the 50 parameters (x 4 dates, total size, approx.8 Gbytes) some are constant (e.g., soil properties, land sea mask, land use . . .), others derived from satellite data (e.g., LAI from MODIS), while others from NWP models available from ECMWF.Since at the frequencies involved the impact of sea surface salinity in the brightness temperatures is negligible, only climatological values are used.
The parameters are the following: Other sources of information required are the International Geomagnetic Reference Field [44], International Reference Ionosphere-IRI-2012 [45] and the Low Frequency Radio Sky Map [42].Some of the above parameters (e.g., land use) are static or quasi static, and do not require temporal interpolation.ECMWF parameters can be linearly interpolated both in time (if close enough) and space to the grid cell points of a Gaussian grid used for the final brightness temperature maps (in 3 resolutions (125 km, 40 km and 15 km), and stored in NetCDF format.Alternatively, data can be downloaded via ftp on demand, as needed, since it is likely that the final BT map has a higher resolution than the input parameters, maps are bi-linearly interpolated, taking into account the presence (or not) of a high resolution coastline, to avoid blurring near the coastline.

Simulation Results and Discussion
In the previous section a detailed description of the models used to create synthetic brightness temperature (BT) images to be used in the simulation of SAIRs was presented.In this section, sample synthetic BT images are presented to illustrate them.
Figure 10 shows the simulated full-polarimetric BTs at 1.4 GHz over the Pacific Ocean in the Earth's surface reference frame, as seen from a LEO orbit corresponding to a SMOS-like instrument with a tilt angle of 32 ˝, and for summer and winter seasons.As it can be appreciated, because of the array tilting, the sky appears in the top of the images.In addition, both the vand h-polarizations are the same at nadir (~90 K), at h-pol it decreases with increasing incidence angle, and at v-pol it increases, but then, close the 90 ˝it sharply decreases, which is clear when approaching the Earth-sky horizon.Variations due to salinity and temperature are much smaller than the incidence angle variation.Note that the LEO (low Earth orbit) BTs are on-ground values, and as such, there is a T3 signal in both the sky and sea, and T4 only for sea reflection.10, but over the Iberian Peninsula.These BT images present physically meaningful values both over land and sea, although some interpolation artifacts close to the coastline are evident due to the coarser resolution of the ocean parameters.BTs are on the Earth's reference frame, and as such, there is a T3 signal in both the sky and sea and T4 only for sea reflection.No model for the third and fourth Stokes elements emission over land exist.A difference between summer and winter BTs is observed, due to the different physical temperatures and soil moisture values over land.Drier soils have a larger emissivity.Note the brightness temperature changes in mountainous regions, such as the Pyrenees, with an increase at h-pol, and a decrease at v-pol.
Figure 12 is similar to the previous ones, but over the Amazon river plume.It is important to note the higher and closer values of the BTs at vand h-polarizations due to the dense vegetation cover, the increase in the BTs very close to the river plume, not due to a variation of the local incidence angle, but to the fresh water discharged by the Amazon river, which has a higher emissivity as compared to the saltier water of the Atlantic Ocean.Lastly, it is worth noting that there is a thin horizontal strip of increased BT that corresponds to the reflected noise coming from the galactic plane.It is more noticeable at h-polarization, than at v-polarization, because of the higher reflectivity at h-polarization, than at v-polarization.
Figure 13 is similar to the previous ones, but over Siberia (snow, sea ice, land, and ocean scenario), and at 36.5 GHz to better illustrate the emission of ice/snow and sea, specially when sea ice melts in the summer season (left column) and it is frozen during the winter (right column).BTs are on the Earth's surface local reference plane, and as such, there is a T3 signature (third Stokes parameter) in both the sky and sea, while T4 (fourth Stokes signature) only for a sea reflection.
Figure 14 is similar to the previous ones, but over Antarctica (snow, sea ice, and ocean scenario), and at 18.7 GHz to better illustrate the emission of sea ice vs sea, specially when sea ice melts in the Austral summer season ("winter", right column), and it is frozen during the Austral winter ("summer", left column).Again, BTs are on the Earth's surface local reference plane and as such, there is a T3 signal in both the sky and sea, and T4 only for sea reflection.Note that the galactic plane ("milky way") is present in the background sky and appears as a thin brighter strip.
Figure 15 shows the simulated full-polarimetric BTs at 1.4 GHz from a nadir pointing instrument (tilt angle = 0 ˝) in a GEO (geostationary orbit) at 0 ˝longitude, and for summer and winter seasons.These BT images present physically meaningful values over the ocean and land (bare soil such as the Saharan desert, or over the Equatorial rain forests).In these images the BTs are computed in the instrument reference plane (antennas aligned with the X and Y directions), and as such in the vertical direction the X-pol image corresponds to the horizontal polarization, and in the horizontal direction, the X-pol image corresponds to the vertical polarization, and vice-versa.In addition, as such, there is a strong T3 signal coming from the geometric rotation (α angle in Equation ( 123)), and a much lower signal in T4.Note that there is T3 in the sky, but no T4 (from the Galaxy Map).The temperatures are in line with the 1.4 GHz expected range.
Figure 16 is similar to Figure 14, but at 6.9 GHz.Again, these images present physical meaningful values in all quantities.At this frequency, the atmosphere is also almost transparent.As in Figure 14, BT images are computed in the instrument reference plane values, and as such, there is a strong T3 signal coming from the geometrical rotation and a much lower signal in T4.
Figure 17 is similar to Figures 14 and 15 but at 53.6 GHz, which corresponds to Channel 5 of AMSU [46], a band with high atmospheric attenuation used for tropospheric temperature measurements.For this reason, the BT images in winter look like the ones in summer, but flipped in the vertical direction, and T3 and T4 values are much lower because of the strong attenuation of the emission coming from the Earth's surface.
Figure 18 is similar to Figures 14-16 but at 60 GHz, in the center of the Oxygen absorption band, which corresponds approximately to channels 20-24 of SSMIS [47].The behavior is very similar to Figure 16, but it senses temperatures at higher altitudes (lower values).T3 and T4 values are also much lower than for surface signals because of the strong atmospheric attenuation.Agency (ESA) by the Universitat Politècnica de Catalunya-BarcelonaTech and DEIMOS Engenharia S.A.This implementation includes a comprehensive data set of geophysical variables, emission models for land/cryosphere, ocean, atmosphere, and sky radiation.Fully polarimetric brightness temperatures (Stokes vector) can be computed at the local Earth's surface reference frame or at the antenna polarization reference frame, in this last case accounting for the geometric and Faraday rotation effects.The radiometric accuracy of the generated brightness temperatures has been validated through an extensive validation campaign.This allows us to use the generated brightness temperature maps (at each polarization) as inputs for the SAIRPS Instrument, Calibration, and Image Reconstruction Modules, described in the second part of this manuscript, and to properly assess the performance of algorithms developed to retrieve geophysical parameters from radiometric observations.
In terms of requirements and performance, the simulator is run on an Intel Core i7 CPU M620 at 2.67 GHz with 8 cores and 8 GBytes of RAM, running Windows 7. A total of 8 GBytes are required to store the 50 input ECMWF variables for 4 dates.Using 1, 4, or 8 cores the computation times of the geometry module are 14 s, 5 s, and 3 s, respectively; and those of the Radiative Transfer Module are 907 s, 244 s, and 142 s, respectively.

Figure 1 .
Figure 1.Generic design of a Synthetic Aperture Interferometric Radiometer Performance Simulator (SAIRPS) with the Radiative Transfer Mode l (RTM) module .

Figure 1 .
Figure 1.Generic design of a Synthetic Aperture Interferometric Radiometer Performance Simulator (SAIRPS) with the Radiative Transfer Model (RTM) module.

Figure 2 .
Figure 2. De sign of the RTM and Brightness Temperature (BT) Maps module .

Figure 2 .
Figure 2. Design of the RTM and Brightness Temperature (BT) Maps module.

Figure 4 .
Figure 4. Sample dielectric constant of dry snow v i = 0.5 (a: real part, b: imaginary part); and wet snow with P s = 0.2 g/cm 3 , m v = 5% (c: real part, d: imaginary part).

Figure 8 .
Figure 8. Sample cloud vapor distribution vs. he ight (left) and Rain distribution vs. he ight (right).

Figure 9 .
Figure 9. Sky Radio Map at 1420 MHz use d in SAIRPS: (left) First Stoke s parame te r, (center) Se cond Stokes parameter, and (right) Third Stoke s parame te r (all in [K]).Data are sampled with a 0.25° × 0.25° re solution in de clination × right asce nsion (e quatorial coordinate s).Se nsitivity de fine d as three time s the rms brightness te mperature noise is 0.05 K. 2.2.4.Down-Welling Temperature Scattered Towards a Space-Borne Radiometer A) Over the land

Figure 9 .
Figure 9. Sky Radio Map at 1420 MHz used in SAIRPS: (left) First Stokes parameter, (center) Second Stokes parameter, and (right) Third Stokes parameter (all in [K]).Data are sampled with a 0.25 ˝ˆ0.25 resolution in declination ˆright ascension (equatorial coordinates).Sensitivity defined as three times the rms brightness temperature noise is 0.05 K.

Figure 11
Figure 11  is similar to Figure10, but over the Iberian Peninsula.These BT images present physically meaningful values both over land and sea, although some interpolation artifacts close to the coastline are evident due to the coarser resolution of the ocean parameters.BTs are on the Earth's reference frame, and as such, there is a T3 signal in both the sky and sea and T4 only for sea reflection.No model for the third and fourth Stokes elements emission over land exist.A difference between summer and winter BTs is observed, due to the different physical temperatures and soil moisture values over land.Drier soils have a larger emissivity.Note the brightness temperature changes in mountainous regions, such as the Pyrenees, with an increase at h-pol, and a decrease at v-pol.Figure12is similar to the previous ones, but over the Amazon river plume.It is important to note the higher and closer values of the BTs at vand h-polarizations due to the dense vegetation cover, the increase in the BTs very close to the river plume, not due to a variation of the local incidence angle, but to the fresh water discharged by the Amazon river, which has a higher emissivity as compared to the saltier water of the Atlantic Ocean.Lastly, it is worth noting that there is a thin horizontal strip of increased BT that corresponds to the reflected noise coming from the galactic plane.It is more noticeable at h-polarization, than at v-polarization, because of the higher reflectivity at h-polarization, than at v-polarization.Figure13is similar to the previous ones, but over Siberia (snow, sea ice, land, and ocean scenario), and at 36.5 GHz to better illustrate the emission of ice/snow and sea, specially when sea ice melts in the summer season (left column) and it is frozen during the winter (right column).BTs are on the Earth's surface local reference plane, and as such, there is a T3 signature (third Stokes parameter) in both the sky and sea, while T4 (fourth Stokes signature) only for a sea reflection.Figure14is similar to the previous ones, but over Antarctica (snow, sea ice, and ocean scenario), and at 18.7 GHz to better illustrate the emission of sea ice vs sea, specially when sea ice melts in the Austral summer season ("winter", right column), and it is frozen during the Austral winter ("summer", left column).Again, BTs are on the Earth's surface local reference plane and as such, there is a T3 signal in both the sky and sea, and T4 only for sea reflection.Note that the galactic plane ("milky way") is present in the background sky and appears as a thin brighter strip.Figure15shows the simulated full-polarimetric BTs at 1.4 GHz from a nadir pointing instrument (tilt angle = 0 ˝) in a GEO (geostationary orbit) at 0 ˝longitude, and for summer and winter seasons.These BT images present physically meaningful values over the ocean and land (bare soil such as the Saharan desert, or over the Equatorial rain forests).In these images the BTs are computed in the instrument reference plane (antennas aligned with the X and Y directions), and as such in the vertical direction the X-pol image corresponds to the horizontal polarization, and in the horizontal direction, the X-pol image corresponds to the vertical polarization, and vice-versa.In addition, as such, there is a strong T3 signal coming from the geometric rotation (α angle in Equation (123)), and a much lower signal in T4.Note that there is T3 in the sky, but no T4 (from the Galaxy Map).The temperatures are in line with the 1.4 GHz expected range.Figure16is similar to Figure14, but at 6.9 GHz.Again, these images present physical meaningful values in all quantities.At this frequency, the atmosphere is also almost transparent.As in Figure14, BT images are computed in the instrument reference plane values, and as such, there is a strong T3 signal coming from the geometrical rotation and a much lower signal in T4.Figure17is similar to Figures 14 and 15 but at 53.6 GHz, which corresponds to Channel 5 of AMSU[46], a band with high atmospheric attenuation used for tropospheric temperature measurements.For this reason, the BT images in winter look like the ones in summer, but flipped in the vertical direction, and T3 and T4 values are much lower because of the strong attenuation of the emission coming from the Earth's surface.Figure18is similar to Figures 14-16 but at 60 GHz, in the center of the Oxygen absorption band, which corresponds approximately to channels 20-24 of SSMIS[47].The behavior is very similar to Figure16, but it senses temperatures at higher altitudes (lower values).T3 and T4 values are also much lower than for surface signals because of the strong atmospheric attenuation.

Figure 19 Figure 11 .
Figure19is similar to Figures 14-17 but at 89 GHz, in the center of a transmission band between the Oxygen absorption band around 60 GHz, and the 118 GHz Oxygen absorption line, which corresponds to channel 15 of AMSU[46], used for cloud top and thin snow monitoring.At 89 GHz the atmosphere is more transparent than in the previous two cases, and there is a measurable signal from the surface.

Figure 11 .Figure 12 .
Figure 11.Simulate d full-polarime tric BTs at 1.4 GHz ove r the Ibe rian Pe ninsula, in the Earth's surface re fe rence frame, as seen from a LEO (low Earth orbit) such as for an SMOS-like instrume nt.Figure 11.Simulated full-polarimetric BTs at 1.4 GHz over the Iberian Peninsula, in the Earth's surface reference frame, as seen from a LEO (low Earth orbit) such as for an SMOS-like instrument.

Figure 12 .Figure 13 .
Figure 12.Simulate d full-polarime tric BTs at 1.4 GHz ove r the Amazon rive r plume , in the Earth's surface re ference frame, as seen from a LEO (low Earth orbit) such as for an SMOS-like instrume nt.Figure 12. Simulated full-polarimetric BTs at 1.4 GHz over the Amazon river plume, in the Earth's surface reference frame, as seen from a LEO (low Earth orbit) such as for an SMOS-like instrument.

Figure 13 .
Figure 13.Simulate d full-polarime tric BTs at 36.5 GHz ove r the Sibe ria (snow, se a ice , oce an, and land sce nario), in the Earth's surface reference frame, as seen from a LEO (low Earth orbit).Figure 13.Simulated full-polarimetric BTs at 36.5 GHz over the Siberia (snow, sea ice, ocean, and land scenario), in the Earth's surface reference frame, as seen from a LEO (low Earth orbit).

Figure 14 .
Figure 14.Simulate d full-polarime tric BTs at 18.7 GHz ove r the Sibe ria (snow, se a ice , and oce an sce nario), in the Earth's surface reference frame, as seen from a LEO (low Earth orbit).Figure 14.Simulated full-polarimetric BTs at 18.7 GHz over the Siberia (snow, sea ice, and ocean scenario), in the Earth's surface reference frame, as seen from a LEO (low Earth orbit).

Figure 14 .Figure 15 .
Figure 14.Simulate d full-polarime tric BTs at 18.7 GHz ove r the Sibe ria (snow, se a ice , and oce an sce nario), in the Earth's surface reference frame, as seen from a LEO (low Earth orbit).Figure 14.Simulated full-polarimetric BTs at 18.7 GHz over the Siberia (snow, sea ice, and ocean scenario), in the Earth's surface reference frame, as seen from a LEO (low Earth orbit).

Figure 15 .Figure 16 .
Figure 15.Simulate d full-polarimetric BTs at 1.4 GHz, in the antenna reference frame, as seen from a GEO (ge ostationary orbit).Figure 15.Simulated full-polarimetric BTs at 1.4 GHz, in the antenna reference frame, as seen from a GEO (geostationary orbit).

Figure 16 .Figure 17 .
Figure16.Simulate d full-polarimetric BTs at 6.9 GHz, in the antenna reference frame, as seen from a GEO (ge ostationary orbit).Figure16.Simulated full-polarimetric BTs at 6.9 GHz, in the antenna reference frame, as seen from a GEO (geostationary orbit).

Figure 17 .
Figure 17.Simulate d full-polarimetric BTs at 53 GHz, in the antenna reference frame, as seen from a GEO (ge ostationary orbit, same as Figure s 16 and 17).Figure 17.Simulated full-polarimetric BTs at 53 GHz, in the antenna reference frame, as seen from a GEO (geostationary orbit, same as Figures 16 and 17).
Figure 17.Simulate d full-polarimetric BTs at 53 GHz, in the antenna reference frame, as seen from a GEO (ge ostationary orbit, same as Figure s 16 and 17).Figure 17.Simulated full-polarimetric BTs at 53 GHz, in the antenna reference frame, as seen from a GEO (geostationary orbit, same as Figures 16 and 17).

Figure 18 .
Figure 18.Simulate d full-polarimetric BTs at 60 GHz, in the antenna reference frame, as seen from a GEO (ge ostationary orbit, same as Figure s 16 and 17).Figure 18. Simulated full-polarimetric BTs at 60 GHz, in the antenna reference frame, as seen from a GEO (geostationary orbit, same as Figures 16 and 17).

Figure 18 .
Figure 18.Simulate d full-polarimetric BTs at 60 GHz, in the antenna reference frame, as seen from a GEO (ge ostationary orbit, same as Figure s 16 and 17).Figure 18. Simulated full-polarimetric BTs at 60 GHz, in the antenna reference frame, as seen from a GEO (geostationary orbit, same as Figures 16 and 17).
Figure 18.Simulate d full-polarimetric BTs at 60 GHz, in the antenna reference frame, as seen from a GEO (ge ostationary orbit, same as Figure s 16 and 17).Figure 18. Simulated full-polarimetric BTs at 60 GHz, in the antenna reference frame, as seen from a GEO (geostationary orbit, same as Figures 16 and 17).

Table 1 .
Fitting coefficients to estimate topography effects at horizontal and vertical polarizations.

Table A2 .
Parameters required to compute the oxygen attenuation 0-1THz.