ABSTRACT

The so-called “typical” values for the Earth’s topography are often used in the literature, such as the mean continental altitude (MCA), the Moho depth for “normal” continental crust, or “typical” depth of mid-oceanic ridges. However, the statistical relevance of those values is hardly discussed. Focussed on data for the global topography, this paper presents statistical analyses regarding various environments. It is shown in particular that the definition of the mid-oceanic ridge is not straightforward, and varies considerably according to what is actually considered: the ridge “inner-rift”, the ridge “crest”, or the “virtual ridge” at spreading centre. This definition is also a function of the spreading rate and has strong implications for the rationale on the age-depth relationship of the sea-floor. In addition, the latter relationship is highly dependent on how the topographic data are corrected from sediment load. The correction itself implies numerous aspects that relies on the precision and associated uncertainties of, in particular, the sediment thickness, sediment porosity, and the mantle, water and sediment densities. In this respect, the analysis carried out here favours a plate cooling model (PCM) for the age- depth dependence of the sea-floor. The topographic elevation at trench proves also to be related to the age of the sea-floor through a different PCM equation. Away from the trench, the oceanic lithosphere is affected by flexuration, for which equations can be defined assuming that the end-load position is not located at trench. On the other hand, the elevation of magmatic arc does not appear to be related to sea-floor age or spreading rate. However, the correlation between the arc-trench distance and the topographic elevation of arc for continental crust seems to be an indicator of slab dip and therefore the existence of slab roll-back processes. Along intra-oceanic magmatic arc, a periodicity in topographic elevation suggests a periodicity in the occurrence of magma chambers, and therefore magmatic processes that need to be further studied. At passive margin, the transition between continental and oceanic crust seems to be relatively sharp in average. Subdivision of the datasets according to the age of the continent-ocean boundary (COB) indicates that rift and passive margin shoulders are found within a couple of degree away from the COB and for ages younger than ca. 20 Ma. Finally, the statistical analysis of continental data assumed to be free of thinning or thickening effects suggests that the MCA should rather be considered in terms of “lowlands” and “highlands”. Relying on model of Moho depth, the “normal” crustal depth might be thinner than commonly accepted. In any case, the filtering of reduced topography can help to determine the impact of dynamic topography.

Subject Areas:

Environmental Sciences, Geology

Keywords:

Statistics, Topography, Tectonic Environments, Earth’s Crust

1. Introduction

The present-day topography is the sole example of topography that we can access to. It is consequently used as reference when addressing various issues about the topography in the geological past. Therefore, it is of prime importance to properly characterize this topography on the global scale.

A number of studies in Geosciences use and/or discuss data about the topography of the Earth at global scale, most of the time with implications for palaeogeographies (e.g. [1] [2] [3] [4] ). Such data may include mean bathymetric values for abyssal plains or mean continental altitude, mean depth of mid-oceanic ridges, or mean distance between trench and volcanic arc, etc. However, the significance of such values is hardly discussed in the literature, and in particular whether the mean value (μ) and associated standard deviation (σ) have a real meaning for those data.

A series of statistical analysis about the present-day topography of the Earth is provided here and discussed for the following environments: mid-oceanic sp- reading centres, depth-age dependence of the oceanic lithosphere, active margin settings and intra-oceanic subduction zone, passive margin settings, and mean continental altitude (MCA).

The aim is to clarify the relationships between physical processes and resulting topography in order to decipher what parameters can be used in assessing topographies in deep time.

2. Data Sources

A statistical analysis is highly dependent of the dataset used. For this study, statistics were carried out on the following datasets:

・ Global Relief Model Etopo1 [5] ; Etopo1 is a 1 arc-minute (ca. 1.85 km horizontal resolution) global relief model. The vertical resolution is not published, but is estimated to be 1% of the cell grid elevation at best (Eakins, pers. com., 2010). For sake of simplicity, the data uncertainty is, herein, arbitrary assumed to correspond to 1% of the given altitude with a minimum uncertainty of 10 m.

・ Global sea-floor age model and associated uncertainties [6] .

・ Global sediment thickness map [7] (updated online versions, 2010; hereafter referred as to [8] ); the map grid has 2 × 2 degrees horizontal resolution, and the vertical resolution is not provided. The sediment thickness of the NOAA ( [9] ; 5 arc-minute (ca. 9.27 km) horizontal resolution) is also used, but it covers only oceanic realm with some data gaps, in particular in polar area. The vertical resolution is arbitrary assumed to equals the difference between the two datasets (see below) with minimum uncertainty of 250 m, or 5% of the sediment thickness given by Laske & Masters [8] where data from the NOAA are missing.

Given the uncertainties associated with the datasets used, all statistics were performed assuming a spherical Earth with radius R_Mean:

$R_{\text{Mean}}=\sqrt[3]{{\left(R_{Eq.}\right)}^2\times R_{Pol.}}=\sqrt[3]{{\left(6378.137\right)}^2\times 6356.752}=6371.001\text{km}$ (1)

where, R_eq. is the equatorial radius and R_pol. the polar radius of the ellipsoid of reference (WGS84). In the following, the Earth’s surface is therefore 5.101 × 10⁸ km² and the volume, 1.083 × 10¹² km³; one arc-degree corresponds to 111.195 km. In addition, because of the various spatial resolution of the aforementioned datasets, resampling with linear interpolation between original grid cells were carried out using a geodetic grid in order not to overestimate polar regions (Figure 1(a)). This geodetic grid contains over 2.6 × 10⁶ data points, which correspond approximately to one data point every 16 km world-wide. The various examined quantities (topographic elevation, sediment thickness, etc.) are therefore taken at the exact same location.

Data are separated between points lying on crust assumed to be continental in nature from those lying on crust assumed to be continental in nature. The definition of the Continent-Ocean Boundaries have been made by hand by following “at best” both prominent features in the first derivative of the map of the Earth gravity field model (Grace GGM02, [10] ) and the second derivative in the magnetic anomaly map EMag2 [11] . The result is however not fundamentally different from other publications as reported by Eagles et al. [12] , and does not significantly affect the statistical outcomes provided below.

As first put forward by Otto Krümmel [13] , the distribution of elevation at global scale is predominantly bimodal (Figure 1(b)) ranging here from −10726 m to +7446 m (Figure 1(c)), so that neither the mean (μ) value of −2432 m nor the median (m) value of −3280 m has a real meaning.

Separately, statistics show that neither the data points assumed to lie on continental crust (brown dots in Figure 1(a)) nor those lying on assumed oceanic crust (blue dots in Figure 1(a)) are Gaussian distributed (Table 1; Figure 1(b)). Although one may consider that the median values (respectively m_oceanic = −4325 m and m_continental = +187 m) are more representative or more robust than the mean values, none of them have a correct meaning either and it might be useful to consider the peak of the data distribution (here, peak_oceanic = 6.268 × 10⁴ data

in the bin [−4300 m; −4400 m] and peak_continental = 1.108 × 10⁵ data in the bin [+0 m; +100 m], respectively). However, oceanic and continental distributions are relatively well-symmetric to first order, and the differences in elevation according to what values are regarded as more representative are relatively small (difference between the mean and median is Δ(μ − m) = 136 m for oceanic dataset and Δ(μ − m) = 30 m for continental dataset; Table 1). Notwithstanding, uncer-

tainties based on standard deviation are useless.

Furthermore, mean values of topographic features may be biased by some outliers related to distinct features. For instance, the determination of the “best” sea-floor age-depth relationship may be considered as biased by the presence of oceanic plateaus or seamounts, deep fracture zones, abandoned arcs, … In order to exclude such potential outliers, a buffer zone has been defined (Figure 2) and includes: main volcanoes and plateaus associated with hot-spot magmatism and a buffer zone of 1.5 around them to account for lithospheric flexuration due to loading, abandoned arcs with buffer of 1, active intra-oceanic arcs with buffer starting at trench and ending 1.5 behind arcs, oceanic area closer than 1000 km (buffer of 9) from any subduction zone to account for lithospheric flexuration (which might reach kilometre scale elevation above reference depth; e.g. [14] ), and area affected by present-day ice loading and/or post-glacial rebound (the latter has been defined as area where post-glacial rebound exceeds ±1 mm/yr after the model of Paulson et al. [15] [16] ).

