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

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.

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

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 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 non-linearly interacting multifractal processes and correspond to a broad, multiplicative generalization of the central limit theorem 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é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. Schmitt (2005a, 2005b) 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

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 (λ = 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 Lovejoy, 1987, 1997). Only two parameters with physical interpretation are needed to define K(q) for conservative fields: -C 1 , the mean intermittency co-dimension, which measures the clustering of the (average) intensity at smaller and smaller scales (C 1 = 0 for a homogeneous field); α, the multifractality index (0 ≤ α ≤ 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, andDouglas andBarros, 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(q) ≈ +∞ 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).

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) = 0 for independent fields. If they are power law combinations related with φ λ = c d λ , then r(q, h) is symmetric in the dp-q plane.

Multiplicative combinations of two fields
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.

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 1,X = 0.3, α Y = 0.8 and C 1,Y = 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.

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 pseudo-UM 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 Figure 1. λ (in red) and X λ (in blue) for a realization of X λ and Y λ with α X = 1.8, C 1,X = 0.3, α Y = 0.8 and C 1,Y = 0.3. Definition of Eq. (5) is used. Various values of a are shown; b is tuned to ensure the same C 1 is retrieved on all the fields. λ 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(q D ) = (q D − 1)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.

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 , C 1,X , α and C 1, to be estimated. Without any loss of generality, we can assume that C 1,Y = C 1,X . Indeed C 1,Y is a rather arbitrary quantity that can be changed while the one that actually matters is C 1,Y b α 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 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 α C 1, = C 1,X a α X α X + C 1,Y b α Y α Y , and that the term C 1,Y b α Y is simply equal to C 1, − C 1,X a α X , which enables the non-linear part of the equation to be removed): can be used to estimate b as follows (noting that Figure 2. Illustration 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.

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 (inscale) multivariate analysis of multifractal fields or measures 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 1,X = 0.3, α Y = 0.8, C 1,Y = 0.3, a = 0.6 and b = 0.2. As a consequence we expect to find α = 1.39 and C 1, = 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 sub-section), we consider the two moments q = h = 0.7. Note that with these values we have ha + q = 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

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 Figure 4. Estimate 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.
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 1,Y = C 1,φ and C 1,Z = C 1, . This enables the following calculations to be simplified.

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 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 1, = 0.4, α φ = 0.8, C 1,φ = 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 = 0.19 obtained with h = q = 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, α = α φ = 2, Eq. (14) becomes meaning that once hq has been removed, a is deterministically obtained once a is set and r φ (q, h) = r φ (q, h) = C 1,φ 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 1, = 0.4, α φ = 0.8, C 1,φ = 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 0 ≤ α Y ≤ 2 must be respected, leading to a ≤ min[( 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.

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 Figure 8. Plot of IC φ and IC φ as a function of a (Eq. 17) for the same data that is presented in Fig. 7.
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 Implementation on rainfall data

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

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 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 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 1,R = 0.14 and α N t = 1.78 and C 1,N t = 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 = 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.

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 (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. 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 (Gires et al., 2020).
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. The authors declare that they have no conflict of interest.