Abstract

The multifractal properties of rain are investigated within the framework of universal multifractals. The database used in this study includes measurements performed over several months in different locations by means of a disdrometer, the dual-beam spectropluviometer (DBS). An assessment of the effect of the rain–no rain intermittency shows that the analysis of rain-rate time series may lead to a spurious break in the scaling and to erroneous parameters. The estimation of rain multifractal parameters is, therefore, performed on an event-by-event basis, and they are found to be significantly different from those proposed in scientific literature. In particular, the parameter H, which has often been estimated to be 0, is more likely to be 0.53, thus meaning that rain is a fractionally integrated flux (FIF). Finally, a new model is proposed that simulates high-resolution rain-rate time series based on these new parameters and on a simple threshold.

1. Introduction

Rain is a complex and highly intermittent phenomenon and therefore difficult to model and even to measure. During the past decades, thanks to advances in the turbulence theory and in the description of multiplicative cascades, various models based on scale invariance and multifractals have been proposed (see, e.g., Veneziano et al. 1996; Over and Gupta 1996; Mazzarella 1999; Deidda 2000). These models open the way to a realistic description of rain, since it is not possible to dissociate rain from atmospheric turbulence. In this paper, we focus on the universal multifractals framework and more specifically on the fractionally integrated flux (FIF) model developed by Schertzer and Lovejoy (1987, 1997). This model was chosen because of its strong links to turbulence theory and because most multiplicative processes tend to resort to it. It has been applied successfully to a wide variety of phenomena, from cloud radiance (Lovejoy and Schertzer 2006) to ocean color (Lovejoy et al. 2001), or to financial assets (Schmitt et al. 1999).

As far as rain is concerned, this model has been tested against experimental data measured by rain gauges, meteorological radars [see Lilley et al. (2006) for a review], or—more recently—remote sensing (Lovejoy et al. 2008). Surprisingly, although this type of data is useful in general climatologic or weather statistical studies, very little is known in scientific literature about high-resolution rain-rate time series retrieved from disdrometer measurements. The latter provide individual drop data, thus leading to a large number of explored time scales and overcoming some limitations of others’ instruments—namely, the uncertainty affecting rain-rate retrieval in satellite and radar measurements and the variability of the sampling frequency of standard rain gauges due to the integrating effect of the tipping bucket (de Lima 1998). In this paper, measurements from the dual-beam spectropluviometer (DBS) developed at the Centre d’étude des Environnements Terrestre et Planétaires (CETP) are analyzed. The datasets include rain-rate time series covering several months from three different locations. The experiments are described in section 4, and the multifractal analysis techniques are presented in section 5.

Rain-rate time series have already been extensively studied within the universal multifractals framework, and many authors have reported similar results [typically H ≈ 0, α ≈ 0.5–0.7, and C1 ≈ 0.4–0.6, using daily accumulations in Tessier et al. (1993, 1996), Hubert et al. (1993), Ladoy et al. (1993), and Fraedrich and Larnder (1993), or using accumulations during about 10 min in de Lima and Grasman (1999), de Lima (1998), and Hubert et al. (2002)]. However, when looking at high-resolution DBS data, a break in the scaling is observed at higher frequencies. This deviation from the model [previously pointed out by Fraedrich and Larnder (1993), Pathirana et al. (2004), and Onof et al. (1996)] leads to the following question: do smaller time scales have different scaling properties from larger time scales? In this paper, we investigate the possibility that this problem may not be due to physical reasons but to an artifact related to rain–no rain intermittency. Rain-rate time series include many no-rain periods with zero rain rates, which universal multifractals are not designed to model, and as suggested by Schmitt et al. (1998), the effect of this rain–no rain intermittency on the estimation of the universal multifractals model parameters has to be taken into account. This effect will indeed be shown to possibly lead to an artificial scaling at low resolutions and thus to biased multifractal parameters (section 3). The data are, therefore, analyzed on an event-by-event basis to avoid this difficulty (section 6). Long rain-rate time series that include rain–no rain intermittency (section 7) will also be analyzed to better characterize the effect of rain–no rain intermittency. Using the new parameters found and on a simple threshold, a new model that simulates high-resolution rain-rate time series will be presented (section 8) as well as the properties of the corresponding synthetic time series (section 9).

2. The fractionally integrated flux model

The FIF model is based on the concept of multiplicative cascades that had been suggested by Richardson (1922) in the context of turbulence and notably developed by Kolmogorov (1941), Yaglom (1966), and Mandelbrot (1974). According to this model, the kinetic energy injected into the system at a large scale is transferred by a conservative multiplicative process to smaller and smaller scales and finally dissipated as heat. More generally, this type of process can be modeled by (i) uniformly distributing a given quantity during a given interval; (ii) splitting this interval into subintervals; (iii) assigning these subintervals the original given quantity multiplied by a random variable whose law does not depend on scale; and (iv) repeating steps (ii) and (iii) iteratively. By construction, such processes exhibit interesting scale invariance properties and to describe them, it is necessary to introduce the concept of fractal dimension.

a. The fractal dimension and codimension

The definition of the fractal dimension based on the box-counting method will be recalled here [see Falconer (2003) for more details and the classical example of Cantor dust and other definitions of the fractal dimension]. Assume an observation period of duration T including rain periods and drought periods. Then split this original interval into contiguous subintervals of length Δt (here and throughout, the scale is given by Δt and the resolution by the number of subintervals λ = Tt). Next, count the number of subintervals nrain(λ) where rain is present. If this number is a power law of resolution, then the rain process is fractal and the power law exponent Df is its fractal dimension (here and throughout, ∼ stands for equality within slowly varying and constant factors):

 
formula

