Coupled general circulation models (GCMs) of the atmosphere and ocean are the most detailed tools for studying the climate and its possible response to manmade changes in environmental conditions [10]. The analysis of the simulations has generally been limited to their lower-order statistical moments (mean and variance). However, there is growing awareness that important features of climate change are scale-dependent (e.g., spectral or scaling analyses of surface air temperature fluctuations [2,6,33]) and these are related to higher-order statistics [4].

The multifractal framework was largely developed for studying turbulent intermittency [8,17,20,26,30,31]; it has been shown to be well adapted for providing a parsimonious description of the statistics at all orders, and over a wide range of scales for rainfall [9,19,26,34] and for other intermittent geophysical fields. An overall scaling regime ranging from several minutes up to several weeks has been obtained from high-resolution time-series of rainfall [7,11,14–16,21,22,32]. Surprisingly, the multifractal formalism has rarely been applied to the analysis of fields simulated by general circulation models, although a few pioneering studies did show promising results for characterizing power-law scaling of surface air temperature fluctuations [2,6,33].

The aim of this communication is to further investigate the scaling behavior of simulated precipitation fields and its evolution with respect to climate change. A coupled atmosphere–ocean GCM (CNRM-CM3 [25]) is used to generate the precipitation fields over the past century and in a scenario of greenhouse gas increase over the 21st century (SRES A2, [10]) for investigating a possible sensitivity to greenhouse warming. For this first attempt, the analysis has been restricted to model grid points over France. The method of analysis is described in Section 2, the simulated precipitation series used in the analysis in Section 3, the results of the scale analysis in Section 4, the implications for the evolution of extremes in Section 5, and finally a synthesis and conclusion in Section 6.

The scale symmetry or ‘scaling’ of a multifractal field R means by definition that the statistical properties – e.g., probability distributions or statistical moments – of the field R at different scales are related by power laws. For instance, the scaling of the probability of the field R_λ exceeding a threshold can be defined with respect to its resolution λ = T/τ (defined as the ratio of the scale T chosen for studying the phenomenon divided by the smallest scale of observation τ) as:

The equivalence between the scaling exponents c(γ) and K(q) reduces to a Legendre transform [8]. Though, in the most general case, an infinite number of parameters may be needed to characterize these exponent functions, three parameters have been shown to be particularly important, namely:

• - the degree H of non-conservation of the field, which measures the scale dependency of the average field (H = 0 for a conservative field);
• - the mean intermittency codimension C₁, which measures how concentrated is the average field, and corresponds to the codimension of the singularity of the average field (C₁ = 0 for an homogeneous field);
• - the multifractality index α (0 ≤ α ≤ 2), which measures how fast the intermittency evolves when considering singularities slightly different from the average field singularity.

In fact, these three exponents are sufficient to statistically characterize ‘universal multifractals’ [27], which have been shown, under fairly general conditions, to be an attracting limit for multifractal cascade processes (this can be seen as a broad generalization of the central limit theorem). In this generic case, the scaling functions are analytically defined as follows:

For a given value of α, an increase of the intermittency C₁ corresponds to precipitation becoming less continuous and more sporadic in time, and thus to an increase of intensity extremes, as a given amount of precipitation is then distributed on a smaller set. Conversely, for a given value of C₁, an increase of α corresponds to a larger variation of the range of precipitation intensity, and thus also of extremes. When C₁ and α vary in the same direction, then their effects on the extremes reinforce each other. However, if they vary in opposite directions, their effects will tend to cancel out, and the resulting response of extremes will depend sensitively on which of the two variations is dominant.

A convenient and robust method for determining the parameters α and C₁ is provided by the Double Trace Moment algorithm (DTM) [13,30]. It is based on the fact that, for conservative universal multifractals, the scaling exponent K(q,η) of Rλ(η) obtained by upscaling to the resolution λ the η-th power of the original field RΛη at the highest empirical resolution Λ (RΛ(η)=RΛη) has a simple power-law relationship with the scaling exponent of R_λ:

The data used in this analysis are the grid point precipitation series simulated by the climate model CNRM-CM3 of the ‘Centre national de recherches météorologiques’ (CNRM) in a coupled simulation performed for the IPCC 4th assessment report [25]. Previous simulations with a former version of the CNRM coupled model have been documented and analyzed [5,24]. In the simulation reported here, the observed global concentrations of the main greenhouse gases (CO₂, CH₄, N₂O, and CFCs) have been specified according to their observed annual concentrations over the period 1860–2000, and then according to the scenario SRES A2 [10] for the 21st century. The accumulated daily precipitation series provided on the model grid with an approximately 2.8 × 2.8° resolution have been analyzed for 12 grid points over France (Fig. 1).

The mean annual precipitation is higher in the eastern, northern and central parts of France, while it decreases eastward (Mediterranean climate) in the southern part (Fig. 2a). While this geographical distribution is rather stable for the 20th century, there is a systematic decrease in mean precipitation during the 21st century, particularly in the South. This reduction in rainfall is associated with an increase in the number of dry days, which is not compensated by the concurrent increase of mean precipitation during wet days (not displayed). The maximum R_max of the daily rainfall (at the model scale, i.e. of the order of 250 km) over each time period tends to increase during the 21st century, particularly in the South (Fig. 2b).

