# A Large-Scale Space-Time Stochastic Simulation Tool

Authors: Nicolas Jeannin; Laurent Féral; Henri Sauvageot; Laurent Castanet; Frédéric Lacoste
Journal: International Journal of Antennas and Propagation
 (2012)

---

## 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.


---

## Body

## 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 fieldsR(x,t) are generated on a domain with size N×N (grid step Δx) during a time interval T (time step Δt), where x=(x1,x2) is the spatial index and t the temporal one. It is assumed that for each spatial location and each time index, the rain rate R (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 t, a centered, reduced, stationary Gaussian random field G(x,t) with spatial correlation structure cG(x,y)=〈G(x)G(y)〉 where x and y are spatial indexes and 〈·〉 the expectancy operator is constructed in the Fourier plane. Second, to account for the temporal evolution that defines G(x,t+Δt), it is assumed that each spatial frequency k=(k1,k2) of the Fourier transform a(k,t) of G(x,t) follows a first-order Markov process in which innovation depends on the spatial frequency. Lastly, the field advection between G(x,t) and G(x,t+Δt) is introduced by a phase shift in the frequency domain.Then, the Gaussian fieldsG(x,t) are turned into rainfall fields R(x,t). For that purpose, the rainfall rate distribution is assumed to be zero with probability 1-P0 and log normal with probability P0, where P0 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 fieldR(x,t) that has an arbitrary Probability Distribution Function (PDF) is to define a monotonic transformation ψ that converts a Gaussian field G(x,t) into R(x,t)=ψ[G(x,t)] [19, 21]. If the rain spatiotemporal correlation function cR(x,y,t) is estimated from realizations (radar observations for example) of the rain field R(x,t), a methodology described in Sections 2.3.3 and 2.3.4 allows to define a correlation cG(x,y,t) such that R(x,t)=ψ[G(x,t)] has the expected correlation cR(x,y,t).


#### 2.1.1. Generation of the Gaussian Field in Space

