Nonlinear Processes in Geophysics Rank-Ordered Multifractal Analysis ( ROMA ) of probability distributions in fluid turbulence

Rank-Ordered Multifractal Analysis (ROMA) was introduced by Chang and Wu (2008) to describe the multifractal characteristic of intermittent events. The procedure provides a natural connection between the rankordered spectrum and the idea of one-parameter scaling for monofractals. This technique has successfully been applied to MHD turbulence simulations and turbulence data observed in various space plasmas. In this paper, the technique is applied to the probability distributions in the inertial range of the turbulent fluid flow, as given in the vast Johns Hopkins University (JHU) turbulence database. In addition, a new way of finding the continuous ROMA spectrum and the scaled probability distribution function (PDF) simultaneously is introduced.


Introduction
It is well-known that fully developed turbulent fluid flows are intermittent and multifractal (Frisch, 1995 and references therein).The intermittent turbulence and associated multifractal characteristics are also evident in the analyses of space plasmas observations (e.g., Sorriso-Valvo et al., 1999;Consolini and Chang, 2001;Bruno et al., 2001Bruno et al., , 2003;;Forman and Burlaga, 2003;Tam et al., 2005;Weygand et al, 2005;and Chang, 2009).Through the analysis of probability distribution functions (PDFs) for field fluctuations, intermittency in turbulence is characterized by a strong non-Gaussian behavior of PDF at small scales.The multifractal characteristics have generally been analyzed with structure functions or singular spectra based on the partition functions of the probability measures (Halsey et al., 1986).
Correspondence to: C. C. Wu (ccwu@ucla.edu)Recently, Chang and Wu (2008) introduced the rankordered multifractal analysis (ROMA) for analyzing intermittent fluctuations to describe the explicit multifractal characteristics and indicated how they were distributed by parametrically separating the fluctuations according to their ranks.The technique retains the spirit of the traditional structure function analysis and combines it with the idea of oneparameter scaling of monofractals.It was first applied to the results of a large-scale two-dimensional magnetohydrodynamic (MHD) turbulence simulation.It has also been successfully applied to in-situ solar wind observations (Chang et al., 2008), the broadband electric field oscillations from the auroral zone (Tam et al., 2010), the magnetosphere cusp turbulence (Lamy et al., 2008), and the AE index over approximately a complete solar cycle (Consolini and Michelis, 2011).A brief review of ROMA and its applications to some of these studies were provided in Chang et al. (2010).
In this paper, to show that ROMA can also be useful in the analysis of the fluid turbulence, we apply the technique to the Johns Hopkins University (JHU) large-scale direct numerical simulation turbulence database (Perlman et al., 2007;and Li et al., 2008).With the huge number of data points, a detailed study can be conducted.This paper is structured as follows: in Sect.2, we first provide a brief description of the JHU turbulence data set and present the PDF of longitudinal velocity fluctuations.The PDFs are shown to be well fitted with a modified version of the Castaing et al. (1990) model.In Sect.3, ROMA is applied to these PDFs and two invariant functions in the ROMA analysis are obtained.In this section, a "new" method of the ROMA analysis is also introduced.The method provides a way of obtaining the multifractal spectrum and the scaled PDF as continuous functions of a rank-ordered local scale invariant.In Sect.4, this new approach is applied to fluctuations of the solar wind turbulence and of 2-D MHD simulations.A summary is given in Sect. 5.

