Journal topic
Nonlin. Processes Geophys., 27, 133–145, 2020
https://doi.org/10.5194/npg-27-133-2020
Nonlin. Processes Geophys., 27, 133–145, 2020
https://doi.org/10.5194/npg-27-133-2020

Research article 19 Mar 2020

Research article | 19 Mar 2020

# Approximate multifractal correlation and products of universal multifractal fields, with application to rainfall data

Approximate multifractal correlation and products of universal multifractal fields, with application to rainfall data
Auguste Gires, Ioulia Tchiguirinskaia, and Daniel Schertzer Auguste Gires et al.
• Hydrologie Météorologie et Complexité, Ecole des Ponts ParisTech, Champs-sur-Marne, France

Correspondence: Auguste Gires (auguste.gires@enpc.fr)

Abstract

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/portfolio-archive/taranis-observatory/, 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 (Nt) along with the mass weighed diameter (Dm), 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 Nt and Dm) is discussed.

1 Introduction

Numerous geophysical fields exhibit intermittent features with sharp fluctuations across all scales, skewed probability distribution and long-range 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 scale-invariant properties of the governing Navier–Stokes equations and hence should exhibit scale invariant features as well. The reader is referred to the reviews by and for more details. In the large class of universal multifractals (UMs), which are the stable and attractive limits of non-linearly interacting multifractal processes and correspond to a broad, multiplicative generalization of the central limit theorem , 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énez-Hornero et al. (2011) focused on the links between wind patterns and surface temperature. Xie et al. (2015) used this framework in a non-geophysical domain to better understand the cross-correlation between stock market indexes and the index of volatilities.

Seuront and Schmitt (, ) suggested a refinement of this framework and introduced a re-normalization 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 log-normal 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 Theoretical framework

## 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 λ:

$\begin{array}{}\text{(1)}& 〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}〉\approx {\mathit{\lambda }}^{K\left(q\right)},\end{array}$

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 . Only two parameters with physical interpretation are needed to define K(q) for conservative fields:

• C1, the mean intermittency co-dimension, which measures the clustering of the (average) intensity at smaller and smaller scales (C1=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

$\begin{array}{}\text{(2)}& K\left(q\right)=\frac{{C}_{\mathrm{1}}}{\mathit{\alpha }-\mathrm{1}}\left({q}^{\mathit{\alpha }}-q\right).\end{array}$

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 . One can also note that UM parameters characterize the first and second derivatives of K(q) near q=1:

$\begin{array}{ll}& {K}^{\prime }\left(\mathrm{1}\right)={C}_{\mathrm{1}},\\ \text{(3)}& & {K}^{\prime \prime }\left(\mathrm{1}\right)={C}_{\mathrm{1}}\mathit{\alpha }.\end{array}$

When doing a multifractal analysis, one should keep in mind that such fields can be affected by multifractal phase transitions . 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 qs are not reliable (see , and , for some examples of implementation). In practice, the empirical curve of K(q) will become linear from qs and hence depart (being below) from the theoretical curve. The second one is trickier and associated with the divergence of moments . The issue was also mentioned in and but they did not address the quantification of the spurious statistical estimates on finite samples and their dependence on their size . 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>qD. 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 refined the initial framework of and suggested performing a joint multifractal analysis as follows:

$\begin{array}{}\text{(4)}& \frac{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}{\mathit{\varphi }}_{\mathit{\lambda }}^{h}〉}{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}〉〈{\mathit{\varphi }}_{\mathit{\lambda }}^{h}〉}\approx {\mathit{\lambda }}^{S\left(q,h\right)-{K}_{\mathit{ϵ}}\left(q\right)-{K}_{\mathit{\varphi }}\left(h\right)}\approx {\mathit{\lambda }}^{r\left(q,h\right)},\end{array}$

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\left(q,h\right)=\mathrm{0}$ for independent fields. If they are power law combinations related with ${\mathit{\varphi }}_{\mathit{\lambda }}=c{\mathit{ϵ}}_{\mathit{\lambda }}^{d}$, then r(q,h) is symmetric in the dpq plane.

3 Multiplicative combinations of two fields

Let us consider two independent UM fields Xλ and Yλ, with their respective characteristic parameters αX, C1,X, αY and C1,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