3. Oceanic Spreading Centres

The often called “mid-oceanic ridges” are important features of the topography of the Earth. They indeed represent a cumulative length of some 65,000 km (e.g. https://en.wikipedia.org/wiki/Mid-ocean_ridge). However, because this feature is fractal (as are length of coast lines for instance), the cumulative length depends on the resolution of the dataset. Using the present resolution and definition (Figure 2), the cumulative length is 1715.955 or 190,805.572 km, among which 61% (116,396.949 km) correspond to “true” spreading centres (or segments) and 39% (74,408.623 km) correspond to transform faults. It must be noted in addition that 7.6% (or 14,459.382 km of ridges, i.e. including both spreading centres and transforms) belong to back-arc basin (Figure 2).

The definition of what is a spreading centre and what is a transform fault is not as straightforward as one may first think. It is merely considered herein that a segment of ridge with direction within ±45 from the orientation of the motion vector between the two adjacent tectonic plates is a transform fault, whereas others are spreading centres (see sketch in Figure 3(a)). The distribution of the orientation of segments of ridge are shown in Figure 3(b), where segments corresponding to spreading centres are centred around 90 and transform fault centred around 0 and 180 relative to motion vectors of the different tectonic plates.

The use of ETopo1 [5] allows to define the location of mid-oceanic ridges quite precisely. The location is most often marked by the presence of an inner rift, surrounded by two ridge crests, and two ridge flanks (Figure 4; see also

[17] ). The definition of the elevation of the ridge is therefore very dependent of the features considered: elevation of the inner rifts of the ridges, elevation of the ridge crests (highest points), or elevation of the “virtual” ridges corresponding to the point where linear fits on ridge flanks cross each other (red lines in Figure 4). The value chosen has important implications, since it controls the entire definition of the sea-floor elevation when using age-depth relationship (e.g. [14] [18] [19] ).

3.1. Mid-Oceanic Ridge Inner Rifts

The mean value of all spreading centres at the location of inner rifts is μ = −3250.780 m. However, even if there is a maximum of data close to the median value (m= −3231 m), statistics show that the data are not Gaussian distributed (Figure 5(a); Table 2). It must be noticed, in addition, that transform faults are,

in general, about 300 m - 350 m deeper than inner rifts of spreading centres (344.5 m deeper according to mean values, and 303 m according to median values; Table 2), so it is relevant to exclude them from statistics on spreading centre topography.

The age of the sea-floor at spreading centres shall be, by definition, 0 million years (Ma). Because the definition of ridges by Müller et al. [6] was coarser than the definition made here (see Figure 4), the age of the sea-floor determined from Müller et al.’s dataset [6] is about 3 Ma (2.93 Ma after the mean value and 3.04 Ma after the median value both determined on the logarithmic distribution of sea-floor ages (using all spreading centres only); yellow histogram in Figure 5(b)). The age uncertainties calculated by Müller et al. [6] accounts for uncertainties on fitting pairs of isochrones. The mean age uncertainty of 0.5 Ma (median = 0.78 Ma) obtained at ridges using Müller et al.’s dataset [6] therefore most likely underestimate the age error when uncertainties on location are taken into account (Figure 5(c)). Now, uncertainties remain low and this should not fundamentally affect results on the age-depth dependence of the sea-floor (see below).

It is commonly accepted that slow spreading ridges have deep inner rifts whereas fast ridges have not. The spreading rate is therefore an important parameter, but as for sea-floor ages, the resolution of the dataset may impact results. The amount of accreted surface of sea-floor per year (in mm²/yr) has been calculated from Euler poles and angles at each spreading centre segment (Figure 5(d)). Knowing every segment lengths, the results have been converted in spreading rates (in mm/yr) and compared with values provided by Müller et al. [6] (Figure 5(e) & Figure 5(f)).

The topographic elevation of ridge inner rifts as function of spreading rates has thus been tested (Figure 6) using both the spreading rates computed by Müller et al. ( [6] ; hereafter termed data. A) and the accretion rates computed here (this study; hereafter termed data. B).

The mid-oceanic ridges are indeed deeper in average with slow accretion rates, but the dispersion of data is maximum for slow accretion rates (even when potential outliers are excluded using the buffer zones defined in Figure 2), and the rise of the inner ridge elevation seems to reach a threshold around 1 mm²/yr (Figure 6(a)) or around 30 mm/yr (Figure 6(b)) whatever the dataset used (i.e. data. A or data. B). Albeit similar relationship using data. A or data. B, the rise of the inner rift elevation is not linear with the spreading rate (Figure 6(c)). The mathematical relationship between elevation and rate is undetermined, although

one may think a power law might rule the system.

3.2. Mid-Oceanic Ridge from Ridge Crests, Ridge Flanks and “Virtual Ridge Axes”

Statistics on the elevation of ridges are considered here as more relevant from analysis of ridge crests, ridge flanks and “virtual ridge axes” than from analysis of ridge inner rifts. For such analysis, data were first selected within the area accreted since chron C3 (ca. 6 Ma; Figure 7(a) & Figure 7(b)). However, the amount of data is rather large, and a good picture of mid-oceanic ridge shape and elevation can be obtained with a dataset limited to 1.333 around ridge axes (Figure 7(c)).

In order to have a clearer picture, the distance to ridge axis has been divided into 0.025 bins. Although the distributions per bin are not Gaussian from a statistical point-of-view, they are relatively symmetric and mean and median values are relatively close (generally below 100 m in difference; Annexe 1 in supplementary data). Mean values per bin are displayed in Figure 8, and show a nice general shape with the inner rift around −3100 m (μ = −3120.387 m), the ridge crest distant of 0.15 (ca. 17 km) away from the axis with elevation around −2650 m (μ = −2661.473 m), and the ridge flank gently dipping away. In the ridge flank, two parts can be distinguished: the nearest part―below 0.80 in distance from ridge axis (yellow dashed line in Figure 8)―with a steeper slope, and the farthest part―beyond 0.75 (yellow dash-dotted line in Figure 8)―with a shallower slope. The “virtual ridge axis” can therefore have different elevation according to the slope chosen to represent the ridge flank. The “virtual ridge” elevation―the point at which the best linear fit on flank data crosses the ridge axis (Figure 8)―is −2601.615 m using the entire ridge flank, −2559.811 m using near section

of the flank, and −2695.038 m using data of the flank from the farthest section, respectively.

However, we previously saw (Figure 6) that the rate of sea-floor accretion is an important parameter in the elevation of the mid-oceanic ridge. The ridge dataset (data within 1.333 away from ridge axis) has been therefore divided into sub-datasets according to the accretion rates (bins of 0.5 mm²/yr). Using running averages (50 data sliding window), the relationship of the elevation of the inner rifts to accretion rates becomes clear (Figure 9(a)). In addition, one can see that the ridge crests are also impacted by the accretion rates within a relatively constant segment length ranging from 0.15 to 0.70 - 0.80 from ridge axis, i.e. the nearest section of the ridge flank, hereafter named “near flank”. Beyond 0.80, elevations and slopes are relatively similar; the cause of the observed variations is undetermined but might be due to effects of the geoid and/or the dynamic topography since data per accretion bin are not evenly distributed around the world. Consequently, linear fits have been established on the three sections (inner rift, near flank, far flank) and normalised to the mean value within the far flank section (at 1.03 away from ridge axis) for comparison purpose (Figure 9(b)). Ridge crests is prominent for accretion rates lower than ca. 1.5 mm²/yr and associated with troughs at ridge axis (inner rift), whereas ridge crest are mostly lower (deeper in bathymetry) than ridge axis elevation when rates are higher. Those relationships between topographic elevations, distances to ridge axis and accretion rates are also well perceptible in Figure 9(c) & Figure 9(d).

If we now consider that the elevation of the ridge axis is not the true elevation affected by processes acting in inner rift but a virtual elevation resulting from the extrapolation of linear fitting on ridge flank (as the red line of inset of Figure 4), then the “virtual ridge elevation” may vary quite considerably, and values are summarised in Table 3.

In conclusion, the depth definition of the often called “mid-oceanic ridges” is subject to caution, since values may vary quite considerably according to what features are taken into account. This issue has a particularly strong implication in the definition of an age-depth dependence of the entire sea-floor.

4. Global Sea-Floor and Age-Depth Dependence

The realization that sea-floor topography (bathymetry) is highest at mid-oceanic ridges and decreases with distance is at least as old as the advent of Plate Tectonics (e.g. [20] [21] [22] [23] ). The question remains open as to whether the age-depth dependence of the sea-floor is best described by a square root mathematical relationship (representing a Half-Space Cooling Model: HSCM) or an exponential relationship (Plate Cooling Model: PCM). Stein & Stein [18] defined the “GDH1” model (“Global Depth and Heat flow” model) which combined the two types of model (HSCM and PCM) on the basis of an arbitrary dichotomy of data. The cut-off values chosen by these authors is however troublesome because they are different for topographic and heat flow data (the distinction between young-aged sea-floor versus old-aged sea-floor is chosen to be 20 Ma for topog-

raphic data and 55 Ma for heat flow data) and their physical meaning remain obscure.

In addition, the equations for the GDH1 model are of the form:

$d\left(t\right)=d_R+f\cdot \sqrt{t}$ , when t < 20 Ma (i.e. young-aged sea-floor) (2.1)

$d\left(t\right)=d_R+\frac{\alpha {\rho }_M\left(T_M-T_0\right)Z_L}{2\left({\rho }_M-{\rho }_W\right)}\left(1-\frac{8}{{\text{π}}^2}\exp \left(-\frac{\kappa {\text{π}}^2}{Z_L^2}t\right)\right)$ , when t ≥ 20 Ma (i.e.

old-aged sea-floor) (2.2)

where,

d(t) is the sea-floor topography (depth in m) as function of age t (in second but usually converted into Ma with a conversion factor τ = 3.15576.10⁺¹³ s per Ma),

d_R is the depth at the mid-oceanic ridge (in m),

f is an ad hoc factor (dimensionless, hereafter shown with the sign Ø),

α is the volume coefficient of thermal expansion (in 1/kelvin or K⁻¹),

κ is the thermal diffusivity (in mm²・s⁻¹),

ρ_M and ρ_W are the mantle and water densities respectively (or mass per unit volume in kg・m⁻³),

Z_L is the asymptotic thermal plate thickness (i.e. at infinite distance from the ridge, in m),

T_M and T₀ are respectively the basal and the top temperature of the thermal plate (in K).

It must be noticed that both equations are dependent upon d_R, which we saw, may vary in definition on the one hand, and varies as function of spreading rates on the other hand.

Furthermore, those equations hold for thermal plates. On oceanic sea-floor however, the crust is covered by sediments, and isostatic calculation shall first be carried out to correct for loading effects. The estimate for sediment thicknesses is subject to caution and the method for correction is not as straightforward as commonly thought.

4.1. Isostatic Correction for Sediment Load

The isostatic correction for sediment load is crucial to define properly the age- depth dependence of the sea-floor. Although the Veining-Meinesz method may be regarded as a more comprehensive technique for accounting for isostatic correction, the method is complex (e.g. [23] for solution in inverse problem) and requires using datasets for which uncertainties are probably as large as the problem we want to address here. The Airy-Pratt method for isostasy is therefore regarded here as more adapted to the present concern. The isostatic correction C for sediment load is simply expressed as follow:

$C=\left(\frac{{\rho }_W-{\rho }_S}{{\rho }_W-{\rho }_M}\right)\cdot h_S$ (3)

The correction thus only depends on the mean sea water density ρ_W, the mean sediment density ρ_S, and the mean mantle density above compensation depth ρ_M.

4.1.1. Sea Water Density ρ_W

The sea water density can be computed from the International Equation of State of Seawater [24] . The one-page-long equation is function of sea water temperature, salinity and pressure. The ranges of sea water temperatures and salinity as function of depth are shown in Figure 10 (given we assume herein a constant atmospheric pressure of 101325 Pa). Using bootstrapping technique (i.e. random resampling within bounds leading to Gaussian distribution of sub-datasets), mean temperature and salinity profiles at global scale have been defined (Figure 10(a) & Figure 10(b)) and used to compute the range and mean values of sea water density as function of depth (purple in Figure 10(c)).

The sea water density ρ_W required for the isostatic correction C corresponds to the mean value over the water column, i.e. corresponds to the integration of sea water densities throughout every water columns. The mean sea water density ρ_W is shown as function of the height of the sea water column in red in Figure 10(c)

within bounds in deep yellow.

4.1.2. Sediment Density ρ_S

The density of the sediment cover at global scale is obviously a much more challenging issue. Multiple factors impact the sediment density, the two most important of which are certainly the lithology and the compaction. Winterbourne et al. [25] have proposed that the mean sediment density ρ_S over the sediment pile can be calculated using:

${\rho }_S=\frac{1}{h_S}{\displaystyle {\int }_0^{h_S}\left(\varphi {\rho }_{igW}+\left(1-\varphi \right){\rho }_g\right)\text{d}z}$ (4.1)

with,

$\varphi ={\varphi }_0{\text{e}}^{\left(-z/\lambda \right)}$ after Athy [26] (4.2)

where,

h_S is the height of the sediment pile or sediment thickness,

ρ_igW is the intergranular water density or water density in sediment pores,

ρ_g is the density of the solid sediment grains,

φ is the porosity as function of height of sediment,

φ₀ is the initial porosity,

z is the compaction decay wavelength,

z is the depth within the sediment pile.

Assuming ρ_igW = 1013 kg・m⁻³, ρ_g = 2650 kg・m⁻³ and ρ_M = 3300 kg・m⁻³, Winterbourne et al. [25] inverted two-way travel times (TWTT) from seismic reflection and wide-angle profiles from the Atlantic Ocean and found best values (i.e. minimum misfit between predicted and observed relationship between TWTT and depth) for initial porosity φ₀ = 0.56 and compaction decay length λ = 4.5 km.

Audet & Fowler [27] in particular, have shown the complexity of the relationship between porosity and sediment thickness. The relationship put forward by Athy [26] is therefore an over-simplification which implies that the mean sediment density is approximated with substantial uncertainty. Moreover, the present-day global mean value for sea water density is ρ_W = 1027.6 kg・m⁻³ (for T = 4 C and S = 34.7‰; P_atm = 101,325 Pa), which is higher than the 1013 kg・m⁻³ used by Winterbourne et al. [25] . Because water located in sediment pore is undoubtedly subject to higher pressure, temperature and salinity, the intergranular water density ρ_igW must certainly be even higher.

Now, for sake of simplicity, the present study nonetheless follows the equations of Winterbourne et al. [25] because the relationship between the TWTT and depth found by those authors is considered here to be sufficiently well-rep- resentative of the change in sediment density for isostatic correction. Bootstrap resampling, however, has been carried out to refine mean values and assess the ranges of uncertainties around them (i.e. within ± 2σ). Doing so, the following parameters were used:

${\phi }_0=0.5616\pm 0.1003$ [Ø]

$\lambda =4578.6\pm 1112.3\left[\text{m}\right]$

In Winterbourne et al. [25] , the density of solid grain is taken to be ρ_g = 2650 ± 250 [kg・m⁻³]. It is found here (see below, inset in Figure 12) that a value of ρ_g = 2300 ± 500 [kg・m⁻³] (“best” value of ρ_g = 2292 [kg・m⁻³]) both improve the correction and is more coherent with the diversity of grains in sediment. The intergranular water density ρ_igW is assumed to merely follow the International Equation of State of Seawater (purple curve and associated pink uncertainties in Figure 10(c)).

4.1.3. Mantle Density ρ_M and Correction Factor C

The mean mantle density above compensation depth is just defined as ρ_M = 3150 ± 300 [kg・m⁻³].

In order to illustrate the amount of correction and associated uncertainties those values correspond to, the sediment density as function of sediment pile and the corresponding correction factor C are defined in Figure 11 assuming a fixed water depth of h_W = −4325 m (mean global ocean depth).

4.2. General Age-Depth Dependence Using Topographic Elevation Data Corrected for Sediment Load

The topographic dataset used for age-depth relationship corresponds to all oceanic data (after ETopo1 [5] , i.e. blue dots in Figure 1) that are not excluded by buffer zones as shown in Figure 2 plotted against ages provided by Müller et al. [6] . In other words, all oceanic data not affected by known main bias were used and corrected from sediment load as described above.

Because the number of data points is large (N = 798,458), age bins of 1 Ma have been used to obtain subdatasets on which statistics have been determined.

Although none of the distributions per bin is properly Gaussian from statistical point-of-view, all are relatively well-symmetrical, “bell-shaped” with median values close to the mean (Figure 12). Actually, most distributions are closer to Laplace distributions implying that the mean values μ are relevant and representative of the distributions whereas the standard deviations σ simply overestimate the dispersions.

Consequently, the mean topographic elevation μ per 1 Ma bin with associated 2σ error bars are plotted against the age of the sea-floor in Figure 13. The sediment load correction appears satisfactory when compared with the average raw topographic elevation from ETopo1 since the old part exhibit a quasi-flat relationship as expected.

It must be noticed that topographic elevations accounting for uncertainties (i.e. uncertainties on the raw datasets augmented by uncertainties on sediment load correction as described above) mainly increase the dispersion around the mean values defined per 1 Ma bin, but the mean values (μ− and μ+) are not themselves drastically shifted up or down (Figure 13). Another method―not carried out here―to define the best parameters for isostatic correction (namely ρ_igW, ρ_g, ρ_M, φ₀, λ) might therefore consist in minimizing the dispersion around the mean in every 1 Ma bin. Nevertheless, the difference between μ− and μ+ has been used to search for the “best” value ρ_g, the density of solid grain in sediment, which is required in Equation (4.1). The cumulative difference between μ− and μ+ is minimal for ρ_g = 2292 kg・m⁻³ (Inset in Figure 12). Because the uncertainty associated with this value would need more proper quantification, the density of solid grain in sediment is, here, just arbitrarily chosen to be ρ_g = 2300 ± 500

kg・m⁻³.

A square-root mathematical equation (HSCM) can fit the relationship between depth (μ) and age within 2σ error bars―and even more easily within bounds when uncertainties are accounted for; i.e. (μ+) + 2σ and (μ−) − 2σ―but the mean (or median) values prove to be badly fitted. The “best” HSCM (shown in pink in Figure 13) has the following equation:

$d\left(t\right)=X+Y\cdot \sqrt{t}=-3425.880-186.594\times \sqrt{t}$ , (5.1)

with coefficient of determination R² = 88.3%.

As emphasized by previous authors, the fit is much better when applied on data for the first 70 Ma (e.g. [18] , although those authors finally chose 20 Ma as limit in their GDH1 model) but the physical meaning of splitting the data between young-aged and old-aged data is unresolved. Using a fit for data younger than 70 Ma only, the HSCM equation (purple in Figure 13) becomes:

$d\left(t\right)=X+Y\cdot \sqrt{t}=-2758.431-30173\times \sqrt{t}$ , (5.2)

with coefficient of determination R² = 99.5% (on data younger than 70 Ma).

The fit using an exponential equation (PCM) is much better. Excluding data younger than 3 Ma because of the presence of the mid-oceanic inner rift (see Section 3.1), the PCM equation (green in Figure 13) is:

$d\left(t\right)=A+B\cdot \exp \left(-C\cdot t\right)=-5632.290+2527.251\times \exp \left(-0.02554\times t\right),$ (6)

with coefficient of determination R² = 98.9%.

4.3. Age-Depth Dependence of the Sea-Floor as Function of Spreading Rate

Instead of dividing the topographic elevation data according to their age (age bins), the dataset can be divided relative to the spreading rates (rate bins) as defined by Müller et al. [6] . PCM equations of the form $d\left(t\right)=A+B\cdot \exp \left(-C\cdot t\right)$ were thus fitted to every sub-datasets and Figure 14 shows the variation of coefficient A, B and C as function of spreading rates (see also Table 4). Polynomial fit of degree 3 (and associated 95% confidence intervals) are merely shown to highlight the first order variation of each coefficient. Beyond 100 - 120 mm・yr⁻¹, however, the number of data per sub-datasets becomes small, and the increase in coefficient B and C suggested by polynomial fit is likely an artefact linked to the natural shape of 3^rd degrees equations. It is suspected that the value of the different coefficients stabilize within the 95% confidence level as spreading rates keep on growing.

Besides the coefficients correspond to:

$A=d_R+\frac{\alpha {\rho }_M\left(T_M-T_0\right)Z_L}{2\left({\rho }_M-{\rho }_W\right)};B=-\frac{\alpha {\rho }_M\left(T_M-T_0\right)Z_L}{2\left({\rho }_M-{\rho }_W\right)}.\frac{8}{{\text{π}}^2};C=-\frac{\kappa {\text{π}}^2}{Z_L^2}$ (7)

in the equation provided by Stein & Stein ( [18] ; see Equation (2)), where d_R corresponds to the “virtual” ridge axis (see analysis of the mid-oceanic ridges above;

Section 3.2) and is given as d_R = 2600 m.

The definition of d_R from the coefficients determined in Figure 14 and Table 4 is simply $d_R=A+B\times \left({\text{π}}^2/8\right)$ . For the other parameters however, the determination is not trivial. Stein & Stein [18] , in particular, noticed the difficulty to quantify uncertainties.

To do so, the genuine equation for Plate Cooling Model (PCM; [14] [22] [28] ) was used and is written as follow:

$d\left(t\right)=\frac{{\rho }_M\alpha \left(T_M-T_0\right)Z_L}{2\left({\rho }_M-{\rho }_W\right)}\left(1-\frac{4}{\text{π}}{\displaystyle {\int }_0^1{\displaystyle {\sum }_{n=1}^{\infty }\frac{1}{n}\exp \left(-\frac{\kappa n^2{\text{π}}^2}{Z_L^2}t\right)\sin \left(n\text{π}Z\right)\text{d}Z}}\right)$ (8.1)

After evaluation of the integral [14] , the equation can be re-written as:

$d\left(t\right)=\frac{{\rho }_M\alpha \left(T_M-T_0\right)Z_L}{2\left({\rho }_M-{\rho }_W\right)}\left(\frac{1}{2}-\frac{4}{{\text{π}}^2}{\displaystyle {\sum }_{m=0}^{\infty }\frac{1}{{\left(1+2m\right)}^2}\exp \left(-\frac{\kappa n^2{\text{π}}^2}{Z_L^2}t\right)}\right)$ (8.2)

The Bayesian-Markov chain Monte-Carlo inversion method employed by Scholer [29] [30] in another context was used here to quantify the unknown parameters, namely ρ_M, α, T_M (T₀ is assumed to be the temperature defined at the sea-floor as in Figure 10), κ, and Z_L (ρ_W being equally defined as above).

This stochastic inversion simply consists in randomly choosing values for each parameter within given bounds, and in only accepting resulting curves from Equation (8.2) that fit within defined limits.

Here, the parameters were chosen with rather large and conservative ranges as follow:

${\rho }_M\in \left[0;5000\right]$ (in kg・m⁻³) ; mean density of the mantle above compensation depth, which value is commonly chosen around 3300 kg・m⁻³ in the literature.

$\alpha \in \left[0;10\times {10}^{-5}\right]$ (in K⁻¹) ; volume coefficient of thermal expansion, which value is commonly chosen around 3 × 10⁻⁵ K⁻¹ in the literature.

$T_M\in \left[0;5000\right]$ (in K); Temperature at compensation depth, which value is commonly chosen around 1625 K in the literature, so that $\left(T_M-T_0\right)\simeq 1350˚\text{C}$ .

$Z_L\in \left[0;200000\right]$ (in m); Thermal plate thickness, which value is commonly chosen around 125 km in the literature.

$\kappa \in \left[0;10\right]$ (in mm²・s⁻¹ or × 10⁻⁶ m²・s⁻¹) ; Thermal diffusivity, which value is commonly chosen around 1 mm².s^-1 in the literature.

Note that in Equation (8.2), the sum is carried out up to m = 10,000.

For sake of clarity, Figure 15 and Figure 16 shows the outcome for all the data together, i.e. not separated by bin of spreading rates.

The acceptance of the resulting curves was defined in two ways: 1) curves that fit within two standard errors around the mean values per bin of age (μ ± 2σ), and 2) curves that fit within 95% of the number of data between the minimum and maximum values for every bin that are the closest to the median (m ± δ_95%). The second method has the advantage to discard the main outliers (5% of the subdatasets) and to be independent of the data distribution. However the range of acceptable curves around the median is larger (Figure 15) than the range around the mean. In order to avoid being limited and biased by minimal values of 2σ and δ_95% respectively, bounds have been smoothed using polynomial functions. The simulation stops after 5000 curves fit within bounds. Because the possibility for simulated curves to fit within bounds is smaller around the mean than around the median, the number of iterations was much larger to obtain the 5000 curves using (μ ± 2σ) than using (m ± δ_95%) (N_iter = 184,603 and N_iter =

26,888 respectively).