The fractal dimension and the fractal codimension are related to each other by the following equation in a one-dimensional embedding space:

 
formula

The fractal codimension is an interesting concept, because it characterizes the probability that a randomly chosen subinterval Iλ contains rain:

 
formula

For the rain-rate process Df ≈ 0.8, Lavergnat and Golé (1998) found 0.82 for time scales ranging from 10−2 to 105 min, and Hubert and Carbonnel (1989) found 0.8 for time scales ranging from 1 day to 45 yr.

b. Singularities and multifractals

The concept of fractal dimension is restricted to ensembles and therefore to the binary case rain or no rain. To study the scale invariance properties of a process whose intensities vary across a whole range of possible values, the multifractal formalism must be introduced. The idea is that the fractal dimension is not unique but rather depends on the intensity that is used to define the ensemble. Moreover, the scale-dependent notion of intensity has to be left aside, and the scale invariant concept of singularity (denoted γ) has to be used instead. The intensity related to a given singularity is the γth power of resolution λ, so that when the scale varies, the intensity is modified but not the singularity. The basic equation of the “codimension multifractal formalism” is as follows [see Schertzer et al. (2002) for a detailed introduction to multifractal processes and multiplicative cascades]:

 
formula

where Φλ is the process at resolution λ and c(γ) the fractal codimension related to the singularity γ. Since the probability density is related to the statistical moments of Φλ by the Mellin transform, Eq. (4) can equivalently be written (Schertzer et al. 2002, section 4.2) as

 
formula

where <·> is the mean operator, q the order of the moment, and K(q) the so-called moment scaling function. Here, K(q) and c(γ) are linked to each other via the Legendre transform, which is a one-to-one relation (Parisi and Frisch 1985):

 
formula

Here, K(q) is a convex function, and it can easily be shown that K(0) = 0 and K(1) = 0. The latter property is a consequence of the assumption of conservation, namely, that the mean of the process is a constant M regardless of resolution,

 
formula

The notion of conservative process is derived from turbulence theory, according to which the energy flux density between smaller and smaller eddies is assumed to be constant down to the dissipation scale. Multiplicative cascades with a constant mean are, therefore, called conservative cascades.

c. Universal multifractals and fractional integration

By densifying the cascade, that is, by increasing the number of multiplicative steps between the inner and the outer scales of the cascade, it has been shown (Schertzer and Lovejoy 1997) that multifractal cascades generally tend toward a universal class of processes called universal multifractals, whose moment scaling function K(q) can be described by only two parameters (α and C1),

 
formula

Here, C1 is the fractal codimension of the level providing the dominant contribution to the mean of the process; α describes how rapidly the fractal codimension varies with the singularity (α varies from 0 to 2; α = 0 is the monofractal case; and α = 2 is the lognormal case).

In general, geophysical processes are not stationary and an additional fractional integration is required, which leads to the following expression for the rain-rate process (see the appendix for a formal definition of fractional integration):

 
formula

where ΔRt) is the rain-rate gradient over time lags of length Δt and where parameter H is related to the degree of smoothness. When Φλ is a conservative cascade, the process is referred to as a fractionally integrated flux [see Gagnon et al. (2006) for a comparison with classical models, such as the fractional Brownian motion or fractional Lévy motion]. Combining Eqs. (9) and (5), the qth-order structure function of the rain-rate process can be written as

 
formula

where ζ(q) is known as the structure function exponent and has the following form:

 
formula

3. The effect of rain–no rain intermittency

By definition, the FIF model cannot generate zero-value periods, because it is a fractional integration of a multiplicative cascade with strictly positive multiplicative increments. Therefore, it is not clear whether it would be possible to perform a multifractal analysis of time series that include rain–no rain intermittency. This section shows that the effect of rain–no rain intermittency may strongly affect the estimation of the universal multifractals model parameters.

As for the estimation of H, the effect of the rain–no rain intermittency can be understood by considering that, for Δt sufficiently long, if R(t) ≠ 0, the probability of having R(t + Δt) = 0 is close to 1 and therefore

 
formula

Thus, for Δt sufficiently long, the first-order structure function does not depend on Δt, and H is found to be 0.

Concerning the estimation parameters of the stationary cascade (α and C1) with the moment scaling function K(q), a simplified case may be considered. Since rain events have a limited extension in space and time, assume that, at a sufficiently long time scale Δt, the successive averaging operations performed to reconstruct the conservative cascade are roughly similar to the degeneracy of a rectangular function (a function equal to zero everywhere except over an interval of duration Δt). Let Θλmax be the rectangular function at resolution λmax = Tt and θλmax the level of the rectangle itself. When Θλmax is degraded to resolution λmax/2 by averaging pairs of contiguous values, since the rectangle is surrounded by zeros, its level falls to

 
formula

If this operation is repeated n times, Eq. (13) can be generalized to

 
formula

Then, the qth-order moments of the rectangular function at resolution λmax/2n are easily calculated because, at all scales, only one value is not zero. This value is the level of the rectangle θmax/2n raised to the qth power and divided by the resolution λmax/2n:

 
formula

Modifying Eq. (15) to make the scaling explicitly appear leads to

 
formula

Then, by taking the logarithm of both sides and since λλmax/2n,

 
formula

Consequently, a scaling is found and the corresponding moment scaling function is

 
formula

