Universal multifractal Martian topography

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.


Introduction
The acquisition of altimetric data from 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 20 can be statistically described with quantitative parameters able to characterize the geologiccal units. Many statistical indicators have been proposed and widely explored in order to study the surface of Mars: RMS height, RMS slope, median slope (Kreslavsky and Head, 2000), autocorrelation length ) (see Shepard et al., 2001 for a review of all those indicators). Useful information has been 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) used 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 a region to the other, revealing a strong intermittency in the statistical variability of the surface of Mars 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, 15 2015), might lead one to reject the idea of a global description of any topographic field at the planetary scale. However, moderns developments in the theory of scale invariant processes might be able to a give full account to 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 20 of fractal dimensions (one for each altitude level). Multifractal simulations performed by 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 DEMs revealing that multifractal parameters may be able to discriminate between different stages in topographic 25 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 to model the variability of topographic fields and have the advantage of being simply 1009 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 satisfying fit of the structure function computed in different area of the globe and allowing 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 5 (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.

Methodology
This section first describes the Martian topography dataset 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 case of the irregularly sampled Martian topography.  Smith et al., 2001) is a laser altimeter on board of 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 measurement every two seconds, 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 dataset 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 millions. 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 acrosstrack fluctuations (measurements extracted from different orbits). This restriction might introduce a bias, 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.

10
The topography is a two-dimensional field providing for each pair of latitude/longitude a elevation value. In order to investigate the scaling properties of such field, 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 on 15 this study is largely based on Lovejoy and Schertzer (2013).
Once fluctuations are defined, a common way to explore 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-20 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 dependancy of ζ (q) is purely linear: ζ (q) = qH where H is called the conservation exponent used in honour of Edwin Hurst (typically a fractional brownian motion). A single coefficient (H) is sufficient to describe the global statistics, 25 that can be easily estimated by computing statistical moments of order 1 and 2 to obtained the slope of ζ (q) = qH. In the multifractal case, the scaling function is corrected by a additional quantity K (q), called moment scaling function: ζ (q) = qH − K (q).
Although three moments of distinct orders are in principle sufficient to establish the 5 curvature of ζ (q) and assess multiscaling properties, a field with general multifractal symmetries would require 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 estimate 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 behaviour in the exponent (and diverging as the sample size increases) . 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 15 three independent multifractal parameters introduced by Schertzer (1997): the mean fluctuation exponent H = ζ (1). 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 measuring the mean 20 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 estimate 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.
Consequence of the above definitions, any estimate of H, α and C 1 will rely on computing non-integers order moments around the mean q = 1. In the most general multifractal case, this set of three parameters only provides a restricted description 5 of statistics in the neighborhood of the mean (1st order moment) and presumes nothing of the variability that may be revealed by higher order moments. Still there exist a class of multifractal called the universal multifractals (Schertzer, 1997) whose associated structure function is fully determined by H, α and C 1 . In otherwords, for a member of this class, a set of only three parameters is enough to characterize all the 10 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 .

Haar fluctuations for irregular signals
Mars is well-known to present huge fluctuations of altitude on its surface despite of its 15 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 20 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 fractal field with H < 0 (a Gaussian white noise has H = −1/2), absolute differences will only characterize fluctuations with high wavenumber instead of the local fluctuations at a given scale. In the case of topography, Introduction 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 numbers 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 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 10 complicated to implement particularly in the case of 1-dimensional irregular series (see Appendix A). For the purpose of this analysis, we computed both 1st 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.

Results
For each scales, we compute moments of several orders q = 0.1, 0.2. . .2. The scaledependence 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 on Fig. 1 is obtained 20 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 on Fig. 1 (right plot) and their resemblance with 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 Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 one order of magnitude. Scaling needs to be studied separately for both of these ranges of scale.