The distributions of potential values for the aforementioned parameters (Figure 16) are neither Gaussian distributed nor clearly symmetrical. It seems therefore that many possibilities exist to combine the values for those parameters and no “best” solution satisfactorily stands out. The conclusion is the same when the dataset is divided per bin of spreading rates.

In order to see the effect of limiting the range of possibilities, the same computation has been carried out using the upper and lower bounds corresponding to the average mean topographic data accounting for uncertainties (μ+ and μ− as in Figure 13, and smoothed using polynomial fits) and using more “realistic” ranges for the parameters:

${\rho }_M=3150\pm 300$ , i.e. ${\rho }_M\in \left[2850;3450\right]$ (in kg・m⁻³); Mean density of the mantle above compensation depth,

$\alpha =3\times {10}^{-5}\pm 4\times {10}^{-5}$ , i.e. $\alpha \in \left[0;7\times {10}^{-5}\right]$ (in K⁻¹); Volume coefficient of thermal expansion,

$\Delta T=1350\pm 350$ , i.e. $\Delta T\in \left[1000;1700\right]$ (in K); Difference in temperature between the base and the top of the thermal plate,

$Z_L=125000\pm 50000$ , i.e. $Z_L\in \left[75000;175000\right]$ (in m); Thermal plate thickness,