$\begin{array}{}\text{(5)}& {\mathit{ϵ}}_{\mathit{\lambda }}=\frac{{X}_{\mathit{\lambda }}^{a}\phantom{\rule{0.25em}{0ex}}{Y}_{\mathit{\lambda }}^{b}}{〈{X}_{\mathit{\lambda }}^{a}\phantom{\rule{0.25em}{0ex}}{Y}_{\mathit{\lambda }}^{b}〉},\end{array}$

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

Figure 1ϵλ (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}$. Definition of Eq. (5) is used. Various values of a are shown; b is tuned to ensure the same C1 is retrieved on all the 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

$\begin{array}{ll}〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}〉={\mathit{\lambda }}^{{K}_{\mathit{ϵ}}\left(q\right)}& =\frac{〈{X}_{\mathit{\lambda }}^{qa}〉〈{Y}_{\mathit{\lambda }}^{qb}〉}{{〈{X}_{\mathit{\lambda }}^{a}〉}^{q}{〈{Y}_{\mathit{\lambda }}^{b}〉}^{q}}\\ \text{(6)}& & ={\mathit{\lambda }}^{{K}_{X}\left(qa\right)-q{K}_{X}\left(a\right)+{K}_{Y}\left(qb\right)-q{K}_{Y}\left(b\right)},\end{array}$

which means we have

$\begin{array}{ll}{K}_{\mathit{ϵ}}\left(q\right)& ={a}^{{\mathit{\alpha }}_{X}}{K}_{X}\left(q\right)+{b}^{{\mathit{\alpha }}_{Y}}{K}_{Y}\left(q\right)\\ & ={a}^{{\mathit{\alpha }}_{X}}\frac{{C}_{\mathrm{1},X}}{{\mathit{\alpha }}_{X}-\mathrm{1}}\left({q}^{{\mathit{\alpha }}_{X}}-q\right)+{b}^{{\mathit{\alpha }}_{Y}}\frac{{C}_{\mathrm{1},Y}}{{\mathit{\alpha }}_{Y}-\mathrm{1}}\left({q}^{{\mathit{\alpha }}_{Y}}-q\right)\\ \text{(7)}& & \approx \frac{{C}_{\mathrm{1},\mathit{ϵ}}}{{\mathit{\alpha }}_{\mathit{ϵ}}-\mathrm{1}}\left({q}^{{\mathit{\alpha }}_{\mathit{ϵ}}}-q\right).\end{array}$

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 pseudo-UM parameters C1,ϵ 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

$\begin{array}{ll}{C}_{\mathrm{1},\mathit{ϵ}}& ={C}_{\mathrm{1},X}{a}^{{\mathit{\alpha }}_{X}}+{C}_{\mathrm{1},Y}{b}^{{\mathit{\alpha }}_{Y}},\\ \text{(8)}& {\mathit{\alpha }}_{\mathit{ϵ}}& =\frac{{C}_{\mathrm{1},X}{a}^{{\mathit{\alpha }}_{X}}{\mathit{\alpha }}_{X}+{C}_{\mathrm{1},Y}{b}^{{\mathit{\alpha }}_{Y}}{\mathit{\alpha }}_{Y}}{{C}_{\mathrm{1},X}{a}^{{\mathit{\alpha }}_{X}}+{C}_{\mathrm{1},Y}{b}^{{\mathit{\alpha }}_{Y}}}.\end{array}$

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 qD are much greater than the moment for which the discrepancies start to be visible. In the cases of Fig. 2, we have qD=5.96 for panel (a), qD=4.58 for panel (b) and qD=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)=\left({q}_{D}-\mathrm{1}\right)D$ using the pseudo-UM 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.

Figure 2Illustration of the scaling moment functions K(q) of Xλ, Yλ and ϵλ, along with the UM approximation for ϵλ (fitted around q=1). Three possible sets of parameters are displayed.

## 3.3 Techniques for retrieving parameters