JHU turbulence flow database and a modified Castaing et al. model
A detailed description of the JHU turbulence database can be found in Perlman et al. (2007) and Li et al. (2008).Briefly, the data is obtained from a direct numerical simulation of forced isotropic turbulence on a 1024 3 periodic grid in a physical domain of 2π × 2π × 2π, using a pseudo-spectral parallel code.Energy is injected by keeping constant the total energy in modes such that their wave-number magnitude is less or equal to 2. After the simulation has reached a statistically stationary state, 1024 frames of data, which includes the 3 components of the velocity vector and the pressure, are generated and stored into the database.The duration of the stored data is about one large-eddy turnover time.The frames are stored at every 10 time-steps of the numerical simulation, for time t between 0 and 2.048.In this duration, the timeaverage Taylor-scale Reynolds number is 433 and the largeeddy turnover time is 2.02.The radial energy spectrum averaged over this duration indicates an inertial range of approximately from the wave number 8 to 60, which corresponds to a spatial range from 17 to 128 , with denoting the grid spacing, which is 2π/1024.Instead of the huge 1024 4 data points, 19 z-planes of data points are used in the present analysis.They are arbitrarily selected at various z locations and various times: (t,z) = (0, 9 ), (0, 499 ), (0.5, 0), (0.5, 499 ), (1, 0), (1, 499 ), (1.5, 0), (1.5, 499 ), (2, 0), (2, 99 ), (2, 499 ), (0.25, 249 ), (0.75, 249 ), (1.25, 249 ), (1.75, 249 ), (0.25, 499 ), (0.75, 499 ), (1.25, 499 ), and (1.75, 499 ).This set of 19 million data points provides sufficient statistics.
First, the fluctuations of longitudinal velocity, δv || , are considered.The fluctuation is defined by δv || (r, δ) = (v(r+δi)−v(r))•i, with i a unit vector, and δ the spatial scale.In the analysis, i is either along the x-axis or the yaxis and δ is in the range of (16 , 160 ).As an example, the PDF of δv || at scale δ = 64 , denoted by P (δv || ,δ), is shown in Fig. 1.In computing the PDF, the range of δv || is divided evenly into 1601 bins, with the bin number extending from −800 to 800.For the bin number i,  (Castaing et al., 1990).The PDF can further be decomposed into a symmetric P + (δv || ) and an asymmetric P − (δv || ) as follows (Chevillard et al., 2006):  1).In this and subsequent figures, the accuracy of the PDF is indicated by the noise level of the plot.For example, in this figure, in the region |δv || | < 300, the noise level is very small -not larger than the thickness of the curves.However, the PDF is approximatively (3.2 ± 0.4) 10 −6 at δv || = 400.The following shows that the PDFs from the JHU database can be fitted quantitatively by a modified Castaing et al. (1990) model.The original model was developed from the classical Kolmogorov's picture of turbulent energy cascade (Kolmogorov, 1941) and took into account that the energy transfer rate in the cascade at different scales is not strictly self-similar.According to the model, for a given scale δ, PDF of fluctuations X, P (X), is a convolution of a Gaussian and a log-normal distribution:

P δv
where a s is a skewness factor to account for the asymmetry, σ 0 is the most probable value of σ , which is the standard deviation of a Gaussian distribution, and λ is the width of the log-normal distribution.The normalization factor It can be shown that if a s > 0,P (X) has a maximum at X < 0. Therefore, in order to account for the fact that P (δv || ) has a maximum at δv || > 0, as shown in Fig. 1, a modified Castaing form is introduced as follows: The additional parameter X 0 gives the location of X where P (X) is a maximum.In this new form, there are four parameters, a s , σ 0 , λ and X 0 .Although there are four parameters to fit the data to Eq. ( 3), the fitting procedure is not complicated.The parameter X 0 is given by the data as mentioned above; σ 0 is specified to fit the PDF data distribution in a small range near X 0 ; λ is then chosen to fit the wings of the PDF; and a s , is adjusted to account for the asymmetry.
In Fig. 2, P (δv || , δ) is shown to be well fitted by Eq. ( 3), with X representing δv || , for several values of δ.The parameters employed for the fits for δ from 12 to 128 are given in Table 1.Also, as shown in Fig. 3, the parameters scale as: X 0 ∼ δ 0.4 , σ 0 ∼ δ 0.4 , a s ∼ constant, and λ 2 ∼ 1/(20+δ).One notes that as δ is increased, λ 2 tends to zero; it means that at large scale, the PDF is dominated by a Gaussian distribution with σ 0 , showing the intermittent characteristic of the turbulence, being Gaussian at large scale and long tails at small scale.Further discussion on these scaling relations will be provided after the presentation of the ROMA analysis.

