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.     ). 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  ; 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  .
・ Global sediment thickness map  (updated online versions, 2010; hereafter referred as to  ); the map grid has 2 × 2 degrees horizontal resolution, and the vertical resolution is not provided. The sediment thickness of the NOAA (  ; 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  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 RMean:
where, Req. is the equatorial radius and Rpol. the polar radius of the ellipsoid of reference (WGS84). In the following, the Earth’s surface is therefore 5.101 × 108 km2 and the volume, 1.083 × 1012 km3; 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 × 106 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,  ) and the second derivative in the magnetic anomaly map EMag2  . The result is however not fundamentally different from other publications as reported by Eagles et al.  , and does not significantly affect the statistical outcomes provided below.
As first put forward by Otto Krümmel  , 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 moceanic = −4325 m and mcontinental = +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, peakoceanic = 6.268 × 104 data
Figure 1. (a) Geodetic data grid with points plotting on assumed continental crust (brown dots) and on assumed oceanic crust (blue dots). Purple lines (and purple names) correspond to tectonic plate boundaries; (b) Histogram showing the distribution of topographic elevation over the globe. Using the full dataset in purple (background), with data assumed to lie on oceanic crust in blue, and data assumed to lie on continental crust in red. “Best fits” using normal distribution (Gaussian) on oceanic dataset (blue curve) and continental dataset (red curve) shown for comparison; (c) Hypsometric curve highlighting that 72.2% of the world is under sea-level versus 27.8% above.
in the bin [−4300 m; −4400 m] and peakcontinental = 1.108 × 105 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-
Table 1. Global topography statistics.
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.  ), 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.   ).
Figure 2. World topography (after ETopo1;  ) and main features discussed from a statistical point-of-view in the text.
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  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
Figure 3. (a) Sketch depicting how ridge segment are defined as spreading centres and transform faults. The motion of a tectonic plate A is defined relative to a tectonic plate B thanks to an Euler pole and angle. The orientation of a ridge segment can thus be determined relative to the motion vector of the plates at the location of the segment. Segments oriented within the quandrant parallel of the motion vector (i.e. ±45˚ of the motion vector orientation) are defined as transforms, the other defined as spreading centres; (b) Distribution of the orientation of all ridge segments (in red in Figure 2). Spreading centres are centred around 90˚, and transform faults are centred around 0˚ and 180˚.
Figure 4. Zoom on a portion of the Central Atlantic Ocean, off the coast of the French Guyana. On the topography map (after ETopo1;  ) are shown in particular the definition of the mid-oceanic ridge with spreading centres (yellow dots) and transform faults (purple dots), defined as in Figure 3; Thick grey dashed line corresponds to the mid-oceanic ridge, and thin grey dashed lines are isochrons after Müller et al. (2008). Inset: example of topographic profile showing the inner rift of the ridge (located at about 1.5˚), the ridge crests (at about 1.2˚ - 1.4˚ and 1.7˚), and the ridge flanks (below 1.2˚ and above 1.7˚). Red lines correspond to linear fitting on the topographic elevation on each flank (see Section 3.2).
 ). 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.    ).
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,
Table 2. Statistics of mid-oceanic ridges (inner rift).
Figure 5. Statistics on mid-oceanic ridges inner rifts; (a) Topographic elevation at inner rift with blue: all data (with blue Gaussian best-fit), yellow: segments corresponding to spreading centres only (with orange Gaussian best-fit), red: spreading centres excluding data covered by the buffer zones defined in Figure 2 (with dark red Gaussian best-fit); (b) Age of the sea-floor at the location of inner rift after Müller et al.’s dataset  . Same colour-coding as in a). Inset: representation using log-scale and Gaussian best-fit, only for data excluding transform segments; (c) Age error of the sea-floor at the location of inner rift after Müller et al.’s dataset  . Same colour-coding as in (a). Inset: representation using log-scale and Gaussian best-fit, only for data excluding transform segments; (d) Accreted surface of sea-floor at the location of inner rift (in mm2/yr) computed in this study. Same colour-coding as in (a). Inset: representation using log-scale and Gaussian best-fit, only for data excluding transform segments; (e) Spreading rates (in mm/yr) computed in this study in green, and from Müller et al.’s dataset  in purple; (f) Spreading rates (in mm/yr) obtained from Müller et al.’s dataset  compared with those computed in this study. Red line (and yellow 95% confidence interval) represents the best linear fit through origin.
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.  was coarser than the definition made here (see Figure 4), the age of the sea-floor determined from Müller et al.’s dataset  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.  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  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 mm2/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.  (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. (  ; 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 mm2/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
Figure 6. (a) Topographic elevation (bathymetry in meter) of mid-oceanic ridges as function of sea-floor accretion rate (in mm2/yr) computed in this study with blue: all data (with blue best linear fit), yellow: segments corresponding to spreading centres only (with orange best linear fit), red: spreading centres excluding data covered by the buffer zones defined in Figure 2 (with dark red best linear fit); (b) Topographic elevation (in m) of spreading centres inner rifts (excluding data from buffer zones as in Figure 2) as function of the spreading rates (in mm/yr) computed by Müller et al. (  ; data. A in purple), and as the spreading rates (in mm/yr) computed by this study (data. B in green); (c) Equations corresponding to the best fit functions shown in (b), i.e. linear fits and power fits with corresponding colour coding.
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)).
Figure 7. (a) World topography (ETopo1;  ) showing selected data (pink dots) belonging to area accreted since chron C3 (ca. 6 Ma) and the mid-oceanic ridge (yellow); (b) Zoom on the Central Atlantic area (similar to Figure 4; Note that the black line correspond to the profile shown in Figure 4) showing the density of data points (pink); (c) Topographic elevation (in meter) of all points (pink in a) as function of the distance (in degrees) to their respective spreading centres. The red line corresponds to the running average (with 101 data running window).
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
Figure 8. Topographic elevation (ETopo1) around (within 1.333˚) the ridge axis, using mean values (black dots and associated 95% confidence intervals in dark red) and standard deviation (1σ as grey bars) per bin of 0.025˚. Green curves represent the best fit for the inner rift section (2nd order polynomial fit from 0˚ to 0.15˚ in distance), the nearest section of the flank (linear fit from 0.15˚ to 0.75˚) and the farthest section (linear fit from 0.75˚ to 1.333˚). Yellow lines represent those linear fit extended to ridge axis (dashed line for the near flank and dash-dotted line for the far flank). The red line corresponds to the linear fit of data for the entire flank (from 0.15˚ to 1.333˚) extended to ridge axis; Inset: Graph showing the variation of number of data per bin (grey), and variation of the 95% confidence interval (CI, dark red), which reflects the precision and dispersion of data around the means.
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 mm2/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 mm2/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.     ). 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  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-
Figure 9. Topographic elevation (ETopo1 in m) around the ridge axis (within 1.333˚) as function of distance to ridge axis (in degrees) and accretion rates (accreted surface since chron C3 in mm2/yr); (a) Running averages (50 data sliding window) using sub-datasets divided according to 0.5 mm2/yr bins of accretion rate (see colour coding); (b) Normalized linear fits per sub-datasets (same colour coding as in (a)) for the three sections: inner rift (0.00˚ - 0.15˚ from ridge axis), near flank (0.15˚ - 0.80˚) and far flank (0.80˚ - 1.333˚); (c) 3D diagram illustrating the relationship between topographic elevation (depth in m), accretion surface (rate in mm2/yr), and distance to ridge axis (in ˚); (d) Same illustration as in (c) but in 2D, with emphasis on the relationship between depth and rate, and colour coding for landmark distances to ridge axis.
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:
, when t < 20 Ma (i.e. young-aged sea-floor) (2.1)
Table 3. Determination of the elevation of “virtual ridge axis” (in metre) as function of accretion rates (in mm2/yr).
, when t ≥ 20 Ma (i.e.
old-aged sea-floor) (2.2)
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+13 s per Ma),
dR 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−1),
κ is the thermal diffusivity (in mm2・s−1),
ρM and ρW are the mantle and water densities respectively (or mass per unit volume in kg・m−3),
ZL is the asymptotic thermal plate thickness (i.e. at infinite distance from the ridge, in m),
TM and T0 are respectively the basal and the top temperature of the thermal plate (in K).
It must be noticed that both equations are dependent upon dR, 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.  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:
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  . 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)
Figure 10. (a) Sea water temperature (in ˚C) profile as function of water depth (in m). In deep yellow are ranges of temperature profiles from polar to equatorial zones. In blue, an example of typical temperature profile from the Atlantic Ocean near the Bahamas (after http://ocp.ldeo.columbia.edu/climatekidscorner/whale_dir.shtml). In red, the global mean temperature profile obtained from bootstrapping method within deep yellow bounds; (b) Sea water salinity (in ‰) profile as function of water depth (in m). Same colour coding as in (a)), i.e. example of salinity profile from the Atlantic Ocean in blue, range of profiles from polar to equatorial zones in deep yellow, and mean values in red; (c) Sea water density (in kg・m−3) profile as function of water depth (in m). In pink are ranges of density profiles and in purple, the mean values calculated from the International Equation of State of Seawater  . In deep yellow are ranges and in red the mean values of the sea water density (in kg・m−3) integrated over the water column and used for isostatic correction.
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.  have proposed that the mean sediment density ρS over the sediment pile can be calculated using:
after Athy  (4.2)
hS 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,
φ0 is the initial porosity,
z is the compaction decay wavelength,
z is the depth within the sediment pile.
Assuming ρigW = 1013 kg・m−3, ρg = 2650 kg・m−3 and ρM = 3300 kg・m−3, Winterbourne et al.  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 = 0.56 and compaction decay length λ = 4.5 km.
Audet & Fowler  in particular, have shown the complexity of the relationship between porosity and sediment thickness. The relationship put forward by Athy  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−3 (for T = 4 C and S = 34.7‰; Patm = 101,325 Pa), which is higher than the 1013 kg・m−3 used by Winterbourne et al.  . 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.  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:
In Winterbourne et al.  , the density of solid grain is taken to be ρg = 2650 ± 250 [kg・m−3]. It is found here (see below, inset in Figure 12) that a value of ρg = 2300 ± 500 [kg・m−3] (“best” value of ρg = 2292 [kg・m−3]) 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−3].
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 hW = −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  , 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.  . 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.
Figure 11. (a) Mean sediment density in red with associated uncertainty bounds in deep yellow as function of height of the sediment pile (or sediment thickness in m) for a fixed water depth of hW = −4325 m; (b) Correction factor C used for the Airy-Pratt isostatic compensation for sediment load as function of height of the sediment pile with fixed water depth (hW = −4325 m). In red, mean values, and in deep yellow, associated uncertainty bounds.
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, φ0, λ) 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−3 (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
Figure 12. Examples of data distribution of the sea-floor elevation (in m) after isostatic correction for sediment load for data of same age within 1 Ma bin; Examples of histogram are shown for sub-datasets every 9 Ma together with their associated “best” Gaussian fit; Because the data distributions are often closer to Laplace distribution, Gaussian curves have lower amplitude and larger width; Inset: Variation of the cumulative difference between μ− and μ+ according to the value chosen for the density of the grain in sediment ρg; “Best” value is obtain for ρg = 2292 kg・m−3.
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:
Figure 13. Sea-floor elevation (in m; raw data from ETopo1 in dark blue, and after isostatic correction (μ ± 2σ in red and grey error bars, and μ− & μ+ in deep yellow) as function of sea-floor age after Müller et al. (2008) averaged over the globe. The HSCM fit (pink) does not well reproduce the mean values (μ) unless the fit relates only to data younger than 70 Ma (purple). The PCM fit (green) better match all data.
with coefficient of determination R2 = 88.3%.
As emphasized by previous authors, the fit is much better when applied on data for the first 70 Ma (e.g.  , 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:
with coefficient of determination R2 = 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:
with coefficient of determination R2 = 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.  . PCM equations of the form 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−1, 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 3rd 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:
in the equation provided by Stein & Stein (  ; see Equation (2)), where dR corresponds to the “virtual” ridge axis (see analysis of the mid-oceanic ridges above;
Figure 14. Variation of the coefficients of PCM equations when fitted to topographic data as function of spreading rates. The PCM equations are of the form ; values for coefficient A are shown in blue (with dashed blue 3rd polynomial fit and associated 95% confidence interval (CI)) and refer to the outer left vertical axis; values for coefficient B are in yellow (with dashed yellow fit and 95% CI) and refer to the inner left vertical axis; values for coefficient C are in red (with dashed red fit and 95% CI) and refer to the inner right vertical axis. Number of data per bin of rate are shown in grey and refer to the short inner vertical axis to the right; the coefficient of determination R2 associated to the fit of PCM equations to every sub-datasets are shown in black and refer to the short outer vertical axis to the right.
Table 4. Coefficient A, B and C as function of spreading rate [mm・yr−1] as provided by Müller et al. (2008).
Section 3.2) and is given as dR = 2600 m.
The definition of dR from the coefficients determined in Figure 14 and Table 4 is simply . For the other parameters however, the determination is not trivial. Stein & Stein  , in particular, noticed the difficulty to quantify uncertainties.
To do so, the genuine equation for Plate Cooling Model (PCM;    ) was used and is written as follow:
After evaluation of the integral  , the equation can be re-written as:
The Bayesian-Markov chain Monte-Carlo inversion method employed by Scholer   in another context was used here to quantify the unknown parameters, namely ρM, α, TM (T0 is assumed to be the temperature defined at the sea-floor as in Figure 10), κ, and ZL (ρ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:
(in kg・m−3) ; mean density of the mantle above compensation depth, which value is commonly chosen around 3300 kg・m−3 in the literature.
(in K−1) ; volume coefficient of thermal expansion, which value is commonly chosen around 3 × 10−5 K−1 in the literature.
(in K); Temperature at compensation depth, which value is commonly chosen around 1625 K in the literature, so that .
(in m); Thermal plate thickness, which value is commonly chosen around 125 km in the literature.
(in mm2・s−1 or × 10−6 m2・s−1) ; Thermal diffusivity, which value is commonly chosen around 1 mm2.s-1 in the literature.
Note that in Equation (8.2), the sum is carried out up to m = 10,000.
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%) (Niter = 184,603 and Niter =
Figure 15. Stochastic inversion concerning the age-depth dependence of the sea-floor; (a) Mean values (red) of the elevation of the sea-floor corrected from sediment load as function of the age (divided in bins of 1 Ma) and associated two standard error (grey bars): μ ± 2σ (same as in fig.13); the errors bounds (2σ) are smoothed by polynomial fits (deep yellow). Median values (dark red) of the elevation of the sea-floor corrected from sediment load as function of the age (divided in bins of 1 Ma) and associated range of the 95% of data the closest to the median (orange): m ± δ95%); the errors bounds (δ95%) are smoothed by polynomial fits (green); (b) 5000 accepted curves (light blue) corresponding to the stochastic inversion within μ ± 2σ (smoothed); (c) 5000 accepted curves (light blue) corresponding to the stochastic inversion within m ± δ95% (smoothed).
Figure 16. Distribution of the possible values for unknown parameters in Equre (8.2); Dark red bins and associated red “best” Gaussian fit for data within 2 standard error around the mean (μ ± 2σ); Grey bins and associated blue “best” Gaussian fit for data within 95% of the number of data the closest to the median (m ± δ95%); (a) Mean mantle density above compensation depth ρM [kg・m−3]; (b) Volume coefficient of thermal expansion α [K−1]; c) Temperature difference [K] ΔT = (TM - T0) between the base and the top of the cooling plate; (d) Thickness of the thermal plate [m]; (e) Thermal diffusivity κ [m2・s−1]; (f) Virtual mid-oceanic ridge elevation [m].
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:
, i.e. (in kg・m−3); Mean density of the mantle above compensation depth,
, i.e. (in K−1); Volume coefficient of thermal expansion,
, i.e. (in K); Difference in temperature between the base and the top of the thermal plate,
, i.e. (in m); Thermal plate thickness,
, i.e. (in mm2・s−1); 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:
Stein & Stein  noticed that their estimates of the different parameters with one standard deviation (1σ) are lower than Parsons & Sclater’s  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.
Figure 17. Stochastic inversion concerning the age-depth dependence of the sea-floor with range limited to μ− and μ+, and associated distribution of the possible values for unknown parameters in Equation (8.2); (a) Mean values (red) of the elevation of the sea-floor corrected from sediment load as function of the age (divided in bins of 1 Ma) (same as μ in Figure 13); the bounds, corresponding to the mean values (per 1 Ma bin) applied upon topographic elevation data accounting for uncertainties (deep yellow; same as μ− & μ+ in Figure 13), are smoothed by polynomial fits (green); (b) Mean mantle density above compensation depth ρM [kg・m−3]; (c) Volume coefficient of thermal expansion α [K-1]; (d) Temperature difference [K] ΔT = (TM - T0) between the base and the top of the cooling plate; (e) Thickness of the thermal plate [m]; (f) Thermal diffusivity κ [m2・s−1]; (g) Virtual mid-oceanic ridge elevation [m].
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  dataset and correspond to the sea-floor of the East Mediterranean Sea. Although many arguments support Palaeozoic ages in this area (e.g.      ), the precise age is questionable and the isochrones provided by Müller et al.  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.
Figure 18. Location of data points defined at trenches and arcs; blue dots are data at trenches; light blue and green dots are data located in exclusion buffer zone; red triangles are data defining the associated magmatic arc; light blue and dark green triangles correspond to data located in exclusion buffer zone.
Figure 19. Histogram of data at trenches; (a) Topographic elevation [m] after ETopo1  ; (b) Age of the sea-floor entering subduction [Ma] after Müller et al.  ; (c) Sediment thickness at trench (deepest point) after Laske & Masters  (dark red) and after the NOAA  (green). For comparison, an exponential fit and log-normal fit have been adjusted on data from Laske & Masters  .
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 dT with age (t) are:
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.   ), 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−1) 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.  ; see also  ). The deflection of the plate w [m] as function
Figure 20. Topographic elevation right at trenches as function of sea-floor age  entering subduction; (a) Raw bathymetric data (after ETopo1;  ) inclusive data from the East Mediterranean Sea with ages older than 150 Ma; Linear (purple) and PCM (blue) fits are shown for comparison; (b) The topographic data younger than 150 Ma have been corrected from sediment load and averaged per 1 Ma bin; mean (red curve; μ) and associated double standard error (grey bars; 2σ); median (dark red; m); Linear (pink) and PCM (light blue) curves fit the median values.
Figure 21. Relationship with subduction rate (in mm/yr); (a) Histogram of subduction rates defined at trenches; (b) Subduction rate [mm・yr−1] as function of the sea-floor age (excluding data older than 150 Ma); (c) Topographic elevation [m] after sediment load correction as function of subduction rate at trenches [mm・yr−1]; black crosses, data younger than 150 Ma; grey crosses, data older than 150 Ma.
of the distance (x) [m] from a loading point is written as follows:
Equation (10.2)― ; wb is the relative elevation
Figure 22. Topographic elevation [m] of data as function of distance [˚] to trench; (a) Data points (black crosses) and running average (deep yellow, with sliding window of 101 data); (b) mean value (μ; red curve) per bin of 0.1˚ in distance to trench and associated two standard error bars (2σ; grey), and median value (m; dark red curve).
of the flexural bulge and w0 is the maximum depression (in m) relative to the plate elevation far from the end load (x → ∞),
Equation (10.3)― ; x0 is the distance at which the bended plate
crosses the elevation of the plate far from the end load (x → ∞),
Equation (10.4)― ; xb is the distance of the bulge to the end load, i.e.
the location of the highest point due to flexuration,
Equation (10.5)― ; φ is the flexure
E: the young modulus in Pascal [Pa].
he: the elastic thickness of the plate in metre [m].
ν: the Poisson’s coefficient [Ø].
g: the gravitational acceleration [m・s−2].
ρM and ρW: the density of the mantle and the density of the water column respectively [kg・m−3].
The equation expressing the depth due to flexuration as function of distance to an end load is thus of the form:
Seeking the “best” fit, the parameters found are:
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 w0 corresponds to the difference between the elevation at trench and the elevation far from the trench, and w0 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 dT(25) = −6092.880 m. The maximum depression would therefore be . The flexural bulge would have a relative elevation of wb = 120.313 m. Assuming:
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
Figure 23. (a) Natural log distribution of sea-floor ages at trench, and “best” Gaussian fit (blue curve); Inset shows the distribution with raw data with “best” exponential fit (blue); (b) mean value (μ; black curve) per bin of 0.1˚ in distance to trench (as per Figure 22(b)); green curve, flexural model when “best” fitting parameters A, B & C; pink curve, flexural model when assuming that the end load is located at trench and the age of the sea-floor at trench is 25 Ma after mean value of the Gaussian fit in a); deep yellow line, “best” linear fit on data (μ) located at distance between 4˚ and 9˚ away from the trench (at 0˚) showing that the plate elevation far from the end load (x → ∞) is about −4600 m; blue curve, flexural model when “best” fitting parameters A, B, C & D, where D is the distance at which the end load is located from the trench.
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:
Using Equation (12), the flexural model (blue in Figure 23(b)) much better fit the data (R2 = 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 wb = 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 xb = 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  as above.
i.e. (in kg・m−3); mean density of the mantle above compensation depth, which value is commonly chosen around
Figure 24. Stochastic inversion concerning the distance-depth dependence of the sea-floor at subduction zones; (a) Mean values (μ; red) of the elevation of the sea-floor corrected from sediment load as function of the distance (divided in bins of 0.1˚); the errors bounds (±2σ, deep yellow) are smoothed by polynomial fits (orange). Median values (m; dark red); the errors bounds (±δ95%; green) are smoothed by polynomial fits (dark green); (b) 5000 accepted curves (light blue) corresponding to the stochastic inversion within μ ± 2δ (smoothed); (c) 5000 accepted curves (light blue) corresponding to the stochastic inversion within m ± δ95% (smoothed).
3300 kg・m−3 in the literature.
g is chosen constant and g = 9.80665 m・s−2.
i.e. (Ø) ; the Poisson’s coefficient value is commonly chosen around 0.25 in the literature.
i.e. (in ×103 m) ; the elastic thickness of oceanic lithosphere is commonly chosen around 50 × 103 m in the literature (see in particular   ).
i.e. (in ×1010 Pa) ; the Young’s modulus value is commonly chosen around 7 × 1010 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:
Figure 25. Distribution of the possible values for parameters in Equation (10); Dark red bins and associated red “best” Gaussian fit for data within 2 standard error around the mean (μ ± 2δ); Grey bins and associated blue “best” Gaussian fit for data within 95% of the number of data the closest to the median (m ± δ95%); a) Young’s modulus E [Pa]; (b) Elastic thickness he of the plate [m]; (c) Poisson’s coefficient ν [Ø]; (d) Mean mantle density above compensation depth ρM [kg・m−3]; (e) Maximal depression w0 [m] below the depth of the plate far from the loading point (x → ∞); (f) Position of the loading point [˚] onto the bending plate relative to the trench (positive values backward, i.e. in the direction of the upper plate).
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 3rd order polynomial fits does not highly depart from linear fits with equations:
; R2 = 41.0% (relative to data); R2 = 97.4% (relative to 3rd polynomial fit);
; R2 = 7.25% (relative to data); R2 = 90.0% (relative to 3rd polynomial fit);
Figure 26. (a) Variation of flexural models as function of the age of the sea-floor (per bins of age; see colour coding); (b) Variation of the corresponding coefficients: Values for coefficient A are shown in blue (with dashed blue 3rd polynomial fit and associated 95% confidence interval (CI)) and refer to the outer left vertical axis; Values for coefficient B are in yellow (with dashed yellow fit and 95% CI) and refer to the inner left vertical axis; Values for coefficient C are in red (with dashed red fit and 95% CI) and refer to the inner right vertical axis; Values for coefficient D are in purple (with dashed purple fit and 95% CI) and refer to the outer right vertical axis. Number of data per bin of age is shown in grey and refers to the short inner vertical axis to the right; the coefficient of determination R2 associated to the fit of flexural model equations to every sub-datasets is shown in black and refer to the short outer vertical axis to the right.
Table 5. Coefficient A, B, C & D in flexural equation as function of sea-floor age [Ma].
; R2 = 8.52% (relative to data); R2 = 33.8% (relative to 3rd polynomial fit).
Coefficient A is most probably related to age with a PCM equation (light blue dashed curve in Figure 23(b)) with equation:
; R2 = 85.6% (relative to data); R2 = 95.8% (relative to 3rd polynomial fit).
6. Topographic Elevation of Arcs and Distance Arc-Trench
Lallemand et al.  and Heuret & Lallemand  , 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
Figure 27. (a) Histogram showing the distribution of distance [˚] between arcs and trenches; the blue curve corresponds to the “best” Gaussian fit; b) Relationship between the arc-trench distance [˚] and the age of the sea-floor entering subduction [Ma] (age at trench  ); the red line corresponds to the “best” linear fit.
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.
Table 6. Statistics on arc-trench distance and arc elevation.
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  ), the power versus frequency diagram (Figure 30) shows several significant peaks. In particular, the multiple-taper spectral analysis  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).
Figure 28. (a) Histogram showing the distribution of the topographic elevation [m] of arcs; the blue distribution corresponds to arcs located on oceanic crust, with “best” Gaussian fit in dark blue; the red distribution corresponds to arcs located on continental crust, with “best” Gaussian fit in dark red; (b) Relationship between the topographic elevation of arcs [m] and the distance between the arc and the trench [˚]; the dark red line corresponds to the “best” linear fit on data for continental arcs (red crosses); the dark blue line corresponds to the “best” linear fit on data for oceanic arcs (blue crosses).
Figure 29. (a) Change of the topographic elevation [m] along the island arc [m] of the Marianas. The green curve is a polynomial fit showing the first order (main) variation; (b) Detrended variation of the topographic elevation [m] on top of which an a priori sinusoidal function (red) of amplitude A = 1000 and phase = 200 is shown.
Using all data for arcs together (including data from cordillera and from island arcs), Fourier transforms, periodograms and other techniques such as
Figure 30. (a) Periodogram corresponding to Figure 29, i.e. power spectral density against frequency (1/distance); grey, simple periodogram; black, multitaper spectrum (5 tapers); red, walch transform (see  ); (b) Same data as function of the log10 of the frequency; (c) Corresponding F-values (purple) and Degree Of Freedom (DOF, green).
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.     ); 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 (R2 = 99.45%) a hyperbolic tangent function of
Figure 31. (a) Topographic data from passive margins of the world as function of distance to their respective COB; data were selected within −9˚ to +9˚ from COB (grey crosses) and running average (red curve); yellow lines fit the segments between −9˚ & −3˚ and +3˚ & +9˚ respectively; the purple curve corresponds to the “best fit” of a hyperbolic tangent function; (b) 3D diagram showing topographic elevation as function of distance to COB and age of COB; (c) Topographic data shown as function of age of the COBs, and divided by bins of 1˚ relative to the distance of the COBs; sub-datasets are fitted using polynomial functions (see colour coding); (d) Slope of the passive margin (smoothed values from running average in blue) on top of which a double Gaussian model is fitted (pink).
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 6th 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  , 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  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  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   . 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  . 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.  with the
Figure 32. (a) Histogram showing the distribution of elevation of data corresponding to domain assumed to be continental in nature (green; geodetic data grid, same as Figure 1(a)), and same data corrected for ice-loading, post-glacial rebound and geoid anomalies (purple); (b) Histogram of the Moho depth for the same dataset as in (a) after the definition of Laske et al.  ; (c) Relationship between the Moho depth [km] (after Crust1.0) and the topographic elevation [m] (after ETopo1); red line corresponds to the “best” linear fit.
Crust.5.1 model (5 × 5 model) now updated to the Crust.1.0 model (  ; 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:
(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:
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.  is too crude to be compared with the ETopo1 model  , 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−3] as above and using this linear relationship (Equation (13)), the mean crustal density ρC shall correspond to the following equation:
where d' is the Moho depth below −4509.878 m, and hC 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−3].
However, the right relationship between topographic elevation and Moho depth is undetermined. When polynomial fits of degree 2 (for which R2 = 56.35%) or degree 6 (for which R2 = 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−3] and ρC = 2691 ± 256.3 [kg・m−3] 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 σ1) and μnC.2 = +310.2 m (±160.3 m;
Figure 33. (a) Histogram showing the distribution of topographic elevation for data assumed to lie on continental crust (corrected for ice loading, post-glacial rebound and geoid); the blue curve represents the “best” Gaussian fit; (b) Histogram for data points assumed to lie on “normal continental crust”, i.e. with topographic elevation arbitrarily limited to −333 m and +777 m; dark blue curve: single term Gaussian fit; blue curve: double term Gaussian fit; light blue curve: six terms Gaussian fit; (c) Global map showing the spatial distribution of data point for continental crust, with blue dots: topographic elevation below −333 m (thinned crust); yellow and orange dots: topographic elevation of the “normal crust” comprised between −333 m and +777 m, where lowlands (below +270 m) are in yellow, and highlands (above +270 m) are in orange; red dots: topographic elevation above +777 m (thickened crust).
1 σ2). 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
Table 7. Statistics on topographic elevation of continental crust.
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.  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.
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,
Figure 34. (a) Diagram representing the topographic elevation [m] (after ETopo1  corrected from ice loading, post-glacial rebound and geoid) as function of the Moho depth [m]  ; data are fitted using polynomial curves of degree 1 (linear): red, of degree 2: green, of degree 6: pink, of degree 10: light blue; (b) Same dataset reduced from the linear trend (polynomial fit of degree 1); (c) Same dataset reduced from the trend defined by the polynomial fit of degree 10.
Figure 35. (a) Global map of topographic elevation after reduction from linear fit (see Figure 34(b)); from blue, below zero to red, above zero (white: zero elevation); (b) Gaussian filter (18˚ × 18˚) applied on data shown in Figure 35(a); (c) Global map of topographic elevation after reduction from polynomial fit of degree 10 (see Figure 34(c)); from blue, below zero to red, above zero (white: zero elevation); (d) Gaussian filter (18˚ × 18˚) applied on data shown in Figure 35(c).
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.
Annexe 1. Statistics of elevation data (ETopo1) for every 0.025˚ bin of distance to ridge axis
 Scotese, C.R., Illich, H., Zumberge, J., Brown, S. and Moore, T. (2011) The Gandolph Project, Year Four Report: Paleogeographic and Paleoclimatic Controls on Hydrocarbon Source Rock Deposition—A Report on the Results of the Paleo-geographic, Paleoclimatic Simulations (FOAM), and Oil/Source Rock Compilation: Conclusions at the End of Year Four: Oligocene (30 Ma), Cretaceous/Tertiary (70 Ma), Permian/Triassic (250 Ma), Silurian/Devonian (400 Ma), Cambrian/Ordovician (480 Ma). GeoMark Research Ltd., Houston, TX, 219 p.
 Vérard, C., Hochard, C., Baumgartner, P.O. and Stampfli, G.M. (2015) Geodynamic Evolution of the Earth over the Phanerozoic: Plate Tectonic Activity and Palaeoclimatic Indicators. Journal of Palaogeography, 4, 167-188.
 Vérard, C., Hochard, C., Baumgartner, P.O. and Stampfli, G.M. (2015) 3D Palaogeograhic Reconstructions of the Phanerozoic versus Sea-Level and Sr-Ratio Variations. Journal of Palaogeography, 4, 64-84.
 Baatsen, M., van Hinsbergen, D.J.J., von der Heydt, A.S., Dijkstra, H.A., Sluijs, A., Abels, H.A. and Bijl, P.K. (2015) A Generalised Approach to Reconstructing Geographical Boundary Conditions for Palaeoclimate Modelling. Climate of the Past, 11, 4917-4942.
 Laske, G. and Masters, G. (1997) A Global Digital Map of Sediment Thickness. EOS Transactions American Geophysical Union, 78, F483.
 Laske, G., Masters, G., Ma, Z. and Pasyanos, M. (2013) Update on Crust 1.0: A 1-Degree Global Model of Earth’s Crust. Geophysical Research Abstracts, 15; Abstract EGU2013-2658.
 Divins, D.L. (2003) Total Sediment Thickness of the World’s Oceans and Marginal Seas. NOAA National Geophysical Data Center, Bouler, Colorado, USA.
 Tapley, B., Ries, J., Bettadpur, S., Chambers, D., Cheng, M., Condi, F., Gunter, B., Kang, Z., Nagel, P., Pastor, R., Pekker, T., Poole, S. and Wang, F. (2005) GGM02— An Improved Earth Gravity Field Model from Grace. Journal of Geodesy, 79, 467-478.
 Maus, S., Barckhausen, U., Bournas, N., Brozena, J., Childers, V., Dostaler, F., Fairhead, J.D., Finn, C., von Frese, R.R.B., Gaina, C., Golynsky, S., Kucks, R., Lühr, H., Milligan, P., Mogren, S., Müller, R.D., Olesen, O., Pilkington, M., Saltus, R., Schreckenberger, B., Thébault, E. and Caratori Tontini, F. (2009) EMAG2: A 2-Arc Min Resolution Earth Magnetic Anomaly Grid Compiled from Satellite, Airborne, and Marine Magnetic Measurements. Geochemistry, Geophysics, Geosystems, 10, Q08005.
 Paulson, A., Zhong, S. and Wahr, J. (2005) Modelling Post-Glacial Rebound with Lateral Viscosity Variations. Geophysical Journal International, 163, 357-371.
 Paulson, A., Zhong, S. and Wahr, J. (2007) Limitations on the Inver-sion for Mantle Viscosity from Post-Glacial Rebound. Geophysical Journal International, 168, 1195- 1209.
 Small, C. (1998) Global Systematics of Mid-Ocean Ridge Morphology. In: Buck, W.R., Delaney, P.T., Karson, J.A. and Lagabrielle, Y., Eds., Faulting and Magmatism at Mid-Oceanic Ridges, Geophysical Monograph 106, American Geophysical Union, Washington DC, 25 p.
 Doin, M.-P. and Fleitout, L. (2000) Flattening of the Oceanic Topography and Geoid: Thermal versus Dynamic Origin. Geophysical Journal International, 143, 582-594. https://doi.org/10.1046/j.1365-246X.2000.00229.x
 Sclater, J.G. and Francheteau, J. (1970) The Implications of Terrestrial Heat Flow Observations on Current Tectonic and Geochemical Models of the Crust and Upper Mantle of the Earth. Geophysical Journal of the Royal Astronomical Society, 20, 509-542.
 Parsons, B. and Sclater, J.G. (1977) An Analysis of the Variation of Ocean Floor Bathymetry and Heat Flow with Age. Journal of Geophysical Research, 82, 803-827.
 Winterbourne, J., Crosby, A. and White, N. (2009) Depth, Age and Dynamic Topography of Oceanic Lithosphere Beneath Heavily Sedimented Atlantic Margins. Earth and Planetary Science Letters, 287, 137-151.
 Scholer, M. (2011) Estimating Vadose Zone Hydraulic Properties Using Ground Penetrating Radar: The Impact of Prior Information. Water Resources Research, 47, W10512.
 Scholer, M. (2012) Bayesian Markov chain Monte Carlo Inversion of Geophysical Data for the Purpose of Estimating the Hydraulic Properties of the Unsaturated Zone. Ph.D. Thesis, The University of Lausanne, R007229980, Lausanne, Switzerland, 135 p.
 Lallemand, S., Heuret, A. and Boutelier, D. (2005) On the Relationships between Slab Dip, Back-Arc Stress, Upper Plate Absolute Motion, and Crustal Nature in Subduction Zones. Geochemistry, Geophysics, Geosystems, 6, Q09006.
 Burov, E. and Diament, M. (1995) The Effective Elastic Thickness (Te) of Continental Lithosphere: What Does It Really Mean? Journal of Geophysical Research, 100, 3905-3927.
 Grimaud, S., Boillot, G., Collette, B.J., Mauffret, A., Miles, P.R. and Roberts, D.B. (1982) Western Extension of the Iberian-European Plate Boundary during the Early Cenozoic (Pyrenean) Convergence: A New Model. Marine Geology, 45, 63-77.
 Boillot, G., Recq, M., Winterer, E.L., Meyer, A.W., Applegate, J., Baltuck, M., Bergen, J.A., Comas, M.C., Davies, T.A., Dunham, K., Evans, C.A., Girardeau, J., Goldberg, G., Haggerty, J., Jansa, L.F., Johnson, J.A., Kasahara, J., Loreau, J.-P., Luna-Sierra, E., Moullade, M., Ogg, J., Sarti, M., Thurow, J. and Williamson, M. (1987) Tectonic Denudation of the Upper Mantle along Passive Margins: A Model Based on Drilling Results (ODP Leg 103, Western Galicia Margin, Spain). Tectonophysics, 132, 335-342.
 Manatschal, G. (2004) New Models for Evolution of Magma-Poor Rifted Margins Based on a Review of Data and Concepts from West Iberia and the Alps. International Journal of Earth Sciences, 93, 432-466.
 Manatschal, G., Sutra, E. and Péron-Pinvidic, G. (2010) The Lesson from the Iberia- Newfoundland Rifted Margins: How Applicable Is It to Other Rifted Margins? Central & North Atlantic Conjugate Margins Conference, Vol. II, Lisbon, 27-37.
 Flament, N. (2006) Evolution à long terme de l’altitude moyenne des continents. Report of the Ecole Normale Supérieure de Lyon (ENS Lyon), France, 39 p.
 Reigber, C., Schmidt, R., Flechtner, F., Konig, R., Meyer, U., Neumayer, K.-H., Schwintzer, P. and Zhu, S.Y. (2005) An Earth Gravity Field Model Complete to Degree and Order 150 from Grace: Eigen-Grace02S. Journal of Geodynamics, 39, 1-10.
 Kaban, M.K., Schwintzer, P. and Reigber, C. (2005) Dynamic Topography as Reflected in the Global Gravity Field. In: Reigber, C., Lühr, H., Schwintzer, P. and Wickert, J., Eds., Earth Observation with Champ: Results from Three Years in Orbit, Springer, Berlin, 199-204.