In this sub-section 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, C1,X, αϵ and C1,ϵ to be estimated. Without any loss of generality, we can assume that ${C}_{\mathrm{1},Y}={C}_{\mathrm{1},X}$. Indeed C1,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}r\left(q,h\right)& ={K}_{X}\left(ha+q\right)-K\left(ha\right)-K\left(q\right)\\ \text{(9)}& & =\frac{{C}_{\mathrm{1},X}}{{\mathit{\alpha }}_{X}-\mathrm{1}}\left(\left(ha+q{\right)}^{{\mathit{\alpha }}_{X}}-\left(ha{\right)}^{{\mathit{\alpha }}_{X}}-\left(q{\right)}^{{\mathit{\alpha }}_{X}}\right).\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 remp(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 qs and qD.

• Step 3: estimating αY. Using Eq. (8), one can easily obtain the following (noting that ${\mathit{\alpha }}_{\mathit{ϵ}}{C}_{\mathrm{1},\mathit{ϵ}}={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{ϵ}}-{C}_{\mathrm{1},X}{a}^{{\mathit{\alpha }}_{X}}$, which enables the non-linear part of the equation to be removed):

$\begin{array}{}\text{(10)}& {\mathit{\alpha }}_{Y}=\frac{\frac{{C}_{\mathrm{1},\mathit{ϵ}}}{{C}_{\mathrm{1},X}}{\mathit{\alpha }}_{\mathit{ϵ}}-{a}^{{\mathit{\alpha }}_{X}}{\mathit{\alpha }}_{X}}{\frac{{C}_{\mathrm{1},\mathit{ϵ}}}{{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{ϵ}}-{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(\frac{{C}_{\mathrm{1},\mathit{ϵ}}}{{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 in-scale cascades. It consists of iteratively repeating a cascade step with a non-infinitesimal 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 (in-scale) multivariate analysis of multifractal fields or measures .

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{ϵ}}=\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 C1 respectively, which is close to the values input in the simulations. The small discrepancy in C1 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 (qs=1.95 for the input UM parameters). It can be noted that for ϵλ we have a greater qs equal to 1.95, while it is even greater for Yλ(=4.5). The values of qD 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 sub-section), 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 qs 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.

Figure 3Results of numerical analysis with α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 input parameters. (a) TM analysis i.e. Eq. (1) in log–log plot, for ϵλ. (b) Scaling moment functions K(q) for ϵλ and Xλ. (c) Joint multifractal analysis (Eq. 4 in log–log plot) for $q=h=\mathrm{0.7}$. (d) Illustration of the estimation of a with the values r(0.7,0.7) computed in  (c).

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.

Figure 4Estimate 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. The blue grid at the constant value of 0.6 corresponds to the value of a inputted in the simulations.

4 Toward an indicator of correlation

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

$\begin{array}{ll}{\mathit{ϵ}}_{\mathit{\lambda }}& =\frac{{\mathit{\varphi }}_{\mathit{\lambda }}^{a}\phantom{\rule{0.25em}{0ex}}{Y}_{\mathit{\lambda }}^{b}}{〈{\mathit{\varphi }}_{\mathit{\lambda }}^{a}\phantom{\rule{0.25em}{0ex}}{Y}_{\mathit{\lambda }}^{b}〉},\\ \text{(12)}& {\mathit{\varphi }}_{\mathit{\lambda }}& =\frac{{\mathit{ϵ}}_{\mathit{\lambda }}^{{a}^{\prime }}{Z}_{\mathit{\lambda }}^{{b}^{\prime }}}{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{{a}^{\prime }}{Z}_{\mathit{\lambda }}^{{b}^{\prime }}〉},\end{array}$

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{ϵ}}$. 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:

$\begin{array}{ll}\frac{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}{\mathit{\varphi }}_{\mathit{\lambda }}^{h}〉}{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}〉〈{\mathit{\varphi }}_{\mathit{\lambda }}^{h}〉}& =\frac{〈{\mathit{\varphi }}_{\mathit{\lambda }}^{aq+h}{Y}_{\mathit{\lambda }}^{bq}〉}{〈{\mathit{\varphi }}_{\mathit{\lambda }}^{ah}〉〈{Y}_{\mathit{\lambda }}^{bq}〉〈{\mathit{\varphi }}_{\mathit{\lambda }}^{q}〉}=\frac{〈{\mathit{\varphi }}_{\mathit{\lambda }}^{aq+h}〉}{〈{\mathit{\varphi }}_{\mathit{\lambda }}^{aq}〉〈{\mathit{\varphi }}_{\mathit{\lambda }}^{h}〉}\\ & ={\mathit{\lambda }}^{{r}_{\mathit{ϵ}\mathit{\varphi }}\left(q,h\right)}={\mathit{\lambda }}^{\frac{{C}_{\mathrm{1},\mathit{\varphi }}}{{\mathit{\alpha }}_{\mathit{\varphi }}-\mathrm{1}}\left[\left(aq+h{\right)}^{{\mathit{\alpha }}_{\mathit{\varphi }}}-\left(aq{\right)}^{{\mathit{\alpha }}_{\mathit{\varphi }}}-\left(h{\right)}^{{\mathit{\alpha }}_{\mathit{\varphi }}}\right]},\\ \frac{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}{\mathit{\varphi }}_{\mathit{\lambda }}^{h}〉}{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}〉〈{\mathit{\varphi }}_{\mathit{\lambda }}^{h}〉}& =\frac{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q+{a}^{\prime }h}{Z}_{\mathit{\lambda }}^{h{b}^{\prime }}〉}{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}〉〈{\mathit{ϵ}}_{\mathit{\lambda }}^{{a}^{\prime }h}〉〈{Z}_{\mathit{\lambda }}^{h{b}^{\prime }}〉}=\frac{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q+{a}^{\prime }h}〉}{〈{\mathit{ϵ}}_{\mathit{\lambda }}^{q}〉〈{\mathit{ϵ}}_{\mathit{\lambda }}^{{a}^{\prime }h}〉}\\ \text{(13)}& & ={\mathit{\lambda }}^{{r}_{\mathit{\varphi }\mathit{ϵ}}\left(q,h\right)}={\mathit{\lambda }}^{\frac{{C}_{\mathrm{1},\mathit{ϵ}}}{{\mathit{\alpha }}_{\mathit{ϵ}}-\mathrm{1}}\left[\left(q+{a}^{\prime }h{\right)}^{{\mathit{\alpha }}_{\mathit{ϵ}}}-\left(q{\right)}^{{\mathit{\alpha }}_{\mathit{ϵ}}}-\left({a}^{\prime }h{\right)}^{{\mathit{\alpha }}_{\mathit{ϵ}}}\right]},\end{array}$