$\kappa =1\pm 5$ , i.e. $\kappa \in \left[0;6\right]$ (in mm²・s⁻¹); Thermal diffusivity.

The number of iteration has to reach 2,728,362 to obtain only 100 accepted curves (Figure 17(a)).

Similarly, the outcome for each parameter does not allow clearly defining favoured values (Figure 17). The following mean μ (parameter) and median m (parameter) values are therefore provided for information only:

$\mu \left({\rho }_M\right)=3144.930\text{kg}\cdot {\text{m}}^{-3}\text{and}m\left({\rho }_M\right)=3154.485\text{kg}\cdot {\text{m}}^{-3},$

$\mu \left(\alpha \right)=3.569\times {10}^{-5}{\text{K}}^{-1}\text{and}m\left(\alpha \right)=3.351\times {10}^{-5}{\text{K}}^{-1},$

$\mu \left(\Delta T\right)=1359.687\text{K}\text{and}m\left(\Delta T\right)=1367.602\text{K},$

$\mu \left(Z_{Plate}\right)=130786.628\text{m}\text{and}m\left(Z_{Plate}\right)=133353.594\text{m},$

$\mu \left(\kappa \right)=1.079\times {10}^{-6}{\text{m}}^{\text{2}}\cdot {\text{s}}^{-1}\text{and}m\left(\kappa \right)=1.074\times {10}^{-6}{\text{m}}^{\text{2}}\cdot {\text{s}}^{-1},$