A Gaussian random fieldG(x) can be interpreted as a collection of Gaussian random variables correlated the ones with each other. G(x) is stationary in space whenever its correlation function cG(x,y) between two points x=(x1,x2) and y=(y1,y2) depends only on the distance d=(x1-y1)2,(x2-y2)2:
(1)cG(x,y)=cG(d).
The stationarity assumption is essential from a computational point of view. Indeed, it allows reducing the size of the correlation matrix from N2×N2 to N×N to use efficient algorithms such as turning band methods [22] or circular embedded correlation matrices [23].Another convenient methodology to generate a stationary Gaussian fieldG(x) is to start from its spatial spectrum s(k)=sk that is linked to the spatial correlation function cG(d) [19] by Fourier transform. The N×N stationary Gaussian random field G(x), can then be expressed for each point x=(x1,x2) of the simulation grid by its Fourier expansion:
(2)G(x)=∑k1=0N-1∑k2=0N-1akexp⁡⁡(i2πkTxN),
where k=(k1,k2) is the spatial frequency, ak=a(k) is the Fourier transform of G(x) and T is the transposition operator. Coefficients ak are shown to be equal to [19]:
(3)ak=skek,
where ek are independent standard centered complex Gaussian random variables. According to (2) and (3), the definition of the Gaussian field spatial spectrum sk—or reciprocally the correlation function cG(d)—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 t and time t+Δt of the spatial spectral components ak (3) of the underlying Gaussian field is driven by [19]
(4)ak,t+Δt=ak,texp⁡⁡(-ikTVΔt)βk+(1-βk2)skek,t+Δt,
where V is the advection vector and βk the autocorrelation lag of the Markov process in the Fourier domain. In (4), the term ak,texp(-ikTVΔt) defines the field spectral component at spatial frequency k and time t with a phase shift that corresponds to the rain advection between t and t+Δt. The term skek,t+Δt corresponds to the process innovation. The variance of ak,t 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]
(5)βk=exp⁡⁡(-Δtτk),
where τk 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 fieldsG(x,t) obtained from (2), (3), (4), and (5) into rain fields R(x,t), the local statistical distribution of rain rate R 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 P0 and the empirical rainfall rate Cumulative Distribution Function (CCDF) P(R≥R*) all over the world [26, 27]. In other respects, many authors ([28] or [29]) have shown that the rain rate conditional CCDF P(R≥R*/R>0) is log normal (mean μR, standard deviation σR) so that the absolute CCDF P(R≥R*) can be modelled by:
(6)P(R≥R*)=P0∫R*+∞12πσR2rexp⁡⁡[-(ln⁡⁡r-μR)22σR2]drifR*>0,P(R≥0)=1,
where the probability of rain P0 is given by Rec. ITU-R P.837, 2007 [25] and (μR,σR) 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.Figure 1
Rain rate absolute CCDFP(R≥R*) derived from ITU-R Rec P.837-5 [25] (crosses) and log normal regression (6) (solid line).From (6), [19] proposes a transformation ψ that turns the Gaussian fields G(x,t) into rainfall rate fields R(x,t):
(7)R(x)=ψ[G(x)]=0ifG(x)<α,R(x)=ψ[G(x)]=exp⁡⁡[μR+2σRerfc-1(erfc⁡[G(x)/2]P0)]ifG(x)≥α,
where α=2σRerfc-1(2P0), erfc is the complementary error function and erfc-1 its inverse.Equation (7) converts the Gaussian random variables G(x,t) into random variables R(x,t) 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 parameterscG(d), βk, and τk 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 cG(d)—that is, the spatial spectrum sk in (3)—and the parameters βk and τk 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 Z-R relation for midlatitudes Z=300×R1.5 [31], where Z is the radar reflectivity factor in mm6 m−3 and R 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 G(x,t) into rain fields R(x,t) aims at preserving the rain rate local distribution. Nevertheless, it introduces an unrealistic feature on the generated rain fields. Particularly, (7) introduces on R(x,t) an unrealistic link between the spatial fraction fR of R(x,t) that is affected by rain and the spatial average of the strictly positive rain rates values 〈R∣R>0〉. This point is illustrated in Figure 2, where (7) defines 〈R∣R>0〉  as a growing function of fR. This is in contradiction with the ergodicity assumption of rain rate fields [26, 32, 33] that implies that 〈R∣R>0〉 does not change with respect to the fractional area fR 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 N of the simulation area reduces with respect to the decorrelation distance driven by cG(d), the random variables R(x,t) start to deviate from a log normal distribution with parameters μR and σR in (6). Therefore, (7) has to be corrected, introducing the correction terms μoff and σoff:
(8)R(x)=ψ[G(x)]=0ifG(x)<α,R(x)=ψ[G(x)]=exp⁡⁡[μR-μoff⁡σoff⁡+2σRσoff⁡erfc-1(erfc⁡(G(x)/2)P0)]ifG(x)≥α,
The values of the numerical correction terms μoff and σoff must be determined from simulations of R(x,t): they are adjusted to insure that the random variables R(x,t) that compose the N×N rain field are log normal, with parameters μR and σR. 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 cG(d): 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.Figure 2
Relation between the fraction of the areafR that is affected by rain and the conditional mean rain rate 〈R∣R>0〉 computed from 35 000 radar scans and from 35 000 simulations using transformation ψ defined by (7) [19] and the correction (8).Finally, (8) guarantees the statistical independence between the fraction of the simulation area fR that is affected by rain and the conditional mean rain rate 〈R∣R>0〉 as illustrated on Figure 2 where 〈R∣R>0〉 is constant whatever fR 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 processR(x,t). In such conditions, (8) implicitly relates the spatial correlation cG(d) of the underlying Gaussian field G(x,t) to the spatial correlation cR(d) of R(x,t). Particularly, let r1 and r2 be random variables representing rainfall rates that are separated by a distance d. Assuming that r1=ψ(g1) and r2=ψ(g2), where ψ is given by (8) and g1 and g2 are standard centered variable with spatial correlation cG(d), the joint Probability Density Function (PDF) p(g1,g2) is given by
(9)p(g1,g2)=12π1-cG2(d)exp⁡⁡[-12g12+g22-2cG(d)g1g21-cG2(d)].
From (7), α=2erfc-1(2P0) and the rain rate cross expectation 〈r1r2〉 can be related to cG(d) by
(10)〈r1r2〉=〈ψ(g1)ψ(g2)〉=∬-∞+∞ψ(g1)ψ(g2)p(g1,g2)dg1dg2=∬α+∞ψ(g1)ψ(g2)12π1-cG2(d)×exp⁡⁡[-12g12+g22-2cG(d)g1g21-cG2(d)]dg1dg2.
As the numerical integration of (10) is inaccurate and requires a significant computation time, the rain field covariance 〈r1r2〉-〈r1〉〈r2〉 is related to the Gaussian field correlation function cG(d) through a method based on Hermit polynomials [21, 34]:
(11)〈r1r2〉-〈r1〉〈r2〉=∑k=1+∞ψk2k!cG(d)k,
where ψk 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 G̃(x,t)=ψ̃-1[R(x,t)]=Ω[G(x,t)], defined from the analytical inversion ψ-1 of (8):
(12)G̃(x,t)=ψ-1[R(x,t)]=G(x,t)ifR(x,t)>0,G̃(x,t)=αifR(x,t)=0.
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:
(13)〈g̃1g̃2〉-〈g̃1〉〈g̃2〉=∑k=1+∞Ωk2k!cG(d)k,
where g̃1,g̃2,Ωk are as in (11) but now refer to G̃(x,t) and transformation Ω.From a practical point of view, that is, to estimatecG(d) from radar observations of rain field, ψ̃-1 is first applied to the radar data. From (12) and (8), the rain fields R(x,t) are thus turned into truncated Gaussian field G̃(x,t). Afterwards, the covariance of the truncated Gaussian field G̃(x,t) is computed for each radar map performing the inverse Fourier transform of the spatial power spectrum of G̃(x,t). Finally, cG(d) is computed inverting numerically (13), from the average covariance function of G̃(x,t) over all the radar maps. A two-dimensional representation of cG(d) (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 cG(d) reduces to cG(d). The isotropic correlation function cG(d) finally obtained from the radar of Bordeaux is presented in Figure 3.Figure 3
Isotropic correlation functioncG(d) of the Gaussian field G(x,t) derived from the radar observations of Bordeaux and analytical regression (14).A regression process concludes thatcG(d) can be appropriately approximated by
(14)cG(d)=0.59exp⁡⁡(-d31)+0.41exp⁡⁡(-d800).
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 G(x,t) through (4) and (5). Here, the definition of the correlation time τk in (5) is required. From radar observations of rain fields in tropical areas, Bell [19] proposes the formulation:
(15)τk=0.24×(π|k|)2/3.
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 G(x,t) 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 R̂(x,t)=ψ⌊Ĝ(x,t)⌋ which average motion is null are defined. This advection correction allows getting rid of the shift phase exp(-ikTVΔt) in (4). Now, the temporal evolution of the underlying Gaussian fields Ĝ(x,t) reduces to
(16)ak,t+Δt=ak,tβk+(1-βk2)skek,t+Δt,
where ak,t are the coefficients of the spatial Fourier expansion of Ĝ(x,t)—that is, of G(x,t)—and where βk is given by (5). Recalling that ek,t is normal, standard, centered, and uncorrelated, it follows from (16) that
(17)〈ak,t+Δtak,t*〉=〈ak,tak,t*〉βk,
where * is the complex conjugate. Considering (5), the correlation time τk is finally given by:
(18)τk=-Δtln⁡(βk)=-Δtln⁡(〈ak,t+Δtak,t*〉)-ln⁡(〈ak,tak,t*〉).
The term 〈ak,t+Δtak,t*〉 is the spatial cross-spectral density of the spatial processes Ĝ(x,t) and Ĝ(x,t+Δt). It is linked by Fourier transform to the spatiotemporal correlation function cĜ(d,Δt) of Ĝ(x,t). The latter is inferred from the methodology developed in Section 2.3.3. Particularly, the rain fields R̂(x,t) corrected for advection are first transformed into truncated Gaussian fields ĜT(x,t)=ψ̃-1[R̂(x,t)] with ψ̃-1 defined in (12). From fields ĜT(x,t), the spatiotemporal power spectrum is computed, and, by inverse Fourier transform, the space-time correlation function cĜT(d,Δt) is obtained. Equation (13) then allows the derivation of cĜ(d,Δt) and, finally, 〈ak,t+Δtak,t*〉 is obtained by Fourier transform.To evaluate the value ofτk in (18) from the radar data of Bordeaux (Δt=5 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 R̂(x,t) of 2 hours corrected for advection have been selected in the database. For each event, the value of τk has been computed according to the methodology described above. As shown in Figure 4, the average dependency of τk with respect to the spatial frequency k can be approximated by
(19)τk=0.06(1|k|)0.94.
Equation (19) was found assuming that τk follows a power-law function of the spatial frequency |k|. Equation (19) is thus an adjustment of (13) [15] for midlatitude areas. As shown in Figure 4, the dependency of τk with respect to |k| suggests, as expected, that the evolution of large features (rain field stratiform component associated to small values of |k|) is slower (higher correlation time τk) than the one of smaller features (rain field convective components with limited spatial extent, associated to highest |k| values).Figure 4
Evolution of parameterτk as a function of the spatial frequency |k|.


### 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 β0 is equal to 1. Consequently, (4) indicates that the spatial average of the Gaussian field MG=a0 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).


## 2.1. Generation of Correlated Gaussian Random Field

A convenient framework to simulate a correlated stationary random fieldR(x,t) that has an arbitrary Probability Distribution Function (PDF) is to define a monotonic transformation ψ that converts a Gaussian field G(x,t) into R(x,t)=ψ[G(x,t)] [19, 21]. If the rain spatiotemporal correlation function cR(x,y,t) is estimated from realizations (radar observations for example) of the rain field R(x,t), a methodology described in Sections 2.3.3 and 2.3.4 allows to define a correlation cG(x,y,t) such that R(x,t)=ψ[G(x,t)] has the expected correlation cR(x,y,t).


### 2.1.1. Generation of the Gaussian Field in Space

A Gaussian random fieldG(x) can be interpreted as a collection of Gaussian random variables correlated the ones with each other. G(x) is stationary in space whenever its correlation function cG(x,y) between two points x=(x1,x2) and y=(y1,y2) depends only on the distance d=(x1-y1)2,(x2-y2)2:
(1)cG(x,y)=cG(d).
The stationarity assumption is essential from a computational point of view. Indeed, it allows reducing the size of the correlation matrix from N2×N2 to N×N to use efficient algorithms such as turning band methods [22] or circular embedded correlation matrices [23].Another convenient methodology to generate a stationary Gaussian fieldG(x) is to start from its spatial spectrum s(k)=sk that is linked to the spatial correlation function cG(d) [19] by Fourier transform. The N×N stationary Gaussian random field G(x), can then be expressed for each point x=(x1,x2) of the simulation grid by its Fourier expansion:
(2)G(x)=∑k1=0N-1∑k2=0N-1akexp⁡⁡(i2πkTxN),
where k=(k1,k2) is the spatial frequency, ak=a(k) is the Fourier transform of G(x) and T is the transposition operator. Coefficients ak are shown to be equal to [19]:
(3)ak=skek,
where ek are independent standard centered complex Gaussian random variables. According to (2) and (3), the definition of the Gaussian field spatial spectrum sk—or reciprocally the correlation function cG(d)—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 t and time t+Δt of the spatial spectral components ak (3) of the underlying Gaussian field is driven by [19]
(4)ak,t+Δt=ak,texp⁡⁡(-ikTVΔt)βk+(1-βk2)skek,t+Δt,
where V is the advection vector and βk the autocorrelation lag of the Markov process in the Fourier domain. In (4), the term ak,texp(-ikTVΔt) defines the field spectral component at spatial frequency k and time t with a phase shift that corresponds to the rain advection between t and t+Δt. The term skek,t+Δt corresponds to the process innovation. The variance of ak,t 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]
(5)βk=exp⁡⁡(-Δtτk),
where τk 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.1.1. Generation of the Gaussian Field in Space

A Gaussian random fieldG(x) can be interpreted as a collection of Gaussian random variables correlated the ones with each other. G(x) is stationary in space whenever its correlation function cG(x,y) between two points x=(x1,x2) and y=(y1,y2) depends only on the distance d=(x1-y1)2,(x2-y2)2:
(1)cG(x,y)=cG(d).
The stationarity assumption is essential from a computational point of view. Indeed, it allows reducing the size of the correlation matrix from N2×N2 to N×N to use efficient algorithms such as turning band methods [22] or circular embedded correlation matrices [23].Another convenient methodology to generate a stationary Gaussian fieldG(x) is to start from its spatial spectrum s(k)=sk that is linked to the spatial correlation function cG(d) [19] by Fourier transform. The N×N stationary Gaussian random field G(x), can then be expressed for each point x=(x1,x2) of the simulation grid by its Fourier expansion:
(2)G(x)=∑k1=0N-1∑k2=0N-1akexp⁡⁡(i2πkTxN),
where k=(k1,k2) is the spatial frequency, ak=a(k) is the Fourier transform of G(x) and T is the transposition operator. Coefficients ak are shown to be equal to [19]:
(3)ak=skek,
where ek are independent standard centered complex Gaussian random variables. According to (2) and (3), the definition of the Gaussian field spatial spectrum sk—or reciprocally the correlation function cG(d)—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 t and time t+Δt of the spatial spectral components ak (3) of the underlying Gaussian field is driven by [19]
(4)ak,t+Δt=ak,texp⁡⁡(-ikTVΔt)βk+(1-βk2)skek,t+Δt,
where V is the advection vector and βk the autocorrelation lag of the Markov process in the Fourier domain. In (4), the term ak,texp(-ikTVΔt) defines the field spectral component at spatial frequency k and time t with a phase shift that corresponds to the rain advection between t and t+Δt. The term skek,t+Δt corresponds to the process innovation. The variance of ak,t 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]
(5)βk=exp⁡⁡(-Δtτk),
where τk 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 fieldsG(x,t) obtained from (2), (3), (4), and (5) into rain fields R(x,t), the local statistical distribution of rain rate R 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 P0 and the empirical rainfall rate Cumulative Distribution Function (CCDF) P(R≥R*) all over the world [26, 27]. In other respects, many authors ([28] or [29]) have shown that the rain rate conditional CCDF P(R≥R*/R>0) is log normal (mean μR, standard deviation σR) so that the absolute CCDF P(R≥R*) can be modelled by:
(6)P(R≥R*)=P0∫R*+∞12πσR2rexp⁡⁡[-(ln⁡⁡r-μR)22σR2]drifR*>0,P(R≥0)=1,
where the probability of rain P0 is given by Rec. ITU-R P.837, 2007 [25] and (μR,σR) 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.Figure 1
Rain rate absolute CCDFP(R≥R*) derived from ITU-R Rec P.837-5 [25] (crosses) and log normal regression (6) (solid line).From (6), [19] proposes a transformation ψ that turns the Gaussian fields G(x,t) into rainfall rate fields R(x,t):
(7)R(x)=ψ[G(x)]=0ifG(x)<α,R(x)=ψ[G(x)]=exp⁡⁡[μR+2σRerfc-1(erfc⁡[G(x)/2]P0)]ifG(x)≥α,
where α=2σRerfc-1(2P0), erfc is the complementary error function and erfc-1 its inverse.Equation (7) converts the Gaussian random variables G(x,t) into random variables R(x,t) 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 parameterscG(d), βk, and τk 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 cG(d)—that is, the spatial spectrum sk in (3)—and the parameters βk and τk 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 Z-R relation for midlatitudes Z=300×R1.5 [31], where Z is the radar reflectivity factor in mm6 m−3 and R 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 G(x,t) into rain fields R(x,t) aims at preserving the rain rate local distribution. Nevertheless, it introduces an unrealistic feature on the generated rain fields. Particularly, (7) introduces on R(x,t) an unrealistic link between the spatial fraction fR of R(x,t) that is affected by rain and the spatial average of the strictly positive rain rates values 〈R∣R>0〉. This point is illustrated in Figure 2, where (7) defines 〈R∣R>0〉  as a growing function of fR. This is in contradiction with the ergodicity assumption of rain rate fields [26, 32, 33] that implies that 〈R∣R>0〉 does not change with respect to the fractional area fR 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 N of the simulation area reduces with respect to the decorrelation distance driven by cG(d), the random variables R(x,t) start to deviate from a log normal distribution with parameters μR and σR in (6). Therefore, (7) has to be corrected, introducing the correction terms μoff and σoff:
(8)R(x)=ψ[G(x)]=0ifG(x)<α,R(x)=ψ[G(x)]=exp⁡⁡[μR-μoff⁡σoff⁡+2σRσoff⁡erfc-1(erfc⁡(G(x)/2)P0)]ifG(x)≥α,
The values of the numerical correction terms μoff and σoff must be determined from simulations of R(x,t): they are adjusted to insure that the random variables R(x,t) that compose the N×N rain field are log normal, with parameters μR and σR. 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 cG(d): 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.Figure 2
Relation between the fraction of the areafR that is affected by rain and the conditional mean rain rate 〈R∣R>0〉 computed from 35 000 radar scans and from 35 000 simulations using transformation ψ defined by (7) [19] and the correction (8).Finally, (8) guarantees the statistical independence between the fraction of the simulation area fR that is affected by rain and the conditional mean rain rate 〈R∣R>0〉 as illustrated on Figure 2 where 〈R∣R>0〉 is constant whatever fR 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 processR(x,t). In such conditions, (8) implicitly relates the spatial correlation cG(d) of the underlying Gaussian field G(x,t) to the spatial correlation cR(d) of R(x,t). Particularly, let r1 and r2 be random variables representing rainfall rates that are separated by a distance d. Assuming that r1=ψ(g1) and r2=ψ(g2), where ψ is given by (8) and g1 and g2 are standard centered variable with spatial correlation cG(d), the joint Probability Density Function (PDF) p(g1,g2) is given by
(9)p(g1,g2)=12π1-cG2(d)exp⁡⁡[-12g12+g22-2cG(d)g1g21-cG2(d)].
From (7), α=2erfc-1(2P0) and the rain rate cross expectation 〈r1r2〉 can be related to cG(d) by
(10)〈r1r2〉=〈ψ(g1)ψ(g2)〉=∬-∞+∞ψ(g1)ψ(g2)p(g1,g2)dg1dg2=∬α+∞ψ(g1)ψ(g2)12π1-cG2(d)×exp⁡⁡[-12g12+g22-2cG(d)g1g21-cG2(d)]dg1dg2.
As the numerical integration of (10) is inaccurate and requires a significant computation time, the rain field covariance 〈r1r2〉-〈r1〉〈r2〉 is related to the Gaussian field correlation function cG(d) through a method based on Hermit polynomials [21, 34]:
(11)〈r1r2〉-〈r1〉〈r2〉=∑k=1+∞ψk2k!cG(d)k,
where ψk 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 G̃(x,t)=ψ̃-1[R(x,t)]=Ω[G(x,t)], defined from the analytical inversion ψ-1 of (8):
(12)G̃(x,t)=ψ-1[R(x,t)]=G(x,t)ifR(x,t)>0,G̃(x,t)=αifR(x,t)=0.
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:
(13)〈g̃1g̃2〉-〈g̃1〉〈g̃2〉=∑k=1+∞Ωk2k!cG(d)k,
where g̃1,g̃2,Ωk are as in (11) but now refer to G̃(x,t) and transformation Ω.From a practical point of view, that is, to estimatecG(d) from radar observations of rain field, ψ̃-1 is first applied to the radar data. From (12) and (8), the rain fields R(x,t) are thus turned into truncated Gaussian field G̃(x,t). Afterwards, the covariance of the truncated Gaussian field G̃(x,t) is computed for each radar map performing the inverse Fourier transform of the spatial power spectrum of G̃(x,t). Finally, cG(d) is computed inverting numerically (13), from the average covariance function of G̃(x,t) over all the radar maps. A two-dimensional representation of cG(d) (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 cG(d) reduces to cG(d). The isotropic correlation function cG(d) finally obtained from the radar of Bordeaux is presented in Figure 3.Figure 3
Isotropic correlation functioncG(d) of the Gaussian field G(x,t) derived from the radar observations of Bordeaux and analytical regression (14).A regression process concludes thatcG(d) can be appropriately approximated by
(14)cG(d)=0.59exp⁡⁡(-d31)+0.41exp⁡⁡(-d800).
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 G(x,t) through (4) and (5). Here, the definition of the correlation time τk in (5) is required. From radar observations of rain fields in tropical areas, Bell [19] proposes the formulation:
(15)τk=0.24×(π|k|)2/3.
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 G(x,t) 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 R̂(x,t)=ψ⌊Ĝ(x,t)⌋ which average motion is null are defined. This advection correction allows getting rid of the shift phase exp(-ikTVΔt) in (4). Now, the temporal evolution of the underlying Gaussian fields Ĝ(x,t) reduces to
(16)ak,t+Δt=ak,tβk+(1-βk2)skek,t+Δt,
where ak,t are the coefficients of the spatial Fourier expansion of Ĝ(x,t)—that is, of G(x,t)—and where βk is given by (5). Recalling that ek,t is normal, standard, centered, and uncorrelated, it follows from (16) that
(17)〈ak,t+Δtak,t*〉=〈ak,tak,t*〉βk,
where * is the complex conjugate. Considering (5), the correlation time τk is finally given by:
(18)τk=-Δtln⁡(βk)=-Δtln⁡(〈ak,t+Δtak,t*〉)-ln⁡(〈ak,tak,t*〉).
The term 〈ak,t+Δtak,t*〉 is the spatial cross-spectral density of the spatial processes Ĝ(x,t) and Ĝ(x,t+Δt). It is linked by Fourier transform to the spatiotemporal correlation function cĜ(d,Δt) of Ĝ(x,t). The latter is inferred from the methodology developed in Section 2.3.3. Particularly, the rain fields R̂(x,t) corrected for advection are first transformed into truncated Gaussian fields ĜT(x,t)=ψ̃-1[R̂(x,t)] with ψ̃-1 defined in (12). From fields ĜT(x,t), the spatiotemporal power spectrum is computed, and, by inverse Fourier transform, the space-time correlation function cĜT(d,Δt) is obtained. Equation (13) then allows the derivation of cĜ(d,Δt) and, finally, 〈ak,t+Δtak,t*〉 is obtained by Fourier transform.To evaluate the value ofτk in (18) from the radar data of Bordeaux (Δt=5 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 R̂(x,t) of 2 hours corrected for advection have been selected in the database. For each event, the value of τk has been computed according to the methodology described above. As shown in Figure 4, the average dependency of τk with respect to the spatial frequency k can be approximated by
(19)τk=0.06(1|k|)0.94.
Equation (19) was found assuming that τk follows a power-law function of the spatial frequency |k|. Equation (19) is thus an adjustment of (13) [15] for midlatitude areas. As shown in Figure 4, the dependency of τk with respect to |k| suggests, as expected, that the evolution of large features (rain field stratiform component associated to small values of |k|) is slower (higher correlation time τk) than the one of smaller features (rain field convective components with limited spatial extent, associated to highest |k| values).Figure 4
Evolution of parameterτk as a function of the spatial frequency |k|.


## 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 Z-R relation for midlatitudes Z=300×R1.5 [31], where Z is the radar reflectivity factor in mm6 m−3 and R 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 G(x,t) into rain fields R(x,t) aims at preserving the rain rate local distribution. Nevertheless, it introduces an unrealistic feature on the generated rain fields. Particularly, (7) introduces on R(x,t) an unrealistic link between the spatial fraction fR of R(x,t) that is affected by rain and the spatial average of the strictly positive rain rates values 〈R∣R>0〉. This point is illustrated in Figure 2, where (7) defines 〈R∣R>0〉  as a growing function of fR. This is in contradiction with the ergodicity assumption of rain rate fields [26, 32, 33] that implies that 〈R∣R>0〉 does not change with respect to the fractional area fR 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 N of the simulation area reduces with respect to the decorrelation distance driven by cG(d), the random variables R(x,t) start to deviate from a log normal distribution with parameters μR and σR in (6). Therefore, (7) has to be corrected, introducing the correction terms μoff and σoff:
(8)R(x)=ψ[G(x)]=0ifG(x)<α,R(x)=ψ[G(x)]=exp⁡⁡[μR-μoff⁡σoff⁡+2σRσoff⁡erfc-1(erfc⁡(G(x)/2)P0)]ifG(x)≥α,
The values of the numerical correction terms μoff and σoff must be determined from simulations of R(x,t): they are adjusted to insure that the random variables R(x,t) that compose the N×N rain field are log normal, with parameters μR and σR. 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 cG(d): 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.Figure 2
Relation between the fraction of the areafR that is affected by rain and the conditional mean rain rate 〈R∣R>0〉 computed from 35 000 radar scans and from 35 000 simulations using transformation ψ defined by (7) [19] and the correction (8).Finally, (8) guarantees the statistical independence between the fraction of the simulation area fR that is affected by rain and the conditional mean rain rate 〈R∣R>0〉 as illustrated on Figure 2 where 〈R∣R>0〉 is constant whatever fR 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 processR(x,t). In such conditions, (8) implicitly relates the spatial correlation cG(d) of the underlying Gaussian field G(x,t) to the spatial correlation cR(d) of R(x,t). Particularly, let r1 and r2 be random variables representing rainfall rates that are separated by a distance d. Assuming that r1=ψ(g1) and r2=ψ(g2), where ψ is given by (8) and g1 and g2 are standard centered variable with spatial correlation cG(d), the joint Probability Density Function (PDF) p(g1,g2) is given by
(9)p(g1,g2)=12π1-cG2(d)exp⁡⁡[-12g12+g22-2cG(d)g1g21-cG2(d)].
From (7), α=2erfc-1(2P0) and the rain rate cross expectation 〈r1r2〉 can be related to cG(d) by
(10)〈r1r2〉=〈ψ(g1)ψ(g2)〉=∬-∞+∞ψ(g1)ψ(g2)p(g1,g2)dg1dg2=∬α+∞ψ(g1)ψ(g2)12π1-cG2(d)×exp⁡⁡[-12g12+g22-2cG(d)g1g21-cG2(d)]dg1dg2.
As the numerical integration of (10) is inaccurate and requires a significant computation time, the rain field covariance 〈r1r2〉-〈r1〉〈r2〉 is related to the Gaussian field correlation function cG(d) through a method based on Hermit polynomials [21, 34]:
(11)〈r1r2〉-〈r1〉〈r2〉=∑k=1+∞ψk2k!cG(d)k,
where ψk 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 G̃(x,t)=ψ̃-1[R(x,t)]=Ω[G(x,t)], defined from the analytical inversion ψ-1 of (8):
(12)G̃(x,t)=ψ-1[R(x,t)]=G(x,t)ifR(x,t)>0,G̃(x,t)=αifR(x,t)=0.
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:
(13)〈g̃1g̃2〉-〈g̃1〉〈g̃2〉=∑k=1+∞Ωk2k!cG(d)k,
where g̃1,g̃2,Ωk are as in (11) but now refer to G̃(x,t) and transformation Ω.From a practical point of view, that is, to estimatecG(d) from radar observations of rain field, ψ̃-1 is first applied to the radar data. From (12) and (8), the rain fields R(x,t) are thus turned into truncated Gaussian field G̃(x,t). Afterwards, the covariance of the truncated Gaussian field G̃(x,t) is computed for each radar map performing the inverse Fourier transform of the spatial power spectrum of G̃(x,t). Finally, cG(d) is computed inverting numerically (13), from the average covariance function of G̃(x,t) over all the radar maps. A two-dimensional representation of cG(d) (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 cG(d) reduces to cG(d). The isotropic correlation function cG(d) finally obtained from the radar of Bordeaux is presented in Figure 3.Figure 3
Isotropic correlation functioncG(d) of the Gaussian field G(x,t) derived from the radar observations of Bordeaux and analytical regression (14).A regression process concludes thatcG(d) can be appropriately approximated by
(14)cG(d)=0.59exp⁡⁡(-d31)+0.41exp⁡⁡(-d800).
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 G(x,t) through (4) and (5). Here, the definition of the correlation time τk in (5) is required. From radar observations of rain fields in tropical areas, Bell [19] proposes the formulation:
(15)τk=0.24×(π|k|)2/3.
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 G(x,t) 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 R̂(x,t)=ψ⌊Ĝ(x,t)⌋ which average motion is null are defined. This advection correction allows getting rid of the shift phase exp(-ikTVΔt) in (4). Now, the temporal evolution of the underlying Gaussian fields Ĝ(x,t) reduces to
(16)ak,t+Δt=ak,tβk+(1-βk2)skek,t+Δt,
where ak,t are the coefficients of the spatial Fourier expansion of Ĝ(x,t)—that is, of G(x,t)—and where βk is given by (5). Recalling that ek,t is normal, standard, centered, and uncorrelated, it follows from (16) that
(17)〈ak,t+Δtak,t*〉=〈ak,tak,t*〉βk,
where * is the complex conjugate. Considering (5), the correlation time τk is finally given by:
(18)τk=-Δtln⁡(βk)=-Δtln⁡(〈ak,t+Δtak,t*〉)-ln⁡(〈ak,tak,t*〉).
The term 〈ak,t+Δtak,t*〉 is the spatial cross-spectral density of the spatial processes Ĝ(x,t) and Ĝ(x,t+Δt). It is linked by Fourier transform to the spatiotemporal correlation function cĜ(d,Δt) of Ĝ(x,t). The latter is inferred from the methodology developed in Section 2.3.3. Particularly, the rain fields R̂(x,t) corrected for advection are first transformed into truncated Gaussian fields ĜT(x,t)=ψ̃-1[R̂(x,t)] with ψ̃-1 defined in (12). From fields ĜT(x,t), the spatiotemporal power spectrum is computed, and, by inverse Fourier transform, the space-time correlation function cĜT(d,Δt) is obtained. Equation (13) then allows the derivation of cĜ(d,Δt) and, finally, 〈ak,t+Δtak,t*〉 is obtained by Fourier transform.To evaluate the value ofτk in (18) from the radar data of Bordeaux (Δt=5 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 R̂(x,t) of 2 hours corrected for advection have been selected in the database. For each event, the value of τk has been computed according to the methodology described above. As shown in Figure 4, the average dependency of τk with respect to the spatial frequency k can be approximated by
(19)τk=0.06(1|k|)0.94.
Equation (19) was found assuming that τk follows a power-law function of the spatial frequency |k|. Equation (19) is thus an adjustment of (13) [15] for midlatitude areas. As shown in Figure 4, the dependency of τk with respect to |k| suggests, as expected, that the evolution of large features (rain field stratiform component associated to small values of |k|) is slower (higher correlation time τk) than the one of smaller features (rain field convective components with limited spatial extent, associated to highest |k| values).Figure 4
Evolution of parameterτk as a function of the spatial frequency |k|.


## 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 β0 is equal to 1. Consequently, (4) indicates that the spatial average of the Gaussian field MG=a0 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 Section3.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 areafR that is affected by rain above a given threshold and the (unconditional) area-averaged rain rate 〈R〉 [40, 41]. As an example, Figure 5 shows the linear dependency of 〈R〉 with respect to the fractional area fR deduced from the radar data of Bordeaux.Figure 5
Area averaged rain rate〈R〉 as a function of the fractional area fR above 0.5 mm·h−1 from the radar data set of Bordeaux. Each cross derives from one radar observation. The correlation coefficient with the linear regression (red line) is 0.85.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) VR which corresponds to a spatially averaged rainfall rate integrated during a duration ΔT and the fraction of the GCM pixel that is under rain fR:
(20)fR=VRr̅ΔT,
where r̅=〈R∣R>0〉 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 r̅ can be obtained considering the log normal local distribution of rainfall rates [28, 32]:
(21)r̅=〈R∣R>0〉≥exp⁡⁡(μR+σR22).
Therefore, from the rain amount time series VR given by ERA-40, the associated time series of the fractional area fR can be inferred worldwide from (20) and (21), every ΔT=6 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 areafR and the simulation parameters. Recently, Authors in [43] have proposed a model to link the fraction fG of one Gaussian field realization that is above a given threshold with the Gaussian field average value MG. 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 N×N is normal, with mean MG and standard deviation 1-σG2, where σG2 is the spatial average of the Gaussian field correlation function over the N×N simulation grid. In such conditions, Authors in [43] have shown that the fraction fG of a simulated Gaussian field that is over a preset threshold αG is given by
(22)fG=12erfc⁡(αG-MG2(1-σG2)).
From the Gaussian modeling of rain fields developed in Section 2, αG is linked to the rain no-rain threshold through the equation αG=2erfc-1(2P0), where P0 is the probability of rain. Obviously fG=fR so that the mean value MG of the Gaussian field can be related to the rain fractional area fR by
(23)MG=2erfc⁡-1(2P0)-2(1-σG2)×erfc⁡-1(2fR).
By definition of the Fourier transform, the 0 frequency term a0 of the Fourier expansion of the Gaussian field G(x,t) is equal to its mean value MG. Therefore, whenever a0 is set to MG, the fraction of the simulation area that undergoes rain is fR, 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ΔT is 6 hours. It is thus mandatory to find an interpolation scheme to get values of the rain fractional area fR at a finer time resolution. This high-resolution time series of rain fractional area denoted by fHR(t) and sampled at Δt(Δt<ΔT) must be continuous and, obviously, its average value over ΔT=6 hours must be equal to the low-resolution rain fractional area derived each 6 hours from ERA-40 and denoted by fLR(t) hereinafter. Therefore, considering that fLR(t) holds in the temporal interval [-T/2;T/2], fHR(t) is arbitrarily defined as a piecewise linear function such as
(24)fHR(t)=a′t+b,t∈[-T2;0],fHR(t)=a′′t+b,t∈[0;T2].
To ensure the continuity at the edges of the interval, we set
(25)fHR(-T2)=fLR(-T)fLR(0),fHR(T2)=fLR(T)fLR(0),
combined with the constraint on the average value on [-T/2;T/2]:
(26)1T∫-T/2T/2fHR(t)dt=1T∫-T/20(a′t+b)dt+1T∫0T/2(a′′t+b)dt=fLR(0).
Equations (24), (25), and (26) define a set of three equations with three unknown coefficients a′, a′′, and b. 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 6
Interpolation at higher temporal resolution of the rain fractional area time series derived from ERA-40 (ΔT=6 h). The black dashed line results from the rain fractional area computed from rain field spatiotemporal simulations with high-resolution parameters (i.e., Δt=6mn).Figure7 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 (ΔT=6hours) time series of rain amount, a low-resolution (ΔT=6hours) time series of rain fractional area is deduced from (20). The latter is interpolated at Δt=6mn (i.e., the radar observation frequency) from the piecewise approach defined above. At this stage, a time series of Gaussian field mean values MG(t) is derived from (23). It is injected in the spatiotemporal model through the Fourier coefficient a0(t)=MG(t). 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.Conversion of the rain amount time series given by ERA-40 into rain fractional areafLR(t) and into spatial average of the Gaussian fieldMG(t) that drives the model spatiotemporal evolution in the Fourier plane through a0(t)=MG(t).
(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 size2N×2N is simulated according to the methodology described in Section 2. Then, for each point x=(x1,x2) of the simulation area, the value of the global Gaussian field G results from the weighted sum of the Gaussian fields G1,G2,G3,G4 generated for each of the ERA-40 adjacent resolution cells:(27)G(x1,x2)=(d-x1)(d-x2)G1(x1,x2)+x1(d-x2)G2(x1,x2)+(d-x1)x2G3(x1,x2)+x1x2G4(x1,x2)d,where d is the distance between two ERA-40 resolution cells (d=N). Figure 8 gives an overview of the interpolation scheme.Figure 8
Spatial extension of the modelling to several resolution cells of the ERA-40 reanalysis database. Each gray square figures an ERA-40 resolution cell. The interpolation procedure (27) is applied inside the white square.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 Gi. 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 a0 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.Figure 9
Example of rain field simulated over Europe on 24/02/1999 at 19:30 UTC using 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 V 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 Δt of 6 mn.


## 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 Section3.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 areafR that is affected by rain above a given threshold and the (unconditional) area-averaged rain rate 〈R〉 [40, 41]. As an example, Figure 5 shows the linear dependency of 〈R〉 with respect to the fractional area fR deduced from the radar data of Bordeaux.Figure 5
Area averaged rain rate〈R〉 as a function of the fractional area fR above 0.5 mm·h−1 from the radar data set of Bordeaux. Each cross derives from one radar observation. The correlation coefficient with the linear regression (red line) is 0.85.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) VR which corresponds to a spatially averaged rainfall rate integrated during a duration ΔT and the fraction of the GCM pixel that is under rain fR:
(20)fR=VRr̅ΔT,
where r̅=〈R∣R>0〉 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 r̅ can be obtained considering the log normal local distribution of rainfall rates [28, 32]:
(21)r̅=〈R∣R>0〉≥exp⁡⁡(μR+σR22).
Therefore, from the rain amount time series VR given by ERA-40, the associated time series of the fractional area fR can be inferred worldwide from (20) and (21), every ΔT=6 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 areafR and the simulation parameters. Recently, Authors in [43] have proposed a model to link the fraction fG of one Gaussian field realization that is above a given threshold with the Gaussian field average value MG. 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 N×N is normal, with mean MG and standard deviation 1-σG2, where σG2 is the spatial average of the Gaussian field correlation function over the N×N simulation grid. In such conditions, Authors in [43] have shown that the fraction fG of a simulated Gaussian field that is over a preset threshold αG is given by
(22)fG=12erfc⁡(αG-MG2(1-σG2)).
From the Gaussian modeling of rain fields developed in Section 2, αG is linked to the rain no-rain threshold through the equation αG=2erfc-1(2P0), where P0 is the probability of rain. Obviously fG=fR so that the mean value MG of the Gaussian field can be related to the rain fractional area fR by
(23)MG=2erfc⁡-1(2P0)-2(1-σG2)×erfc⁡-1(2fR).
By definition of the Fourier transform, the 0 frequency term a0 of the Fourier expansion of the Gaussian field G(x,t) is equal to its mean value MG. Therefore, whenever a0 is set to MG, the fraction of the simulation area that undergoes rain is fR, 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ΔT is 6 hours. It is thus mandatory to find an interpolation scheme to get values of the rain fractional area fR at a finer time resolution. This high-resolution time series of rain fractional area denoted by fHR(t) and sampled at Δt(Δt<ΔT) must be continuous and, obviously, its average value over ΔT=6 hours must be equal to the low-resolution rain fractional area derived each 6 hours from ERA-40 and denoted by fLR(t) hereinafter. Therefore, considering that fLR(t) holds in the temporal interval [-T/2;T/2], fHR(t) is arbitrarily defined as a piecewise linear function such as
(24)fHR(t)=a′t+b,t∈[-T2;0],fHR(t)=a′′t+b,t∈[0;T2].
To ensure the continuity at the edges of the interval, we set
(25)fHR(-T2)=fLR(-T)fLR(0),fHR(T2)=fLR(T)fLR(0),
combined with the constraint on the average value on [-T/2;T/2]:
(26)1T∫-T/2T/2fHR(t)dt=1T∫-T/20(a′t+b)dt+1T∫0T/2(a′′t+b)dt=fLR(0).
Equations (24), (25), and (26) define a set of three equations with three unknown coefficients a′, a′′, and b. 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 6
Interpolation at higher temporal resolution of the rain fractional area time series derived from ERA-40 (ΔT=6 h). The black dashed line results from the rain fractional area computed from rain field spatiotemporal simulations with high-resolution parameters (i.e., Δt=6mn).Figure7 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 (ΔT=6hours) time series of rain amount, a low-resolution (ΔT=6hours) time series of rain fractional area is deduced from (20). The latter is interpolated at Δt=6mn (i.e., the radar observation frequency) from the piecewise approach defined above. At this stage, a time series of Gaussian field mean values MG(t) is derived from (23). It is injected in the spatiotemporal model through the Fourier coefficient a0(t)=MG(t). 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.Conversion of the rain amount time series given by ERA-40 into rain fractional areafLR(t) and into spatial average of the Gaussian fieldMG(t) that drives the model spatiotemporal evolution in the Fourier plane through a0(t)=MG(t).
(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 size2N×2N is simulated according to the methodology described in Section 2. Then, for each point x=(x1,x2) of the simulation area, the value of the global Gaussian field G results from the weighted sum of the Gaussian fields G1,G2,G3,G4 generated for each of the ERA-40 adjacent resolution cells:(27)G(x1,x2)=(d-x1)(d-x2)G1(x1,x2)+x1(d-x2)G2(x1,x2)+(d-x1)x2G3(x1,x2)+x1x2G4(x1,x2)d,where d is the distance between two ERA-40 resolution cells (d=N). Figure 8 gives an overview of the interpolation scheme.Figure 8
Spatial extension of the modelling to several resolution cells of the ERA-40 reanalysis database. Each gray square figures an ERA-40 resolution cell. The interpolation procedure (27) is applied inside the white square.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 Gi. 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 a0 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.Figure 9
Example of rain field simulated over Europe on 24/02/1999 at 19:30 UTC using the concurrent ERA-40 rain amount data.


## 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 areafR that is affected by rain above a given threshold and the (unconditional) area-averaged rain rate 〈R〉 [40, 41]. As an example, Figure 5 shows the linear dependency of 〈R〉 with respect to the fractional area fR deduced from the radar data of Bordeaux.Figure 5
Area averaged rain rate〈R〉 as a function of the fractional area fR above 0.5 mm·h−1 from the radar data set of Bordeaux. Each cross derives from one radar observation. The correlation coefficient with the linear regression (red line) is 0.85.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) VR which corresponds to a spatially averaged rainfall rate integrated during a duration ΔT and the fraction of the GCM pixel that is under rain fR:
(20)fR=VRr̅ΔT,
where r̅=〈R∣R>0〉 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 r̅ can be obtained considering the log normal local distribution of rainfall rates [28, 32]:
(21)r̅=〈R∣R>0〉≥exp⁡⁡(μR+σR22).
Therefore, from the rain amount time series VR given by ERA-40, the associated time series of the fractional area fR can be inferred worldwide from (20) and (21), every ΔT=6 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 areafR and the simulation parameters. Recently, Authors in [43] have proposed a model to link the fraction fG of one Gaussian field realization that is above a given threshold with the Gaussian field average value MG. 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 N×N is normal, with mean MG and standard deviation 1-σG2, where σG2 is the spatial average of the Gaussian field correlation function over the N×N simulation grid. In such conditions, Authors in [43] have shown that the fraction fG of a simulated Gaussian field that is over a preset threshold αG is given by
(22)fG=12erfc⁡(αG-MG2(1-σG2)).
From the Gaussian modeling of rain fields developed in Section 2, αG is linked to the rain no-rain threshold through the equation αG=2erfc-1(2P0), where P0 is the probability of rain. Obviously fG=fR so that the mean value MG of the Gaussian field can be related to the rain fractional area fR by
(23)MG=2erfc⁡-1(2P0)-2(1-σG2)×erfc⁡-1(2fR).
By definition of the Fourier transform, the 0 frequency term a0 of the Fourier expansion of the Gaussian field G(x,t) is equal to its mean value MG. Therefore, whenever a0 is set to MG, the fraction of the simulation area that undergoes rain is fR, 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ΔT is 6 hours. It is thus mandatory to find an interpolation scheme to get values of the rain fractional area fR at a finer time resolution. This high-resolution time series of rain fractional area denoted by fHR(t) and sampled at Δt(Δt<ΔT) must be continuous and, obviously, its average value over ΔT=6 hours must be equal to the low-resolution rain fractional area derived each 6 hours from ERA-40 and denoted by fLR(t) hereinafter. Therefore, considering that fLR(t) holds in the temporal interval [-T/2;T/2], fHR(t) is arbitrarily defined as a piecewise linear function such as
(24)fHR(t)=a′t+b,t∈[-T2;0],fHR(t)=a′′t+b,t∈[0;T2].
To ensure the continuity at the edges of the interval, we set
(25)fHR(-T2)=fLR(-T)fLR(0),fHR(T2)=fLR(T)fLR(0),
combined with the constraint on the average value on [-T/2;T/2]:
(26)1T∫-T/2T/2fHR(t)dt=1T∫-T/20(a′t+b)dt+1T∫0T/2(a′′t+b)dt=fLR(0).
Equations (24), (25), and (26) define a set of three equations with three unknown coefficients a′, a′′, and b. 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 6
Interpolation at higher temporal resolution of the rain fractional area time series derived from ERA-40 (ΔT=6 h). The black dashed line results from the rain fractional area computed from rain field spatiotemporal simulations with high-resolution parameters (i.e., Δt=6mn).Figure7 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 (ΔT=6hours) time series of rain amount, a low-resolution (ΔT=6hours) time series of rain fractional area is deduced from (20). The latter is interpolated at Δt=6mn (i.e., the radar observation frequency) from the piecewise approach defined above. At this stage, a time series of Gaussian field mean values MG(t) is derived from (23). It is injected in the spatiotemporal model through the Fourier coefficient a0(t)=MG(t). 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.Conversion of the rain amount time series given by ERA-40 into rain fractional areafLR(t) and into spatial average of the Gaussian fieldMG(t) that drives the model spatiotemporal evolution in the Fourier plane through a0(t)=MG(t).
(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 size2N×2N is simulated according to the methodology described in Section 2. Then, for each point x=(x1,x2) of the simulation area, the value of the global Gaussian field G results from the weighted sum of the Gaussian fields G1,G2,G3,G4 generated for each of the ERA-40 adjacent resolution cells:(27)G(x1,x2)=(d-x1)(d-x2)G1(x1,x2)+x1(d-x2)G2(x1,x2)+(d-x1)x2G3(x1,x2)+x1x2G4(x1,x2)d,where d is the distance between two ERA-40 resolution cells (d=N). Figure 8 gives an overview of the interpolation scheme.Figure 8
Spatial extension of the modelling to several resolution cells of the ERA-40 reanalysis database. Each gray square figures an ERA-40 resolution cell. The interpolation procedure (27) is applied inside the white square.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 Gi. 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 a0 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.Figure 9
Example of rain field simulated over Europe on 24/02/1999 at 19:30 UTC using 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 V 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 Δt 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(28)A=∫0LskRα(x)dx,
where R(x) is the rain rate intensity (mm h−1) at x and Ls is the slant path (km) through rain. k 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 N×N 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.Figure 10
Rain height for Bordeaux (France) from Rec. ITU-R P.839 [43] and derived from ERA-40 temperature profiles.


### 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 Section3 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.Figure 11
1 Hz attenuation time series simulated on two sites located 7 km apart, at 20 GHz. The green circles represent the values of the initial time series extracted from the spatiotemporal model (Δt=6mn).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).


## 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(28)A=∫0LskRα(x)dx,
where R(x) is the rain rate intensity (mm h−1) at x and Ls is the slant path (km) through rain. k 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 N×N 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.Figure 10
Rain height for Bordeaux (France) from Rec. ITU-R P.839 [43] and derived from ERA-40 temperature profiles.


## 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 Section3 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.Figure 11
1 Hz attenuation time series simulated on two sites located 7 km apart, at 20 GHz. The green circles represent the values of the initial time series extracted from the spatiotemporal model (Δt=6mn).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 Figure12.6 h rain amount derived from ERA-40 database on 12/24/1999 at 08:00 UTC (a) and corresponding 6 h rain amount computed from 6 mn rain field spatiotemporal simulations (b).
(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.Rain attenuation CCDF and attenuation CCDF at 20 GHz computed from one year of simulations compared to ITU-R P. Rec 837-5 and ITU-R P. Rec 618-10.
(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.Figure 14
Attenuation CCDF at 20 GHz measured during the OLYMPUS campaign (1992) in two sites near Paris (La-Folie-Bessin and Gometz-la-ville, color lines) and attenuation CCDF derived from simulated attenuation time series (black lines). The two sites are located 7 km apart.According to Figure14, 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.Figure 15
Diversity gain as a function of distance for different single site attenuation values derived from radar data (blue curves) and derived from simulations (red curves).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.Figure 16
Joint attenuation CCDF derived from simulations and derived from OLYMPUS measurements at 20 GHz, at La-Folie-Bessin and Gometz-La-Ville (7 km apart).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 Figure17, 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.Figure 17
Temporal correlation function estimated from radar-derived attenuation time series and from simulations.


## 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 Figure12.6 h rain amount derived from ERA-40 database on 12/24/1999 at 08:00 UTC (a) and corresponding 6 h rain amount computed from 6 mn rain field spatiotemporal simulations (b).
(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.Rain attenuation CCDF and attenuation CCDF at 20 GHz computed from one year of simulations compared to ITU-R P. Rec 837-5 and ITU-R P. Rec 618-10.
(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.Figure 14
Attenuation CCDF at 20 GHz measured during the OLYMPUS campaign (1992) in two sites near Paris (La-Folie-Bessin and Gometz-la-ville, color lines) and attenuation CCDF derived from simulated attenuation time series (black lines). The two sites are located 7 km apart.According to Figure14, 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.Figure 15
Diversity gain as a function of distance for different single site attenuation values derived from radar data (blue curves) and derived from simulations (red curves).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.Figure 16
Joint attenuation CCDF derived from simulations and derived from OLYMPUS measurements at 20 GHz, at La-Folie-Bessin and Gometz-La-Ville (7 km apart).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 Figure17, 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.Figure 17
Temporal correlation function estimated from radar-derived attenuation time series and from simulations.


## 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.



---

*Source: 749829-2012-05-07.xml*