$\begin{array}{ll}\forall \phantom{\rule{0.33em}{0ex}}h,q& \frac{{C}_{\mathrm{1},\mathit{\varphi }}}{{\mathit{\alpha }}_{\mathit{\varphi }}-\mathrm{1}}\left[\left(qa+h{\right)}^{{\mathit{\alpha }}_{\mathit{\varphi }}}-\left(qa{\right)}^{{\mathit{\alpha }}_{\mathit{\varphi }}}-\left(h{\right)}^{{\mathit{\alpha }}_{\mathit{\varphi }}}\right]\\ \text{(14)}& & =\frac{{C}_{\mathrm{1},\mathit{ϵ}}}{{\mathit{\alpha }}_{\mathit{ϵ}}-\mathrm{1}}\left[\left(q+{a}^{\prime }h{\right)}^{{\mathit{\alpha }}_{\mathit{ϵ}}}-\left(q{\right)}^{{\mathit{\alpha }}_{\mathit{ϵ}}}-\left({a}^{\prime }h{\right)}^{{\mathit{\alpha }}_{\mathit{ϵ}}}\right].\end{array}$

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 left-hand 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{ϵ}}=\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{ϵ}}={\mathit{\alpha }}_{\mathit{\varphi }}=\mathrm{2}$, Eq. (14) becomes

$\begin{array}{}\text{(15)}& \forall \phantom{\rule{0.33em}{0ex}}h,q\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}{C}_{\mathrm{1},\mathit{\varphi }}ahq={C}_{\mathrm{1},\mathit{ϵ}}{a}^{\prime }hq,\end{array}$

meaning that once hq has been removed, a is deterministically obtained once a is set and ${r}_{\mathit{ϵ}\mathit{\varphi }}\left(q,h\right)={r}_{\mathit{ϵ}\mathit{\varphi }}\left(q,h\right)={C}_{\mathrm{1},\mathit{\varphi }}ahq$ is linear with regard to h and q.

