Research article 19 Mar 2020
Research article  19 Mar 2020
Approximate multifractal correlation and products of universal multifractal fields, with application to rainfall data
 Hydrologie Météorologie et Complexité, Ecole des Ponts ParisTech, ChampssurMarne, France
 Hydrologie Météorologie et Complexité, Ecole des Ponts ParisTech, ChampssurMarne, France
Correspondence: Auguste Gires (auguste.gires@enpc.fr)
Hide author detailsCorrespondence: Auguste Gires (auguste.gires@enpc.fr)
Universal multifractals (UMs) have been widely used to simulate and characterize, with the help of only two physically meaningful parameters, geophysical fields that are extremely variable across a wide range of scales. Such a framework relies on the assumption that the underlying field is generated through a multiplicative cascade process. Derived analysis techniques have been extended to study correlations between two fields not only at a single scale and for a single statistical moment as with the covariance, but across scales and for all moments. Such a framework of joint multifractal analysis is used here as a starting point to develop and test an approach enabling correlations between UM fields to be analysed and approximately simulated.
First, the behaviour of two fields consisting of renormalized multiplicative power law combinations of two UM fields is studied. It appears that in the general case the resulting fields can be well approximated by UM fields with known parameters. Limits of this approximation will be quantified and discussed. Techniques to retrieve the UM parameters of the underlying fields as well as the exponents of the combination have been developed and successfully tested on numerical simulations. In a second step tentative correlation indicators are suggested.
Finally the suggested approach is implemented to study correlation across scales of detailed rainfall data collected with the help of disdrometers of the Fresnel platform of Ecole des Ponts ParisTech (see available data at https://hmco.enpc.fr/portfolioarchive/taranisobservatory/, last access: 12 March 2020). More precisely, four quantities are used: the rain rate (R), the liquid water content (LWC) and the total drop concentration (N_{t}) along with the mass weighed diameter (D_{m}), which are commonly used to characterize the drop size distribution. Correlations across scales are quantified. Their relative strength (very strong between R and LWC, strong between DSD features and R or LWC, almost null between N_{t} and D_{m}) is discussed.
Numerous geophysical fields exhibit intermittent features with sharp fluctuations across all scales, skewed probability distribution and longrange correlations. A common framework to analyse and simulate such fields is multifractals. The underlying idea of this framework is that these fields are the result of an underlying multiplicative cascade process. It is physically based in the sense that it is assumed the fields inherit the scaleinvariant properties of the governing Navier–Stokes equations and hence should exhibit scale invariant features as well. The reader is referred to the reviews by Schertzer and Lovejoy (2011) and Schertzer and Tchiguirinskaia (2017) for more details. In the large class of universal multifractals (UMs), which are the stable and attractive limits of nonlinearly interacting multifractal processes and correspond to a broad, multiplicative generalization of the central limit theorem (Schertzer and Lovejoy, 1987, 1997), a conservative field is fully described with the help of only two parameters with a physical interpretation. The UM framework was initially developed to address wind fluctuations and has also been implemented on numerous other geophysical fields ranging from rainfall, discharge, temperature or humidity to soil properties and phytoplankton concentration, for example.
Much less work has been devoted to the analysis of the correlations and couplings between two fields exhibiting multifractal properties. A framework was originally presented by Meneveau et al. (1990), who suggested studying the properties of joint moments of two multifractal fields (i.e. the product of the two fields raised to two different powers) across scales. The behaviour of the scaling exponent as a function of the two moments provides information on the correlations between the two fields. They tested their framework on velocity and temperature as well as velocity and vorticity. Such a framework has been implemented in many other contexts. Bertol et al. (2017) used it to extract information on the tillage technique by joint analysis of water and soil losses. Siqueira et al. (2018) studied the correlations between soil properties (e.g. pH, organic carbon, exchangeable cations and acidity) and altitude. Wang et al. (2011) focused on joint properties of soil water retention parameters and soil texture, while JiménezHornero et al. (2011) focused on the links between wind patterns and surface temperature. Xie et al. (2015) used this framework in a nongeophysical domain to better understand the crosscorrelation between stock market indexes and the index of volatilities.
Seuront and Schmitt (2005a, 2005b) suggested a refinement of this framework and introduced a renormalization of these joint moments to define an exponent called “generalized correlation function” and used the properties of this function to better understand the coupling between fluorescence (which is related to phytoplankton concentration) and temperature for various levels of turbulence. A similar formalism is used by Calif and Schmitt (2014) to study the coupling between wind fluctuations and the aggregate power output from a wind farm. The generalized correlation function is found to be symmetrical with regard to the chosen moments for the two studied fields, suggesting a simple relation of proportionality between the two quantities.
Actually the previously discussed frameworks have only been implemented for lognormal cascades, for which computations basically boil down to a single parameter and correlation functions are represented by linear ones. Furthermore only two specific cases have been primarily studied, either a proportional or a power law relation between the two studied fields. In this paper, we suggest relying on this theoretical framework and extending its use to UMs and to relations between fields consisting of multiplicative power law combinations.
In Sect. 2, the theoretical framework of UM and joint multifractal analysis is presented. Its theoretical consequences on the analysis of multiplicative power law combinations of UM fields are explored in Sect. 3. Numerical simulations are used to confirm the validity of the suggested analysis techniques. A new indicator of correlation is presented in Sect. 4 and its limitations discussed. Finally the framework is implemented on rainfall data to study the correlation between rain rate, liquid water content and quantities characterizing the drop size distribution.
2.1 Universal multifractals
The goal is to represent the behaviour of a field ϵ_{λ} across scales. The resolution λ is defined as the ratio between the outer scale L (i.e. the duration or size of studied event) and the observation scale l ($\mathit{\lambda}=L/l$). In practice, the field at resolution λ is computed by averaging over adjacent time steps or pixels of the field measured or simulated at a maximum resolution (λ_{max}). Multifractal fields exhibit a power law relation between their statistical moment of order q and the resolution λ:
where K(q) is the scaling moment function that fully characterizes the variability across scales of the field. UMs are a specific case towards which multiplicative cascade processes converge (Schertzer and Lovejoy, 1987, 1997). Only two parameters with physical interpretation are needed to define K(q) for conservative fields:

C_{1}, the mean intermittency codimension, which measures the clustering of the (average) intensity at smaller and smaller scales (C_{1}=0 for a homogeneous field);

α, the multifractality index ($\mathrm{0}\le \mathit{\alpha}\le \mathrm{2}$), which measures the clustering variability with regard to the intensity level.
For UM, we have
K(q) is computed through trace moment (TM) analysis which basically consists of plotting Eq. (1) in log–log and estimating the slope of the retrieved straight line. Double trace moment (DTM) analysis, specifically designed for UM fields, is commonly used to estimate UM parameters (Lavallée et al., 1993). One can also note that UM parameters characterize the first and second derivatives of K(q) near q=1:
When doing a multifractal analysis, one should keep in mind that such fields can be affected by multifractal phase transitions (Schertzer and Lovejoy, 1992). One is associated with sampling limitations. It results from the fact that due to the limited size of studied samples, estimates of statistical moments greater than a given moment q_{s} are not reliable (see Hubert et al., 1993, and Douglas and Barros, 2003, for some examples of implementation). In practice, the empirical curve of K(q) will become linear from q_{s} and hence depart (being below) from the theoretical curve. The second one is trickier and associated with the divergence of moments (Schertzer and Lovejoy, 1987). The issue was also mentioned in Mandelbrot (1974) and Kahane (1985) but they did not address the quantification of the spurious statistical estimates on finite samples and their dependence on their size (Schertzer and Lovejoy, 1992). This is due to the fact the field generated by a cascade process can become so concentrated that its average over a given area can diverge. This results in $K\left(q\right)\approx +\mathrm{\infty}$ for q>q_{D}. In practice the K(q) will obviously be computed but its value will be an overestimation of the theoretical K(q) (hence it will be greater).
2.2 Joint multifractal analysis
Let us consider two fields, ϕ_{λ} and ϵ_{λ}, that exhibit multifractal properties. In order to study the correlation across scales Seuront and Schmitt (2005a) refined the initial framework of Meneveau et al. (1990) and suggested performing a joint multifractal analysis as follows:
where r(q,h) is a “generalized correlation exponent”. If ϕ_{λ} and ϵ_{λ} are lognormal multifractal processes (i.e. α=2), then r(q,h) is linear with regard to both h and q. $r(q,h)=\mathrm{0}$ for independent fields. If they are power law combinations related with ${\mathit{\varphi}}_{\mathit{\lambda}}=c{\mathit{\u03f5}}_{\mathit{\lambda}}^{d}$, then r(q,h) is symmetric in the dp–q plane.
Let us consider two independent UM fields X_{λ} and Y_{λ}, with their respective characteristic parameters α_{X}, C_{1,X}, α_{Y} and C_{1,Y}. The goal of this section is to understand the behaviour of a field ϵ_{λ} consisting of renormalized multiplicative power law combinations of X_{λ} and Y_{λ}. ϵ_{λ} is then defined by
where a and b are exponents characterizing the relative weight of X_{λ} and Y_{λ} in the combination.
3.1 Intuitive understanding of a and b
Let us first discuss intuitively the influence of the parameters a and b. Figure 1 displays the fields ϵ_{λ} (in red) and X_{λ} (in blue) for a realization of X_{λ} and Y_{λ} with α_{X}=1.8, ${C}_{\mathrm{1},X}=\mathrm{0.3}$, α_{Y}=0.8 and ${C}_{\mathrm{1},Y}=\mathrm{0.3}$ (Eq. 5 is used). Values of a ranging from 1 to 0 are shown. b was tuned to ensure the same C_{1} is retrieved on all the fields. For a=1 and b=0 (upper left), the two fields are obviously equal and hence superposed. The opposite case is a=0 and b=1 (lower right), for which ϵ_{λ} is simply equal to Y_{λ}, and hence fully independent of X_{λ}. In the intermediate cases, the progressive decorrelation between the two fields is visible with decreasing values of a. In that sense the parameters a and b characterize the level of correlation between the two fields.
3.2 Theoretical expectations
In order to evaluate the expected multifractal behaviour of ϵ_{λ}, its statistical moments of order q are computed to evaluate K_{ϵ}(q). Given that X_{λ} and Y_{λ} are independent, it yields
which means we have
The exact computation of K_{ϵ}(q) is written in the second line of Eq. (7). The third line is not exact and corresponds the form K_{ϵ}(q) would have if ϵ_{λ} was actually a UM. This is not true in the general case. In order to assess pseudoUM parameters C_{1,ϵ} and α_{ϵ}, we suggest to use the properties of Eq. (3) and equalize the first and second derivatives of the two last lines of Eq. (7) for q=1. This yields
It should be noted that in the specific case of α_{X}=α_{Y}, α_{ϵ} is also equal to this value and ϵ_{λ} is actually an exact UM field.
Figure 2 displays the scaling moment functions of the previously discussed fields for various sets of parameters. Similar results are found for other sets of UM parameters and combinations of a and b exponents. In Fig. 2a, the same α is used for both X_{λ} and Y_{λ}, and the expected exact UM behaviour is correctly retrieved. When α_{X}≠α_{Y}, ϵ_{λ} is not exactly a UM. As it is illustrated in Fig. 2b and c, the smaller the differences, the better the UM approximation for ϵ_{λ}. In the extreme case when α_{Y}=0 (Fig. 2c), the approximation remains valid only for q ranging from ∼0.6 to 1.6. This range is much wider when the α values are closer. It should be noted that for great moments, some discrepancies are visible, with the exact value of K_{ϵ}(q) always being greater than its UM approximation. This could wrongly be interpreted suggesting that a multifractal phase transition associated with the divergence of moments is occurring, whereas it is merely an illustration of the limits of validity of the approximation of ϵ_{λ} as a UM field. Indeed, the values of q_{D} are much greater than the moment for which the discrepancies start to be visible. In the cases of Fig. 2, we have q_{D}=5.96 for panel (a), q_{D}=4.58 for panel (b) and q_{D}=119 for panel (c), for which the approximation as a UM field is less valid. These values are obtained by looking for the solution >1 to the equation $K\left({q}_{D}\right)=({q}_{D}\mathrm{1})D$ using the pseudoUM parameters of ϵ_{λ} (D is the dimension of the embedding space and is equal to 1 for time series). When confronted with such behaviour, keeping in mind this sort of interpretation could be interesting.
3.3 Techniques for retrieving parameters
In this subsection an empirical technique to estimate the UM parameters of Y_{λ} and the exponents a and b from a joint multifractal analysis of X_{λ} and ϵ_{λ} is presented. The following steps should be implemented:

Step 1: performing a UM analysis of each field X_{λ} and ϵ_{λ} independently. This enables the quality of the scaling behaviour to be confirmed and α_{X}, C_{1,X}, α_{ϵ} and C_{1,ϵ} to be estimated. Without any loss of generality, we can assume that ${C}_{\mathrm{1},Y}={C}_{\mathrm{1},X}$. Indeed C_{1,Y} is a rather arbitrary quantity that can be changed while the one that actually matters is ${C}_{\mathrm{1},Y}{b}^{{\mathit{\alpha}}_{Y}}$.

Step 2: estimating a. This is actually the trickiest portion of the process and requires a joint multifractal analysis. More precisely Eq. (4) is implemented with X_{λ} and ϵ_{λ}. In that case, it turns out that the ratio does not depend any more on Y_{λ} but only on X_{λ}. One obtains
$$\begin{array}{ll}{\displaystyle}r(q,h)& {\displaystyle}={K}_{X}(ha+q)K\left(ha\right)K\left(q\right)\\ \text{(9)}& {\displaystyle}& {\displaystyle}={\displaystyle \frac{{C}_{\mathrm{1},X}}{{\mathit{\alpha}}_{X}\mathrm{1}}}\left(\right(ha+q{)}^{{\mathit{\alpha}}_{X}}(ha{)}^{{\mathit{\alpha}}_{X}}(q{)}^{{\mathit{\alpha}}_{X}}).\end{array}$$Hence, for a given value of h and q, r(q,h) is an increasing function of a. This property is used to compute an estimate of a. The simplest approach is to set h and q, compute an empirical value of r_{emp}(q,h), and find the a that yields this value. When implementing this technique, one should keep in mind that empirical fields are subject to multifractal phase transitions affecting their scaling behaviour. This means that ha+q, ha and q should remain within the range of values for which the estimations of the scaling moment functions remain reliable, i.e. smaller that the corresponding q_{s} and q_{D}.

Step 3: estimating α_{Y}. Using Eq. (8), one can easily obtain the following (noting that ${\mathit{\alpha}}_{\mathit{\u03f5}}{C}_{\mathrm{1},\mathit{\u03f5}}={C}_{\mathrm{1},X}{a}^{{\mathit{\alpha}}_{X}}{\mathit{\alpha}}_{X}+{C}_{\mathrm{1},Y}{b}^{{\mathit{\alpha}}_{Y}}{\mathit{\alpha}}_{Y}$, and that the term ${C}_{\mathrm{1},Y}{b}^{{\mathit{\alpha}}_{Y}}$ is simply equal to ${C}_{\mathrm{1},\mathit{\u03f5}}{C}_{\mathrm{1},X}{a}^{{\mathit{\alpha}}_{X}}$, which enables the nonlinear part of the equation to be removed):
$$\begin{array}{}\text{(10)}& {\mathit{\alpha}}_{Y}={\displaystyle \frac{\frac{{C}_{\mathrm{1},\mathit{\u03f5}}}{{C}_{\mathrm{1},X}}{\mathit{\alpha}}_{\mathit{\u03f5}}{a}^{{\mathit{\alpha}}_{X}}{\mathit{\alpha}}_{X}}{\frac{{C}_{\mathrm{1},\mathit{\u03f5}}}{{C}_{\mathrm{1},X}}{a}^{{\mathit{\alpha}}_{X}}}}.\end{array}$$ 
Step 4: computing b. Once α_{Y} is known, Eq. (8) (top) can be used to estimate b as follows (noting that ${C}_{\mathrm{1},Y}{b}^{{\mathit{\alpha}}_{Y}}={C}_{\mathrm{1},\mathit{\u03f5}}{C}_{\mathrm{1},X}{a}^{{\mathit{\alpha}}_{X}}$ and that we have ${C}_{\mathrm{1},Y}={C}_{\mathrm{1},X}$):
$$\begin{array}{}\text{(11)}& b={\left({\displaystyle \frac{{C}_{\mathrm{1},\mathit{\u03f5}}}{{C}_{\mathrm{1},X}}}{a}^{{\mathit{\alpha}}_{X}}\right)}^{\mathrm{1}/{\mathit{\alpha}}_{Y}}.\end{array}$$
3.4 Implementation on numerical simulations (discrete UM)
The approach presented above is tested on numerical simulations obtained with discrete inscale cascades. It consists of iteratively repeating a cascade step with a noninfinitesimal scale ratio in which a “parent” structure is divided into “daughter” structures whose affected value is the one of the “parent” structure multiplied by a random factor, ensuring that Eqs. (1) and (2) remain valid. Such a simple field generation process is sufficient for the purposes of this paper. The recent introduction of multifractal operators and vectors paves the way for physically based, continuous (inscale) multivariate analysis of multifractal fields or measures (Schertzer and Tchiguirinskaia, 2015, 2020).
A set of 10 000 realizations of 512 long 1D discrete cascades is used, and analyses are carried out on ensemble averages.
Before starting, let us clarify the objective of this section. X_{λ} and Y_{λ} are first simulated and then ϵ_{λ} is built with some values of a and b. The purpose is to retrieve the values of a, b and α_{Y} by simply analysing X_{λ} and ϵ_{λ}, which are assumed to be known.
The parameters used for these simulations are α_{X}=1.8, ${C}_{\mathrm{1},X}=\mathrm{0.3}$, α_{Y}=0.8, ${C}_{\mathrm{1},Y}=\mathrm{0.3}$, a=0.6 and b=0.2. As a consequence we expect to find α_{ϵ}=1.39 and ${C}_{\mathrm{1},\mathit{\u03f5}}=\mathrm{0.20}$. Other sets of parameters have been tested and yield similar results.
Results of this analysis are displayed in Fig. 3. As expected, the scaling behaviour observed on both X_{λ} and ϵ_{λ} is excellent. TM analysis, i.e. Eq. (1) in a log–log plot, for ϵ_{λ} is shown in Fig. 3a and all the coefficients of determination of the straight lines used to compute K(q) are greater than 0.99. With regard to the estimates of UM parameters retrieved via the DTM technique, for X_{λ} they are equal to 1.79 and 0.27 for α and C_{1} respectively, which is close to the values input in the simulations. The small discrepancy in C_{1} has already been noted with such discrete simulations. The respective estimates for ϵ_{λ} are 1.35 and 0.18, which are in agreement with the theoretical expectations. These small differences are visible in Fig. 3b, which displays the empirical and theoretical fitting of K(q). For X_{λ}, it can be noted that the empirical estimate of K(q) is smaller that its theoretical value (using UM estimates retrieved from the DTM analysis) for q greater than ∼1.7. This is consistent with a behaviour affected by the multifractal phase transition associated with sampling limitations (q_{s}=1.95 for the input UM parameters). It can be noted that for ϵ_{λ} we have a greater q_{s} equal to 1.95, while it is even greater for Y_{λ}(=4.5). The values of q_{D} are greater in all cases, meaning that the multifractal phase transition associated with divergence of moment will not bias our analysis.
In order to estimate a (step 2 of the process described in the previous subsection), we consider the two moments $q=h=\mathrm{0.7}$. Note that with these values we have $ha+q=\mathrm{1.12}$, which is much smaller than the minimum q_{s} for the chosen values of UM parameters. This means that the estimates should not be affected by expected biases associated with multifractal phase transitions. Figure 3c shows the output of the joint multifractal (Eq. 4 in log–log plot). It appears that the scaling is excellent and the slope gives an estimate of r(0.7,0.7). It is then used to estimate a by adjusting the value of a so that r(0.7,0.7)(a) equals the computed empirical value (Fig. 3d). This yields a=0.59. Finally (Eqs. 10 and 11) we obtain an estimate of b equal to 0.20 and an estimate of α_{Y} equal to 0.77. These values are very close to those input in the simulations. In summary, there is a very good agreement between theoretical expectations and numerical simulations, which confirms the validity of the framework presented in this section.
Finally, let us discuss the uncertainties in the estimates of a. Fig. 4 displays the estimates of a on the simulated fields (see Fig. 3) as a function of the moment orders q and h used in the joint multifractal analysis. It appears that as long as the studied moments remain within the range of reliability of the multifractal analysis (i.e. $ha+q<{q}_{s}$ as previously discussed), the estimates are rather stable. For greater values, there is an underestimation of a.
Let us consider two fields ϵ_{λ} and ϕ_{λ}. It is assumed that they both exhibit UM properties, with known UM parameters. The purpose of this section is to present a framework to study the correlations across scales between the two fields. This relies on the joint multifractal analysis presented in Sect. 2.1, with the suggestion of a simplified indicator. It furthermore opens the path to numerical simulations of one field from the other.
More precisely, the consequences of describing each field as a multiplicative power law combination of the other and an independent one will be explored. The notations are
where a, b, a^{′} and b^{′} characterize the level of correlation between the two fields, while Y_{λ} and Z_{λ} are independent random UM fields. As shown in the previous section, without any loss of generality it can be assumed that ${C}_{\mathrm{1},Y}={C}_{\mathrm{1},\mathit{\varphi}}$ and ${C}_{\mathrm{1},Z}={C}_{\mathrm{1},\mathit{\u03f5}}$. This enables the following calculations to be simplified.
4.1 Limitations of this symmetric framework
If both lines of Eq. (12) were to be correct, then the joint multifractal correlation of ϵ_{λ} and ϕ_{λ} could be computed in two equivalent ways:
leading to
In the general case, Eq. (14) is not valid for any q and h. To better understand this, let us consider a given level of correlation by setting the parameters a and b. The goal is to compute a^{′} and b^{′} from the available parameters. The lefthand side of Eq. (14) is known, and after setting given values of q and h it is possible to implement the same process as in Sect. 3.3 to determine a^{′}, b^{′} and α_{Z}. Figure 5 displays the outcome of this analysis, according to the values of h and q used, for a=0.2 in the case α_{ϵ}=0.8, ${C}_{\mathrm{1},\mathit{\u03f5}}=\mathrm{0.4}$, α_{ϕ}=0.8, ${C}_{\mathrm{1},\mathit{\varphi}}=\mathrm{0.2}$ (meaning that b=0.30 and α_{Y}=0.68). As can be seen, the estimates of a^{′} exhibit a dependency on q and h. The dependency is stronger on q than on h and estimates remain rather stable as long as q<0.8. Both sides of Eq. (14) are plotted in Fig. 6 for this set of UM parameters with a=0.2 and estimates of ${a}^{\prime}=\mathrm{0.19}$ obtained with $h=q=\mathrm{0.7}$. Expected differences are visible for the larger values of h and q. It should be mentioned that these results are presented for a bad case with strong differences between α_{ϵ} and α_{ϕ}. They are actually much smaller if both values are closer to 2. For the specific case, ${\mathit{\alpha}}_{\mathit{\u03f5}}={\mathit{\alpha}}_{\mathit{\varphi}}=\mathrm{2}$, Eq. (14) becomes
meaning that once hq has been removed, a^{′} is deterministically obtained once a is set and ${r}_{\mathit{\u03f5}\mathit{\varphi}}(q,h)={r}_{\mathit{\u03f5}\mathit{\varphi}}(q,h)={C}_{\mathrm{1},\mathit{\varphi}}ahq$ is linear with regard to h and q.
Figure 7 illustrates the relation between the parameters retrieved by setting different values of a in the same case α_{ϵ}=0.8, ${C}_{\mathrm{1},\mathit{\u03f5}}=\mathrm{0.4}$, α_{ϕ}=0.8, ${C}_{\mathrm{1},\mathit{\varphi}}=\mathrm{0.2}$. First it should be mentioned that for a given set of UM parameters, not all values of a are possible. Indeed the inequality $\mathrm{0}\le {\mathit{\alpha}}_{Y}\le \mathrm{2}$ must be respected, leading to $a\le min\left[\right(\frac{{C}_{\mathrm{1},\mathit{\u03f5}}{\mathit{\alpha}}_{\mathit{\u03f5}}}{{C}_{\mathrm{1},\mathit{\varphi}}{\mathit{\alpha}}_{\mathit{\varphi}}}{)}^{\mathrm{1}/{\mathit{\alpha}}_{\mathit{\varphi}}},\left(\frac{{C}_{\mathrm{1},\mathit{\u03f5}}(\mathrm{2}{\mathit{\alpha}}_{\mathit{\u03f5}})}{{C}_{\mathrm{1},\mathit{\varphi}}(\mathrm{2}{\mathit{\alpha}}_{\mathit{\varphi}})}{)}^{\mathrm{1}/{\mathit{\alpha}}_{\mathit{\varphi}}}\right]$. In this case we must have a≤0.43. We retrieved the expected behaviour and are able to quantify it: b decreases with increasing a (Fig. 7a), a^{′} increases with increasing a (Fig. 7b), α_{Y} decreases with increasing a (Fig. 7c), and similar behaviour is found in terms of dependency in a^{′} for the symmetric case.
4.2 A simplified indicator
In Sect. 4.1, limitations of this fully symmetric framework are highlighted. However, it is possible to suggest a rather intuitive indicator enabling most of the information obtained from the joint multifractal correlation analysis (i.e. the computation of r(q,h)) to be extracted. This corresponds to the portion of intermittency C_{1} of one field explained by the other:
Both “indicators of correlation” (ICs) are displayed in Fig. 8 for the data corresponding to Fig. 7. Both curves are close, and this symmetric behaviour is what is wanted for such an indicator of correlation. Again, much closer curves are obtained with greater values of α and identical curves are retrieved when the α of ϵ and ϕ are both equal to 2. Forcing IC_{ϵϕ}=IC_{ϕϵ} can actually be a way to find an estimate of a^{′} once a is known without having to implement the process described above. It yields
Equation (17) is actually plotted in a dashed line in Fig. 7b and provides very good estimates. Hence this IC appears to be a good candidate for characterizing in a simple way the correlations across scales between two fields. One should keep in mind that it is mainly relevant in the case where the studied fields do not exhibit values of α that are too small (typically <0.8).
5.1 Presentation of the data
The rainfall data used in this paper were collected by a OTT Parsivel^{2} disdrometer (Battaglia et al., 2010; OTT, 2014) located on the roof of the Carnot building of the Ecole des Ponts ParisTech campus near Paris between 15 January 2018 and 9 December 2018. It is part of the TARANIS observatory of the Fresnel Platform of École des Ponts ParisTech (https://hmco.enpc.fr/portfolioarchive/fresnelplatform/, last access: 12 March 2020). Data are collected with 30 s time steps. Data will only be briefly presented in this paper and interested readers are referred to Gires et al. (2018b), who discusses available database in detail along with some data samples for a similar measurement campaign.
In this paper four quantities are studied:

R, the rain rate (mm h^{−1});

LWC, the liquid water content (g m^{−3});

N_{t}, the total drop concentration (m^{−3});

D_{m}, the mass weight diameter (mm).
N_{t} and D_{m} are used to characterize the drop size distribution (DSD, N(D), ${\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{mm}}^{\mathrm{1}}$) of the rainfall. N(D)dD is the number of drops per unit volume (in m^{−3}) with an equivolumic diameter between D and D+dD (in mm). DSDs are commonly written in the form N(D)=N_{t}f(D_{m}), with D_{m} being an indicator of the shape of the DSD and N_{t} an indicator of the total intensity. They can be computed from the DSD as follows (Leinonen et al., 2012; Jaffrain and Berne, 2012):
It should be noted that the disdrometer provides data binned per class of equivolumic diameter and fall velocity, from which a discrete DSD is computed and then used to evaluate the integrals of Eqs. (18) and (19) (see Gires et al., 2018b, for more details).
Multifractal analyses are carried out on ensemble analyses, i.e. on average over various samples. Once rainfall events (an event is defined as a rainy period during which more than 1 mm is collected and that is separated by more than 15 min of dry conditions before and after) have been selected within the longer time series, a process similar to in Gires et al. (2016) and Gires et al. (2018a) is implemented to extract the various samples of data: “for each event (i) a sample size is chosen (a power of two, if possible); (ii) the maximum number of samples for this event is computed; (iii) the portion of the event of length equal to the sample size multiplied by the number of samples found in (ii) with the greatest cumulative depth is extracted; (iv) the extracted series is cut into various samples.” Since D_{m} is not defined when there is no rain, only samples with no zeros are used.
Dyadic sample sizes are simpler to use for multifractal analysis, which results in some data not being used. With the process described above, 63, 52, 38 and 22 % of the data is actually not used for sample sizes of 32, 64, 128 and 256 respectively. A size of 32 time steps, corresponding to 16 min, is used, to maximize the amount of data used while keeping an acceptable length for the studied time series. An example of a sample for the four studied quantities during a rainfall event that occurred on 15 January 2018 is shown in Fig. 9. A total of 491 such samples are used in the analysis.
5.2 Joint analysis and discussion
Let us first discuss the results of the joint multifractal analysis carried out between N_{t} and R. The purpose is to check if the scaleinvariant analysis of correlations is relevant for these fields and then to quantify their correlations in this framework (i.e. write the fields as in Eq. 12 (top) and estimate a, b and α_{Y} from the two fields only).
The main curves are shown in Fig. 10, with ϵ_{λ} being the fluctuations of N_{t} and ϕ_{λ} being the fluctuations of R. The analysis directly on the field showed that they were nonconservative, meaning that the TM and DTM analysis would be biased. Hence multifractal analysis was carried out on an approximation of the underlying conservative fields consisting of their fluctuations (Lavallée et al., 1993). Numerical values of the various parameters of the analysis are in Tables 1 and 2. R exhibits a very good scaling behaviour on the whole range of scales taken into account, as shown by the TM analysis where the coefficients of correlation r^{2} of the linear regressions for q around 1 are all greater than 0.98 (Fig. 10b). Similar scaling behaviour was found on a previous campaign with the same devices (Gires et al., 2016). The scaling for N_{t} is worse, with r^{2} values only slightly greater than 0.9, but it remains acceptable (Fig. 10a). We find α_{R}=1.86 and ${C}_{\mathrm{1},R}=\mathrm{0.14}$ and ${\mathit{\alpha}}_{{N}_{\mathrm{t}}}=\mathrm{1.78}$ and ${C}_{\mathrm{1},{N}_{\mathrm{t}}}=\mathrm{0.10}$. The values of UM parameters observed mean that we are in the domain of highest relevance of the framework developed in the previous section. For R, and to a lesser extent N_{t}, there is a clear departure of the fitted K(q) from the empirical K(q), with much greater values for the fitted curve. Furthermore the empirical K(q) values exhibit a linear behaviour for q greater than approximately 1.5 (Fig. 10c). Such behaviour is consistent with the expected one when a multifractal phase transition associated with sampling limitations occurs.
The joint multifractal analysis (Eq. 4 in log–log) for $q=h=\mathrm{0.7}$ of the two studied fields is displayed in Fig. 10d. The scaling is good, with a value of r^{2}=0.97 for the linear fit. It enables the exponents a and b to be estimated at 0.33 and 0.75 respectively (Fig. 10e). The corresponding IC is equal to 0.18. In addition to quantifying the level of correlations between the two fields, it suggests how to simulate one from the other. More precisely, once a time series of fluctuations of R is available, it is possible to simulate a realistic corresponding time series of fluctuations of N_{t}, by raising to the power a=0.33 the R series and multiplying it with an independent random field Y with α=1.76 and C_{1}=0.14 raised to the power b=0.75, and renormalizing the ensemble. Formally it suggests the fluctuations of N_{t} can be written as
Such relations open the path for techniques to simulate fluctuations of N_{t} knowing only the temporal evolution of the rain rate.
Similar qualitative results are found for the other combinations, and numerical values are reported in Table 1. Both LWC and D_{m} exhibit a good scaling behaviour and their UM parameters are in Table 1. As expected given the observed values of α, the ICs computed in one way or the other (i.e. inverting the role of ϵ_{λ} and ϕ_{λ}) are very similar. Furthermore the values of a^{′} found using Eq. (17) (not shown) are very close to those obtained by inverting the role of the two fields. This confirms the relevancy of the framework of Sect. 4 in this case. It appears that the correlation found between R and LWC is much stronger than between R and N_{t} or D_{m}. There is no correlation between N_{t} or D_{m} which suggests but is not proof of independence (it would be proof for Gaussian variables). Note that the very bad scaling for the joint analysis of these two quantities is partially due to the very small values found for r(q,h) which are basically equal to zero. R exhibits a slightly greater correlation with D_{m} (IC=0.26) than with N_{t} (IC=0.18). It is the inverse for LWC with values of IC equal to 0.15 and 0.27 respectively.
In this paper, we used the framework of joint multifractal analysis to characterize the correlation across scales between two multifractal fields. We extended the existing framework to universal multifractals and also to analyse the correlations between two fields consisting of renormalized multiplicative power law combinations of two known UM fields. In general, the resulting fields can be well approximated by UM fields. Estimates of the corresponding pseudoUM parameters can be theoretically computed by focusing on the behaviour for moments close to one. These estimates remain valid for a range of moments between ∼ 0.6 and ∼ 1.6 in the worst case. The closer the two α of the initial fields, the better the approximation. When both α values are equal, the approximation is exact. An analysis technique to estimate the properties of the underlying fields (UM parameters and power law exponents used in the combination) was developed and validated with the help of numerical simulations.
In a second step, this analysis was used to develop an innovative framework to investigate the correlations between two UM fields. It basically consists of looking at the best parameters, enabling one field to be written as a power law multiplicative combination of the other field and a random one. In this context, a good candidate for a simple indicator of the strength of the correlation (called IC) is the proportion of intermittency of a field explained by the other one. In the general case, this framework is not symmetric, which is a limitation. However when the α values are typically greater than ∼0.8, it is approximately symmetric, meaning that it is relevant to extract some information on the correlations between two fields.
Finally this was implemented on rainfall data collected by a disdrometer installed on the roof the Ecole des Ponts ParisTech. More precisely the correlations between R and LWC and DSD features (N_{t} and D_{m}) are investigated. First it should be mentioned that the scaling behaviour of both R and LWC is excellent, while that of the DSD features is only good. The α values are rather similar and greater than 1.7, meaning that it is a favourable context in which to use the newly developed approach. It appears that the correlation between R and LWC is as expected very strong, the one between R or LWC and the DSD features is medium, and the one between N_{t} and D_{m} is basically null. Besides quantifying these correlations, the developed framework suggests a simple technique to simulate one field from the other. Indeed, it is sufficient to compute a power law multiplicative combination between one field and a random one to obtain the other. The characteristic parameters of the random field as well as the power law exponents of the relation can be obtained through a joint multifractal analysis of the two studied fields.
Further investigations on other fields in various contexts should be carried out to confirm the ability of this framework to both characterize and simulate correlations across scales between two multifractal fields. In future work, this framework should also be extended to more than two fields.
The data and python scripts which have been used for this paper have been made available on a public repository. It can be accessed at https://doi.org/10.5281/zenodo.3707904 (Gires et al., 2020).
AG designed the study and wrote the paper. The joint analysis by the authors of the obtained results shaped the paper into its actual form.
The authors declare that they have no conflict of interest.
Authors gratefully acknowledge partial financial support from the chair “Hydrology for Resilient Cities” (endowed by Veolia) of Ecole des Ponts ParisTech and the ÎledeFrance region RadX@IdF Project.
This paper was edited by Stéphane Vannitsem and reviewed by two anonymous referees.
Battaglia, A., Rustemeier, E., Tokay, A., Blahak, U., and Simmer, C.: PARSIVEL Snow Observations: A Critical Assessment, J. Atmos. Ocean. Tech., 27, 333–344, https://doi.org/10.1175/2009JTECHA1332.1, 2010. a
Bertol, I., Schick, J., Bandeira, D. H., PazFerreiro, J., and Vázquez, E. V.: Multifractal and joint multifractal analysis of water and soil losses from erosion plots: A case study under subtropical conditions in Santa Catarina highlands, Brazil, Geoderma, 287, 116–125, https://doi.org/10.1016/j.geoderma.2016.08.008, 2017. a, b
Calif, R. and Schmitt, F. G.: Multiscaling and joint multiscaling description of the atmospheric wind speed and the aggregate power output from a wind farm, Nonlinear Process. Geophys., 21, 379–392, https://doi.org/10.5194/npg213792014, 2014. a, b
Douglas, E. M. and Barros, A. P.: Probable maximum precipitation estimation using multifractals: Application in the eastern United States, J. Hydrometeorol., 4, 1012–1024, 2003. a
Gires, A., Tchiguirinskaia, I., and Schertzer, D.: Multifractal comparison of the outputs of two optical disdrometers, Hydrol. Sci. J., 61, 1641–1651, https://doi.org/10.1080/02626667.2015.1055270, 2016. a, b
Gires, A., Tchiguirinskaia, I., and Schertzer, D.: Pseudoradar algorithms with two extremely wet months of disdrometer data in the Paris area, Atmos. Res., 203, 216–230, https://doi.org/10.1016/j.atmosres.2017.12.011, 2018a. a
Gires, A., Tchiguirinskaia, I., and Schertzer, D.: Two months of disdrometer data in the Paris area, Earth Syst. Sci. Data, 10, 941–950, https://doi.org/10.5194/essd109412018, 2018b. a, b
Gires, A., Tchiguirinskaia, I., and Schertzer, D.: Data for “Approximate multifractal correlation and products of universal multifractal fields, with application to rainfall data” by Auguste Gires, Ioulia Tchiguirinskaia, and Daniel Schertzer, NPG 2020 [Data set], Zenodo, https://doi.org/10.5281/zenodo.3707904, 2020. a
Hubert, P., Tessier, Y., Ladoy, P., Lovejoy, S., Schertzer, D., Carbonnel, J. P., Violette, S., Desurosne, I., and Schmitt, F.: Multifractals and extreme rainfall events, Geophys. Res. Lett., 20, 931–934, 1993. a
Jaffrain, J. and Berne, A.: Influence of the Subgrid Variability of the Raindrop Size Distribution on Radar Rainfall Estimators, J. Appl. Meteorol. Clim., 51, 780–785, https://doi.org/10.1175/JAMCD110185.1, 2012. a
JiménezHornero, F., PavónDomínguez, P., de Ravé, E. G., and ArizaVillaverde, A.: Joint multifractal description of the relationship between wind patterns and land surface air temperature, Atmos. Res., 99, 366–376, https://doi.org/10.1016/j.atmosres.2010.11.009, 2011. a, b
Kahane, J.: Sur le Chaos Multiplicatif, Ann. Sci. Math. Que., 9, 435–444, 1985. a
Lavallée, D., Lovejoy, S., and Ladoy, P.: Nonlinear variability and landscape topography: analysis and simulation, in: Fractas in geography, edited by: de Cola, L. and Lam, N., 171–205, Englewood Cliffs, PrenticeHall, 1993. a, b
Leinonen, J., Moisseev, D., Leskinen, M., and Petersen, W. A.: A Climatology of Disdrometer Measurements of Rainfall in Finland over Five Years with Implications for Global Radar Observations, J. Appl. Meteorol. Clim., 51, 392–404, https://doi.org/10.1175/JAMCD11056.1, 2012. a
Mandelbrot, B. B.: Intermittent turbulence in selfsimilar cascades: divergence of high moments and dimension of the carrier, J. Fluid Mech., 62, 331–358, https://doi.org/10.1017/S0022112074000711, 1974. a
Meneveau, C., Sreenivasan, K. R., Kailasnath, P., and Fan, M. S.: Joint multifractal measures: Theory and applications to turbulence, Phys. Rev. A, 41, 894–913, https://doi.org/10.1103/PhysRevA.41.894, 1990. a, b, c
OTT: Operating instructions, Present Weather Sensor OTT Parsivel2, 2014. a
Schertzer, D. and Lovejoy, S.: Physical modelling and analysis of rain and clouds by anisotropic scaling and multiplicative processes, J. Geophys. Res., 92, 9693–9714, 1987. a, b, c
Schertzer, D. and Lovejoy, S.: Hard and soft multifractal processes, Physica A, 185, 187–194, 1992. a, b
Schertzer, D. and Lovejoy, S.: Universal multifractals do exist!: Comments, J. Appl. Meteorol., 36, 1296–1303, 1997. a, b
Schertzer, D. and Lovejoy, S.: Multifractals, generalized scale invariance and complexity in geophysics, Int. J. Bifurc. Chaos, 21, 3417–3456, 2011. a
Schertzer, D. and Tchiguirinskaia, I.: Multifractal vector fields and stochastic Clifford algebra, Chaos, 25, 123127, https://doi.org/10.1063/1.4937364, 2015. a
Schertzer, D. and Tchiguirinskaia, I.: An Introduction to Multifractals and Scale Symmetry Groups, in: Fractals: Concepts and Applications in Geosciences, edited by: Ghanbarian, B. and Hunt, A., 1–28, CRC Press, 2017. a
Schertzer, D. and Tchiguirinskaia, I.: A century of turbulent cascades and the emergence of multifractal operators, Earth Space Sci., 7, e2019EA000608, https://doi.org/10.1029/2019EA000608, 2020. a
Seuront, L. and Schmitt, F. G.: Multiscaling statistical procedures for the exploration of biophysical couplings in intermittent turbulence. Part I. Theory, DeepSea Res. Pt. II, 52, 1308–1324, https://doi.org/10.1016/j.dsr2.2005.01.006, 2005a. a, b, c
Seuront, L. and Schmitt, F. G.: Multiscaling statistical procedures for the exploration of biophysical couplings in intermittent turbulence. Part II. Applications, DeepSea Res. Pt. II, 52, 1325–1343, https://doi.org/10.1016/j.dsr2.2005.01.005, 2005b. a
Siqueira, G. M., Silva, Ê F. F., VidalVázquez, E., and PazGonzález, A.: Multifractal and joint multifractal analysis of general soil properties and altitude along a transect, Biosyst. Eng., 168, 105–120, https://doi.org/10.1016/j.biosystemseng.2017.08.024, 2018. a, b
Wang, Z.Y., Shu, Q.S., Xie, L.Y., Liu, Z.X., and Si, B.: Joint Multifractal Analysis of Scaling Relationships Between Soil WaterRetention Parameters and Soil Texture, Pedosphere, 21, 373–379, https://doi.org/10.1016/S10020160(11)601380, 2011. a, b
Xie, W.J., Jiang, Z.Q., Gu, G.F., Xiong, X., and Zhou, W.X.: Joint multifractal analysis based on the partition function approach: analytical analysis, numerical simulation and empirical application, New J. Phys., 17, 103020, https://doi.org/10.1088/13672630/17/10/103020, 2015. a, b