Abstract
The design and optimization of propagation impairment techniques for space telecommunication systems operating at frequencies above 20 GHz require a precise knowledge of the propagation channel both in space and time. For that purpose, space-time channel models have to be developed. In this paper the description of a model for the simulation of long-term rain attenuation time series correlated both in space and time is described. It relies on the definition of a stochastic rain field simulator constrained by the rain amount outputs of the ERA-40 reanalysis meteorological database. With this methodology, realistic propagation conditions can be generated at the scale of satellite coverage (i.e., over Europe or USA) for many years. To increase the temporal resolution, a stochastic interpolation algorithm is used to generate spatially correlated time series sampled at 1 Hz, providing that way valuable inputs for the study of the performances of propagation impairment techniques required for adaptive SatCom systems operating at Ka band and above.
1. Introduction
With the congestion of conventional frequency bands such as C (4–6 GHz) or Ku (11–14 GHz) band and the need to convey higher data rates for multimedia services, new SatCom systems are progressively pushed towards the use of higher-frequency bands such as Ka (20–30 GHz) or Q/V band (40–50 GHz) where larger bandwidth is available. However, at those frequencies, strong tropospheric impairments occur on Earth-space links with a significant impact on the system quality of service. The latter are due to rain, clouds, water vapor, and oxygen. The resulting attenuation cannot be mitigated simply acting on the static margin taken on the overall link budget as it will result in a major loss of efficiency. In order to perform a good trade-off between quality of service, system capacity, and link availability, adaptive resource allocation mechanisms using Fade Mitigation Techniques (FMTs) have to be implemented [1].
FMT implementation leads also to consider the system resource allocation and therefore upper layer issues that can be studied through protocol simulations. Network simulations for SatCom systems operating at Ka and Q/V bands have to take into account the influence of the propagation channel, not only in terms of dynamics but also in terms of the spatial variations [2]. Therefore, in order to assess FMT efficiency or the resulting system availability, the propagation channel has to be simulated both from a temporal and from a spatial point of view.
To emulate channel dynamic, attenuation time series synthesizers [3–7] have been developed. They are crucial to design and optimize Uplink Power Control (ULPC) systems or data rate reduction mechanisms [8]. Nevertheless, the knowledge of the temporal dynamic of the propagation channel is not sufficient to implement all the FMTs. For instance, the optimization of site diversity systems [9] or adaptive TDMA (Time Division Multiple Access) satellite systems requires also the knowledge of the spatial variability of the propagation impairments.
Several models and methods exist to depict the spatial variability of tropospheric impairments, especially the attenuation due to rain. Indeed, attenuation fields and the underlying rain fields can be modeled by rain cells [10–16], by fractals [17] or random fields [5, 18]. Even if the spatiotemporal description of the propagation channel is of great interest for the design and optimization of the previously mentioned FMTs, it is not adapted to solve the issues raised by onboard radio resource management for which a description both in space and in time of the propagation channel over the whole satellite coverage is required.
Up to now, few models exist to combine the temporal and the spatial description of the propagation channel. Authors in [5] have adapted the model of rain fields correlated in space and in time developed by [19], while those in [20] have proposed to correlate scaled time series from measurements made during the ITALSAT campaign. The main limitation of these models lies in their temporal or spatial range of validity. Indeed, the model proposed by [5] is realistic for short durations (some hours) and small areas (less than ~300 × 300 km²) due to the stationarity assumption necessary to construct the random field in the Fourier plane. Nevertheless, this approach is unable to reproduce the alternation of clear sky and rainy periods, preventing that way the long-term simulation of propagation conditions. As for the model of correlated time series presented in [19], it is limited to duration of one hour. However, to be statistically reliable, network simulations require attenuation data continuously acquired during several years and over the whole satellite coverage. Therefore, propagation inputs correlated in space and in time for large areas and long durations are needed.
The aim of this paper is to extend the range of validity of the stochastic rainfall field simulator described in [19] and adapted for rain attenuation in [5] from the midscale (~300 × 300 km²) to the continental scale (Europe, USA) and from some hours to many years. To do it, rain amount time series provided by general circulation models are used. Particularly, a methodology to constrain the model with inputs from the low-resolution reanalysis database ERA-40 (ECMWF reanalysis for 40 years) provided by the ECMWF (European Centre for Medium-range Weather Forecast) is presented. The use of ERA-40 data ensures a realistic long term evolution of the rain fields as they reproduce large-scale conditions from past meteorological situations.
The basics and limitations of the stochastic modelling of rain rate field presented in [19] are recalled in the first part of the paper. A set of parameters more suitable for temperate areas is then determined from a radar dataset acquired by the meteorological radar of Bordeaux (South-Western France, midlatitude oceanic climate). In a second part of the paper, a methodology to parameterize the model with inputs from the ECMWF ERA-40 reanalysis database is developed. Some statistical properties of the simulations are lastly compared to statistics derived from weather radar datasets, ITU-R (International Telecommunication Union section Radiocommunications) Recommendations and to statistics obtained from beacon measurements.
2. Spatiotemporal Modelling of Rain Fields
In this section, the basics of the rain field space-time stochastic modelling firstly presented in [19] are detailed. This approach was originally developed to study the performances of satellite rainfall sensors (such as the precipitation radar of the Tropical Rainfall Measuring Mission TRMM) from simulated rain fields at midscale (~300 × 300 km²) in which statistical characteristics (spatiotemporal correlation structure, local climatologic Complementary Cumulative Distribution Function (CCDF) of rain rates) must be as close as possible to the ones of real rain fields.
The rain fields are generated on a domain with size (grid step ) during a time interval (time step ), where is the spatial index and the temporal one. It is assumed that for each spatial location and each time index, the rain rate (i.e., the amount of water in mm crossing a horizontal cross section of unit area over 1 hour) is the realization of a random variable. The simulation is then conducted in three main steps. First, at time , a centered, reduced, stationary Gaussian random field with spatial correlation structure where and are spatial indexes and the expectancy operator is constructed in the Fourier plane. Second, to account for the temporal evolution that defines , it is assumed that each spatial frequency of the Fourier transform of follows a first-order Markov process in which innovation depends on the spatial frequency. Lastly, the field advection between and is introduced by a phase shift in the frequency domain.
Then, the Gaussian fields are turned into rainfall fields . For that purpose, the rainfall rate distribution is assumed to be zero with probability and log normal with probability , where denotes the local probability of rain.
The different steps of this modelling are detailed in this section. Then, a methodology to derive the correlation parameters from weather radar observations is presented. It is later applied to the radar observations of Bordeaux.
2.1. Generation of Correlated Gaussian Random Field
A convenient framework to simulate a correlated stationary random field that has an arbitrary Probability Distribution Function (PDF) is to define a monotonic transformation that converts a Gaussian field into [19, 21]. If the rain spatiotemporal correlation function is estimated from realizations (radar observations for example) of the rain field , a methodology described in Sections 2.3.3 and 2.3.4 allows to define a correlation such that has the expected correlation .
2.1.1. Generation of the Gaussian Field in Space
A Gaussian random field can be interpreted as a collection of Gaussian random variables correlated the ones with each other. is stationary in space whenever its correlation function between two points and depends only on the distance : The stationarity assumption is essential from a computational point of view. Indeed, it allows reducing the size of the correlation matrix from to to use efficient algorithms such as turning band methods [22] or circular embedded correlation matrices [23].
Another convenient methodology to generate a stationary Gaussian field is to start from its spatial spectrum that is linked to the spatial correlation function [19] by Fourier transform. The stationary Gaussian random field , can then be expressed for each point of the simulation grid by its Fourier expansion: where is the spatial frequency, is the Fourier transform of and is the transposition operator. Coefficients are shown to be equal to [19]: where are independent standard centered complex Gaussian random variables. According to (2) and (3), the definition of the Gaussian field spatial spectrum —or reciprocally the correlation function —is enough to generate a spatially correlated Gaussian field.
2.1.2. Addition of the Temporal Dimension
The rain field temporal evolution can be decomposed into two components. The first one is related to the rain field advection due to air streams. The second one corresponds to the rain cell evolution along their path. In [19], both terms are taken into account separately in the simulation algorithm. The rain field temporal evolution is introduced in the Fourier domain. Considering that each spatial frequency evolves according to a first-order Markov process, the temporal evolution between time and time of the spatial spectral components (3) of the underlying Gaussian field is driven by [19] where is the advection vector and the autocorrelation lag of the Markov process in the Fourier domain. In (4), the term defines the field spectral component at spatial frequency and time with a phase shift that corresponds to the rain advection between and . The term corresponds to the process innovation. The variance of is assumed to be constant in time, preserving that way the spatial correlation of the Gaussian field when performing the Fourier expansion (2). The parameterization proposed in [19] where is the correlation time in the Fourier plane, complies with the definition of a Markov process. Equation (5) induces a non-separable spatiotemporal correlation function as the latter cannot be expressed as the product of a spatial correlation function by a temporal one [24]. The choice and parameterization of (5) is discussed in Section 2.3.4 from the radar data of Bordeaux.
2.2. Conversion into Rain Rate Fields
In order to convert the Gaussian random fields obtained from (2), (3), (4), and (5) into rain fields , the local statistical distribution of rain rate has to be known. The latter is highly dependent on the geographical area. For that purpose, Rec. ITU-R P.837, 2007 [25] gives the probability of rain and the empirical rainfall rate Cumulative Distribution Function (CCDF) all over the world [26, 27]. In other respects, many authors ([28] or [29]) have shown that the rain rate conditional CCDF is log normal (mean , standard deviation ) so that the absolute CCDF can be modelled by: where the probability of rain is given by Rec. ITU-R P.837, 2007 [25] and () are derived by least square fitting with Rec. ITU-R P.837, 2007 [25] using the procedure given in Rec. ITU-R P.1057. The result is shown in Figure 1. A good agreement can be observed between (6) and the absolute CCDF derived from Rec. ITU-R P.837, 2007 [25], at least for the temperate and terrestrial areas considered.
From (6), [19] proposes a transformation that turns the Gaussian fields into rainfall rate fields : where , is the complementary error function and its inverse.
Equation (7) converts the Gaussian random variables into random variables that follow the CCDF defined by (6). However, from the radar data of Bordeaux, it is shown in Section 2.3.2 that (7) introduces a numerical artifact that needs to be corrected.
2.3. Parameterization and Refinement of the Modeling
The spatial and temporal parameters , , and in the model proposed by [19] are the key parameters that need to be accurately assessed in order to obtain a realistic description of the statistical dependence of rain rates in space and in time. In [19], the model parameterization derives from radar data sets collected over tropical oceans during the experimental campaign GATE (GARP Atlantic Tropical Experiment) from the GARP (Global Atmospheric Research Program) [30]. As the spatial and temporal characteristics of the rain fields over oceanic tropical areas are likely to be different from the ones over temperate areas, a parameter retrieval specific to midlatitude areas is undertaken. For that purpose, the spatial correlation function —that is, the spatial spectrum in (3)—and the parameters and that drive the temporal evolution in (4) and (5) are determined from a meteorological radar dataset collected over Bordeaux (France) that is first presented. The latter is also used to illustrate the necessity of a slight modification of (7) to match features observed on the radar data.
2.3.1. Description of the Dataset
The data come from the weather radar located at Bordeaux (44.5°N, −0.5°E), south western France, in the vicinity of the Atlantic Ocean, in a very flat region. The radar is an S-band radar which is part of the French operational radar network managed by Météo France. The polar scans acquired each 5 mn are projected into a Cartesian grid with a pixel size of 1 × 1 km2. The dataset contains 35286 scans acquired from January to December 1996. Images including ground clutter or melting layer echoes were removed from the data set so that the used radar data only refer to rain fields [12–14]. The conversion of the reflectivity into rain rate is made using the standard relation for midlatitudes [31], where is the radar reflectivity factor in mm6 m−3 and the rainfall rate in mm h−1.
2.3.2. Refinement of the Conversion of Gaussian Fields into Rain Fields
Transformation (7) proposed by [19] to convert Gaussian fields into rain fields aims at preserving the rain rate local distribution. Nevertheless, it introduces an unrealistic feature on the generated rain fields. Particularly, (7) introduces on an unrealistic link between the spatial fraction of that is affected by rain and the spatial average of the strictly positive rain rates values . This point is illustrated in Figure 2, where (7) defines as a growing function of . This is in contradiction with the ergodicity assumption of rain rate fields [26, 32, 33] that implies that does not change with respect to the fractional area affected by rain. This point is confirmed in Figure 2 by the radar observations of Bordeaux. The artifact introduced by (7) is purely numerical: whenever the size of the simulation area reduces with respect to the decorrelation distance driven by , the random variables start to deviate from a log normal distribution with parameters and in (6). Therefore, (7) has to be corrected, introducing the correction terms and : The values of the numerical correction terms and must be determined from simulations of : they are adjusted to insure that the random variables that compose the rain field are log normal, with parameters and . This leads to assume that considering a sufficiently large area the spatial CCDF can be approximated by the temporal one as already underlined by [21, 28]. Therefore, their values are specific to the simulation area considered and to the definition of : they must be evaluated once for all whenever the simulation configuration has been defined. That is the reason why their values are not specified in the present paper.
Finally, (8) guarantees the statistical independence between the fraction of the simulation area that is affected by rain and the conditional mean rain rate as illustrated on Figure 2 where is constant whatever in compliance with the radar observations.
2.3.3. Parameterization of the Spatial Correlation
The spatial correlation of the Gaussian field can be determined from the analysis of radar data. As part of a stochastic approach, each radar image is a realization of the random process . In such conditions, (8) implicitly relates the spatial correlation of the underlying Gaussian field to the spatial correlation of . Particularly, let and be random variables representing rainfall rates that are separated by a distance . Assuming that and , where is given by (8) and and are standard centered variable with spatial correlation , the joint Probability Density Function (PDF) is given by From (7), and the rain rate cross expectation can be related to by As the numerical integration of (10) is inaccurate and requires a significant computation time, the rain field covariance is related to the Gaussian field correlation function through a method based on Hermit polynomials [21, 34]: where are the coefficients of the Hermit polynomial expansion of . To reduce inaccuracies due to the finite number of terms in the sum (11) and the slow convergence rate [21], it is preferable to evaluate the covariance of the truncated field , defined from the analytical inversion of (8): The main advantage of this transformation is that the terms of the Hermite expansion of transformation are bounded so that (13) converges rapidly, allowing a fast numerical computation: where are as in (11) but now refer to and transformation .
From a practical point of view, that is, to estimate from radar observations of rain field, is first applied to the radar data. From (12) and (8), the rain fields are thus turned into truncated Gaussian field . Afterwards, the covariance of the truncated Gaussian field is computed for each radar map performing the inverse Fourier transform of the spatial power spectrum of . Finally, is computed inverting numerically (13), from the average covariance function of over all the radar maps. A two-dimensional representation of (not shown) suggests slightly larger correlation values along a north south axis indicating a possible local anisotropy. Nevertheless, this trend is weak, probably due to the flatness of the region of Bordeaux, where orographic forcing is almost absent. Therefore, the correlation is assumed to be isotropic in what follows, so that reduces to . The isotropic correlation function finally obtained from the radar of Bordeaux is presented in Figure 3.
A regression process concludes that can be appropriately approximated by Equation (14) indicates a slowly decaying part, denoting the existence of a possible long range dependence of rain rates. Obviously, (14) is specific to the region of Bordeaux since it is derived from local radar observations. Nevertheless, (14) is very close to the formulation proposed by Barbaliscia et al. in [35] from rain gauge measurement across Italy. Therefore, (14) should be appropriated to model rain fields in midlatitude areas.
2.3.4. Parameterization of the Temporal Evolution
The rain field temporal evolution is decomposed into one term related to the advection and the other related to the rain rate field dynamics along their tracks. Bell in [19] introduces the temporal evolution of the underlying Gaussian field through (4) and (5). Here, the definition of the correlation time in (5) is required. From radar observations of rain fields in tropical areas, Bell [19] proposes the formulation: Equation (15) corresponds to the rain field correlation time and does not refer to the temporal evolution of the underlying Gaussian fields in (4) or (5). Moreover, the effect of the advection was not isolated, biasing the analysis. Lastly, a parameterization for midlatitude areas is required. Consequently, an attempt toward a more rigorous derivation of the temporal evolution of the Gaussian field from (4) is conducted from the midlatitude radar dataset of Bordeaux.
First, the rain fields are corrected for advection. To do it, the average advection vector is computed from two successive radar observations of the rain field, by maximizing their spatial cross-correlation [36]. Proceeding that way, successive rain fields which average motion is null are defined. This advection correction allows getting rid of the shift phase in (4). Now, the temporal evolution of the underlying Gaussian fields reduces to where are the coefficients of the spatial Fourier expansion of —that is, of —and where is given by (5). Recalling that is normal, standard, centered, and uncorrelated, it follows from (16) that where * is the complex conjugate. Considering (5), the correlation time is finally given by: The term is the spatial cross-spectral density of the spatial processes and . It is linked by Fourier transform to the spatiotemporal correlation function of . The latter is inferred from the methodology developed in Section 2.3.3. Particularly, the rain fields corrected for advection are first transformed into truncated Gaussian fields with defined in (12). From fields , the spatiotemporal power spectrum is computed, and, by inverse Fourier transform, the space-time correlation function is obtained. Equation (13) then allows the derivation of and, finally, is obtained by Fourier transform.
To evaluate the value of in (18) from the radar data of Bordeaux mn, successive rainy maps corrected for advection have been used. As the radar spatial coverage is limited, the correction for advection on the successive radar maps has been restricted to areas of 100 × 100 km² while only durations less than 2 hours have been considered. Finally, 25 rainy events of 2 hours corrected for advection have been selected in the database. For each event, the value of has been computed according to the methodology described above. As shown in Figure 4, the average dependency of with respect to the spatial frequency can be approximated by Equation (19) was found assuming that follows a power-law function of the spatial frequency . Equation (19) is thus an adjustment of (13) [15] for midlatitude areas. As shown in Figure 4, the dependency of with respect to suggests, as expected, that the evolution of large features (rain field stratiform component associated to small values of ) is slower (higher correlation time ) than the one of smaller features (rain field convective components with limited spatial extent, associated to highest values).
2.4. Limitations of the Presented Methodology
Several hypotheses limit the validity range of the space-time modeling developed above. The most significant one is the stationarity assumption of the random field. Even if Ferraris et al. [37] have shown that the rain field stationarity is realistic for areas ~200 × 200 km², this assumption is not likely to hold for larger areas where orography or climatology related effects may introduce unstationarity in the rain field. For that reason, the model presented above should be limited to midscale areas (~200 × 200 km²). Therefore, its use for areas of the size of a typical satellite regional coverage (such as Europe or USA) is not appropriate. The same kind of restriction applies for the simulation duration. Indeed, considering (5) and (19), the 0 frequency term of the temporal evolution coefficient is equal to 1. Consequently, (4) indicates that the spatial average of the Gaussian field remains constant, preventing that way the temporal evolution of integrated values during the simulation process such as the average rain amount or the fractional area affected by rain. Particularly, the alternation of clear sky and rainy periods cannot be reproduced by the simulated rain fields, the rain amount being approximately constant during the simulation process. In order to overcome these limitations, a coupling with a meteorological reanalysis database is investigated. It is described in Section 3. It allows benefiting from the statistical robustness of the stochastic model for midscale features while ensuring a realistic repartition of the generated rain amount over large areas (continental scale) and long durations (many years).
3. Extension of the Model Using Inputs from Meteorological Reanalysis Database
To cope with the problem of the temporal evolution of the fractional area affected by rain and hence being able to generate realistic long-term correlated time series of rain rate or rain attenuation, the use of GCM (General Circulation Model) outputs is proposed. For that purpose, the ECMWF ERA-40 reanalysis database [38] is chosen as it offers a worldwide estimation of rainfall amounts for durations long enough to be statistically representative. We underline that the purpose of this work is not to forecast propagation conditions from outputs of numerical weather prediction as proposed in [39] but rather to generate long term propagation scenarios with high spatiotemporal resolutions from realistic past meteorological situations.
3.1. Description of the ECMWF ERA-40 Dataset
The dataset used in this study is a subset from the ECMWF ERA-40 reanalysis database freely available on the ECMWF website. This database strives at reproducing the state of the atmosphere in the past (the considered period ranges from 1957 to 2002), constraining a numerical weather forecast model with observations. The model internal variables are pressure, temperature, wind, and humidity profiles from which several other parameters such as rain amounts are deduced. This database is provided with a coarse spatial resolution of 2.5° × 2.5° (i.e., ~250 × 250 km2), with one sample each 6 hours.
3.2. Rain Amount Parameterization
A relation between the rain amount of one ERA-40 resolution cell with the fraction of this resolution cell that is affected by rain is presented in Section 3.2.1. Then, this information related to the rain fractional area is introduced in the stochastic approach developed in Section 2. Proceeding that way, a realistic (since driven by ERA-40 outputs) spatiotemporal evolution of the simulated rain fields is insured at large-scale (continental scale) and for long durations (several years).
3.2.1. Link between the Rain Amount from a GCM and the Fractional Area Affected by Rain
Many studies on rainfall remote sensing have put forward that if a sufficiently large number of rain cells at different stages of their life cycle are observed over an area (~200 × 200 km²), a tight link exists between the fractional area that is affected by rain above a given threshold and the (unconditional) area-averaged rain rate [40, 41]. As an example, Figure 5 shows the linear dependency of with respect to the fractional area deduced from the radar data of Bordeaux.
This relation is all the more tight because a large number of rain cells are observed simultaneously in the rain field, giving an exhaustive representation (in a statistical sense) of the climatological rain rate distribution [42].
Besides, Eltahir and Bras [32] have shown that considering rainfall outputs from low-resolution GCM, the same kind of relation holds between the average rain accumulation (equivalent height of water) which corresponds to a spatially averaged rainfall rate integrated during a duration and the fraction of the GCM pixel that is under rain : where is the conditional mean rain rate (i.e., knowing that it is raining) on the GCM considered area. From (6), the conditional mean rain rate can be obtained considering the log normal local distribution of rainfall rates [28, 32]: Therefore, from the rain amount time series given by ERA-40, the associated time series of the fractional area can be inferred worldwide from (20) and (21), every hours, over cells with size 2.5° × 2.5°.
3.2.2. Parameterization of the Mean Value of the Gaussian Field
In order to include the information derived from the reanalysis database in the high-resolution stochastic rain model, a link has to be found between the rain fractional area and the simulation parameters. Recently, Authors in [43] have proposed a model to link the fraction of one Gaussian field realization that is above a given threshold with the Gaussian field average value . Under some assumptions on the shape of the correlation function, the authors have shown that the spatial distribution of the samples of one Gaussian field realization simulated on a finite grid with size is normal, with mean and standard deviation , where is the spatial average of the Gaussian field correlation function over the simulation grid. In such conditions, Authors in [43] have shown that the fraction of a simulated Gaussian field that is over a preset threshold is given by From the Gaussian modeling of rain fields developed in Section 2, is linked to the rain no-rain threshold through the equation , where is the probability of rain. Obviously so that the mean value of the Gaussian field can be related to the rain fractional area by By definition of the Fourier transform, the 0 frequency term of the Fourier expansion of the Gaussian field is equal to its mean value . Therefore, whenever is set to , the fraction of the simulation area that undergoes rain is , as expected.
In the foregoing, the spatiotemporal evolution of the simulated rain fields is derived from ERA-40 rain amount time series which temporal resolution is 6 hours. It is thus mandatory to find an interpolation scheme to get values of the rain fractional area at a finer time resolution. This high-resolution time series of rain fractional area denoted by and sampled at must be continuous and, obviously, its average value over hours must be equal to the low-resolution rain fractional area derived each 6 hours from ERA-40 and denoted by hereinafter. Therefore, considering that holds in the temporal interval , is arbitrarily defined as a piecewise linear function such as To ensure the continuity at the edges of the interval, we set combined with the constraint on the average value on : Equations (24), (25), and (26) define a set of three equations with three unknown coefficients , , and . The effect of this interpolation procedure is illustrated in Figure 6, where it is shown to reproduce accurately the rain fractional area derived from ERA-40. Nevertheless, the effect on the simulations of this piecewise linear interpolation is not well assessed and constitute. It will be shown further that it has not a significant impact on the realism of the simulations in terms of temporal correlation.
Figure 7 sums up the step-by-step methodology developed to parameterize the stochastic model from ERA-40 rain amount time series. Particularly, from the ERA-40 low-resolution () time series of rain amount, a low-resolution () time series of rain fractional area is deduced from (20). The latter is interpolated at (i.e., the radar observation frequency) from the piecewise approach defined above. At this stage, a time series of Gaussian field mean values is derived from (23). It is injected in the spatiotemporal model through the Fourier coefficient . Finally, the model generates realistic rain fields with long-term evolution at the scale of one ERA-40 resolution cell (2.5° × 2.5°, i.e., ~250 × 250 km2) with a spatial resolution of 1 × 1 km² and with a temporal resolution of 6 mn.
(a)
(b)
(c)
3.2.3. Extension to Several Cells of the GCM
The extension of the simulation domain to areas larger than 2.5° × 2.5° requires taking into account more than one ERA-40 cell. It is not possible to juxtapose directly the 2.5° × 2.5° rain fields obtained from the methodology described above for one ERA-40 resolution cell because the rain field continuity between two adjoining cells is not ensured. For that purpose, an interpolation procedure that preserves the statistical features of each 2.5° × 2.5° subfield while ensuring the continuity of the global rain field has been developed. Particularly, for each ERA-40 cell contained in the simulation area, a Gaussian field with size is simulated according to the methodology described in Section 2. Then, for each point of the simulation area, the value of the global Gaussian field results from the weighted sum of the Gaussian fields , , , generated for each of the ERA-40 adjacent resolution cells:where is the distance between two ERA-40 resolution cells (). Figure 8 gives an overview of the interpolation scheme.
Common interpolation methodologies such as bilinear or cubic spline interpolation are not suitable in our case because they do not preserve the Gaussian field statistical features in terms of variance and correlation. On the contrary, (27) does not change the variance and the correlation of the Gaussian subfields . Moreover, by construction, the model reproduces both the rain amount given by ERA-40 for each resolution cell (as a laplacian filtering is applied to the coefficient of each subfield) and the local rain rate CCDF given as input parameters. Figure 9 gives an example of rain rate field simulated at large-scale, over Europe, on 24/02/1999 at 19:30 UTC from the concurrent ERA-40 rain amount data.
3.3. Advection
Whenever the rain field spatiotemporal simulation lasts a few hours, a common and practical approach is to consider the advection vector as constant both in speed and direction. However, this frozen storm hypothesis is not satisfying if long-durations or large-scale simulations are considered. Indeed, the wind field and the resulting rain field motion may evolve considerably with distance or time. In such conditions, wind data from the ERA-40 database can provide realistic inputs for this parameter. Indeed, comparing radar-derived advection with wind outputs from meteorological model [36, 44], have shown that the rain field motion can be accurately derived from the ERA-40 wind data at the 700 hPa pressure level. Consequently, for each ERA-40 resolution cell making part of the large-scale simulation area, the advection vector defined in (4) is driven by the ERA-40 concurrent wind data at a pressure level of 700 hPa. Lastly, in order to get a smooth temporal evolution of the rain advection, both components of the winds are interpolated by cubic splines with a time step of 6 mn.
4. Conversion into Attenuation
4.1. Conversion of Rain Rate Fields into Attenuation Fields
Rain field simulations are converted into attenuation by integration of the specific attenuation along the link path. The attenuation (dB) endured by the electromagnetic wave during its propagation through rain is given by where is the rain rate intensity (mm h−1) at and is the slant path (km) through rain. and are coefficients, function of the elevation angle, the frequency, and the polarization of the electromagnetic wave. Their values are given by Rec. ITU-R P.838 [45]. The definition of the geometric (latitude, longitude, altitude of the satellite) and radiowave (frequency, polarization) characteristics of the telecommunication link allows the conversion of the rainfall field into an attenuation field. Here, the altitude of the points composing the simulation grid is derived from Rec. ITU-R P.1511 [46]. The rain height is 0.36 km above the 0°C isotherm height given by the ERA-40 temperature profiles. The use of ERA-40 data, sampled each 6 hours, to define the 0°C isotherm height allows a better description of the rain height with respect to Rec. ITU-R P.839 [47] that provides only the yearly average of the rain height as shown in Figure 10.
4.2. Over-Sampling of the Attenuation Time Series
For SatCom systems using adaptive coding, adaptive modulation, or other FMTs, it is crucial to get an estimate of the propagation channel evolution with a temporal resolution of about 1 second. However, the temporal resolution of the model developed in Section 3 is 6 mn. Therefore, it is mandatory to simulate the fast dynamic fluctuations of the propagation channel to increase the model temporal resolution. For that purpose, the on-demand time series synthesizer described in [6] is used. The latter makes a 1 s stochastic interpolation of the 6 mn attenuation time series extracted from the space-time model, generating that way high-frequency (1 Hz) attenuation time series as illustrated in Figure 11.
The high-resolution interpolation [6] holds whenever the attenuation time series are Markov and follow a log normal distribution which parameters can be derived from ITU-R P.618 [48]. The parameter describing the dynamic fluctuations in [6] is taken equal to the advocated value of 2.10−4 s−1. Coupled with the spatiotemporal model, the stochastic interpolation enables to simulate 1 Hz spatially correlated attenuation time series for thousands of links disposed across the satellite coverage, for long durations (several years).
5. Preliminary Validations, Limitations
5.1. First-Order Statistics
Preliminary analysis shows that the rain amount from the reanalysis database is reproduced by the simulations for each 6 h period with an RMS error of about 15% as illustrated in Figure 12.
(a)
(b)
By construction, the model reproduces the rain rate log normal approximation (6) of Rec. ITU-R P.837 [25] on each point of the simulation grid. In other respects, long-term attenuation time series reproduce the attenuation CCDFs given by ITU-R P.618 [48], as illustrated in Figure 13.
(a)
(b)
(c)
(d)
In addition, from OLYMPUS satellite measurements at 20 GHz led in 1992 near Paris (France), the model ability to reproduce attenuation statistics derived from beacon measurements has been investigated. Here, high-resolution (1 Hz) attenuation time series have been synthesised from ERA-40 data concurrent to OLYMPUS measurement period, assuming a link with the OLYMPUS geostationary satellite (19°W) at 20 GHz [49]. Figure 14 shows the results in terms of attenuation CCDF derived from simulations and from OLYMPUS experiment.
According to Figure 14, statistics derived from spatiotemporal simulations of attenuation fields compare satisfactorily with the experimental ones up to small time percentages. Nevertheless, it is noticeable that the variability of the experimental statistics is higher than the one derived from the simulations.
5.2. Second-Order Statistics
In order to assess the model ability to reproduce the spatial variability of the attenuation due to rain, the diversity gain has been computed first from the radar data of Bordeaux collected in 1996 every 5 mn [36]. As the latter refer to rain fields, the radar data have been converted to attenuation fields assuming a 30 GHz link with OLYMPUS. On the other hand, the same exercise has been conducted from the rain fields simulated by the spatiotemporal model each 6 mn using ERA-40 data of 1996. The results are presented in Figure 15, for single site attenuation values ranging from 4 dB to 32 dB. Figure 15 shows that the spatiotemporal model reproduces satisfactorily the diversity gain derived from radar data, even if a slight trend toward underestimation for the highest single site attenuation values and for low distances can be observed.
A similar test has been carried out with the 1 Hz over-sampled time series. Indeed, the joint attenuation CCDF has been computed first from OLYMPUS [49] measurements at La-Folie-Bessin and Gometz-La-Ville that are 7 km apart and, second, from the 1 Hz time series derived at both locations from the spatiotemporal model. According to Figure 16, the experimental and model-based distributions match satisfactorily, showing thus the model ability to account for the rain attenuation spatial variability.
To evaluate the relevance of the model temporal parameterization, the temporal correlation function of attenuation time series has been evaluated first from yearly attenuation time series extracted from the radar observations of Bordeaux in 1996 and, second, from the 6 mn attenuation time series simulated over Bordeaux in 1996.
The results, reported in Figure 17, show that the temporal correlation function computed from simulations is close to the one derived from radar. Nevertheless, the temporal validation of the model has to be deepened, notably to assess the ability of the simulated time series to reproduce fade slopes, fade duration, or interfade durations statistics derived from beacon measurements. The autocorrelation from function from time series simulated at one should also be compared with autocorrelation from measurements made during ITALSAT or OLYMPUS campaigns.
6. Conclusion
A model to generate spatially and temporally correlated rain fields or rain attenuation fields for propagation studies has been presented. It lies on a nonlinear transformation of random Gaussian fields constructed in the Fourier plane.
After a description of the mathematical framework, a method to derive the spatiotemporal correlation parameters from radar data has been proposed. It has been applied to the meteorological radar data of Bordeaux (France, midlatitude area). As the initial approach proposed by Bell in [19] was not intended to model rain rate fields at large-scale and for long durations, a theoretical framework to enlarge the validity domain of the model has been developed. A methodology to constrain the model with rain outputs from the ERA-40 reanalysis database has then been presented. Particularly, the rain amount generated each 6 h by the spatiotemporal model is constrained to be the one given by ERA-40. The rain advection is also modelled using ERA-40 wind data. Coupled with a large-scale interpolation scheme, the spatiotemporal model is thus able to generate realistic rain fields for long durations (several years) and large areas (i.e., over Europe or USA, the size of a typical SatCom system coverage) with a spatial resolution of 1 × 1 km² and a temporal resolution of 6 mn.
Then, defining the geometric and radio electric characteristics of a link, the rain rate fields are turned into attenuation fields considering the wave path through rain. Here, the stochastic interpolation scheme described in [6] allows reducing the temporal resolution up to 1 s. The model then shows its ability to reproduce several statistical characteristics observed on different and independent sources of data. In addition to the first-order statistics of rain and rain attenuation, the model is able to reproduce the spatial repartition of rain attenuation in terms of diversity gain or joint attenuation CDF. From the temporal point of view, first validations in terms of temporal correlation function are promising but further enquiries are needed to consolidate the results.
Reliable inputs for FMTs and RRM (Radio Resource Management) design and optimization may be obtained from the presented model, for several configurations, since the simulation area can reach the size of a typical satellite coverage with a temporal resolution that can be reduced down to 1 s for several years. Particularly, considering complex network simulations that involve thousands of links dispersed over the whole coverage, the spatiotemporal model is able to provide realistic propagation conditions, both correlated in space and in time.
However, the use of the ERA-40 reanalysis database may lead to some bias in the simulations, especially in coastal or mountainous areas, that is, in places where the quality of the reanalysis database is questionable. Nevertheless, the use of reanalysis data turns out to be promising as it allows simulating test cases from specific past weather conditions. For instance, this can provide a good opportunity to study diurnal cycles, seasonal cycles, or interannual variability but also worst cases identified on the reanalysis data.
A logical following to this work is the inclusion, using the same approach, of the gaseous and liquid water attenuation as they become critical for SatCom system operating at Q or V band or even at Ka band for other applications such as radar altimetry.
Acknowledgments
The authors are very grateful to Météo-France for providing the radar data from Bordeaux. The authors also acknowledge the ECMWF for providing ERA-40 data.