Figure 5Estimates of a as a function of h and q using Eq. (14) and the process described in Sect. 3.3. Computations are carried out with αϵ=0.8, ${C}_{\mathrm{1},\mathit{ϵ}}=\mathrm{0.4}$, αϕ=0.8, ${C}_{\mathrm{1},\mathit{\varphi }}=\mathrm{0.2}$ and a=0.2. The blue horizontal grid corresponds to the value obtained with Eq. (17).

Figure 6Both sides of Eq. (14) for αϵ=0.8, ${C}_{\mathrm{1},\mathit{ϵ}}=\mathrm{0.4}$, αϕ=0.8 and ${C}_{\mathrm{1},\mathit{\varphi }}=\mathrm{0.2}$ in the case a=0.2 and ${a}^{\prime }=\mathrm{0.19}$. rϵϕ(q,h) is in red and corresponds to the left-hand side of Eq. (14) while rϕϵ(q,h) is the right-hand side and is in blue. Two views of the same figure are provided to improve visualization.

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{ϵ}}=\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[\left(\frac{{C}_{\mathrm{1},\mathit{ϵ}}{\mathit{\alpha }}_{\mathit{ϵ}}}{{C}_{\mathrm{1},\mathit{\varphi }}{\mathit{\alpha }}_{\mathit{\varphi }}}{\right)}^{\mathrm{1}/{\mathit{\alpha }}_{\mathit{\varphi }}},\left(\frac{{C}_{\mathrm{1},\mathit{ϵ}}\left(\mathrm{2}-{\mathit{\alpha }}_{\mathit{ϵ}}\right)}{{C}_{\mathrm{1},\mathit{\varphi }}\left(\mathrm{2}-{\mathit{\alpha }}_{\mathit{\varphi }}\right)}{\right)}^{\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.

Figure 7Illustration of the relations between the various parameters characterizing the correlation across scales between two UM fields in the case αϵ=0.8, ${C}_{\mathrm{1},\mathit{ϵ}}=\mathrm{0.4}$, αϕ=0.8 and ${C}_{\mathrm{1},\mathit{\varphi }}=\mathrm{0.2}$. The dash line in (b) corresponds to the relation obtained by implementing Eq. (17).

## 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 C1 of one field explained by the other:

$\begin{array}{ll}{\text{IC}}_{\mathit{ϵ}\mathit{\varphi }}& =\frac{{C}_{\mathrm{1},\mathit{\varphi }}{a}^{{\mathit{\alpha }}_{\mathit{\varphi }}}}{{C}_{\mathrm{1},\mathit{ϵ}}},\\ \text{(16)}& {\text{IC}}_{\mathit{\varphi }\mathit{ϵ}}& =\frac{{C}_{\mathrm{1},\mathit{ϵ}}{a}^{\prime {\mathit{\alpha }}_{\mathit{ϵ}}}}{{C}_{\mathrm{1},\mathit{\varphi }}}.\end{array}$

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

$\begin{array}{}\text{(17)}& {a}^{\prime }={\left(\frac{{C}_{\mathrm{1},\mathit{\varphi }}}{{C}_{\mathrm{1},\mathit{ϵ}}}\right)}^{\mathrm{2}/{\mathit{\alpha }}_{\mathit{ϵ}}}{a}^{{\mathit{\alpha }}_{\mathit{\varphi }}/{\mathit{\alpha }}_{\mathit{ϵ}}}.\end{array}$

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

Figure 8Plot of ICϵϕ and ICϕϵ as a function of a (Eq. 17) for the same data that is presented in Fig. 7.

5 Implementation on rainfall data

## 5.1 Presentation of the data

The rainfall data used in this paper were collected by a OTT Parsivel2 disdrometer 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/portfolio-archive/fresnel-platform/, 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 , 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);

• Nt, the total drop concentration (m−3);

• Dm, the mass weight diameter (mm).

Nt and Dm 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)=Ntf(Dm), with Dm being an indicator of the shape of the DSD and Nt an indicator of the total intensity. They can be computed from the DSD as follows :