The DTM technique as described above was applied separately on the daily precipitation series at each one of the 12 grid points over different time spans (50, 25 and 10 years in length). As a preliminary step, it has been checked with the help of spectral analyses that the spectra E(f) display a scaling behavior as a function of frequency f, which corresponds to E(f) ≈ f^–β, and appears as a linear slope when in a log–log plot (Fig. 3). Such scaling behavior is found over the range of two days to about 16 days for most points with a spectral slope β < 1, and this confirms that DTM can be directly applied to those points. However, the points in the southern part of the domain tend to have a smaller scaling range and spectral slopes β close to 1, especially for point Pt32, which has a slope slightly above 1. The results shown in Fig. 3 are for the time span 1950–1999; we have checked that the fluctuations of the spectral slopes over the other time spans are usually smaller than the variation of the slopes between the different points (Fig. 4). A spectral slope β = 1 corresponds to so-called ‘1/f noise’, which has been found to be representative of tropical convective variability [36]. Though the application of the DTM for these points is less reliable, we have kept them in the analysis for comparison with the other points. We will discuss below that a value β ≈ 0.8 corresponds to a conservative rain field (H ≈ 0), as often hypothesized and/or empirically observed.

The first step in the DTM method is illustrated in Fig. 5: a log–log plot of the Double Trace Moment (of order q = 1.5 and power η = 0.93) vs. its time resolution (here shown for grid point Pt13) confirms that the scaling holds from one day up to 16 days. This scaling range is rather in good agreement with what has been found in most of the analyses of observed rainfall series [14]. In particular, Tessier et al. [35] found a change of slope in rainfall spectra at about 14 days. Such a time span is roughly consistent with the characteristic duration of weather regimes; it is roughly the ‘eddy turn-over time’ (lifetime) of planetary-scale structures. Using a straightforward linear regression, the slope of the fitted line yields an estimate of K(q,η) with an excellent accuracy (the smallest coefficient of determination R² is > 0.97).

To study the non-stationarity and temporal evolution of the multifractal parameters, the original 240-year series was split into: 5 subseries of 50 years (actually only 40 years before 1900), 10 subseries of 25 years (only 20 years before 1900), and 24 subseries of 10 years. The DTM technique was applied to each subseries at each grid point. Since similar results have been found for the shorter subseries, only the results for the 50-year series are presented here for clarity (and lack of space). Fig. 6a and b display the α and C₁ estimates respectively for the different time slices. It can be seen that they tend to vary in opposite directions between the different grid points. The geographical patterns remains rather stable, but a systematic evolution can be seen between the successive time spans.

The parameter estimates from these simulations are systematically different from those obtained from analysis of local precipitation measured at observing stations, i.e. α ∈ 0.5–0.8 and C₁ ≈ 0.4 [1,3,12,34], but are closer to the parameters obtained from observed daily river flows for large basins, e.g. α ≈ 1.45, C₁ ≈ 0.2 ([33]) or α ≈ 1.7, C₁ ≈ 0.12 ([23]). This seems consistent with the fact that precipitation at a model grid point represents a quantity integrated over a very large surface (∼70,000 km² for the current model mesh) more similar to that of a river basin than to that of a rain gauge. Also relevant are the recent accurate estimates of the spatial parameters (from 4.3 to 20,000 km) from satellite radar reflectivity of precipitation: α ≈ 1.5, C₁ = 0.63 ± 0.02, H = 0.00 ± 0.01, which correspond roughly to α ≈ 1.5, C₁ = 0.3, H = 0.00 ± 0.01 for the rain rate [16]. Another difference with the observed characteristics of precipitation series can come from the fact that the parameterizations of precipitation processes (large-scale condensation, convection) that operate at the scale of the GCM grid are unable to reproduce the high spatial and temporal variability of the real processes.

Comparison of the different periods shows a systematic evolution of the parameters with time, with a decrease of the multifractality index α during the simulation, and on the contrary an increase of the mean intermittency codimension C₁. The trends in the multifractal parameters seem to follow those of the greenhouse gas concentrations used to force the coupled model during the simulation, and of the resulting temperature and precipitation response, with a small change during the 20th century and an accelerating trend during the 21st century. The increase of the fractal codimension C₁ (decrease of the fractal dimension 1–C₁) of rainfall occurrences means a higher frequency of dry periods at different scales. It therefore quantitatively confirms and generalizes the increase of the number of dry days over the 21st century that is actually found in this simulation.

We also checked that, independently of the time evolution of both α and C₁, the resulting moment scaling function K(2) remains rather stable around the value K(2) ≈ 0.2, confirming that the rainfall remains conservative (H ≈ 0) due to the relationship: H = (β – 1 + K(2))/2 and the estimate obtained above of the spectral slope β ≈ 0.8 for the points in the middle of the domain. The points in the northern part have slightly negative values (−0.2 < H < 0), while those in the southern part have positive values (0 < H < 0.2).