$\mu \left(Z_{Ridge}\right)=2839.611\text{m}\text{and}m\left(Z_{Ridge}\right)=2834.526\text{m}.$

Stein & Stein [18] noticed that their estimates of the different parameters with one standard deviation (1σ) are lower than Parsons & Sclater’s [22] estimate. Although the values given just here are closer to those of the latter authors, one can see that the range of possible values is large and the use of standard deviation has no real meaning.

5. Topographic Elevation at Subduction Zones

Data related to subduction zones are shown as blue hatched zones in Figure 2. Data points for trenches correspond to the border of these zones, at the boundary between tectonic plates (compare Figure 1 and Figure 2). Those points defined at trenches (blue in Figure 18) are also related to data points in arcs (red in Figure 18). However, data points affected by present-day ice loading and/or post-glacial rebound (light blue symbols) or main volcanoes and plateaus associated with hot-spot magmatism (green symbols) correspond to the exclusion buffer zones defined in Figure 2 and are discarded in the following analysis.

5.1. Topographic Elevation at Trenches

Trenches are defined in Figure 1 and Figure 2 as the deepest topographic points where one tectonic plate subducts beneath another. Elevation, however, varies considerably (Figure 19(a)). Most of the sea-floor entering subduction is rather young in age (age < 50 Ma; Figure 19(b)) and the amount of sediment at trenches (deepest point) is not very large (Figure 19(c)).

Old ages of the sea-floor (age > 180 Ma) exists in Müller et al.’s [6] dataset and correspond to the sea-floor of the East Mediterranean Sea. Although many arguments support Palaeozoic ages in this area (e.g. [31] [32] [33] [34] [35] ), the precise age is questionable and the isochrones provided by Müller et al. [6] are highly speculative. The sediment thickness in this region is equally subject to discussion. Consequently, data older than 150 Ma at trenches have merely been discarded herein.

After correction from sediment load (as detailed above), and when data are cleared from main disturbing features (use of the exclusion buffer zones other than subduction zone buffer themseleves, as shown in Figure 2), a relationship between sea-floor depth and age at trenches arises (Figure 20). Within two standard errors around the mean (μ ± 2σ), a linear mathematical relationship cannot be formally ruled out but an exponential function (PCM) may be viewed as more coherent. Numerically, the equations relating the depth at trench d_T with age (t) are:

$d_T\left(t\right)=A\times t+B=-22.219\times t-5615.355;R^2=58.43\%$ (9.1)

$d_T\left(t\right)=A+B\cdot \exp \left(-C\cdot t\right)=-8619.514+4079.148\times \exp \left(-0.01916\times t\right);R^2=67.63\%$ (9.2)

On the contrary, any relationship between sea-floor depth and the rate at which the lower tectonic plate subducts beneath the upper one cannot be determined (Figure 21). As found by other authors (e.g. [36] [37] ), there is no apparent relationship between the subduction rate (i.e. the sum of the upper plate motion and lower plate motion defined from Euler poles at every point location at trenches; in mm・yr⁻¹) and the age of the sea-floor entering subduction (in Ma).

5.2. Flexuration of Subducting Lithospheric Plates

Looking at data according to their distance to the trench (maximum distance of 9 corresponding to data within blue buffer in Figure 2, but discarding data belonging to exclusion buffers), the lithospheric plate clearly shows a flexuration (Figure 22). The distance of the bulge―the most elevated point before trench― is located at 1.09 and 1.00 from trench using polynomial smoothing on median and mean values respectively.

The flexuration can be modelled using equations for the bending of an elastic lithosphere (e.g. [14] ; see also [38] ). The deflection of the plate w [m] as function

of the distance (x) [m] from a loading point is written as follows:

$w\left(x\right)=w_b\sqrt{2}\cdot \exp \left(\frac{\text{π}}{4}\right).\exp \left(-\frac{\text{π}}{4}\left(\frac{x-x_0}{x_b-x_0}\right)\right)\cdot \sin \left(\frac{\text{π}}{4}\left(\frac{x-x_0}{x_b-x_0}\right)\right)$ (10.1)

where,

Equation (10.2)― $w_b=-w_0\cdot \exp \left(\frac{-3\text{π}}{4}\right)\cdot \cos \left(\frac{3\text{π}}{4}\right)$ ; w_b is the relative elevation

of the flexural bulge and w₀ is the maximum depression (in m) relative to the plate elevation far from the end load (x → ∞),