$\begin{array}{}\text{(18)}& & {N}_{\mathrm{t}}=\underset{{D}_{min}}{\overset{{D}_{max}}{\int }}N\left(D\right)\mathrm{d}D,\text{(19)}& & {D}_{\mathrm{m}}=\frac{{\int }_{{D}_{min}}^{{D}_{max}}N\left(D\right){D}^{\mathrm{4}}\mathrm{d}D}{{\int }_{{D}_{min}}^{{D}_{max}}N\left(D\right){D}^{\mathrm{3}}\mathrm{d}D}.\end{array}$

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 and 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 Dm 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.

Figure 9Illustration of the four studied rainfall quantities corresponding to a 32 min sample (i.e. 64 time steps) that occurred on 15 January 2018.

Figure 10Results of joint multifractal analysis for ϵλ being the fluctuations of Nt and ϕλ being the fluctuations of R. (a) TM analysis, i.e. Eq. (1) in log–log plot, for ϵλ. (b) Same as in (a) for ϕλ. (c) Scaling moment functions K(q) for ϵλ and ϕλ. (d) Joint multifractal analysis (Eq. 4 in log–log) for $q=h=\mathrm{0.7}$. (e) Illustration of the estimation of a with the values r(0.7,0.7) computed in (d).

## 5.2 Joint analysis and discussion

Let us first discuss the results of the joint multifractal analysis carried out between Nt and R. The purpose is to check if the scale-invariant 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 Nt and ϕλ being the fluctuations of R. The analysis directly on the field showed that they were non-conservative, 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 . 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 r2 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 . The scaling for Nt is worse, with r2 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 Nt, 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.

Table 1UM parameters for the studied fields.

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 r2=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 Nt, by raising to the power a=0.33 the R series and multiplying it with an independent random field Y with α=1.76 and C1=0.14 raised to the power b=0.75, and renormalizing the ensemble. Formally it suggests the fluctuations of Nt can be written as

$\frac{{R}_{\text{fluctuations}}^{\mathrm{0.33}}{Y}^{\mathrm{0.75}}}{〈{R}_{\text{fluctuations}}^{\mathrm{0.33}}{Y}^{\mathrm{0.75}}〉}.$

Such relations open the path for techniques to simulate fluctuations of Nt knowing only the temporal evolution of the rain rate.

Table 2Numerical output of the joint multifractal analysis of the four studied fields. For each box, using the notations of Eq. (12) ϵλ corresponds to the field of the column and ϕλ to the row.

Similar qualitative results are found for the other combinations, and numerical values are reported in Table 1. Both LWC and Dm 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 Nt or Dm. There is no correlation between Nt or Dm 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 Dm (IC=0.26) than with Nt (IC=0.18). It is the inverse for LWC with values of IC equal to 0.15 and 0.27 respectively.

6 Conclusions

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 pseudo-UM 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 (Nt and Dm) 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 Nt and Dm 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.

Data availability
Data availability.

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 ().

Author contributions
Author contributions.

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.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

Authors gratefully acknowledge partial financial support from the chair “Hydrology for Resilient Cities” (endowed by Veolia) of Ecole des Ponts ParisTech and the Île-de-France region RadX@IdF Project.

Review statement
Review statement.

This paper was edited by Stéphane Vannitsem and reviewed by two anonymous referees.

References

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., Paz-Ferreiro, 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/npg-21-379-2014, 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.: Pseudo-radar 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/essd-10-941-2018, 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/JAMC-D-11-0185.1, 2012. a

Jiménez-Hornero, F., Pavón-Domínguez, P., de Ravé, E. G., and Ariza-Villaverde, 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, Prentice-Hall, 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/JAMC-D-11-056.1, 2012. a

Mandelbrot, B. B.: Intermittent turbulence in self-similar 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, Deep-Sea 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, Deep-Sea Res. Pt. II, 52, 1325–1343, https://doi.org/10.1016/j.dsr2.2005.01.005, 2005b. a

Siqueira, G. M., Silva, Ê F. F., Vidal-Vázquez, E., and Paz-Gonzá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 Water-Retention Parameters and Soil Texture, Pedosphere, 21, 373–379, https://doi.org/10.1016/S1002-0160(11)60138-0, 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/1367-2630/17/10/103020, 2015. a, b