ROMA analysis of the JHU turbulence flow
In this section, we apply the ROMA analysis to the PDFs obtained in Sect. 2 for the turbulent flows in the JHU database.
We start with a summary of the idea of ROMA.In order to understand the scaling behavior among the family of PDFs of the turbulent fluctuations, the idea of fractals is employed.
If the form of PDFs, P (X,δ), is invariant as the scale, δ, changes, then the following interesting one parameter scaling form is obtained (Hnat et al., 2002;Chang et al. 2004): It means that the PDFs would collapse onto one scaled PDF P s (Y ) where Y = X/δ s is a global scale invariant for all ranges of X for a constant scaling exponent s.The scaling δ Fig. 3. Scaling of the parameters vs. the scale in the modified Castaing model used for the fits of P (δv || ,δ) shown in Fig. 2 and given in Table 1.Circles denote the values used; dash curves are: σ 0 = 17δ 0.4 ; X 0 = 2.55δ 0.4 ; and λ 2 = 3./(20 + δ), respectively.In addition, a s ∼ constant, see Table 1.exponent s may be interpreted as a single fractal (monofractal) measure that characterizes the fluctuations of all scales through the scaling relation (Eq.4).In practical applications, the scaling law (Eq.4) is sometimes satisfied for the full range of the scaling variable and sometimes only for some restricted ranges of the variable Y (Chang et al., 2004, Podesta et al., 2006, Hnat et al., 2002).When scaling law (Eq.4) is satisfied over the full range of Y , the scaling is called self-similar or monofractal.Otherwise, the fluctuating phenomenon is multifractal.
For multifractal fluctuations, Chang and Wu (2008) proposed the ansatz of ROMA with the following scaling: Here, the scaling exponent s(Y ) is now a function of Y , meaning that the scaling exponent depends on the rank of the local scale invariant Y .According to the monofractal scaling (Eq.4), a function P (X, δ) of two variables is described by a function of a single variable plus a single number.For the case of ROMA scaling relations (Eq.5), P (X, δ) of two variables is replaced by two functions, s(Y ) and P s (Y ), of a single variable Y and the multifractal nature of the fluctuations are provided by these two functions.
One can rewrite the ROMA scaling relations (Eq.5) as: These relations suggest a new way of finding s(Y ) and P s (Y ).By using the spectrum, s(Y ), and the scaled PDF function, P s (Y ), the distributions P (X,δ) can be calculated, and can then be compared directly with PDFs from observations or simulation results.Suppose that the ROMA scaling is applicable at larger scale, a Gaussian P (δv || , δ) may result as shown in Fig. 8 for δ = 640 .However, it should be noted that, at much larger δ, the scaling relation may break down; because X may no longer increase monotonically with Y .Figure 9  As was discussed in Sect.2, the PDFs based on the JHU data can be fitted with a modified Castaing model, Eq. ( 3), with four parameters.Furthermore, the parameters scale in the following way: X 0 ∼ δ 0.4 , σ 0 ∼ δ 0.4 , a s ∼ constant, and λ 2 ∼ 1/(20 + δ).This can be understood because s(Y ) is monofractal in a small range of Y .It can be shown that if PDFs satisfy the monofractal scaling (Eq.4), then X 0 ∼ δ s , σ 0 ∼ δ s , a s ∼ constant, and λ 2 ∼ constant.The first three relations are consistent with the fit in Sect. 2. The last relation for the parameter λ, however, is not.This is because the parameter λ affects mostly at the wings of the PDFs, which correspond to the region of large Y , where s(Y ) is multifractal.The trend of λ tending to zero at large scale (intermittency) shown in Sect. 2 is due to the multifractal nature of s(Y ).| is in the units of bin size, which is set to be |δB 2 | max /800.The PDFs from the simulation are from Chang and Wu (2008).
In addition to the fluctuations of the longitudinal velocities, the PDFs of the fluctuations of the kinetic energy density, v 2 , as shown in Fig. 10, also satisfy the ROMA scaling relations.

ROMA results for 2-D MHD simulations and solar wind data
ROMA was first applied by Chang and Wu (2008) to the results of a two-dimensional MHD simulation for homogeneous turbulence (Chang et al., 2004).It has also been applied to the analyses of turbulence data observed in solar wind (Chang et al., 2008) and other regimes of the space environment.In the following, the results of the 2-D MHD simulations and the analysis of the set of solar wind data are revisited.In previous ROMA studies, s(Y ) and P s (Y ) were obtained by using the ranked-ordered structure functions S m (X,δ) for several small ranges of Y (see Eq. 6); thus, because of limitation of the statistics, s(Y ) and P s (Y ) were not continuous on Y .Now, with the new approach following Eq.( 7) and the discussion in Sect. 3 for the JHU turbulence database, continuous functions of s(Y ) and P s (Y ) can be obtained.The details of the MHD simulation and the description of the set of solar wind data can be found in the references cited above and are not included here.Figure 11 shows the results for the MHD simulation and Fig. 12 is for the solar wind data.

Summary
Traditional structure and partition function methods of multifractal analysis determine the fractal properties of various moment orders based on the entire set of fluctuations.Since most of the observed or simulated intermittent fluctuations are dominated by those of small amplitudes, the subdominant fractal characteristics of the minority fluctuations with larger amplitudes are easily masked by those characterized by the dominant population.Furthermore, these procedures are not reversible; i.e., it is not possible to recover from the structure functions and the singular spectra based on the partition functions the original scale dependent PDFs.ROMA, on the other hand, provides a physically meaningful description of intermittency and is quantitatively accurate because of the cleanliness of the procedure of statistical sampling.By using the ROMA spectrum and the scaled PDF, the distributions P (X,δ) of intermittency can be calculated and compared with the PDFs from observations or simulation results.Thus, the results are unique.In this paper, we used data from the JHU turbulence database to show that for turbulence in www.nonlin-processes-geophys.net/18/261/2011/ Nonlin.Processes Geophys., 18, 261-268, 2011 incompressible fluids the PDFs of longitudinal velocity fluctuations throughout the inertial range can be described using the ROMA technique.We have also used the model of Castaing et al. (1990) to describe the intermittency of the flow.The relationship between the intermittency and the multifractal nature of the turbulence is discussed.A new way of finding ROMA spectrum is devised; s(Y ) and P s (Y ) can now be obtained as continuous functions of the local invariant Y .The ROMA spectra for both fluid flows and MHD simulation results are qualitatively similar: they both have a decreasing trend with Y and are antipersistent with s(Y ) < 0.5.

