Universal multifractal Martian topography

<jats:p>Abstract. In the present study, we investigate the scaling properties of the topography of Mars. Planetary topographic fields are well known to roughly exhibit (mono)fractal behavior. Indeed, the fractal formalism reproduces much of the variability observed in topography. Still, a single fractal dimension is not enough to explain the huge variability and intermittency. Previous studies have claimed that fractal dimensions might be different from one region to another, excluding a general description at the planetary scale. In this article, we analyze the Martian topographic data with a multifractal formalism to study the scaling intermittency. In the multifractal paradigm, the apparent local variation of the fractal dimension is interpreted as a statistical property of multifractal fields. We analyze the topography measured with the Mars Orbiter Laser altimeter (MOLA) at 300 m horizontal resolution, 1 m vertical resolution. We adapted the Haar fluctuation method to the irregularly sampled signal. The results suggest a multifractal behavior from the planetary scale down to 10 km. From 10 to 300 m, the topography seems to be simple monofractal. This transition indicates a significant change in the geological processes governing the Red Planet's surface.                     </jats:p>

section_content: The acquisition of altimetric data from the Mars Orbiter Laser altimeter (MOLA) has motivated numerous analyses of the Martian topography, each one aiming to properly characterize the surface roughness. A possible approach is to assume that topography can be statistically described with quantitative parameters able to characterize the geological units. Many statistical indicators have been proposed and widely explored in order to study the surface of Mars: root mean square (RMS) height, RMS slope, median slope (Kreslavsky and Head, 2000), and autocorrelation length (Aharonson et al., 2001) (see Shepard et al., 2001 for a review of all those indicators). Useful information has been obtained by the use of those indicators, but they have the disadvantage of being defined at a given scale. By construction, they do not directly take into account the scale symmetry that often occurs in the case of natural surfaces (Mandelbrot, 1967; Turcotte, 1997; Gilbert, 1989; Pommerol et al., 2012). 