A comparison between Eq. (18) and the theoretical form of K(q) given by Eq. (8) leads to α = 0 and C1 = 1. This simplified example shows that, at low resolutions, rain–no rain intermittency may result in an artificial scaling, with α close to 0 and C1 close to 1. To avoid this artifact, the multifractal analysis of the spectropluviometer data will be performed on an event-by-event basis in section 6.

4. Measurements

The rain rate has been measured by means of a DBS (Delahaye et al. 2006). The DBS is more accurate than a previously used instrument, the optical spectropluviometer (OSP; Hauser et al. 1984), because it comprises two closely spaced flat optical beams and an increased catchment surface (100 against 50 cm2 for commercial instruments). With this instrument, false detections due to air turbulence can be reduced and raindrops as small as 0.3 mm can be detected unambiguously by verifying that each drop actually crosses both beams. In the laboratory, the instrument’s accuracy for particle diameters is 1% for bias and 3% for standard deviation. This accuracy has also been tested in the outside turbulent atmosphere with beads larger than 1 mm. For smaller diameters, calibration can only be performed indoors. The DBS also provides an accurate estimation of the vertical velocity of raindrops. Its outputs are raindrop diameters, velocities, and times of arrival (in milliseconds). These outputs were transformed into rain-rate time series sampled at 1 s. However, the retrieved rain rate is not accurate with a 1-s integration time, and the time series have to be resampled at a longer sampling period, typically 60 s (see sections 5a and 6 for more details). The standard deviations of the rain-rate estimation error are 31% and 4% for sampling periods of 1 and 60 s, respectively, but the error is much larger for smaller rain rates (Mallet and Barthès 2009).

The experimental dataset (Table 1) comprises three subsets collected in different climatic areas: Paris [Site Instrumental de Recherche par Télédétection Atmosphérique (SIRTA) experimental site, France], Iowa City [Disdrometer Evaluation Experiment (DEVEX) campaign, United States; Krajewski et al. 2006], and Djougou [African Monsoon Multidisciplinary Analyses (AMMA) campaign, Benin]. The latter corresponds to the African monsoon period and is characterized by strong rain events.

Table 1.

Experimental datasets.

Experimental datasets.
Experimental datasets.

5. Multifractal analysis techniques

This section describes the analysis tools needed to perform the multifractal analysis. Each technique is illustrated based on an example derived from a strong rain event observed on 12 July 2002 during the DEVEX experiment in Iowa (Fig. 1).

Fig. 1.

Rain event observed in Iowa on 12 Jul 2002. (top) Original rain-rate time series sampled at 1 Hz, (middle) power spectrum with corresponding fit and (bottom) rain rate averaged over 32-s time lags.

Fig. 1.

Rain event observed in Iowa on 12 Jul 2002. (top) Original rain-rate time series sampled at 1 Hz, (middle) power spectrum with corresponding fit and (bottom) rain rate averaged over 32-s time lags.

a. Power spectrum analysis

The power spectrum is a useful tool to decide whether the process is scaling or not. If scaling is present, the power spectrum, plotted in a log–log graph, should be a straight line. Then, it has to be decided whether the process is a stationary cascade or a FIF. This can be easily done by considering the spectral slope β: if the process is a stationary cascade, then, as a consequence of the Wiener–Khintchine theorem, the theoretical slope is

 
formula

Therefore, if β >1, the process cannot be stationary. An additional fractional integration is required, which leads to the FIF model with the following theoretical slope:

 
formula

The power spectrum of the 12 July 2002 event is presented in Fig. 1. For this particular event, β is found to be 1.72, and the process is therefore a FIF.

The limits of the power spectrum scaling are also important, because they indicate the maximum and minimum scales over which it is meaningful to perform a multifractal analysis. For the 12 July 2002 event, the spectrum is clearly leveling off at frequencies above 1/64 Hz (Fig. 1) and therefore the time series has to be averaged over 32-s time lags to keep only the scaling part of the power spectrum. This is a fundamental feature because, in the case of a FIF, the multifractal analysis depends on the gradient of the process at the finest time scale [see Eq. (21)], which therefore would have no physical meaning with a very noisy process.

b. Structure function analysis

Parameter H can be derived by means of Eq. (20). However, generally, the power spectrum’s slope cannot be determined with sufficient accuracy. Another possibility is to use the first-order structure function: Eqs. (11) and (8) show that for q = 1, the structure function exponent simplifies to ζ(1) = H. The first-order structure function as a function of the time lag Δt (in a log–log graph) is, therefore, a straight line with a slope of H. Using this method, H is found to be 0.53 for the 12 July 2002 event (Fig. 2).

Fig. 2.

First-order structure function of the rain event observed in Iowa on 12 Jul 2002 with corresponding fit.

Fig. 2.

First-order structure function of the rain event observed in Iowa on 12 Jul 2002 with corresponding fit.

c. Moment scaling analysis

To check the scaling of the statistical moments of Φλ [Eq. (5)] and to estimate the values of α and C1, the multiplicative cascade Φλ has to be retrieved over a wide range of scales. This is done first by deducing Φλ at the finest scale (denoted Φλmax) from the rain-rate time series, and then by reconstituting the multiplicative cascade at longer time scales [see Gagnon et al. (2006) for an instructive example]. Theoretically, the first step requires an H-order fractional differentiation of the rain-rate process. However, since H ≈ 0.5, an adequate approximation (Lavallée et al. 1993) is to take the absolute value of the gradient at the finest time scale Δtmin (i.e., the highest resolution λmax), which is normalized to have M = 1:

 
formula