Equation (10.3)― $x_0=\frac{\text{π}}{2}\phi$ ; x₀ is the distance at which the bended plate

crosses the elevation of the plate far from the end load (x → ∞),

Equation (10.4)― $x_b=\frac{3\text{π}}{4}\phi$ ; x_b is the distance of the bulge to the end load, i.e.

the location of the highest point due to flexuration,

Equation (10.5)― $\phi =\left(\frac{1}{1000}\right)\times {\left(\frac{4E{h}_e^3}{12\left(1-{\nu }^2\right)}/\left({\rho }_M-{\rho }_W\right)\cdot g\right)}^{\frac{1}{4}}$ ; φ is the flexure

parameter, with:

E: the young modulus in Pascal [Pa].

h_e: the elastic thickness of the plate in metre [m].

ν: the Poisson’s coefficient [Ø].

g: the gravitational acceleration [m・s⁻²].

ρ_M and ρ_W: the density of the mantle and the density of the water column respectively [kg・m⁻³].

The equation expressing the depth due to flexuration as function of distance to an end load is thus of the form:

$\begin{array}{c}z\left(x\right)=A+B\times \exp \left(-C\times x-\left(\text{π}/2\right)\right)\times \sin \left(C\times x-\left(\text{π}/2\right)\right)\\ =A-B\times \exp \left(-C\times x-\left(\text{π}/2\right)\right)\times \cos \left(C\times x\right)\end{array}$ (11)

Seeking the “best” fit, the parameters found are:

$A=-4528.251;B=5078.168;C=2.98073;$

And the resulting curve is shown in green in Figure 23(b). The match between data (mean values μ per bin) and model is quite poor.

In theory, if the end load is applied at trench, the maximal depression w₀ corresponds to the difference between the elevation at trench and the elevation far from the trench, and w₀ can be estimated from Equation (6) and Equation (9.2).

As shown in Figure 19(b), the age of the sea-floor at trenches varies considerably. Assuming however, these ages roughly follow an exponential decay (inset in Figure 23(a)), the mean age defined from natural log distribution (Figure 23(a)) is about 25 Ma. The PCM equation (Equation (6)) provides an estimation of the depth of the sea-floor of d(25) = −4297.689 m. The depth at trench is given by Equation (9.2) and corresponds to d_T(25) = −6092.880 m. The maximum depression would therefore be $w_0=d\left(25\right)-d_T\left(25\right)=1795.190\text{m}$ . The flexural bulge would have a relative elevation of w_b = 120.313 m. Assuming:

$A=d\left(25\right)=-4297.689;B=w_b\cdot {\left(2\right)}^{1/2}.\exp \left(\text{π}/4\right)=373.183;C=2.98073$ as before;

The flexuration model is shown as pink curve in Figure 23(b). However, such calculation does not fit the data at all.

The reasons are: 1) the peak in natural log distribution of age is pk ≈ 3.55 (instead of the value of 3.184 with the mean), which corresponds to a sea-floor age of roughly 35 Ma and sea-floor depth from PCM equation of d(35) = −4598.500 m. This value is much more consistent with the mean value of d = −4599.630 m

determined from a linear fit on data located between 4 and 9 away from trench (deep yellow line in Figure 23(b)); 2) It does not seem correct to regard the trench as the position of the end load.

Assuming the position of the force responsible for the bending of the elastic plate is not located at trench (i.e. end load position ≠ trench position), the best fitting parameters can be found with eq.11 re-written as:

$z\left(x\right)=A+B\times \exp \left(-C\times \left(x-D\right)-\left(\text{π}/2\right)\right)\times \sin \left(C\times \left(x-D\right)-\left(\text{π}/2\right)\right),$ (12)

with:

$A=-4583.442;B=25358.096;C=1.12803;D=\mathrm{0.9503.}$

Using Equation (12), the flexural model (blue in Figure 23(b)) much better fit the data (R² = 80.26%). The end load is located −0.95 or 105.671 km in the back of the trench (i.e. in the direction of the upper plate). The elevation of the bulge is w_b = 353.290 m above the plate elevation far from the end load (x → ∞) found to be A = −4583.442 m. The bulge location is defined at x_b = 1.138 from the trench.

The Bayesian-Markov chain Monte-Carlo inversion method has been used here as well (Figure 24) to try to determine the “best” parameters of Equation (10) (Equations (10.1) to (10.5)) within bounds defined as the two standard errors around the mean (μ ± 2σ) and 95% of the data between the minimum and maximum values per bin around the median (m ± δ_95%).

The parameters were chosen within the following ranges:

ρ_W is defined from the International Equation of State of Seawater [24] as above.

${\rho }_M=3150\pm 300$ i.e. ${\rho }_M\in \left[2850;3450\right]$ (in kg・m⁻³); mean density of the mantle above compensation depth, which value is commonly chosen around

3300 kg・m⁻³ in the literature.

g is chosen constant and g = 9.80665 m・s⁻².

$v=0.25\pm 0.15$ i.e. $v\in \left[0.1;0.4\right]$ (Ø) ; the Poisson’s coefficient value is commonly chosen around 0.25 in the literature.

$h_e=50\pm 50$ i.e. $h_e\in \left[0;100\right]$ (in ×10³ m) ; the elastic thickness of oceanic lithosphere is commonly chosen around 50 × 10³ m in the literature (see in particular [39] [40] ).

$E=10\pm 5$ i.e. $E\in \left[5;15\right]$ (in ×10¹⁰ Pa) ; the Young’s modulus value is commonly chosen around 7 × 10¹⁰ Pa in the literature.

The outcome for the parameters used in Equation (10) does not allow clearly defining favoured values (Figure 25). The following mean μ (parameter) and median m (parameter) values are therefore provided for information only:

$\mu \left({\rho }_M\right)=3152.340\text{kg}\cdot {\text{m}}^{-3}\text{and}m\left({\rho }_M\right)=3150.950\text{kg}\cdot {\text{m}}^{-3},$

$\mu \left(\nu \right)=0.25070\text{and}m\left(\nu \right)=0.24874,$

$\mu \left(h_e\right)=40938\text{m}andm\left(h_e\right)=44255\text{m},$

$\mu \left(E\right)=9.88249\times {10}^{10}\text{Pa}\text{and}m\left(E\right)=9.88255\times {10}^{10}\text{Pa},$

$\mu \left(w_0\right)=-1289.030\text{m}\text{and}m\left(w_0\right)=-2590.540\text{m},$

$\mu \left(EndLd\right)={2.484}^{\circ }\text{and}m\left(EndLd\right)={2.5113}^{\circ }.$

5.3. Flexuration of Subducting Lithospheric Plates as Function of Sea-Floor Age

The set of data that belong to the subduction zones (i.e. within 9 from the trench) have been divided according to the age of the sea-floor (per bins of 10 Ma). For every subdataset, a curve using Equation (12) was fitted (Figure 26(a)) and the variation of coefficients A, B, C & D as function of age is shown in Figure 26(b) (see also Table 5). Polynomial fit of degree 3 (and associated 95% confidence intervals) are merely shown to highlight the first order variation of each coefficient. Concerning the coefficients B, C & D, the variation of the 3^rd order polynomial fits does not highly depart from linear fits with equations:

$B=542.191\times \left(\text{Age}\right)+1055.494$ ; R² = 41.0% (relative to data); R² = 97.4% (relative to 3^rd polynomial fit);

$C=-0.00401\times \left(\text{Age}\right)+1.466$ ; R² = 7.25% (relative to data); R² = 90.0% (relative to 3^rd polynomial fit);

$D=-0.01276\times \left(\text{Age}\right)-0.506$ ; R² = 8.52% (relative to data); R² = 33.8% (relative to 3^rd polynomial fit).

Coefficient A is most probably related to age with a PCM equation (light blue dashed curve in Figure 23(b)) with equation:

$A=-5035.361+2410.361\times \exp \left(-0.04050\times \left(\text{Age}\right)\right)$ ; R² = 85.6% (relative to data); R² = 95.8% (relative to 3^rd polynomial fit).

6. Topographic Elevation of Arcs and Distance Arc-Trench