On Mars, different authors (Nikora and Goring, 2004 ) have explored the scaling properties of topography by the use of scale-invariant parameters. Malamud and Turcotte (2001) has performed a wavelet analysis to study the polar topography of Mars. Other studies performed more classical wavelet analyses in the spatial domain, following the procedure proposed by Shepard et al. (2001) ; Orosei et al. (2003) using a quantity called the RMS deviation computed at different scales and exhibiting a power-law dependance on scales with a scaling exponent H (for Gaussian processes, the Hurst exponent). The value of H, evaluated at different locations, was found to be different from one region to the other, revealing a strong intermittency in the statistical variability of the surface of Mars (Aharonson et al., 2001; Orosei et al., 2003). The observed intermittency, also found in the case of Earth (Baldassarri et al., 2008) and less pronounced in the case of the Moon (Rosenburg et al., 2011 (Rosenburg et al.,, 2015)), might lead one to reject the idea of a global description of any topographic field at the planetary scale. However, modern developments in the theory of scale-invariant processes might be able to a give full account of the observed variability and intermittency. As proposed by Lavallee et al. (1993), it is possible to extend the fractal interpretation of topography to a multifractal object requiring an infinite number of fractal dimensions (one for each altitude level Gagnon et al. (2006) revealed that synthetic multifractal fields tend to reproduce very convincingly the variability and intermittency of natural surfaces. A multifractal study has also been performed on artificial topographic digital elevation models (DEMs), revealing that multifractal parameters may be able to discriminate between different stages in topographic evolution (Vidal Vazquez et al., 2008). Moreover, a particular class of multifractal called the universal multifractal has been proposed by Schertzer and Lovejoy (1987) (see Sect. 2. 2). Stable and attractive, universal multifractals are good candidates for modeling the variability of topographic fields and have the advantage of being simply characterized by three parameters (see Sect. 2. 2). In the case of Earth, this formalism has proven to be relevant (Gagnon et al., 2006), providing a satisfying fit of the structure function computed in different areas of the globe and allowing one to quantify multifractal behavior at the planetary scale through the measurement of three parameters: H (the conservation exponent or fluctuation exponent), α (the degree of multifractality) and C 1 (the codimension of the mean). In this study, we aim to measure the global value of H, α and C 1 in the case of Mars. After a careful choice in the definition of fluctuations computed from all the available data in the MOLA database, statistical moments are computed on a large range of scales in order to estimate the universal multifractal parameters.
section_name: Introduction

section_content: This section first describes the Martian topography data set used in this study and the MOLA instrument. The second part contains some elements of the theory of universal multifractality. The third part contains the description of the Haar fluctuation tool we adapted in the case of the irregularly sampled Martian topography. Smith et al., 2001 ) is a laser altimeter on board the Mars Global Surveyor (NASA) spacecraft following a quasi-polar and circular orbit. MOLA has recorded a huge amount of topographic data providing a well-detailed global mapping of the planet. With 20 measurements every 2 s, a vertical accuracy of 1 m, a surface spot size of 168 m and approximatively 300 m between two consecutive measurements, the MOLA database stored in PDS (Planetary Data System, http://pds-geosciences. wustl. edu) constitutes an ideal data set to study the global properties of Mars' topography down to the kilometer scale. Topography is calculated using the gravity field from Goddard Mars Model 3 (mgm1025, Lemoine et al., 2001). The global areoid error through degree 60 × 60 for the mgm1025 model is 1. 8 m. The total number of individual topographic measurements used for the purpose of this study is close to 600 mil-lion. For methodological reasons, we have decided to limit the scope of our analysis to the along-track series, meaning that we exclude all kind of across-track fluctuations (measurements extracted from different orbits). This restriction might introduce a bias, a consequence of the single preferred direction (close to north-south) of MOLA's orbit. Hence, anisotropic statistical properties if they occur will not be detected by our method. However, the assumption that topography is isotropic at the planetary scale seems reasonable.
section_name: Methodology

section_name: MOLA

section_name: MOLA (Mars Orbiter Laser altimeter

section_content: The topography is a two-dimensional field providing, for each pair of latitude-longitude, an elevation value. In order to investigate the scaling properties of such fields, we must study the distribution of slopes at different scales or, equivalently, the fluctuations of elevation from one point to the other. The simplest way to define fluctuations is to compute the first difference for each couple of elevation data. The formalism used in this study is largely based on Lovejoy and Schertzer (2013). 

Once fluctuations are defined, a common way to explore the scaling properties of any geophysical field is to compute statistical moments of several orders and at different scales. If the field is scaling and if the fluctuations are defined properly, statistical moments will appear as straight lines on log-log plots, meaning they follow a power-law dependency on scales. This property is expressed by Eq. ( 1 ), where h stands for fluctuation of elevation and x stands for scale. In the case of simple fractality (i. e., monofractality), the q dependency of ζ (q) is purely linear: ζ (q) = qH, where H is called the conservation exponent used in honor of Edwin Hurst (typically a fractional Brownian motion). A single coefficient (H ) is sufficient to describe the global statistics that can be easily estimated by computing statistical moments of order 1 and 2 to obtain the slope of ζ (q) = qH. 

In the multifractal case, the scaling function is corrected by a additional quantity K(q), called the moment scaling function:

Although three moments of distinct orders are in principle sufficient to establish the curvature of ζ (q) and assess multiscaling properties, a field with general multifractal symmetries would require one to determine all statistical moments (non-integer orders included) to fully describe the statistics. Practically, the estimations of moments of multiple orders can be misleading. Indeed, the estimates of high-order moments (typically moment of order > 2) are often biased by multifractal phase transition in the effective ζ (q): beyond critical values of q, statistical moments may all depend on the largest value in the sample, spuriously leading to linear behavior in the exponent (and diverging as the sample size increases) (Schertzer and Lovejoy, 1993). Therefore, we will cautiously focus on moments of order < 2 that can still provide meaningful information on the shape of ζ (q). In order to characterize multiscaling, we will use three independent multifractal parameters introduced by Schertzer (1997). 

-The mean fluctuation exponent H = ζ (1). A consequence of Eq. ( 1 ), the value H = 0 is obtained for a field that is strictly scale-invariant (conservative), whereas H = 0 rather corresponds to a fractional integration of a conservative field. 

-The codimension of the mean field H -C 1 = dξ(q) dq q=1

measuring the mean intermittency. If C 1 = 0, the field is homogenous. A non-zero value of C 1 indicated intermittency. In the case of topography, C 1 is expected to be close to 0. 1 (see the case of Earth Gagnon et al., 2006). 

that estimates the curvature of ζ (q) near q = 1. It provides information about the relative variation of intermittency around the mean. If α = 0, ζ (q) is linear around q = 1, indicating a monofractal behavior. 

As a consequence of the above definitions, any estimate of H, α and C 1 will rely on computing non-integer-order moments around the mean q = 1. In the most general multifractal case, this set of three parameters only provides a restricted description of statistics in the neighborhood of the mean (first-order moment) and presumes nothing of the variability that may be revealed by higher-order moments. Still, there exists a class of multifractal called the universal multifractals (Schertzer, 1997) whose associated structure function is fully determined by H, α and C 1. In other words, for a member of this class, a set of only three parameters is enough to characterize all the statistical moments. Equation (3) indicates the mathematical form expected in the case of universal multifractality with a moment scaling function depending only on α and C 1. 

(3)
section_name: Universal multifractality

section_content: Mars is well known for presenting huge fluctuations of altitude on its surface despite its relatively small radius in comparison to Earth (R Mars = 3390 km). Indeed, Olympus Mons culminates at 22. 5 km above the adjacent lowlands. In this context, we have to define fluctuations in order to provide an accurate statistical study. It is well established that data analysis strongly relies on defining the fluctuations at a given scale and location. Of course, the chosen definition must adequately characterize the fluctuations and has to be picked among all the possible definitions depending on the context. As mentioned earlier, the simplest definition is the absolute difference that might not be accurate in some cases. For instance, in the case of a fractal field with H < 0 (a Gaussian white noise has H = -1/2), absolute differences will only characterize fluctuations with high wave numbers instead of the local fluctuations at a given scale. In the case of topography, H is expected to be > 0, so absolute differences might be meaningful. Still, it is possible to consider a more general definition of fluctuations that can apply for both cases (-1 < H < 1): the Haar fluctuations. The definition is given by Eq. (A1). For each even number 2n of along-track consecutive points, we compute the average A of the first n points ( h i 1, i = 1. . . n) and the average B of the n last points ( h i 2, i = n + 1. . . 2n). The fluctuation at scale L is defined by the absolute difference | B -A |, L being the spherical distance between the first point and the last point. On synthetic multifractal series obtained by simulation (Lovejoy and Schertzer, 2012), Haar fluctuations have proven to be strong estimators of the input multifractal parameters. However, Haar fluctuations unlike first differences are more complicated to implement, particularly in the case of one-dimensional irregular series (see Appendix A). For the purpose of this analysis, we computed both first-difference and Haar fluctuations. Although the two methods produced similar results (as expected since H > 0), we noticed a better convergence in the case of Haar fluctuations. Hence, in the result Sect. 4, we will focus on the results obtained by Haar fluctuations.
section_name: Haar fluctuations for irregular signals

section_content: For each scale, we compute moments of several orders q = 0. 1, 0. 2. . . 2. The scale dependence of every moment is expected to be a straight line on a log-log plot with a slope predicted by Eq. (2). Figure 1 shows the results (plot on the left). In order to compare the experimental results to simulations, the right plot in Fig. 1 is obtained by performing the same analysis on synthetic series with multifractal parameters H = 0. 5, α = 1. 7 and C 1 = 0. 1. The global aspect of curves in Fig. 1 (right plot) and their resemblance to synthetic curves over a significant range of scales strongly confirm the scaling behavior of topography with a value of H close to 0. 5. Still, some features need to be discussed. We will focus on scales > 600 m where reasonably straight lines occur (due to the definition of Haar fluctuations, no relevant data can be found outside that range). In Fig. 1 (right plot), statistical moments seem to exhibit two distinct scaling regimes. The first regime starts at 10 km to the planet scale 20 000 km covering nearly 4 orders of magnitude. The second one is identified in the range of scale (600 m to 10 km) covering 1 order of magnitude. Scaling needs to be studied separately for both of these ranges of scale. In both plots, red diamonds (or red triangles) stand for moment of order 1 (or order 2). Integer-ordered moments are also computed in order to test the multifractal formalism. They are plotted in shades of red (left plot) and blue (right plot).
section_name: Results

section_content: Before going any further, we discuss in this section possible artifacts that may cause the observed transition around 10 km. 

-Uncertainty due to the accuracy of measurement at small scales: the characteristics presented in Sect. 2. 1 indicate that the error in the topography measurement is a few meters (areoid error + instrument error), lower than the mean fluctuations at the lowest relevant scale (∼ 10 m in Fig. 1 ). Hence, if there is a bias at small scales, it can not be responsible for the transition identified around 10 km, the mean fluctuation at that scale being close to 40 m. 

-Due to the quasi-polar orbit of MOLA, the density of measurements depends on latitudes. The minimum is obtained near the Equator, whereas the maximum is observed near the poles. This could introduce a bias into our results. Indeed, a global analysis might only reflect the statistics of high-latitude regions where the signal is oversampled in comparison to low latitudes. To check that hypothesis, we performed the following: for each contributing fluctuation, we apply a factor cos(θ ), θ being the latitude of the center of the fluctuations in order to artificially restore the homogeneity. This statistical correction reduces the weight of the contributions over a range of latitudes to compensate for the oversampling. It turned out that the final result is unaffected by this correction, rejecting the hypothesis of a bias due to an inhomogeneous density of measures. 

-By construction, all fluctuations have a privileged orientation (north-south) due to the orientation of alongtrack series. We can not exclude the possibility that a similar analysis would produce a different result with a different along-track direction. This would occur in the case of a strong anisotropy of topography at the planetary scale, outside the scope of this analysis.
section_name: Possible artifacts

section_content: Over nearly 4 orders of magnitude, beginning from around 10 km up to the planetary scale, the linear behavior of statistical moments up to order 2 is quite accurate, as shown in the right plot from Fig. 1 (the area between 4 and 7 in log10 units). The global scaling behavior is established over that range of scales, but the multiscaling properties (if they occur) still need to be tested. In order to apply the universal multifractal formalism, we evaluate the slope of each curve in Fig. 1 by computing linear fits. Figure 2 illustrates the results obtained on each curve. As one can see, the linear regression is rather satisfying. For scales > 10 km, each fit provides a scaling exponent that can be studied as a function ζ (q) (q being the statistical moment order) in order to test the validity of Eqs. ( 1 ) and ( 2 ). In Fig. 3 (red points), we have plotted the experimental structure function ζ (q). The apparent curvature of ζ (q) in Fig. 3 indicates multiscaling. Indeed, in the case of simple scaling behavior, the structure function is expected to be linear: ζ (q) = qH. Instead, the 21 red points in Fig. 3 can not be adjusted by a simple linear fit. 

The multifractal formalism may be tested by adjusting the experimental structure function ξ(q) with Eq. ( 2 ). In Fig. 3, the red lines show that the experimental structure is accurately adjusted, indicating that the universal multifractal for- malism is well suited to fitting the data with only three parameters H, α and C 1. Moreover, it allows us to estimate the values of the degree of multifractality α and the codimension of the mean C 1. Those estimates are presented in Table 1.
section_name: Scales > 10 km

section_content: Over the range of scales covering only 1 order of magnitude, the behavior is clearly different. Still, topography seems to exhibit scaling behavior over that range. Figure 2 shows the linear fit obtained, indicating that the log-log variations of moments according to scale is again satisfyingly linear. However, slopes are significantly steeper: fluctuations decrease faster when the scale decreases, indicating a larger H. The regression slope is computed for each statistical moment in order to study the experimental structure function. Figure 3 (blue points) shows the result. As one can see, the behavior differs from the one analyzed earlier. This time, the q dependence of the structure function is clearly linear, indicating that topography is not multiscaling in that range of scale.
section_name: Scale < 10 km

section_content: Our goal was to validate the accuracy of the universal multifractal formalism to describe the global scaling properties of the Martian surface in order to test whether the roughness intermittency is scaling. From our results, multiscaling seems to occur over a large but restricted range of scale (above 10 km). At a smaller scale, the topography is still scaling, but the symmetry is only monofractal, with a parameter H = 0. 75. 

This result is consistent with Aharonson et al. (2001), who obtained similar results by studying the power spectra of two Table 1. Estimates of the fractal and multifractal parameters for scales below and above 10 km. The last two columns were taken from Aharonson et al. (2001). Region 1 is a heavily cratered southern terrain. Region 2 is an area in the northern lowland. large significant regions of Mars (similar scale break). The power spectrum of a topographic transect is related to the second-order moment (square of Fourier transform) and is therefore expected to obey a scaling law if the global topography is scaling. Let β be the scaling exponent for the power spectrum. It can be shown that β is related to the structure function for q = 2. This relation is given by the following equation: Aharonson et al. (2001) found two scaling regimes similar to ours with a transition around 10 km. Our estimates of β for the small-and large-scale regimes are given in Table 1. Corresponding results from Aharonson et al. (2001) are also reproduced. For scales greater than 10 km, we found β = 1. 85 to be compared to the estimate of Aharonson et al. (2001) in region 1 (heavily cratered southern terrain) where β = 1. 4 and region 2 (area in the northern lowland) where β = 2. 0. Our global estimate of β = 1. 85 stands between the values obtained by the authors. As we found multiscaling properties in that range of scale (> 10 km), we can interpret this difference in terms of intermittency, i. e., local variations of β.
section_name: Conclusions

section_content: Indeed, intermittency is an expected feature of multifractal fields. For scales less than 10 km, our estimate of β differs clearly from the one obtained by Aharonson et al. (2001). Still, it might be interesting to notice that the authors found a similar value of β in both of these regions, indicating a low intermittency in that range of scale. This is consistent with our previous conclusion about the monofractal nature of scaling in that range of scale. At low scales, the measured value of H (0. 75) is close to the value obtained by Schmittbuhl (1995) in the case of fractured topography (H = 0. 8), indicating that this process may play a role in that range of scale. 

We demonstrate that a change in processes governing the Martian topography occurs at 10 km. The generic multifractal process is a multiplicative cascade. Such processes can reproduce the statistical behavior at scales greater than 10 km, but a simpler monofractal scaling process occurs at smaller scales. Craterization is well known to be a fractal process with a single fractal dimension (Rosenburg et al., 2015). We propose that the low scales are dominated by craterization processes, at the origin of the monofractal scaling law. This process has already been proposed by Aharonson et al. (2001). Future investigations are required to understand the multiplicative cascade processes at a large scale. 

Appendix A: Haar fluctuations for irregular signals
section_name: F. Landais et al.: Universal multifractal Martian topography

section_content: The above definition for Haar fluctuations implicitly requires that along-track points are regularly spaced. However, the actual MOLA data sets exhibit different kinds of irregularities, mainly due to the occurrence of clouds of various lengths and random dysfunction of the instrument that might introduce a bias into our analysis. Figure A1 provides typical examples of such irregularities in the MOLA along-track series and Fig. A2 shows the distribution of distance between two along-track points computed among the entire MOLA data set. It clearly indicates that most of the points are ∼ 300 m spaced, as expected, but a subsequent numbers of irregularities occur (one point missing, two points missing, etc. ). 

In order to take this into account, we define an adjustable quality criterion for the Haar fluctuations (see Lovejoy, 2014) : for a given fluctuation at scale L, we compute the spherical length of the first half L 1 and the spherical length of the second half L 2. In the case of perfect regularity, L 1 is exactly equal to L 2. If the measures are non-equally spaced, L 1 might be inferior or superior to L 2. The degree of irregularity may be estimated by a ratio R defined by Eq. (A2), equal to 0 if the measures are regular, and close to 1 2 in the case of extreme irregularity. 

The next step is to define a threshold to be applied to the ratio so that irregular fluctuations can be excluded from the global analysis. The threshold has to be chosen carefully: a very restrictive threshold might be damageable by excluding a lot of fluctuations. Indeed, a huge amount of fluctuations is needed for a statistical purpose. On the other hand, a nonrestrictive threshold might involve irrelevant fluctuations in the calculation. That may also introduce a bias.
section_name: A1 Irregularity

section_content: In this paragraph, we try to quantify how the choice of a threshold might affect the analysis. For that purpose, we have analyzed multifractal simulations after having manually introduced holes of different kinds into the data. To produce relevant statistics, 1000 universal multifractal series have been generated using the simulation technique called the fractionally integrated flux (FIF) developed by Schertzer and Lovejoy (1987) (http://www. physics. mcgill. ca/). Each series is composed of 100 000 equally spaced points. The universal scaling parameters used for this simulation are H = 0. 5, α = 1. 7 and C 1 = 0. 1. That particular choice relies on the assumption that scaling laws on Mars might be similar to the ones observed on Earth (see Gagnon et al., 2006). The first step is to introduce irregularities in order to reproduce the MOLA distribution of distance between consecutive points plotted in Fig. A2. This is obtained by suppressing data from the synthetic series. The accurate number of artificial holes of different kinds (one point missing, two points missing, etc. ) is adjusted by comparison with the actual distribution of holes plotted in Fig. A2. This figure compares the effective irregularity of MOLA (blue line) and the synthetic irregularity distribution manually introduced in simulations (red points). The kind and distribution of irregularities are realistic if red points coincide with maxima of the blue lines. As one can see in Fig. A2, the correspondence between simulated and actual data is satisfying. 

We can now compute Haar fluctuations on both biased and unbiased data and see how the choice of a threshold impacts the analysis of the biased series. We evaluate statistical moments of order q = 0. 5, q = 1, q = 1. 5 and q = 2 at different scales with different values of threshold R. Figure A3 shows how the biased moments deviate from the unbiased ones depending on the threshold. As one can see, the difference between blue dots and red dots in Fig. A3 is minimal, regardless of the threshold. Our conclusion is that in the case of a huge data set and given a relatively small density of missing points, the statistical multifractal properties are not significantly affected by high-ratio fluctuations. Hence, we are left to chose the least restrictive ratio in order to involve as many fluctuations as possible in the calculation and optimize the convergence. In the following analysis, the ratio is fixed at R = 0. 5 (completely unrestrictive). The total number of fluctuations computed to perform this analysis is over a billion that are spread into 74 bins of scales, equally spaced logarithmically from the largest scale possible (10 000 km, the half-circumference of the planet) down to the meter scale. Figure A4 shows the distribution of Haar fluctuations according to scale. Under 600 m (∼ 10 2. 7 ), the number of fluctuations decreases really fast, a consequence of the definition of Haar fluctuations.
section_name: A2 Choice of a threshold

section_content: . Statistical moments of order q = 0. 5, q = 1, q = 1. 5 and q = 2 computed for biased and unbiased synthetic data at different scales. From left to right, the threshold is R = 0, 0. 25, and 0. 5 from the more to the less restrictive. Blue points stand for the regular reference data. Red points represent irregular series (all ratios). Red points are satisfyingly centered inside blue points, regardless of the threshold. Therefore, the kind of irregularity that occurs in real data does not significantly impact the quality of Haar analysis.
section_name: Figure A3

section_content: www. nonlin-processes-geophys. net/22/713/2015/ F. Landais et al. : Universal multifractal Martian topography Acknowledgements. We acknowledge support from the Institut National des Sciences de l'Univers (INSU), the Centre National de la Recherche Scientifique (CNRS), the Centre National d'Etude Spatiale (CNES) and the Programme National de Planetologie. 

Edited by: L. Telesca Reviewed by: two anonymous referees