The multiplicative cascade is then reconstituted at longer time scales by consecutive averagings of pairs of contiguous values. Finally, the statistical moments are calculated at each scale and plotted in a log–log graph as functions of the scale. According to Eq. (5), these functions should be straight lines with slopes equal to K(q), where q is the order of the moment. Such straight lines are indeed found for the 12 July 2002 event (Fig. 3, top), so that the corresponding moment scaling function K(q) can be estimated (Fig. 3, bottom).

Fig. 3.

(top) Scaling of moments for the rain event observed in Iowa on 12 Jul 2002 with q = (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and (bottom) corresponding moment scaling function.

Fig. 3.

(top) Scaling of moments for the rain event observed in Iowa on 12 Jul 2002 with q = (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and (bottom) corresponding moment scaling function.

d. Estimation of α and C1

The multifractal parameters α an C1 can be estimated using the moment scaling function K(q) given in Eq. (8). However, the analysis of synthetic time series created with known multifractal parameters shows that α is not retrieved correctly by directly fitting K(q) (the simulation technique is described in the appendix). A more sophisticated technique known as double trace moment (DTM) can improve the accuracy of the estimation (see Lavallée et al. 1993 for details). For that purpose, the ηth power of Φλmax is first computed and the multiplicative cascade is then reconstituted at larger time scales. The statistical moments of such a process, Mλ(η,q), are scaling, and their scaling exponent is denoted K(η,q):

 
formula

Exponents K(η,q) have the following property in the framework of universal multifractals:

 
formula

Therefore, when q has a fixed value (different from the special case 0 or 1), a log–log plot of K(η, q) as a function of η should give a straight line with a slope of α. The scaling of the double trace moments and the corresponding scaling function K(η,q) are shown in Fig. 4 for the 12 July 2002 event. The fit of this scaling function to its linear part leads to α = 1.69. The classical moment scaling function K(q) is then fitted according to Eq. (8) to retrieve C1, which is found to be 0.12.

Fig. 4.

(top) Scaling of DTM for the rain event observed in Iowa on 12 Jul 2002 with η varying (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and q set to 2. (bottom) Scaling function of the DTM.

Fig. 4.

(top) Scaling of DTM for the rain event observed in Iowa on 12 Jul 2002 with η varying (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and q set to 2. (bottom) Scaling function of the DTM.

6. Event-by-event analysis

Since the presence of rain–no rain intermittency may bias the estimated parameters, it was chosen to avoid this difficulty by performing the multifractal analysis of spectropluviometer data on an event-by-event basis. In this section, 30 rain events were analyzed separately (10 for each of the three datasets, see Tables 2 –4). These rain events were extracted in such a way that all zero values were excluded. The rain-rate process is thus said to be apparently uninterrupted, because although there are no zero values at the resolution at which they are analyzed, such values would appear if the resolution was increased.

Table 2.

Characteristics of the events extracted from the SIRTA dataset.

Characteristics of the events extracted from the SIRTA dataset.
Characteristics of the events extracted from the SIRTA dataset.
Table 4.

Same as Table 2 but for the AMMA dataset.

Same as Table 2 but for the AMMA dataset.
Same as Table 2 but for the AMMA dataset.

The power spectrum was calculated for each rain event, and the highest resolution at which the scaling is valid was determined (so as to avoid flattening the power spectrum at the highest frequencies). For most events, a 64-s averaging time was used. However, for certain events, a 32-s averaging time was sufficient. The leveling off of the power spectrum occurring at higher frequencies could be because the sampling frequency exceeds the inner scale of the multiplicative cascade (i.e., the dissipative scale in the turbulence framework). However, Lovejoy and Schertzer (2008) have shown that rain becomes decoupled from atmospheric turbulence at a spatial scale of 1 m. Since most of the drops are falling faster than 1 m s−1 (Gunn and Kinzer 1949), the rain-rate process should at least scale up to the 1-s time scale. A more realistic explanation could be that, when the rain rate is calculated with an integration time that is too short, the number of raindrops detected by the instrument is too small and the measurement is not representative of the process (although the DBS has been designed with a relatively large catchment surface of 100 cm2 to counter this limitation). Another explanation could be that the vertical velocity of the drops is affected at higher frequencies by the turbulence due to the proximity to the earth’s surface. However this point remains to be investigated.

The event-by-event multifractal analysis led to the multifractal parameters presented in Table 5. Instead of the values H ≈ 0, α ≈ 0.5–0.7, and C1 ≈ 0.4–0.6 generally found in scientific literature (see references cited in section 1), significantly different values were found, namely, H ≈ 0.53, α ≈ 1.7, and C1 ≈ 0.13. The fact that H is not 0 indicates that the process is not a stationary cascade but rather a nonstationary FIF. Although this result is consistent with those of Harris et al. (1996) who used high-resolution rain-rate time series (15-s integration time), this contradicts most other studies in which the cascade is found to be stationary. This difference can be explained by the fact that these studies used low-resolution data and therefore were restricted to the flat part of the spectrum observed at the lower frequencies. The standard deviation of H is large (27% of its mean value). Therefore, it is not possible to conclude that H is constant for all rain events. However, H was difficult to estimate because of its sensitivity to the noise level and the sample size, which is small in an event-by-event analysis. Thus, the variability of H might be due to instrumental reasons and, to decide whether H is constant or not, it would be interesting to increase the catchment surface so that finer time scales could be explored.

Table 5.

Multifractal parameter derived from the event-by-event analysis.

Multifractal parameter derived from the event-by-event analysis.
Multifractal parameter derived from the event-by-event analysis.

As to the other parameters, α is large, which means that the multiplicative cascade is close to the lognormal case. As expected, the value of C1 is smaller than the one found in scientific literature, since no rain intermittency artificially increases the mean level of intermittency. Note that the values of α and C1 are remarkably similar regardless of the dataset and that their variability is small, although the data have been collected in very different climatic areas (the climate is continental for Iowa, temperate–oceanic for Paris, and tropical with monsoon for Djougou).

7. Analysis of long rain-rate time series

In this section, a multifractal analysis of the three rain-rate time series is performed without extracting individual events. Since these time series last several months, they inherently include zero-value periods, so that rain–no rain intermittency should have an effect.

Figure 5 shows the power spectrum of each time series. At time scales longer than 1 min, a scaling is observed up to the 1-h scale (corresponding to 1/7200 Hz in Fig. 5) and then the spectrum becomes flat. Between the 1-min and 1-h time scales, the spectral slope was found to be 1.76 for AMMA data, 1.79 for DEVEX data, and 1.77 for SIRTA data. Figure 6 shows the first-order structure function for each dataset and corresponding fits that were performed only for time lags shorter than 15 min. It was found that H = 0.55 for AMMA data, H = 0.52 for DEVEX data, and H = 0.5 for SIRTA data. At larger time intervals, a flattening is also observed. These results are consistent with the assessment of the effect of rain–no rain intermittency performed in section 3, because above a time interval that roughly corresponds to the typical duration of a rain event (with data sampled every minute), the first-order structure function flattens and yields H ≈ 0.

Fig. 5.

Power spectrum of the experimental rain-rate time series and corresponding fits between 1-min and 1-h time scales (corresponding to 1/120 and 1/7200 Hz, respectively). The curves have been vertically shifted, (top) SIRTA, (middle) AMMA, and (bottom) DEVEX data.

Fig. 5.

Power spectrum of the experimental rain-rate time series and corresponding fits between 1-min and 1-h time scales (corresponding to 1/120 and 1/7200 Hz, respectively). The curves have been vertically shifted, (top) SIRTA, (middle) AMMA, and (bottom) DEVEX data.

Fig. 6.

First-order structure function of the experimental rain-rate time series and corresponding fits [the curves have been vertically shifted, (top) SIRTA, (middle) AMMA, and (bottom) DEVEX data].

Fig. 6.

First-order structure function of the experimental rain-rate time series and corresponding fits [the curves have been vertically shifted, (top) SIRTA, (middle) AMMA, and (bottom) DEVEX data].

Figure 7 shows the scaling of moments and the corresponding moment scaling function K(q) for the AMMA dataset. This scaling also exhibits a break and, if only low resolutions are considered (time scales longer than 1 h), the parameters are found to be α = 0.24 and C1= 0.63 [note that these values are the same whether the process is assumed to be a stationary cascade (the multiplicative cascade is reconstructed directly from the rain rate) or a nonstationary FIF (the multiplicative cascade is reconstructed from the absolute value of the gradient)]. These parameters are close to those suggested by other authors who used low-resolution data (α ≈ 0.5–0.7; C1 ≈ 0.4–0.6) and consistent with section 3, where it is shown that rain–no rain intermittency may lower α and increase C1.

Fig. 7.

Moment scaling function analysis of the AMMA time series. (top) Scaling of the moment with q = (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and (bottom) moment scaling function K(q).

Fig. 7.

Moment scaling function analysis of the AMMA time series. (top) Scaling of the moment with q = (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and (bottom) moment scaling function K(q).

8. Simulation of the rain–no rain intermittency

A new model called “thresholded FIF” is proposed, which simulates high-resolution rain-rate time series based on the new parameters and a simple threshold. A FIF process is first generated with H ≈ 0.53, α ≈ 1.7, and C1 ≈ 0.13. Then lower values lying below a given threshold are set to 0. Those above the threshold are shifted downward by subtracting the threshold value. The threshold is chosen so that rain–no rain intermittency reproduces the fractal dimension of the rain “support.” Contrary to stationary cascades, the values of FIF processes do not have an absolute meaning. Therefore, it is possible to shift them up or down without breaking the scaling (this operation does not modify the differentiated process). Note that in this model, the threshold corresponds to an intensity of the fractionally integrated cascade, which means that it does not correspond to a unique singularity of the cascade itself.

9. Analysis of synthetic rain-rate time series

The length of the synthetic time series is set at 32 768 samples. Since the finest scale is 1 min, these time series simulate the rain-rate process over a period of three weeks, which may be considered to be the average lifetime of large-scale structures in the atmosphere. The parameters are set in accordance with the results given in section 6 (H = 0.53; α = 1.7; C1 = 0.13). The threshold is set to obtain a fractal dimension of the rain “support” equal to 0.82 [this value has been reported by Lavergnat and Golé (1998), who also used a spectropluviometer]. It may be noticed that, with these parameters and this fractal dimension, only 5% of the process is not set to 0, which corresponds to the rain time percentage observed in the experimental data (at a 1-min time scale). Figure 8 shows the effect of thresholding on the power spectrum and on the first-order structure function. The effect of rain–no rain intermittency is clearly to flatten both curves at low resolutions. Figure 9 shows that the effect of the rain–no rain intermittency is to straighten out K(q) and to cause an artificial break in the scaling of moments. When K(q) is only estimated on the basis of the scaling at the lower resolutions, this leads to 0.1 < α < 0.5 and 0.4 < C1 < 0.6, depending on the synthetic time series realization. This is consistent with the parameters found for long time series (section 7), with the parameters previously proposed in scientific literature (section 1) and with the expected effect of rain–no rain intermittency (section 3).

Fig. 8.

(top) Power spectrum and (bottom) first-order structure function of (left) a FIF generated with H = 0.53, α = 1.7, and C1 = 0.13 and (right) the same FIF thresholded such that Df = 0.82.

Fig. 8.

(top) Power spectrum and (bottom) first-order structure function of (left) a FIF generated with H = 0.53, α = 1.7, and C1 = 0.13 and (right) the same FIF thresholded such that Df = 0.82.

Fig. 9.

(top) Scaling of the moments with q = (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and (bottom) moment scaling function of (left) a FIF generated with H = 0.53, α = 1.7, C 1 = 0.13 and (right) the same FIF thresholded such that Df = 0.82.

Fig. 9.

(top) Scaling of the moments with q = (0, 0.3, 0.6, 0.9, 1.2, 1.5, and 1.8) and (bottom) moment scaling function of (left) a FIF generated with H = 0.53, α = 1.7, C 1 = 0.13 and (right) the same FIF thresholded such that Df = 0.82.

Moreover, the scale at which the break in the scaling is occurring can be roughly predicted. If there is no rain–no rain intermittency, the number of intervals during which rain is present decreases as 2−n, where n is the number of consecutive time averagings of pairs of contiguous intervals (0 corresponds to the highest resolution). If the rain support is assumed to have a fractal dimension Df , then this number decreases as 2−nDf. Therefore, 2−nDf − 2−n gives the approximate value of the number of time intervals that will be averaged with zero values at step n during the reconstitution of the cascade. If it is assumed that the effect of rain–no rain intermittency becomes the major contributor to the scaling when the number of time intervals that are averaged with zero values is equal to the theoretical number of time intervals without rain–no rain intermittency, then the scaling should break at step n such that (2−nDf − 2−n) = 2−n. This condition is actually equivalent to n = 1/Cf , which provides an estimation of the scale at which the break is expected to occur. In the case of rain, if we assume Df to be 0.8, this leads to n = 5. Since the finest time scale is 1 min, five consecutive time averagings correspond to a time scale of 25 = 32 min. It is in agreement with the breaks observed in Figs. 7 and 9 that occur between 32- and 64-min time scales.

The alternative model based on the parameters previously proposed in scientific literature (H ≈ 0; α ≈ 0.5–0.7; C1 ≈ 0.4–0.6) has also been investigated through simulations. It consists in generating a stationary bounded cascade (H = 0; α < 1) and then in setting to 0 small values that are below a given singularity. It was found that with α = 0.6 and C1=0.5, the singularity corresponding to Df = 0.8 was very small and that setting the values below it to 0 had almost no effect on the scaling. This explains why the effect of rain–no rain intermittency was generally not considered. Conversely, in the model proposed here based on high-resolution data, the effect on rain–no rain intermittency is significant because, as previously mentioned, the threshold does not correspond to a unique small singularity of the cascade but to a given intensity of the fractionally integrated cascade. In other words, a differentiated “thresholded FIF” is close to a stationary cascade multiplied by an independent fractal support. This means that all kinds of singularities may be set to 0 by the threshold, not only the smallest singularities.

10. Conclusions

The presence of zero-value intervals in analyzed time series may cause a spurious break in the scaling and a change in the universal multifractals model parameters. The multifractal analysis was, therefore, performed on an event-by-event basis, so that the analyzed process is apparently uninterrupted. This approach led to new parameter estimates (H ≈ 0.53; α ≈ 1.7; C1 ≈ 0.13) significantly different from those previously suggested in scientific literature (H ≈ 0; α ≈ 0.5–0.7; C1 ≈ 0.4–0.6). In particular, the rain process is found to be a FIF, meaning that the multifractal analysis should not be performed directly on the rain rate but on the absolute value of its gradient.

A new model has been proposed to simulate high-resolution rain-rate time series including rain-no rain intermittency based on these new parameters and on a simple threshold. The resulting time series correctly reproduce the parameters that were found for individual events as well as the break of the scaling due to rain–no rain intermittency and the fractal dimension of the rain “support.”

This study underlines the need for improved analysis techniques that would be capable of modeling the degeneracy of the multiplicative cascade due to the averaging with zero values during the reconstitution. Moreover, although the parameters proposed in this study were estimated from uninterrupted rain periods, they are also affected by rain–no rain intermittency because the process is only apparently uninterrupted—that is, the measurement is actually an average of rain periods and no rain periods. As a consequence, future analysis techniques should also take into account the degeneracy of the fractionally integrated process itself due to the measurement and relate it to the resolution of the data.

Table 3.

Same as Table 2 but for the DEVEX dataset.

Same as Table 2 but for the DEVEX dataset.
Same as Table 2 but for the DEVEX dataset.

Acknowledgments

The authors thank the French DGA and especially Thierry Marsault of CELAR for funding this study.

REFERENCES

REFERENCES
Deidda
,
R.
,
2000
:
Rainfall downscaling in a space–time multifractal framework.
Water Resour. Res.
,
36
,
1779
1794
.
Delahaye
,
J-Y.
,
L.
Barthès
,
P.
Golé
,
J.
Lavergnat
, and
J. P.
Vinson
,
2006
:
A dual-beam spectropluviometer concept.
J. Hydrol.
,
328
,
110
120
.
de Lima
,
M. I. P.
,
1998
:
Multifractals and the temporal structure of rainfall. Ph.D. thesis, Wageningen Argicultural University, 229 pp
.
de Lima
,
M. I. P.
, and
J.
Grasman
,
1999
:
Multifractal analysis of 15-min and daily rainfall from a semi-arid region in Portugal.
J. Hydrol.
,
220
,
1
11
.
Falconer
,
K.
,
2003
:
Fractal Geometry: Mathematical Foundations and Applications.
2nd ed. John Wiley and Sons, 366 pp
.
Fraedrich
,
K.
, and
C.
Larnder
,
1993
:
Scaling regimes of composite rainfall time series.
Tellus
,
45A
,
289
298
.
Gagnon
,
J-S.
,
S.
Lovejoy
, and
D.
Schertzer
,
2006
:
Multifractal earth topography.
Nonlinear Processes Geophys.
,
13
,
541
570
.
Gunn
,
R.
, and
G. D.
Kinzer
,
1949
:
The terminal velocity of fall for water droplets in stagnant air.
J. Meteor.
,
6
,
243
248
.
Harris
,
D.
,
M.
Menabde
,
A.
Seed
, and
G.
Austin
,
1996
:
Multifractal characterization of rain fields with a strong orographic influence.
J. Geophys. Res.
,
101
,
26405
26414
.
Hauser
,
D.
,
P.
Amayenc
,
B.
Nutten
, and
P.
Waldteufel
,
1984
:
A new optical instrument for simultaneous measurement of raindrop diameter and fall speed distributions.
J. Atmos. Oceanic Technol.
,
1
,
256
269
.
Hubert
,
P.
, and
J. P.
Carbonnel
,
1989
:
Fractal dimensions of rain occurrence in Soudano-Sahelian climate (in French).
Hydrol. Cont.
,
4
,
3
10
.
Hubert
,
P.
,
Y.
Tessier
,
S.
Lovejoy
,
D.
Schertzer
,
F.
Schmitt
,
P.
Ladoy
,
J. P.
Carbonnel
, and
S.
Violette
,
1993
:
Multifractals and extreme rainfall events.
Geophys. Res. Lett.
,
20
,
931
934
.
Hubert
,
P.
,
A.
Biaou
, and
D.
Schertzer
,
2002
:
De la méso-echelle à la micro-echelle: Désagrégation/agrégation multifractale et spatio-temporelle des precipitations (From meso-scale to micro-scale: Multifractal and spatio-temporal disaggregation/aggregation of precipitation).
Armines-EDF Report, 30 pp
.
Jennane
,
R.
,
R.
Harba
, and
G.
Jacquet
,
2001
:
Analysis methods for fractional Bronian motion: Thoery and comparative results (in French).
Trait. Signal
,
18
,
419
436
.
Kolmogorov
,
A. N.
,
1941
:
The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers.
Dokl. Akad. Nauk SSSR
,
30
,
299
303
.
Krajewski
,
W. F.
, and
Coauthors
,
2006
:
DEVEX—Disdrometer evaluation experiment: Basic results and implications for hydrologic studies.
Adv. Water Resour.
,
29
,
311
325
.
Ladoy
,
P.
,
F.
Schmitt
,
D.
Schertzer
, and
S.
Lovejoy
,
1993
:
The multifractal temporal variability of Nîmes rainfall data (in French).
C.R. Acad. Sci. Ser. 2
,
317
,
775
782
.
Lavallée
,
D.
,
S.
Lovejoy
,
D.
Schertzer
, and
P.
Ladoy
,
1993
:
Nonlinear variability and landscape topography: Analysis and simulation.
Fractals in Geography, N. Lam and L. De Cola, Eds., Prentice-Hall, 171–205
.
Lavergnat
,
J.
, and
P.
Golé
,
1998
:
A stochastic raindrop time distribution model.
J. Appl. Meteor.
,
37
,
805
818
.
Lilley
,
M.
,
S.
Lovejoy
,
N.
Desaulniers-Soucy
, and
D.
Schertzer
,
2006
:
Multifractal large number of drops limit in rain.
J. Hydrol.
,
328
,
20
37
.
Lovejoy
,
S.
, and
D.
Schertzer
,
2006
:
Multifractals, cloud radiances and rain.
J. Hydrol.
,
322
,
59
88
.
Lovejoy
,
S.
, and
D.
Schertzer
,
2008
:
Turbulence, rain drops and the l1/2 number density law.
New J. Phys.
,
10
,
075017
.
doi:10.1088/1367-2630/10/7/075017
.
Lovejoy
,
S.
,
D.
Schertzer
,
Y.
Tessier
, and
H.
Gaonac’h
,
2001
:
Multifractals and resolution-independent remote sensing algorithms: The example of ocean colour.
Int. J. Remote Sens.
,
22
,
1191
1234
.
Lovejoy
,
S.
,
D.
Schertzer
, and
V.
Allaire
,
2008
:
The remarkable wide range spatial scaling of TRMM precipitation.
Atmos. Res.
,
90
,
10
32
.
Mallet
,
C.
, and
L.
Barthès
,
2009
:
Estimation of gamma raindrop size distribution parameters: Statistical fluctuations and estimation errors.
J. Atmos. Oceanic Technol.
,
in press
.
Mandelbrot
,
B.
,
1974
:
Intermittent turbulence in self-similar cascades: Divergence of high moments and dimension of the carrier.
J. Fluid Mech.
,
62
,
331
350
.
Mazzarella
,
A.
,
1999
:
Multifractal dynamic rainfall processes in Italy.
Theor. Appl. Climatol.
,
63
,
73
78
.
Onof
,
C.
,
P.
Northrop
,
H. S.
Wheater
, and
V.
Isham
,
1996
:
Spatiotemporal storm structure and scaling property analysis for modeling.
J. Geophys. Res.
,
101
,
26415
26425
.
Over
,
T. M.
, and
V. K.
Gupta
,
1996
:
A space–time theory of mesoscale rainfall using random cascades.
J. Geophys. Res.
,
101
,
26319
26331
.
Parisi
,
G.
, and
U.
Frisch
,
1985
:
A multifractal model of intermittency.
Turbulence and Predictability in Geophysical Fluid Dynamics and Climate Dynamics, M. Ghil, R. Benzi, and G. Parisi, Eds., North-Holland, 84–88
.
Pathirana
,
A.
,
S.
Herath
, and
T.
Yamada
,
2004
:
Estimating rainfall distributions at high temporal resolutions using a multifractal model.
Hydrol. Earth Syst. Sci.
,
7
,
668
679
.
Pecknold
,
S.
,
S.
Lovejoy
,
D.
Schertzer
,
C.
Hooge
, and
J. F.
Malouin
,
1993
:
The simulation of universal multifractals.
Cellular Automata: Prospects in Astrophysical Applications, J. M. Perdang and A. Lejeune, Eds., World Scientific, 228–267
.
Richardson
,
L. F.
,
1922
:
Weather Prediction by Numerical Process.
Cambridge University Press, 236 pp
.
Schertzer
,
D.
, and
S.
Lovejoy
,
1987
:
Physically based rain and cloud modeling by anisotropic, multiplicative turbulent cascades.
J. Geophys. Res.
,
92
,
9692
9714
.
Schertzer
,
D.
, and
S.
Lovejoy
,
1997
:
Universal multifractals do exist!: Comments on “A statistical analysis of mesoscale rainfall as a random cascade.”.
J. Appl. Meteor.
,
36
,
1296
1303
.
Schertzer
,
D.
,
S.
Lovejoy
, and
P.
Hubert
,
2002
:
An introduction to stochastic multifractal fields.
Mathematical Problems in Environmental Science and Engineering, A. Ern and L. Weiping, Eds., Series in Contemporary Applied Mathematics, Vol. 4, Higher Education Press, 106–179
.
Schmitt
,
F.
,
S.
Vannitsem
, and
A.
Barbosa
,
1998
:
Modeling of rainfall time series using two-state renewal processes and multifractals.
J. Geophys. Res.
,
103
,
23181
23193
.
Schmitt
,
F.
,
D.
Schertzer
, and
S.
Lovejoy
,
1999
:
Multifractal analysis of foreign exchange data.
Appl. Stochastic Models Data Anal.
,
15
,
29
53
.
Tessier
,
Y.
,
S.
Lovejoy
, and
D.
Schertzer
,
1993
:
Universal multifractals in rain and clouds: Theory and observations.
J. Appl. Meteor.
,
32
,
223
250
.
Tessier
,
Y.
,
S.
Lovejoy
,
P.
Hubert
,
D.
Schertzer
, and
S.
Pecknold
,
1996
:
Multifractal analysis and modeling of rainfall and river flows and scaling, causal transfer functions.
J. Geophys. Res.
,
101
,
26427
26440
.
Veneziano
,
D.
,
R. L.
Bras
, and
J. D.
Niemann
,
1996
:
Nonlinearity and self-similarity of rainfall in time and a stochastic model.
J. Geophys. Res.
,
101
,
26371
26392
.
Yaglom
,
A. M.
,
1966
:
The influence on the fluctuation in energy dissipation on the shape of turbulent characteristics in the inertial interval.
Sov. Phys.-Dokl.
,
2
,
26
30
.

APPENDIX

FIF Simulation Technique

FIF processes are simulated following the method presented in Pecknold et al. (1993). The main steps of the algorithm are (i) generation of a noise with a Levy maximally skewed α-stable law; (ii) adequate normalization; (iii) (1 − 1/α)-order fractional integration; (iv) exponentiation; and (v) H-order fractional integration. The Mathematica code is available in Schertzer et al. (2002). Note that this code is designed to generate 2D fields and has to be slightly modified to generate 1D time series. It appeared that the time series generated with this code were not entirely satisfactory, because the scaling was broken at the highest frequencies (finite-size effect). However, averaging over 50 sample intervals was found to be sufficient to remove this break. The resulting degradation of the synthetic time series has no influence on the multifractal properties of the process, because scaling ensures that neither upscaling nor downscaling modifies the multifractal parameters. Nevertheless, particular attention must be paid to the mean absolute value of the gradient [M, defined in Eq. (7)]. Although it does not depend on the time lag Δt over which the gradient is calculated, it is, however, dependent on the scale at which the rain rate is considered. It can be shown that the following property is a consequence of Eq. (9) (Jennane et al. 2001):

 
formula

where n, k are integers, and ΔRnt) is the rain-rate gradient derived over a time lag of length Δt from the rain-rate averaged over n samples. Therefore, if a given value of M is required, the synthetic time series first has to be generated with M′ = M/50H and then averaged over 50 samples.

Another remark is that the code used here is based on the left–right symmetric “Riemann–Liouville” fractional integration, which is an unrealistic feature for temporal time series:

 
formula

The correct way to simulate temporal processes (Tessier et al. 1996) is to use causal “Liouville” fractional integrals:

 
formula

However, this does not change the conclusions of this study.

Footnotes

Corresponding author address: Louis de Montera, LATMOS; 10-12 avenue de l’Europe, 78140 Vélizy-Villacoublay, France. Email: louis.demontera@cetp.ipsl.fr