The opposite trends noted in the temporal evolution of C₁ (increase) and α (decrease) are rather intriguing, and might explain the difficulties encountered in assessments of the response of rainfall extremes to global change. Indeed, they both measure different aspects of the variability of the field. As discussed above, the mean intermittency codimension C₁ measures the fragmentation of the wet periods into smaller and smaller events, therefore producing more and more intense events if the total large-scale amount does not vary too much. The multifractality index α measures the range of variation of precipitation amounts. Therefore, physically speaking, extremes increase with both C₁ and α. Indeed, it is easy to check that the scaling moment function K(q) (Eq. (3)) is strictly increasing with both C₁ and α. Therefore, the opposite trends of C₁ and α actually lead to a partial mutual compensation for the time evolution of high moments and rainfall extremes, but, as it will be shown below, the C₁ trend is ultimately dominant.

The notion of the maximal probable singularity γ_s that can be observed on a unique sample [4,11,28] can be useful to refine the analysis for the extreme fluctuations. This singularity is defined by the maximum value of the singularity γ_s for which the fractal co-dimension of the support of the field c(γ_s) as defined in Eq. (1) becomes equal to the dimension d of the embedding space (d = 1 for time series as analyzed here):

The main interest of this notion is that it is scale independent, whereas the corresponding rainfall is sensitively scale dependent, since γ_s is its scale exponent. It could be used to analytically infer the maximal probable rainfall R_s(λ) downscaled from a rainfall R at a large scale T (for example, one month) down to a much smaller scale τ = T λ (for example, one day) if both belong to the same scaling range. Indeed, γ_s yields the following downscaling rainfall ratio:

Combining Eqs. (4) and (6) yield the following analytical expression of γ_s for universal multifractals:

Inserting the values of α and C₁ obtained with the DTM into Eq. (8) yields the corresponding values of the maximum probable singularity γ_s displayed in Fig. 7. We notice that this singularity has only a very small variation during the 20th century, but increases during the 21st century, especially for the points in the southern part of France. With the help of Eq. (7), this evolution generalizes at various scales the evolution of the maximum daily precipitation at the model grid scale displayed in Fig. 2b. Actually, if we somewhat extrapolate the scaling relation (7) outside of the empirical scaling range, the maximum probable singularity γ_s should be proportional to the logarithm of the ratio R_max/R_mean (the maximum daily rainfall R_max to the mean rainfall R_mean). The relationship found between these two quantities, displayed in Fig. 8, confirms empirically the relevance of the maximum probable singularity γ_s to infer extremes of the rainfall distribution. One of its main interests is that, since γ_s is a combination of the two parameters characterizing universal multifractals, it makes it possible to decompose the evolution in the extreme behavior in terms of the contributions due to changes in intermittency (C₁) and in multifractality (α). The climate scenarios show that the increase in the extremes in precipitation, such as the maximum daily rainfall, is mainly due to the increase in the intermittency, which is only partly compensated by the decrease in the multifractality. An illustration of this is provided in Fig. 9, where the relative changes of all these parameters computed over successive 10-year time spans are shown as a time series for the point Pt32 (similar results are found for the other points, but we have chosen to show this point since it is the one for which DTM analysis would be less reliable because of the smaller scaling range and the spectral slope close to 1). It can be seen that the multifractal parameters evolve quite smoothly in time, contrary to R_max. The evolution of K(2) is very close to that of C₁, while that of γ_s is somewhat influenced by that of α, but nevertheless follows the C₁ trend. The maximum daily rainfall R_max is less stable and has more sampling fluctuations than the multifractal parameters, since it is actually based on the single maximum value of the daily precipitation over each 10 year span, but its long-term evolution tends to follow that of γ_s, although its trend is less obvious.

The standard way of characterizing climate variability including the extremes is to concentrate on a single – usually the highest available – scale/resolution. However, climate signals are generally variable over a wide range of scales, so that the available resolution has no special climatological significance, and a resolution-independent characterization is more physically relevant. Here, we study the evolution of multiscale variability of precipitation in long-term climate time-series showing that a scaling characterization of the maxima – the highest order singularity present in the sample – can be used in the chosen climate scenario to demonstrate a contrasted but effective increase in precipitation extremes in France over the next century. Other multifractal characterizations were used to provide complementary information with respect to the more usual statistical quantities, and have also the additional advantage of being based on an extensive theoretical basis, combining fractal geometry, statistical physics, nonlinear processes which has been developed and tested in different applications ranging from turbulence to large-scale geophysical fields. The multifractal parameters determined by the DTM technique appear rather robust as they are determined by fitting over the entire scaling range, and are thus more stable than empirical extremes determined from only a few isolated values in the time-series. This study, limited to a small geographical region, has shown promising results of the multifractal formalism as applied to the analysis of rainfall in climate change simulations. It will be particularly interesting to extend it further to different types of climate in other regions and to simulations with different models.