Lallemand et al. [36] and Heuret & Lallemand [37] , in particular, discussed the bending of subducting lithosphere according to the nature of the upper crust and the stress field. The distance of the arc relative to the trench is clearly related to the angle at which the lower plate plunges in the mantle. Moreover, the elevation of the arc is also clearly dependent upon the nature of the upper plate (continental or oceanic), the stress and strain fields in the upper plate (e.g. extension leads to lower elevation) and the magmatic and erosion activities.

Dealing with topographic elevation only, the results shown herein are therefore a global overview. Further analysis combining at least the aforementioned components would be necessary to decipher the relationship between those processes and better understand the observed topographic result. Nevertheless, the statistical outcomes provided below are not reported in the literature.

6.1. Distance between Trench and Arc

The distance between arc and trench is broadly distributed around a mean value of ca. 215 km (1.937 ± 1.949; μ ± 2σ being aware that the distribution is not Gaussian, Figure 27(a); Table 6). There is no relationship between this distance

and the age of the sea-floor entering subduction (Figure 27(b); the linear fit is near flat with equation y = 0.00159 × x + 1.89714 and mean y value = 1.958 ).

6.2. Topographic Elevation of Arcs

The arc elevation obviously differs as function of the nature of the upper plate. The arc related to intra-oceanic subduction zones are largely under-water and the distribution of elevation is quasi-Gaussian with a mean value close to -1300 m (Figure 28(a); Table 6). The distribution of arc elevation in continental crust is sharper and lesser symmetrical so that the mean and median values differs significantly. In addition, a number of data corresponding to the region of the Altiplano in South America marks a cluster of elevation between 4000 and 5000 m (Figure 28(a) & Figure 28(b)) and shift the mean value.

No relationship can be determined between the arc elevation and the arc ? trench distance for intra-oceanic subduction zones (quasi-flat blue linear regression in Figure 28(b)). For active margins however, arcs closer to their trenches seem to predominantly exhibit lower elevation than arcs located farther away. The reason seems to be straightforward: When the arc is close to the trench, the subducting slab is dipping with high angle usually in agreement with slab roll- back processes. The upper plate undergoes extension leading to the lowering of the elevation of the arc. On the contrary, under compression, the upper plate is squeezed and the arc uplifted.

6.3. Variation of Topography along Arc

The topographic elevation along arcs does not appear to be randomly distributed.

Along the Marianas, for example, the topographic elevation varies as in Figure 29. At first glance, one may see an irregular longwave variation, in the order of 1000 km (outlined with polynomial fit in green in Figure 29(a)), and a shorter, with variation of the order of 200 km (outlined with a sinusoidal function in purple in Figure 29(a)).

Using periodograms (created with Past [41] ), the power versus frequency diagram (Figure 30) shows several significant peaks. In particular, the multiple-taper spectral analysis [42] highlights peaks at −2.3336, −1.8720, and −1.6581 on the log10 of the frequency.

While the longwave variations might be caused by numerous factors, one can speculate that the short-wave variations are linked to volcanoes, and therefore to a rather periodic distribution of magma chambers beneath arcs. Adjusting Gaussian fit, the first peak of the periodogram corresponds to a phase of 215.278 km (±4.5 km at the 95% confidence level).

Using all data for arcs together (including data from cordillera and from island arcs), Fourier transforms, periodograms and other techniques such as

wavelets analysis fail to provide clear picture of periodic signal. However, it might be useful to carry out further analysis with refined datasets in order to decipher typical signal of magma-crust interaction versus modified signal related to other processes.

7. Topographic Elevation at Passive Margins

Passive margins are no plate limit but encompass the transition between continental crust and oceanic crust. The definition of continent-ocean boundaries (COBs) is however not straightforward for three main reasons: 1) COBs are buried under a large amount of sediment and/or under volcanic material, which make the COB difficult to identify accurately even with powerful geophysical tools; 2) fragments of tilted block of continental nature (or “extensional allochthons”) may be left apart and separated from the main continent (e.g. [43] [44] [45] [46] ); 3) the tearing of the continental crust may not necessarily lead to the creation of a proper oceanic crust but rather to denuded and often altered mantle; the nature of those rocks lead some authors to use the term “transitional crust”.

Because the thickness of the crust is highly varying and therefore subject to high uncertainty, the thickness of sediment is equally subject to caution, the presence, thickness and extension of magmatic rocks is difficult to assess, and the nature of the underlying mantle itself may be problematic, no correction for sediment load was attempted for data from passive margins.

Using all topographic data within a distance of ±9 away from COBs (as defined in section.2; Figure 1), a clear dichotomy is observed between continental and oceanic crust (Figure 31). The running average (running window of 1001 data; red curve in Figure 31(a)) shows that the change is relatively abrupt and occurs mainly between −3 and +3 from COB in average.

If the general shape resembles (R² = 99.45%) a hyperbolic tangent function of

equation f(x) = 2236.076 × tanh(−0.77391 × x) − 1936.493 (purple curve in Figure 31(a)), the general slope of the passive margins is even steeper at COB than the latter equation. The slope is maximal between −0.2 and +0.2 from COB and reaches nearly 4% in average (s = 3.974% which is equivalent to α = 2.275 ).

Data have then been divided by bins of 1 relative to the distance to COBs and plotted against the age of the COB (age at which the two continental domains are separated whatever the nature of the rocks at sea-floor; Figure 31(b) & Figure 31(c)). Every sub-dataset is then fitted with a 6^th order polynomial curve (Figure 31(c)). Although the fitting technique creates artefacts, positive curvatures visible at distance of few degrees from COB and for ages younger than ca. 20 Ma are attributed to rift and passive margin shoulder formation.

8. Topographic Elevation of Continents

Although largely neglected, the topographic elevation of continents has major implications for many processes of the Earth evolution, including eustatism, climate, erosion-sedimentation-sediment fluxes, etc. Flament [47] , for instance, studied the long term evolution of the Mean Continental Altitude (MCA) and showed the importance of MCA upon the Earth cooling history, magmatic evolution and the growth of the continental crust, and the implications for eustatism or sea-water chemistry (see also [48] for instance).

The global statistics for topographic elevation of crust of continental nature are provided in Figure 1 and Table 1. Now, the data stem from the ETopo1 model [5] which provides the topographic elevation of the bedrock. A large number of data from Antarctica and Greenland are affected by ice loading. Those region and beyond are also impacted by the post-glacial rebound associated with the last glacial maximum [15] [16] . The topographic data must therefore be corrected from those influences. Although an order of magnitude lower than the post-glacial rebound, the data were additionally corrected from the geoid using the Earth gravity field model Eigen-Grace02S [49] . To first order, however, those corrections do not fundamentally impact the global distribution of topographic elevation (Figure 32(a)). Now, the residual long-wave variations might be useful to infer information upon dynamic topography.

8.1. Relationship between Continental Elevation and Moho Depth

The thickness of the crust has been investigated by Mooney et al. [50] with the