5
Before going any further, we discuss in this section possible artefacts 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 on the topography measurement is of a few meters (areoid error + instrument error), lower than the 10 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 of 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 15 observed near the poles. This could introduced a bias in our results. Indeed, a global analysis might only reflect the statistics of high latitudes 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 20 the homogeneity. This statistical correction, reduce 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 inhomogenous density of measures.
-By construction, all fluctuations have a privileged orientation (Nord-South) due to Introduction would produce a different result with a different along-track direction. This would occur in case of a strong anisotropy of topography at the planetary scale, outside the scope of this analysis.

Scales < 10 km
Over nearly 4 orders of magnitude, beginning from around 10 km up to the planetary 5 scale, the linear behavior of statistical moments up to order 2 is quite accurate, as shown on the right plot from Fig. 1 (area between 4 and 7 on 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). On Fig. 3 (red points), we have plotted the experimental structure function ζ (q). The apparent curvature of ζ (q) on 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 on 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). On Fig. 3, the red lines shows that the experimental structure 20 is accurately adjusted, indicating that the universal multifractal formalism is well suited to fit 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

Scale > 10 km
Over the range of scales covering only one 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: 5 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 topographyis not multiscaling on 10 that range of scale. Indeed, the estimation of the universal multifractal parameters (see Table 1) leads to a value of C 1 close to 0 (Tabel 1).

Conclusions
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 if the roughness 15 intermittency is scaling. From our results, multiscaling seems to occur over a large but restricted range of scale (above 10 km) with. At 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 large significant regions of Mars (similare scale 20 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:

) 1017
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 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 > 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 5 (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 that range of scale (> 10 km), we can interpret this difference in terms of intermittency, i.e. local variations of β. 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 10 Aharonson et al. (2001). Still it might be interesting to notice that the authors found similar value of β in both of these regions, indicating a low intermittency on that range of scale. This is consistent with our previous conclusion about the monofractal nature of scaling on that range of scale.
We demonstrate that a change of processes governing the Martian topography 15 occurs at 10 km. The generic multifractal process is a multiplicative cascade. Such process can reproduce the statistical behaviour at scales > 10 km but a simpler monofractal scaling process is occurring a smaller scales. Craterisation 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 craterisation processes, at the origin of 20 the monofractal scaling law. This process has already been proposed by Aharonson et al. (2001). Future investigation are required to understand the multiplicative cascade processes at large scale.

10
In order to take this into account, we define a adjustable quality criterium on the Haar fluctuations (see Lovejoy, 2014): for a given fluctuation at scale L we compute the spherical length of 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 de L 2 . If the measures are nonequally spaced, L 1 might be inferior or superior to L 2 . The degree of irregularity may 15 be estimated by a ratio R defined by Eq. (A2), equal to 0 if the measures are regular and close to 1 2 in case of extreme irregularity.
The next step is to define a threshold to be applied on the ratio so that irregular fluctuations can be excluded of the global analysis. The threshold has to to be chosen

A2 Choice of a threshold
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 kind in the data. To produce relevant statistics, 1000 series of universal multifractal series have been generated using the simulation 5 technique called the Fractionally Integrated Flux (FIF) developed by Schertzer and Lovejoy (1987)  how the choice of a threshold impacts the analysis on the biased series. We evaluate statistical moment of order q = 0.5, q = 1, q = 1.5 and q = 2 at different scales with different values of threshold R. Figure 6 shows how the biased moment deviate from the unbiased ones depending on the threshold. As one can seen, the difference between blue dots and red dots on Fig. 6 Fig. 2. Red points (resp. blue points) correspond to the range of scales larger (resp. smaller) than 10 km.  Clouds ~50km Proportion of occurrences regular distribution Figure 5. Comparison between effective irregularities in the MOLA database and artificial irregularities on simulations. Irregularities are expressed as a proportion of occurrences from the maximum possible in log10 space (0 means all fluctuations are observed). The blue line represents the effective irregularity of actual data. Red points are obtained from the synthetic irregularities on multifractal simulations. The good agreement between maxima of the blue curve and red points indicates that effective irregularity is accurately reproduced by simulations.