Fig. 1 .
Fig. 1.PDF of δv || at δ = 64 , with denoting the grid spacing, and δv || in the units of the bin size, see text.In (a) and (b), P (δv || ,δ) in solid blue curve and P (−δv || ,δ) in dash red curve to show the asymmetry of the PDF with respect to δv || = 0. (b) is an enlarged version of (a) for a smaller range of δv || .The asymmetry is evident in (c) and (d) using P + and P decomposition of the PDF according to Eq. (1).In this and subsequent figures, the accuracy of the PDF is indicated by the noise level of the plot.For example, in this figure, in the region |δv || | < 300, the noise level is very small -not larger than the thickness of the curves.However, the PDF is approximatively (3.2 ± 0.4) 10 −6 at δv || = 400.

Fig. 5 .
Fig. 5. ROMA spectrum s(Y ) and scaled PDF function P s (Y ) for the fluctuations of longitudinal velocities in JHU database.(a) s(Y ), which is symmetric with respect to Y = 0, shows monofractal behavior at |Y | < 25 and decreases monotonically at larger |Y |. (bc) P s (Y ), which is asymmetric about Y = 0. (d-e) P + (Y ) and P − (Y ) of P s (Y ), similarly defined as in Eq. (1).In (b-e), open circles show P s (Y ); solid red curves show a fit with a modified Castaing form of Eq. (3) with X replaced by Y , with σ 0 = 18.25, Y 0 = 2.55, λ = 0.425, a s = 0.4.

Fig. 7 .
Fig. 7. Plots of P (δv || ,δ): solid curves are from JHU data and markers are from ROMA scaling relations: red (circles) for δ = 12 ; green (squares) for δ = 64 ; blue (triangles) for δ = 140 .(a) and (b) are for positive δv || ; (b) is an enlarged plot for δv || ≤ 300.(c) and (d) are for negative δv || ; (d) is an enlarged plot for δv || ≥ −300.Note the agreement between the calculated values from ROMA and data at δ = 12 and 140 is not as good as the one in Fig. 6 for δ from 24 to 128 .This probably is because δ = 12 is close to the dissipation range and δ = 140 close to the injection range.

Fig. 9 .
Fig. 9. P (δv || ,δ) at δ = 48 from ROMA and monofractal scaling relations.Blue curve with squares is based on multifractal ROMA spectrum s(Y ) shown in Fig. 5a; red circles are based on the monofractal scaling with s = 0.4.Both use the scaled PDF P s (Y ) given in Fig. 5b.There are little differences between the two for |δv || | < 200.However, at larger |δv || |, the value from the ROMA scaling is smaller than the value from the monofractal scaling.For instance, at δv || = 400, the value based on the ROMA scaling, which agrees with the data as shown in Fig. 6a, is smaller than the value based on the monofractal scaling by a factor of 15.
shows the difference in the PDF between monofractal and multifractal scaling relations.Since s(Y ) ∼ 0.4 at |Y | < 25, the resulting P (δv || , δ) shows little difference for |δv || | < 200 at δ = 48 .At larger |δv || |,P (δv || , δ) based on the multifractal ROMA spectrum is smaller than the values based on the monofractal scaling; this is because s(Y ) in ROMA decreases monotonically at large Y and is less than s = 0.4 assumed in the monofractal scaling.

Fig. 12 .
Fig. 12. Fluctuations of B 2 based on solar wind data.(a) and (b): ROMA spectrum s(Y ) and scaled PDF P s (Y ).(c): plots of P (δB 2 ,τ ): computed PDFs from the ROMA relations with s(Y ) and P s (Y ) are shown in solid curves; PDFs from solar wind observations are given by markers.Green (o): τ = 1000 s; blue (x): 96 s;and red (+): 9 s.δB 2 is in the units of (nT) 2 .The PDFs from solar wind data are fromChang et al. (2008).

Table 1 .
Values of four parameters in the modified Castaing model, Eq. (2), used in the fit of P (δv || ,δ) based on the JHU turbulence data for the scale δ between 12 and 128 .