Crust.5.1 model (5 × 5 model) now updated to the Crust.1.0 model ( [8] ; 1 × 1 model available online: http://igppweb.ucsd.edu/~gabi/crust1.html). The distribution of the Moho depth for data assumed to lie on continental crust mainly exhibits two peaks (Figure 32(b)). The first is located at a depth below 40 km and the second below 30 km (−37.00 km and −29.36 km after multiple terms Gaussian fit, respectively). The linear relationship between the Moho depth d [km] and the topographic elevation e [m] is given by the equation:

$d=-0.006206095756\times e-32.49861243;R^2=56.27\%$ (13)

(red line in Figure 32(c)).

Note that, in fact, eq.13 does not significantly depart from the linear fit obtained from data corrected from ice loading-post-glacial rebound-geoid since the equation is:

$d=-0.006218679627\times e-32.14407886;R^2=56.78\%.$ (14)

From the histogram of the Moho depth (Figure 32(b)), two main peaks stand out: the most prominent at a depth of 36.998 km and the second, at 29.359 km. From the linear relationship (Equation (13)), and even from the polynomial fits of higher degrees, the corresponding two main topographic elevations should be +725.047 m and −505.971 m (+741.088 m & −491.721 m with a polynomial fit of degree 2; +642.537 m and −422.962 m with degree 6), respectively. Those results do not compare with the histogram of the topographic elevation (Figure 32(a)).

It can be inferred from this that the Moho depth defined by Laske et al. [8] is too crude to be compared with the ETopo1 model [5] , because a large number of data for the Moho does not properly correspond to the histogram of the topographic elevation. However, the general trend between topography and Moho depth is clear and Equation (13) can be regarded as a good first approximation.

8.2. Implications for the Mean Density of Continental Crust

According to Equation (13), the continental crust thickness is null at a depth of −4509.878 m. Assuming a mantle density ρ_M = 3150 ± 300 [kg・m⁻³] as above and using this linear relationship (Equation (13)), the mean crustal density ρ_C shall correspond to the following equation:

${\rho }_C=\left({\rho }_M\times d^{\prime }\right)/h_C,$ (15)

where d' is the Moho depth below −4509.878 m, and h_C is the thickness of the crust (i.e. (e + d)).

As a consequence, the mean density of the continental crust is ρ_C = 2712.870 ± 258.369 [kg・m⁻³].

However, the right relationship between topographic elevation and Moho depth is undetermined. When polynomial fits of degree 2 (for which R² = 56.35%) or degree 6 (for which R² = 58.86%) are used for instance, the mean densities of continental crust does not directly correspond to a single value (as per eq.15), but can be inferred from the data distribution. Hence, the peak values are ρ_C = 2678 ± 254.0 [kg・m⁻³] and ρ_C = 2691 ± 256.3 [kg・m⁻³] for polynomial fits of degree 2 and degree 6, respectively.

8.3. A “Normal” Continental Crust

The literature often refers to a “normal” or “typical” continental crust, id est a crust that has not been affected by thinning or thickening processes. It is therefore commonly accepted to use the value of ca. 250 m for the topographic elevation and 35 - 40 km for the Moho depth.

Equation.13 places the Moho at a depth of 34.050 km if the elevation is 250 m (33.699 km using Equation (14)), and conversely, the topography elevation shall rise to an elevation comprised between +403.053 m and +1208.712 m if the “typical” Moho depth is considered. Note that the Moho is even set at a depth of 33.917 km and 34.376 km if the polynomial fits of degree 2 and degree 6 are respectively considered. As previously noticed, those values are not quite in agreement with those generally accepted.

Focussing herein on topography, the definition of a “normal” continental crust will be tempted from the present-day global topography, corrected from ice, post-glacial rebound, and geoid.

From the distribution of the global topography (Figure 33(a); same as per Figure 32(a)), the mean altitude of continent is μ_C = 273.607 m. Although the distribution is clearly not Gaussian, the peak is well-marked and the distribution is broadly symmetrical, so that the median value is relatively close to the mean (m = 235.208 m; Δ(μ − m) = 38.399 m; Table 7). This observation probably explains why the standard value of +250 m is commonly adopted.

Now, the shape of the data distribution (Figure 33(a)) clearly indicates that continental crust with topographic elevation lower than 0 m to −1000 m are altered by thinning effects. In parallel, a bulge in data distribution around ca. +1000 m suggests the effects of crustal thickening. The range of values for “normal crust” certainly lies between those bounds but cannot be precisely determined from this dataset only. Consequently, the bounds for “normal crust” (i.e. crust considered to be not significantly affected by thinning or thickening processes) have been arbitrarily chosen to be −333 m and +777 m. The corresponding distribution within those limits is shown in Figure 33(b), and the spatial distribution is shown in Figure 33(c) (blue dots below -333 m, yellow and orange dots comprised between −333 m and +777 m, and red dots above +777 m).

Using the limited subdataset (i.e. [−333; +777]), Figure 33(b)―which actually corresponds to a zoom on the peak of Figure 33(a)―discloses the irregular shape of the data distribution. A single term Gaussian fit (dark blue curve in Figure 33(b)) proves the mismatch; the mean value for “normal crust” is μ_nC = +185.6 m (± 408.1 m; 1σ) and is not statistically representative. Although a Gaussian fit with multiple terms obviously better fit the distribution (example with a six term fit in light blue in Figure 33(b)), the signification of the location and shape parameters of each term is obscure. At first, nevertheless, the data distribution displays two main peaks. A Gaussian fit with two terms (blue curve in Figure 33(b)) can bring out those peaks for which the location and shape parameters are: μ_nC.1 = +24.35 m (± 97.02 m; 1 σ₁) and μ_nC.2 = +310.2 m (±160.3 m;

1 σ₂). The boundary between the two peaks can be located at ca. +270 m. It seems to separate what can be termed the “lowlands” (with μ_nC.1 = +24.35 m) from the “highlands” (with μ_nC.2 = +310.2 m). And the latter largely corresponds

to cratonic areas (Figure 33(c)).

Interestingly however, those peak values (μ_nC.1 = 24.35 m & μ_nC.2 = 310.2 m) are quite far from the standard of ca. +250 m. And, according to the definition taken herein, “normal” or “typical” continental crust represents ca. 65.66% (Table 7) of continental crust in the world at present-day.

Besides, albeit the misgivings set out above (Section 8.1), the corresponding Moho depth shall be―to first order (Equation (13))―of 32.650 km under lowlands, and 34.424 km under highlands.

8.4. Map of Dynamic Topography

Kaban et al. [51] described the dynamic topography as “the long-wavelength part of non-isostatic topography [which] is supposed to be generated by mantle dynamics and is up to now not well studied”. Those authors produced a global map of dynamic topography from gravity field data.

It can also be considered that if continental crust was of same thickness everywhere, the long-wave length variation of topography, once corrected from ice loading-post-glacial rebound-geoid, would be due to dynamic topography.

Considering the trend between topographic elevation and Moho depth as manifest, an attempt is made here to correct the effect of that trend in order to highlight the dynamic topography component. It is not determined, however, how the trend should be underlined. Polynomial fits of various degrees (from degree 1 (i.e. linear) to degree 10) are presented in Figure 34(a), and Figure 34(b) and −34.c show the correction that the two end-members (degree 1 and degree 10) imply. In other words, Figure 34(b) and Figure 34(c) show how the topographic elevation would reduce if the continental crust was in general constant. Spatially, the topographic elevation is shown in Figure 35(a) and Figure 35(c) after reduction from polynomial fit of degree 1 (linear) and degree 10 respectively. The use of a Gaussian filter (Figure 35(b) & Figure 35(d)) consequently highlights area where the dynamic topography presumably acts by uplift (red) or subsidence (blue). The amplitude of the obtained dynamic topography is relative because highly dependent on the filter applied. Now, most of the amplitude topography is comprised with ± 1000 m, although stronger negative values are obtained close to subduction zones.

9. Conclusion

The statistics on the Earth’s topography presented herein shows that most of the values commonly chosen in the literature as “typical” are subject to caution. The so-called “typical” MCA of +250 m, “typical” Moho depth of ca. 35 km - 40 km,

or “typical” depth of mid-oceanic ridges of −2600 m are values that largely mismatch the values obtained here.

In general, the use of mean values (μ) is inappropriate. Depending on the precision required, the median value (m) might be more relevant because it is more robust even if the data distribution is relatively symmetrical. In addition, the use of standard deviation (μ) is statistically incorrect because a large part of the data considered herein is not Gaussian distributed.

As the present-day topography is the sole example of topography that we have, more caution should be taken regarding the use of statistical values of the topography, in particular for palaeotopography issues. Indeed, many physical processes are at work behind those “mean values”, and it is important to further study them in order to better understand what parameters are predominant in these processes.

Now, much more should be done using the statistics of the Earth’s topography and also using other datasets (heat flux, magnetics, gravimetry, etc.) in order to decipher the role of the various processes. I hope this paper will pave the way for further studies in this direction.

Cite this paper

Vérard, C. (2017) Statistics of the Earth’s Topography. Open Access Library Journal, 4: e3398. https://doi.org/10.4236/oalib.1103398

References

Annexe

Annexe 1. Statistics of elevation data (ETopo1) for every 0.025˚ bin of distance to ridge axis