Time scales and gaps, Haar fluctuations and multifractal geochronologies

TL;DR

This paper introduces a new method to analyze geochronology data using Haar fluctuation analysis, revealing how measurement density affects the interpretation of Earth's history.

Contribution

The paper introduces geochronology measurement densities as a new paleoindicator and uses Haar fluctuation analysis to characterize their scaling properties.

Findings

Measurement density is a new paleoindicator that is typically correlated with primary indicators.

Haar fluctuation analysis reveals scaling regimes and exponents spanning nine orders of magnitude in time.

Measurement density characteristics are essential for unbiased statistical interpretations of geochronological data.

Abstract

Outcrops and cores are primary sources of information about the Earth’s past. Quantitative analyses rely on geochronologies that take into account highly variable sedimentation and erosion rates as well as gaps from missing strata. Using 23 geochronologies from the Holocene, Quaternary, Phanerozoic and Precambrian, we apply Haar fluctuation analysis to statistically characterize the number of measurements per unit time - the measurement densities. The analysis determines the densities’ (multifractal) scaling regimes and exponents; collectively, the analyses span over nine orders of magnitude in time scale. The measurement density is a new paleoindicator that we show is typically correlated with the primary paleoindicator, biasing and complicating its statistical interpretation. We also analyze the distribution of gaps linking the latter’s (probability) scaling with series incompleteness and the length Sadler effect. The density characteristics are needed to unbias spectra and other statistical characterizations.

Geochronology measurement densities (i.e., the number per time) are a new paleo indicator, as statistically characterized using Haar fluctuation analysis, determining the densities’ scaling regimes, exponents and their correlations with the primary paleoindicators.

Figures

Temperature data recorded in a paleotemperature reconstruction^[56](#CR56)^ (blue) and the corresponding measurement density as a function of time (red).Large gaps in the data correspond to small measurement densities and vice versa. This dataset was a benthic stack at 10kyr nominal resolution, such that data from well-defined layers were used with many data points assigned to each 10kyr interval. There were 11874 data points for 3950 different 10kyr intervals over 470 Myr. The mean density was 25.3/Myr, see also Fig. [2](#Fig2) and S2.

Results of the Haar analysis of representative datasets.Top: The normalized mean Haar fluctuation <Δ $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho \left(\triangle t\right)$$\end{document}$ >, the logarithmic slope of which gives an estimate of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepack

A comparison of the root mean square (RMS) fluctuations of the Isson^[33](#CR33)^ δ^18^O in Carbonates (red, roughly calibrated in K) with the normalized measurement density (blue) and the correlation between the two (*R*(∆*t*), brown, linear vertical scale at right).There were 25399 data points, the oldest being 3.504 Gyrs, so that this fluctuation analysis takes us close to the origin of the Earth. As described in the text, the Isson δ^18^O were roughly calibrated so that the results are expressed in units of *K*. The density fluctuations are dimensionless as is *R*. The dashed reference lin

The probability distribution *Pr*(*τ*′ > *τ*) of intervals *τ* between consecutive datapoints.The black dotted lines show a fitted Gaussian distribution, whereas the straight black lines show the best-fit linear curve at the extremes of the data. The exponents $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{D}$$\end{document}$ are shown for the linear log-log fits, where $\documentclass[12pt]{minimal} \usepac

Full text

Introduction

Much of our knowledge of Earth’s past is gleaned from interpreting the stratigraphy of sedimentary, igneous, and metamorphic sequences, whether exposed in outcrop, recovered from cores through lake and ocean sediments, or preserved in crystalline crust. These were laid down under highly variable conditions, and their stratigraphy generally includes temporal gaps of all sizes^1^. Since erosion, sedimentation, aggradation, denudation and progradation rates of sediment sources and sinks operate over wide ranges of scale^2^, the simplest hypothesis is that they are scaling in accord with scaling Holocene and Quaternary climate records^3–13^ and analyses of Phanerozoic benthic stacks^14^ and macroevolution^15^ physically reflecting interconnections across various temporal scales and planetary geochemical reservoirs. In this study, we analyze chronologies and proxies from records spanning the Holocene, Quaternary, Phanerozoic and (much of) the Precambrian, focusing on the measurement density—i.e., the number of measurements per unit time—that we quantitatively characterize in a multifractal framework. This density is a new and distinct multi-scale proxy that can be used to understand the evolution of Earth’s system history. Implicitly, it includes information on the absence of evidence, thus— paradoxically—providing positive information about the structure of non-recorded history, a significant phenomenon that can be quantified in its own right. We also analyze the scaling of the probability distributions of the inter-measurement time intervals and show that completeness of the records decreases with the number of measurements, thus quantitatively explaining the length Sadler effect^2,16–19^.

Cores and outcrops are typically sampled regularly in depth; the resolutions are determined both by available technology and by their physical characteristics, such as the effects of bioturbation (various oceanic and terrestrial sediments), molecular diffusion (isotopes in ice cores), and material redistribution at the surface (e.g., by winds or currents). However, even when regularly sampled in depth, geochronology is typically non-uniform. Although adequate for identifying and dating strong transitions (“jumps”), or intense events (“spikes”), this can pose serious problems for statistical analysis (including for Lomb periodograms^20^).

Today, series with hundreds of thousands of measurements are increasingly available; therefore, statistical analysis is mandatory. How should we handle highly variable temporal resolution? Either one adapts the method to the data or the data to the method. An example of the former is Haar fluctuation analysis^21^ (see below) that characterizes how fluctuation amplitudes change with time scale and can be conveniently implemented on nonuniform chronologies^14^. Alternatively, one may create a regular timeline, interpolate the data, and then use conventional techniques, such as Fast Fourier Transforms. Although convenient and commonly employed, interpolation can lead to uncontrolled biases^22^. The problem is that the data are typically “rough”: fractal, multifractal, and nondifferentiable, whereas standard interpolation functions are smooth (differentiable of order 1 or more), which leads to biases^8,20^.

Regardless of the analysis strategy, it behooves us to study, characterize, and understand the statistics of the chronologies themselves. This knowledge is required, if only to allow us to assess and statistically correct the biases. Therefore, initial characterizations of inhomogeneities should attempt to characterize the statistical properties of chronologies over a wide range of scales, which requires scaling techniques. For records with only a few hundred measurements, it may be sufficient to theorize inhomogeneity in terms of geometric fractal sets. If time is discretized at the finest resolution, then at any location on the time axis, there is a binary choice: either there is a measurement or none. Either the data or their complements (holes and gaps) may form a fractal set. Such a black/white binary reduction is the temporal analog of treating spatial distributions of meteorological measuring stations as a fractal set, which implies that the measurements have not only limited spatial but also limited dimensional resolutions^23,24^. However, with big data, a more complete characterization is needed, which requires us to follow the spatial example^25^ that instead considers the fundamental quantity to be the measurement density ρ(t) which is a mathematical function with a value at each point on the time axis (more precisely, if it is scaling, ρ(t) is in fact the density of a singular multifractal measure—it is a “generalized function”). Figure 1 shows an example with 11,874 points, showing that ρ(t) may have huge variability, displaying typical multifractal “spikes” (see the “spike plots” in ref. ^26^, see also Fig. S1). The relationship between the fractal set and multifractal density descriptions was clarified in ref. ^27,28^, and discussed in the climate context^29^.Fig. 1. Temperature data recorded in a paleotemperature reconstruction^56^ (blue) and the corresponding measurement density as a function of time (red).Large gaps in the data correspond to small measurement densities and vice versa. This dataset was a benthic stack at 10kyr nominal resolution, such that data from well-defined layers were used with many data points assigned to each 10kyr interval. There were 11874 data points for 3950 different 10kyr intervals over 470 Myr. The mean density was 25.3/Myr, see also Fig. 2 and S2.

This paper is a collaboration of specialists in Holocene, Quaternary, Phanerozoic and Precambrian time scales as well as nonlinear geophysics, presenting an empirical study of geochronology measurement densities, discussed here over nine orders of magnitude in scale (Table 1 and SUPP). By analyzing ρ(t), we show how chronologies—including the gaps—can be understood in a scaling multifractal framework. We also determine the scale-by-scale correlation of fluctuations Δρ(Δt) with fluctuations in the paleo indicators themselves (ΔT(Δt)). If Δρ(Δt) and ΔT(Δt) are statistically independent, then statistically correcting the biases is easier. Alternatively, correlations or similarities in scaling may indicate some ontological causal connections between the time series and the cross-scale unevenness of their measurement densities.Table 1. The datasets that were analysed in this investigationTime PeriodDataProxy typeShort name, ReferenceLocationTime rangeResolution(nominal, in years)Data(°N,°E)pointsHoloceneTemperatureSpeleothemEurSpain^41^Oregon^42^2025-11-10 11:04:00 AM(43, 4)(42.1, −123.41)3972–23 yr8038–247 yr118752678PollenGonghai^43^HOB^43^Kettle^43^, Moosent^43^(38.9, 112.23)(47.69, 9.01)(48.61, −103.6)(47.02, 7.48)15–0 kyr13–0 kyr13–0 kyr19–0 kyr6666777860551513QuaternaryDustIce CoresEDC^44^Talos Dome^45^,NGRIPRECAP^46^(−75.1, 123.35)(−72.82, 159.18)(75, −42.3)(71.3, −26.72)801–0 kyr150–4 kyr108–10 kyr121–5 kyr1143029820510339727742317Marine SedimentCentral PacificCentral South Pacific^47^(4.68, -160.05)(-54.22, -125.43)141–7 kyr474– 1 kyr60402892384LoessXifeng^48^(35.7, 107.6)801–1 kyr200722TemperatureIce CoresFuji Dome^49^EDC^50^EDML^51^WAIS^52^, NEEM^53^, 11/10/25 11:04:00 AM(−77.3, 39.3)(−75.1, 123.35)(−75, 0.07)(−79, −112)(77.45, 51.06)340–0 kyr801–0 kyr150–1 kyrv68–0 kyrv129–2 kyr705710.1118957872303637524434LakeTanganyika^54^(−6.65, 29.8)59–1 kyr60210MarineAlkenone^54^(37, −123)161–3 kyr300284Temperature AnomalyIce CoresVostok^55^(75.1, −42.3)423–0 kyr93309Phanerozoic,PrecambrianTemperatureBenthic StackGrossman^56^-470–0 Myr6 k11874Sea LevelMarine SedimentHaq Sea Level^31^-128 Myr20 k617δ^18^O CarbonatesAssemblageIsson^33^-3.504 Gyr-01 k25399The sample datasets used in the main body of the study were: Moosent, NEEM, South Pacific, EDC, Grossman, and Haq Sea Level. With the exception of the (long-time) benthic stack (Grossman), sea levels (Haq) and assemblage Isson (bottom rows), the shorter time scale series were all from single cores. The nominal resolution of the Isson assemblage is highly depend on the age, see Fig. S2; 1 kyr used in Fig. 3 is a compromise.

The science of stratigraphy has recognized for a long time the significance of temporal gaps in diverse records and across time scales. In concert with the lithological and facies succession features, the genetic discontinuities between neighboring rock formations are used in identifying the presence of gaps in stratigraphical records, and thus defining sequences in a correlational framework of the sequence stratigraphy, which searches for the best ways to objectively distinguish space-time “packages” of the geological record^1^. A precursor to sequence stratigraphy—allostratigraphy—was and still is explicitly concerned with revealing the ranges of discontinuities and gaps in the geological record, which are later used in correlation between distant areas^30^. Gaps appear at all time scales, and depending on their magnitude can be graded from diastems (relatively short interruptions on geological time scales) to disconformities in relatively continuous and lithologically homogenous successions, to angular unconformities which signify not only gaps in the record but also structural and tectonic deformations which happened in the unobserved time, and to complete nonconformities when sedimentary rocks or sediments are found to be overlaying much older igneous and metamorphic rocks, which usually signal gaps with durations reaching 100 Myrs to Gyrs, Therefore the classical stratigraphy already implicitly recognized the scaling and scale free nature of the gaps in the record and chronologies. The present study quantifies this structure of the gaps as a function of time scales over which they occur, while also exploring their correlations with proxies, which show clear scale-dependence of these correlations. The presence of correlations between proxies and the gaps in their records suggests two things: i) proxy record is systematically biased in non-trivial ways with the scale-dependent magnitudes of distortion; ii) there should be a common cause, implicit in both allo- and sequence stratigraphical approaches, which drives paleoclimatic signals and the formation of gaps in their records. Therefore, paraphrasing the famous aphorism: “the absence of evidence is also the evidence of absence”. It means that the time scale-dependent (scaling) structure of gaps in records of proxies can be used as a new and mathematically tractable source of information in revealing the Earth system processes.

Results

Comparisons using the ensemble of data over scales from years to Giga years

We selected 22 paleoclimate series from individual cores covering the Holocene (past 12 kyr), upper and middle Quaternary (past 0.8 Myr). In addition, for the longer time scales, we chose a benthic stack (compiled from numerous ocean cores past 470 Myrs, Fig. 1) as well as a Precambrian δ^18^O carbonate (past 3.054 Gyrs, Fig. 3 and SUPP). These longer time scales were completed with a sea level series based on cores and outcrops^31^ (Table 1, full analyses in the SUPP). Unlike the single-core data, these longer series were therefore able to fill some of the gaps, although many still remain.

To estimate the density ρ(t), we extracted the chronologies and discretized the time axis into bins at the highest nominal data resolutions (see Fig. S2 and discussion), at which scale ρ(t) is the (binary) indicator function (zero in bins with no data, 1 in bins with data). These densities were then analyzed using Haar fluctuation analysis^21^ (based on the Haar wavelet). The Haar fluctuation of ρ(t) at resolution Δt is simply the difference between the average over the first and second halves of the interval: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \rho \left(\Delta t\right)=\overline{{\rho }_{\left[t,t-\Delta t/2\right]}\,}-\overline{\,{\rho }_{\left[t-\Delta t/2,t-\Delta t\right]}}$$\end{document} (overbars for temporal averages, below - except for the correlations – we consider absolute fluctuations).

In a scaling regime, the mean fluctuation is a power law: <Δρ(Δt)> ≈ Δt^H^ where H is the fluctuation exponent and “<>” indicates ensemble averaging, here estimated by averaging over all the available disjoint intervals at scale Δt. H characterizes how the mean fluctuations decay (H < 0) or grow (H > 0) with scale; it is one of an infinite (multifractal) hierarchy of exponents. However, the hierarchy itself can typically be characterized by only two additional exponents: the intermittency exponent C1 and the multifractal index α (Methods). C1 characterizes spikiness and its sharp transitions (e.g., Fig. 1); Gaussian processes are nonintermittent with C1 = 0. C1 was estimated as the logarithmic slope of the function F(Δt) (see Methods), which is related to the ratio of the mean to the RMS fluctuation: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \left( < \Delta \rho \left(\Delta t\right) > / < \Delta \rho {\left(\Delta t\right)}^{2}{ > }^{1/2}\right)\,=a\log F\left(\Delta t\right)$$\end{document} where ½≤a ≤ 1 is a weak function of α (Methods). Finally, the value of the multifractal index α is theoretically restricted^32^. to the range 0 ≤  α  ≤ 2 when H = 0, α = 0 corresponds to a binary monofractal “β model”, and α = 2 to a pure “lognormal” multifractal.

The exponents may be physically interpreted using multifractal sedimentation and erosion models of ρ(t). For these models, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H$$\end{document} controls the basic scaling of the sedimentation rate, whereas the intermittency exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{1}$$\end{document} controls sharp gradients, “spikes” in ρ(t). Combined with α, they also control the statistics of the low ρ(t) regions (where there are long intervals between measurements). However, an additional explicit erosional model may be required to model these gaps, as shown below.

The results of the Haar analysis for normalized mean fluctuations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ < \Delta {{\rm{\rho }}}\left(\Delta t\right) > /\bar{{{\rm{\rho }}}}$$\end{document} are shown in Fig. 2, top. To simplify the figure, only representative datasets are shown; for the remaining series, see the SUPP. Despite often disparate origins (sediments, ice, and speleothems), the series shows two regimes. At higher frequencies, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H\approx -0.5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{1}\approx 0$$\end{document} (the parameters of Gaussian white noise, note that mathematically, 1-D series must have 0 ≤ C1 ≤ 1). H < 0 implies that when averaged over progressively longer timescales, ρ(t) tends to converge to well-defined values. Maximum convergence—often at scales Δt a hundred or more times larger than the finest resolution—occurs at a transition scale τ_c_ beyond which H is positive, so that fluctuations now increase with timescale: the mean measurement density itself is no longer well-defined. At lower frequencies (Fig. 2 middle), we see that this new “wandering” (H > 0) regime (with H typically ≈ 0.2–0.3, Fig. 2 top) has C1 ≳0.1 (Table 2). C1 quantifies the rate at which intermittency builds up with scale, and the empirical values (often ≈ 0.2) are substantial. In comparison, atmospheric turbulence typically^6^ has C1 ≈ 0.05– 0.1, yet in spite of these apparently small values, at human scales, the atmosphere appears highly intermittent because the variability starts at planetary scales, ten million times larger. Finally (Table 2, SUPP), we find that the multifractal index α also transitions at Δt ≈ τ_c_ (from α ≈ 0 to α ≈ 1.0 - 1.5).Fig. 2. Results of the Haar analysis of representative datasets.Top: The normalized mean Haar fluctuation <Δ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho \left(\triangle t\right)$$\end{document} >, the logarithmic slope of which gives an estimate of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H$$\end{document} . The two lower gray curves show the fluctuations in temperature for the EDC (left) and Grossman (right) datasets, showing that macroweather, climate, macroclimate, and megaclimate regimes are quite similar for ρ and for T. Middle: plot of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\log }_{10}F\left(\triangle t\right)$$\end{document} , as defined in the methods section, the slope of which gives an estimate for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{1}$$\end{document} . This plot differs from the log-log plot of the ratio between the fluctuation in measurement density and its root mean square by a constant coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a$$\end{document} for constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} , which is shown for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =1.3$$\end{document} on the right axis. Bottom: the correlation coefficient between fluctuations in measurement density and fluctuations in the measured quantity. The solid lines are absolute correlations, and the dotted lines show negative correlations where relevant. The vertical gray regions in all three plots show the transition timescales between the scaling regimes in the climate system.Table 2. Estimates for the critical timescale \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{C}$$\end{document} , the fluctuation exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H$$\end{document} , the intermittency exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{1}$$\end{document} and the multifractal index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} for all datasets analysedTimePeriodDataProxy typeShort name, Reference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\boldsymbol{\tau }}}}_{{{\boldsymbol{C}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{H}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{\boldsymbol{C}}}}_{{{\bf{1}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\alpha }}}$$\end{document} (yr) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle {{\boldsymbol{t}}} < {{{\boldsymbol{\tau }}}}_{{{\boldsymbol{C}}}}\,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle t > {{{\boldsymbol{\tau }}}}_{{{\boldsymbol{C}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle {{\boldsymbol{t}}} < {{{\boldsymbol{\tau }}}}_{{{\boldsymbol{C}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle t > {{{\boldsymbol{\tau }}}}_{{{\boldsymbol{C}}}}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle {{\boldsymbol{t}}} > {{{\boldsymbol{\tau }}}}_{{{\boldsymbol{C}}}}$$\end{document} HoloceneTemperatureSpeleothemEurSpain^41^Oregon^54^6040−0.5−0.30.30.3000.10.31.11.2PollenGonghai^43^HOB^43^Kettle^43^, Moosent^43^2504001600400−0.5−0.4−0.2−0.80.40.40.30.40.10000.20.20.30.21.11.11.11.0QuaternaryDustIce CoresEDC^44^Talos Dome^45^,NGRIPRECAP^46^,40032050006300−0.5−0.2−0.10.30.30.20.4−0.300000.10.20.50.21.51.21.31.0Marine SedimentCentral Pacific Central South Pacific^47^40004000−0.3−0.50.40.40.200.20.21.21.2LoessXifeng^48^8000−0.30.30.10.11.3TemperatureIce CoresFuji Dome^49^EDC^50^EDML^51^WAIS^52^, NEEM^53^,16,00025001600400110−0.5−0.5−0.5−0.5−0.50.50.50.60.60.3−0.300000.50.20.10.10.30.81.41.51.41.5LakeTanganyika^54^-−0.10.1-MarineAlkenone^54^-−0.10.1-Temperature AnomalyIce CoreVostok^55^3200−0.50.500.41.4Phanerozoic,PrecambrianTemperatureBenthic StackGrossman^56^12 M−0.50.200.21.3Sea LevelMarine sedimentHaq Sea Level^31^20 M−0.20.20.10.21.3δ^18^O Carbonates(see the SUPP)AssemblageIsson^33^3–10 M−0.1−0.050.20.5–1.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{1}$$\end{document} estimates are for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle t < {\tau }_{C}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle t > {\tau }_{C}$$\end{document} , whereas only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle t > {\tau }_{C}$$\end{document} is shown because for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\triangle t < {\tau }_{C}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} transitions from the monofractal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \approx 0$$\end{document} to a multifractal value. Neither the Tanganyika nor the Alkenone datasets undergo a scaling regime transition, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} is observed to increase across all timescales for both these datasets, hence, an estimate for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} is not included for these two datasets.

Figure 2 also shows (gray line) the fluctuation analysis - including the qualitative H < 0 and H > 0 transitions - in paleotemperatures^14^ T(t). This figure demonstrates that the temperature and analogous ρ(t) regimes are qualitatively close; in any given regime, they have the same sign of H, such that T(t) and ρ(t) are both convergent or divergent together.

Correlations between measurements and measurement densities

Each proxy record has two potentially independent pieces of information: the paleoindicator T(t) and the density at which it is estimated ρ(t). Figure 2 bottom, shows the normalized correlation coefficients (R(Δt)) between fluctuations Δρ(Δt) and ΔT(Δt). At small Δt, R(Δt) of each dataset is typically low, indicating statistically independent processes. However, at Δt > τ_c_, correlations generally grow stronger, often roughly as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\left(\Delta t\right)\propto \log \Delta t$$\end{document} linear on the log-linear plots shown: (Fig. 2 bottom).

These correlations make the statistical unbiasing and interpretation of paleoindicators more difficult. As R(Δt) increases, Δρ and ΔT tend to vary together, such that when both are large, paleoindicators tend to be very frequently sampled, and when both are small, paleoindicators tend to be rarely sampled. Both effects lead to biases. This situation arises, for example, in ice cores: the precipitation and hence sedimentation rate are temperature-related, implying Δρ(Δt)-ΔT(Δt) correlations at these time scales (see e.g. the NEEM curve in Fig. 2 bottom). For most series, the sign of R(Δt) was the same at all Δt values, and because the qualitative effects of R ≈ 1 and R ≈ –1 are the same, Fig. 2 bottom shows the absolute correlation (solid). However, in some cases, the sign of R(Δt) did change with Δt, indicated by the dashed lines.

A Precambrian example covering most of the Earth’s history

Figure 2 was limited to Phanerozoic data, hence with a maximum Δt of several hundred Myrs. Let us now briefly consider the Isson^33^ δ^18^O Carbonate series that is significant because it spans most of the Earth’s history (3.045 Gyrs with 25399 measurements, Fig. S1). Being an assemblage from different authors, locations, it is much more heterogeneous than the other series considered here (for a discussion of some of the epoch specific issues, see ref. ^34^). For example, over the most recent 50 Myrs there are many segments with nominal resolutions of 10 kyr and even 100 yr (Fig. S3), yet the parts older than 600 Myrs (roughly the Phanerozoic), generally have resolutions of 1 Myr. This implies that the smaller Δt statistics come exclusively from the Phanerozoic. Since these represent 81% of the measurements, to a good approximation, the statistics for Δt ≈ <600 Myrs characterize only the Phanerozoic (see Figs. S1–7).

With these caveats in mind, Fig. 3 shows the RMS fluctuations of the δ^18^O values (red, multiplied by 4, see below), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ < \Delta \rho \left(\Delta t\right) > /\bar{\rho }$$\end{document} (blue) and the correlation R(Δt) (brown). Significantly, the δ^18^O fluctuations are scaling over the entire megaclimate regime from ≈ 0.5 Myrs <Δt <≈ 0.5 Gyrs, and the density fluctuations (slope −0.05) are also scaling over the same range. The sudden transition at Δt ≈ 0.5 Gyrs agrees with the upper limit of the megaclimate inferred from analysis of the (biostrata-based) Geological Time Scale^29^.Fig. 3A comparison of the root mean square (RMS) fluctuations of the Isson^33^ δ^18^O in Carbonates (red, roughly calibrated in K) with the normalized measurement density (blue) and the correlation between the two (R(∆t), brown, linear vertical scale at right).There were 25399 data points, the oldest being 3.504 Gyrs, so that this fluctuation analysis takes us close to the origin of the Earth. As described in the text, the Isson δ^18^O were roughly calibrated so that the results are expressed in units of K. The density fluctuations are dimensionless as is R. The dashed reference lines with slopes 0.3 and -0.05 are close to the RMS temperature fluctuations and normalized density fluctuations and cover the megaclimate ( ≈ 0.5 Myrs to ≈0.5 Gyrs).

It is not clear how to interpret the larger Δt decrease in variability; Fig. S4 shows that the Phanerozoic and Precambrian variabilities had the same exponent up until about 40 Myrs, but differ at the larger Δt, with the Precambrian having slowly decreasing variability at larger Δt. This indicates that these longest time scales are marginally stable in contrast to the H > 0 (unstable) Phanerozoic part. Even considering the largest Δt (2.5 Gyrs corresponding to 2 (largely) disjoint intervals), this regime’s range of scales (factor of 2.5/0.5 ≈ 5) is probably too short to be meaningfully interpreted as a new dynamical regime.

Returning to the megaclimate range, its exponent H ≈ 0.3 is nearly identical to those from both the Veizer^35^ and Zachos^36^ benthic stack paleotemperatures (analyzed in ref. ^14^), so that their statistics differ from the Isson series by a multiplicative constant. The paleotemperature stacks can thus be used to statistically calibrate the Isson assemblage. In Fig. 3, the latter were multiplied by 4 to be close to the Veizer paleotemperatures (for Zachos, the factor is ≈ 2).

Considering the measurement densities, they show as an inflection at scales ≈ 3 Myrs. While this is much less pronounced than for those in Fig. 2 (top), even flat (H ≈ 0) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ < \Delta \rho \left(\Delta t\right) > /\bar{\rho }$$\end{document} curves imply quite nontrivial statistics: long-range statistical dependencies, and hierarchical clustering of density spikes (e.g., Figs. 1 and S1). We may also note that the density and the values are essentially uncorrelated until there is a fairly sharp rise in R at roughly the age of the Phanerozoic (≈ 0.5 Gyrs), indicating that the interpretation of the longer Δt statistics has some biases. Also note that since there are a few large Δt segments, the correlations are high at the very extreme Δt’s (see the discussion in the SUPP).

Measurement intervals, measurement gaps and the length Sadler effect

There is a well-known tendency for longer and longer paleo records to become less and less complete due to missing strata, we’ll call this the length “Sadler effect”^16^ to distinguish it from the original and related “resolution Sadler effect” in which completeness changes with the resolution which is controlled by a different exponent (see ref. ^29^). To investigate this, we determined the probability of a random intermeasurement interval τ’ exceeding a fixed τ: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Pr \left(\tau \mbox{'} > \tau \right)$$\end{document} (Fig. 4, S7 for Isson). We see that at small τ, the distribution is (often) approximately Gaussian (dashed black lines). However, at large τ, there is a transition to a scaling distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Pr \left(\tau \mbox{'} > \tau \right)\propto {\tau }^{-{q}_{D}}$$\end{document} , indicating a hierarchy of extreme intervals - gaps with small qD implying the existence of huge gaps. For example, for the NEEM, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{D}\le 1$$\end{document} so that the longest gap in a chronology is of the same order as the sum of all shorter gaps.Fig. 4. The probability distribution Pr(τ′ > τ) of intervals τ between consecutive datapoints.The black dotted lines show a fitted Gaussian distribution, whereas the straight black lines show the best-fit linear curve at the extremes of the data. The exponents \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{D}$$\end{document} are shown for the linear log-log fits, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Pr \left({\tau }^{{\prime} } > \tau \right)\propto {\tau }^{-{q}_{D}}$$\end{document} . Dashed red reference lines have slopes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{D}=1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{D}=2$$\end{document} corresponding respectively to diverging means and variances. The NEEM data were particularly extreme with qD ≈ 0.6, and since the Gaussian fits were obtained from the means and standard deviations, in this case, the Gaussians are particularly far from the data. It is worth noting that the two Phanerozoic datasets were taken from multiple locations: the Grossman series is a benthic stack, and the Haq Sea Level curve was assembled from different outcrops in various basins. This effectively fills in holes such that they are more complete and have fewer gaps. See Fig. S7 for analysis of the Isson assemblage (Fig. 3) with qD ≈ 1.5.

To see how the power law gaps explain the length Sadler effect, consider a series of N measurements. If qD > 1, the mean interval will converge so that the overall length is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{N}\propto N$$\end{document} . The extreme gap will have the minimum probability Prmin ≈ 1/N, so that the extreme (maximum) gap length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{\max }\propto {\Pr }_{\min }^{-1/{q}_{D}}\approx {N}^{1/{q}_{D}}{\propto \tau }_{N}^{1/{q}_{D}}$$\end{document} – and therefore incompleteness (maximum gap) – grows with length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{N}$$\end{document} , hence the length Sadler effect (if qD ≤ 1, the effect is even stronger).

Discussion and conclusions

Much of our knowledge of the evolution of life and the planet is inferred from stratigraphic records that are marred by hugely varying sedimentation rates and gaps where swaths of data are missing. The emergence of “big” paleodata makes it particularly necessary to characterize this variability, not only to statistically correct biases but also for understanding and modeling, including the underlying dynamical regimes (defined by the scaling properties) and that constitute the “continuum” variability. For this purpose, we exploit a new paleoindicator, the temporal measurement density ρ(t). To isolate the statistics at different timescales Δt, we used Haar fluctuation analysis to statistically characterize fluctuations Δρ(Δt) from 24 paleoseries, collectively spanning a timescale range Δt from years to 2.5 Gyr. Although the chronologies are products of both objective (and complex) physical processes as well as subjective technical data extraction issues, it is important to objectively determine their statistics, here, through their fluctuation statistics.

Analysis revealed that the series shared various scaling characteristics. Partial exceptions were the Haq paleo sea level and Isson ^18^O Carbonate assemblage, which used data from multiple locations to help fill in gaps. The series could be divided into distinct high- and low-frequency scaling regimes separated by a transition time scale τc. In each regime, we quantified the statistics using the three basic multifractal parameters H, C1, and α. At high frequencies, Δt ≤ τc, we found H ≈ −0.5, C_1_ ≈ 0, α ≈ 0 (Gaussian white noise), whereas at low frequencies H ≈ 0.2–0.3, C_1_ ≈ 0.1–0.25, α ≈ 1–1.5 (Figs. 2 and S4). We note that in certain cases (notably ice cores), the layers are increasingly compacted with depth, so that the chronologies are not statistically homogeneous. However, as shown in the SUPP, this is a low-frequency effect that can be largely removed by polynomial detrending, so that this sediment compression does not impact the results of this study. The benthic stack assemblage and sea level series that combined data from multiple locations to help fill in gaps had similar qualitative behavior, but without evidence for a Gaussian regime.

Overall, four ρ(t) regimes were identified (Fig. 2), which were in close agreement with the paleotemperature (T(t)) dynamical regimes identified in ref. ^14^ on the basis of the sign of H: macroweather (several weeks to centuries to millennia), climate (up to 50–100 kyr), macroclimate (up to 500 kyrs) and megaclimate (from ≈ 0.5–3 Myr; up to ≈ 0.5 Gyrs (Figs. 2, 3 and SUPP). Beyond Δt ≈ 0.5 Gyrs, the Isson analysis showed that fluctuations in T(t) transition from increasing to decreasing (stable) behavior (the ρ(t) fluctuations are roughly constant). The general agreement between the ρ(t) and T(t) transition scales (and hence dynamical regimes) suggests that these have a deep significance, each corresponding to distinct nonlinear scaling processes governing both the climate system and erosion and sedimentation processes. While these scaling processes are presumably nonlinearly connected with each other, there is no implication that one controls the other.

H < 0 implies that as Δt increases, ρ(t) approaches a well-defined value: stable behavior, whereas at lower frequencies, H > 0, implying that fluctuations grow with scale, ρ(t) seems to “wander”, becoming strongly scale dependent. In this regime, the intermittency is strong (C1 ≳0.1–0.2), so that ρ(t) has strong transitions (“jumps”, “spikes”).

We also determined the scale-by-scale correlations R(Δt) between Δρ(Δt) and ΔT(Δt), which are important for unbiasing the paleostatistics. It was found that in the H < 0 regime, R(Δt) was generally small (statistical independence); however, for larger Δt (H > 0), R(Δt) increased roughly logarithmically with Δt.

At low ρ, the discrete nature of the chronologies is important, and it is necessary to quantify the gaps where the measurement interval τ is particularly large. For this, we determined the probability distribution of τ, finding that there were typically abrupt transitions separating small and large τ, distinguishing “ordinary” intervals from those associated with missing strata that had power law probability “tails” with exponent qD. These distributions imply that the length of the extreme gaps increases as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\tau }_{N}^{1/{q}_{D}}$$\end{document} where τ_N_ is the series length, implying that longer records have larger gaps and are less complete: the “length Sadler effect”.

Our Haar fluctuation analysis directly handles inhomogeneous data (for periodic structures, see ref. ^37^) and has implications for statistical analyses. For example, if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ < \Delta T > \propto \Delta {t}^{H}$$\end{document} , then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\, < \Delta T > /\Delta t\propto \Delta {t}^{H-1}$$\end{document} so that the first derivative diverges when H < 1 (as we find here). More generally, H is the highest order of (fractional) differentiation; therefore, linear (H = 1) or polynomial (H ≥ 2) interpolations will be biased since the resulting interpolated series will combine some sections with H < 1 with others H ≥ 1, and this, in an uncontrolled manner. Non-interpolative techniques such as Lomb-Scargle spectral analysis are generally strongly biased when applied to scaling data^22^. This is because scaling spectra are power laws: ω^−β,^ and when β>0, the spectra are dominated by low frequencies. Ignoring intermittency corrections β = 1 + 2H, so that problems arise whenever H > −1/2 (as found here). Other methods, such as the MultiTaper Method (MTM), are only theoretically justified for uniformly spaced data and are often unhelpful^20^. Future work is thus needed both to model ρ(t) from more fundamental sedimentation and erosion processes, and to find optimum techniques for statistically correcting biases.

In the future, Haar fluctuation and trace moment analyses of the fossil record and the stratigraphical successions of environmental variables could be made in a wide variety of other systems, time periods and time scales. Ocean cores provide especially rich records of data and their gaps. In addition, the exploration of extremely deep time composite numerical records of hundreds of geological sections (e.g., Riedman and Sadler, 2018) could give new insights into the processes as well as biases, especially with respect to Precambrian biospheres.

In his seminal work “On the Origin of Species,” Charles Darwin (1859) characterized the geological record as a book where “…history of the world imperfectly kept… entombed in our consecutive, but widely separated, formations”. Here, invoking the power law scaling of extreme intervals (especially in the deep-time megaclimate regime), whose gaps become more extreme with longer records, we can add precision to his words: the gaps in paleo-proxy records become more pronounced at the longest time scales, where the most significant and dramatic changes in biota and the physical Earth become the most apparent.

Correlation analyses suggested a causal connection between the frequency of gaps and indicators of different Earth system parameters, such as surface or deep-water temperatures. This relationship can be seen as a bias in the estimation of numerical patterns (as it really is for many purposes), but also as an indication of the common cause, as was interpreted in the connection of the uneven formation of fossiliferous sedimentary packages and macroevolutionary processes^38^. Deep-sea records suggest that the amount of preserved sediments and all kinds of information recorded in them are dependent on climatic states^37^. In a sense, this is not surprising because a major branch of stratigraphy is based on the correlation of sedimentary packages bounded by sedimentary gaps—the sequence stratigraphy—which is implicitly based on the assumption of spatio-temporal coherency of climate effects (through the changes in the sea level and sedimentary facies) on the gaps in the sedimentary record^39^. Moreover, it was recognized that the gaps in the fossil record (and other stratigraphic records) can be seen not only as a bias of a quantity we want to characterize as accurately as possible, but also as a structure and feature worth exploration and explanation in its own right^40^. Therefore, the appearance of similar scaling regimes and correlations between paleoindicators and the densities of gaps in their records, but also the systemically changing nature of these correlations, reveals this new measure ρ(t) as a new and largely unexplored indicator of past geodynamic, paleoclimatic, paleoecological, and evolutionary changes. At least two quantities need to be statistically characterized—the paleoproxy value itself, and also the temporal densities of values per unit time, which can provide additional and possibly non-redundant information about the paleogeography, and can be further analysed in numerous ways.

Methods

Climate science has become increasingly specialized so that distinct (if overlapping) scientific communities specialize in the Holocene period (past 12 kyr), Quaternary epoch (past 2.6 Myrs) and Phanerozoic eon (past 540 Myrs), beyond the Precambrian. Each period, epoch, and eon exploits different paleoindicators, and each type of record has its own issues, problems, and limitations. The data types chosen here (details, Table 1, full analyses in the SUPP) are typical for the geological period covered, and the individual series were chosen for their quality and length, a total of 24. Despite the diversity, scaling analyses of their statistics (Figs. 2, 3, 4 and the figures in the SUPP) reveal that they can be placed in a common theoretical framework and, as summarized in Table 2, that they share basic statistical features over a wide range of scales.

Analysis methods

Multifractal processes

We may represent geochronologies by the points on the time axis that delimit successive measurements, or—as here—by the temporal density of such points. The former representation is of a geometric set of points, the latter of a function ρ(t). If the geochronology is scaling, then the points form a fractal set while ρ(t) is the (singular) density of a multifractal measure. The relationship between the two descriptions was clarified in the context of deterministic chaos’ strange (fractal) attractors^27,28^. This is discussed in the context of stratigraphy in ref. ^29^ as well as the reasons for preferring the multifractal description that we use throughout.

In scaling processes, fluctuations can be decomposed as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \rho \left(\Delta t\right)={\varphi }_{\lambda }{\Delta t}^{H};\lambda =\frac{{\tau }_{o}}{\Delta t};\Delta t\le {\tau }_{o}$$\end{document}

Where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varphi }_{\lambda }$$\end{document} is the underlying intermittent driving process and τ_ο_ is its outer scale. The statistics of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\varphi }_{\lambda }$$\end{document} may be defined by its qth order statistical moments: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle {\varphi }_{\lambda }^{q}\rangle={\lambda }^{K(q)}$$\end{document} with K(1) = 0 (corresponding to <φ_λ_> = 1), so that Δρ is related to φ_λ_ with the fluctuation exponent H. To characterize the fluctuation statistics at resolution Δt, we can determine the (generalized) qth order structure function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle {\Delta \rho \left(\Delta t\right)}^{q}\right\rangle$$\end{document} . Taking q^th^ powers and averaging the above equation, we obtain:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle {\Delta \rho \left(\Delta t\right)}^{q}\right\rangle =\left\langle {\varphi }_{\lambda }^{q}\right\rangle {\Delta t}^{{qH}}={\tau }_{o}^{K\left(q\right)}{\Delta t}^{\xi \left(q\right)};\xi \left(q\right)={qH}-K\left(q\right)$$\end{document}

so that the structure function can be characterized by its exponent ξ(q). Since K(1) = 0 we have ξ(1) = H.

The convex function K(q) characterizes the multifractal, intermittent part, and due to the existence of stable, attractive multifractal processes^32^ it is often assumed that K(q) is of the “universal” form:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K\left(q\right)=\frac{{C}_{1}}{\alpha -1}\left({q}^{{{\rm{\alpha }}}}-q\right)$$\end{document}

So that K(q) is determined by two fundamental parameters: the codimension of the mean C1 and the intermittency index 0 ≤ α ≤ 2. From this we can verify that K(1) = 0, K’(1) = C1, α = K”(1)/K’(1). A full characterization of the statistics of ρ(t) over a given scaling regime is then determined by the triplet H, C1, α.

The fluctuation exponent H

As outlined above, H characterizes the scaling of the mean (q = 1) with H > 0 indicating that fluctuations tend to grow with lag Δt so that the series tends to “wander”, and when H < 0, the fluctuations decrease with Δt, the fluctuations tend to cancel each other out.

The intermittency exponent C1

The role of C1 is less intuitive; for quasi-Gaussian processes, C1 = 0 so that K(q) = 0 and only H is important. In this special case and when 0 < H < 1, the process is a fractional Brownian motion (with H = 1/2, it is classical Brownian motion), and when –1 < H < 0, the process is fractional Gaussian noise with the special case H = –1/2 corresponding to standard Gaussian white noise. C1 > 0 characterizes the tendency for the series to make “jumps”, transitions: “intermittency”. For ρ(t) this implies that there are occasional high density “spikes” (see Fig. 1) that are significantly larger than the mean, and these spikes form a fractal set with codimension C1 (hence on the d = 1 time axis, a fractal dimension 1-C1), as well as occasional long epochs where ρ(t) is substantially below the average corresponding to long periods with few measurements. A more intuitive interpretation of C1 was proposed in ref. ^26^ to characterize the intermittency of the basic dynamical regimes: it characterizes the rate of divergence of the ratio of mean to root mean square (RMS) fluctuation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle \Delta \rho \left(\Delta t\right)\right\rangle /{\left\langle {\Delta \rho \left(\Delta t\right)}^{2}\right\rangle }^{1/2}={\left(\frac{\Delta t}{{\tau }_{o}}\right)}^{a{C}_{1}}$$\end{document}

where (from Eq. 3 above), we find \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=\frac{\left({2}^{\alpha }-2\right)}{\alpha -1}$$\end{document} and since 0 ≤ α ≤ 2, the constant satisfies ½<a < 1 and depends only (weakly) on α (in the chronologies, we found α ≈ 0 for Δt <  τ_c, α ≈ 1.3 (a ≈ 0.77) for Δt >  τc_, see Fig. S4.

The multifractal index α

The interpretation of α—variously called “the multifractal index”, “Levy index of the generator”, and “degree of multifractality”—is more subtle. When H = 0, α = 0, C1 > 0, the process is the turbulent on/off “beta model”, the process is non-zero over a fractal set with codimension C1. When α = 2 (the maximum), it is a “log-normal” multifractal with fluctuations following log-normal distributions (except for the extremes that have power law probabilities).

Haar fluctuations

Classical fluctuations are differences: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \rho {\left(\Delta t\right)}_{{dif}}=\rho \left(t\right)-\rho \left(t-\Delta t\right)$$\end{document} , and these can easily be estimated from irregular chronologies (e.g., pairwise differences). Another commonly used fluctuation is the “anomaly fluctuation” that is given by averaging segments over durations Δt of series that had their overall means ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\rho }$$\end{document} ) removed: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \rho {\left(\Delta t\right)}_{{anom}}=\overline{{\rho }_{\left[t,t-\Delta t\right]}}-\bar{\rho }$$\end{document} where the overbar indicates the average over the interval indicated in the subscript. These common fluctuations have major—but complementary— shortcomings: for physical processes that decorrelate in time, mean difference fluctuations cannot decrease with lag Δt, whereas, on the contrary, mean anomaly fluctuations cannot increase with Δt. Haar fluctuations (based on the Haar wavelet) overcome both these limitations by combining averaging with differencing. The Haar fluctuation of ρ(t) at resolution Δt is simply the difference between the average over the first half and the second half of the interval:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \rho \left(\Delta t\right)=\overline{{\rho }_{\left[t,t-\Delta t/2\right]}}-\overline{{\rho }_{\left[t-\Delta t/2,t-\Delta t\right]}}$$\end{document} . Like the difference fluctuation, it can easily be estimated from irregular chronologies (see methods and ref. ^14^, appendix B). For scaling processes where the average fluctuation varies as: <Δρ(Δt)> ≈ Δt^H^, differences are adequate when 0 < H < 1 and anomalies are adequate when –1 < H < 0. In comparison, Haar fluctuations are useful over the combined range –1 < H < 1 and this includes most Earth Sciences processes (for Gaussian processes, this corresponds to spectral exponents β in the range –1 < β < 3). Beyond its ease of application and its broad generality, Haar fluctuations are also relatively easy to interpret, being nearly equal to either the difference or anomaly fluctuations (depending on the sign of H).

Haar fluctuation analysis was carried out on each dataset in Table 1. Empirical series have finite resolutions. However, when the chronologies are nonuniform, two slightly different strategies may be used to estimate the fluctuations. The first is outlined in ref. ^14^ here, we used a slightly different method, necessary because we are interested in the correlations between the measurement density and the measured values. We therefore simply discretized the time axis at the finest resolution of the data and produced a uniform resolution series with values Ti or a 0 according to whether or not there was data at the (uniformly spaced) time interval. For the Haar fluctuation of the data values between t and t-Δt, we then took averages of the data values (i.e., excluding zeroes). The approximation being that where there were zeroes in the data, the values were close to the mean of the others in the segment duration Δt/2. For Haar fluctuations of the measurement density, we included the zeroes in the average. Any fluctuation involving segments with no data at all was rejected from the statistics. In numerous cases, we checked that this method gave very close estimates of the statistics as the previous method (we checked using various statistical moments as functions of scale). We also checked the process with various multifractal simulations of regular data with scaling “holes” with various statistics.

Finally, the fluctuations were multiplied by a “canonical” factor of 2 to make the Haar fluctuation closer to the difference fluctuation when H > 0 and closer to anomaly fluctuations when H < 0, see ref. ^8^, for the theory behind this. Note that when calculating the statistical moments, we averaged the absolute values of the fluctuations to the q^th^ power over all available disjoint intervals.

The intermittency: estimating C1, F(Δt)

The fluctuations in measurement density were normalized with respect to the mean measurement density over each dataset, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\rho }$$\end{document} , for comparison purposes. Using this simplified model, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi \left(1\right)=H$$\end{document} , therefore \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\left\langle \Delta \rho \left(\Delta t\right)\right\rangle }{\bar{\rho }}\sim \Delta {t}^{H}$$\end{document} . The slope of the graph of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \frac{\left\langle \Delta \rho \left(\Delta t\right)\right\rangle }{\bar{\rho }}$$\end{document} against \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \Delta t$$\end{document} yields an estimate for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H$$\end{document} (Fig. 2 top). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{1}$$\end{document} can be estimated from the function^6^ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\left(\Delta t\right)$$\end{document} , defined as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log F\left(\Delta t\right)=\log \left\langle \Delta \rho \right\rangle -\frac{\left\langle \Delta \rho \log \Delta \rho \right\rangle }{\left\langle \Delta \rho \right\rangle }$$\end{document}

In a scaling regime:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F\left(\Delta t\right)={\left(\frac{\Delta t}{\tau }\right)}^{{C}_{1}}$$\end{document}

So that plotting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \left(F\left(\Delta t\right)\right)$$\end{document} against \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \left(\Delta t\right)$$\end{document} yields an estimate of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{1}$$\end{document} as the slope of the graph (Fig. 2 middle). Combined with the relation for the mean to RMS ratio (eq. 4), we obtain:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a\log F\left(\Delta T\right)=\log \left[\left\langle \Delta \rho \left(\Delta t\right)\right\rangle /{\left\langle {\Delta \rho \left(\Delta t\right)}^{2}\right\rangle }^{\frac{1}{2}}\right]$$\end{document}

The plot of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \left(F\left(\Delta t\right)\right)$$\end{document} against \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \left(\Delta t\right)$$\end{document} therefore differs from the plot of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log [\langle \Delta \rho (\Delta t)\rangle /{\langle {\Delta \rho (\Delta t)}^{2}\rangle }^{\frac{1}{2}}]$$\end{document} against \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log \left(\Delta t\right)$$\end{document} by a constant factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a$$\end{document} . Taken in Fig. 2 middle to be approximately constant with value a = 0.83 corresponding to α =1.3 (close to the empirical value, see Fig. S4).

Correlations between fluctuations

If the dynamical regimes governing the T(t) and ρ(t) processes are the same, then it is possible that the ΔT(t) and Δρ(t) fluctuations are strongly correlated. At each lag Δt, we therefore estimate the correlation coefficient R(Δt). This was done using:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R\left(\Delta t\right)=\frac{\left\langle \Delta \rho \left(\Delta t\right)\Delta T\left(\Delta t\right)\right\rangle }{{\left\langle \Delta {\rho }^{2}\left(\Delta t\right)\right\rangle }^{\frac{1}{2}}{\left\langle \Delta {T}^{2}\left(\Delta t\right)\right\rangle }^{\frac{1}{2}}}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R$$\end{document} is the correlation coefficient at a given timescale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho$$\end{document} is the measurement density, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T$$\end{document} is the measured quantity. This was done to explore potential systematic biases in the datasets. Note that the fluctuations are independent of any additive constant. However, strong trends could lead to nonzero mean fluctuations if present; these should be removed in the preprocessing (e.g., overall linear or quadratic trends, see Fig. S30).

If ρ(t) is statistically independent of T(t) (here, paleotemperature or dust flux), then high ρ(t) is not associated with particularly large or particularly small T(t) fluctuations, so that chronology corrections are easier. In contrast, strong correlations (R(Δt) ≈ 1) imply that Δρ(Δt) and ΔT(Δt) vary together. Intervals with large Δρ(Δt) have frequent measurements, so that R(Δt) ≈ 1 implies that this tends to occur when fluctuations in ΔT(Δt) are also large. Conversely, there are relatively few measurements (very negative Δρ(Δt)) when the fluctuations also have particularly low values (negative ΔT(Δt)). When R is close to 1, in intervals Δt where ΔT is large, T(t) is sampled very often, and conversely, in intervals where ΔT has low values, T(t) is rarely sampled. Reversing the sign of T(t) changes the signs of the fluctuations (it swaps “low” and “high” values of ΔT (Δt)), and it changes the sign of R(Δt), so that strong negative correlations R(Δt) ≈ –1 are analogous to R(Δt) ≈ 1 but for low values of the indicator.

Time interval probability distributions and the length-Sadler effect

While typical inter-measurement intervals depend on sedimentation processes and rates, individual cores typically also have “gaps” where by—as a consequence of erosion—whole swaths of sediment have completely disappeared from the record. Irrespective of whether or not a time interval is a consequence of a sedimentary or erosion process, at time t, the typical interval between successive measurements will be 1/ρ(t). Therefore, when ρ(t) is small enough, ρ(t) can no longer be treated as a continuous stochastic process. Mathematically, it suffices to instead treat the measurement process as a compound Poisson process controlled by ρ(t), which is then a multifractal subordinator^29^. This model provides realistic simulations of the Geological Time Scale—i.e., the spacings of the “golden spikes” that define the main geological epochs, eras etc^29^. Full models of scaling sedimentation that follow the observed statistics are still needed.

These modeling issues are important, but are outside our present scope, which is the empirical statistical characterization of the measurement process. For this, it is enough to determine the probability of a random interval τ’ exceeding a fixed τ: Pr(τ’ >  τ) ( = 1–the (usual) Cumulative Distribution Function, CDF). Although at small τ, the distribution is (often) roughly Gaussian (Fig. 4, dashed black lines, figures in the SUPP). However, at large τ, there is a transition to a scaling distribution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Pr \left(\tau \mbox{'} > \tau \right)\approx {\tau }^{-{q}_{D}}$$\end{document} , (“power law tails”), indicating a hierarchy of extreme intervals - gaps. Small values of qD imply the existence of huge gaps; for example, a value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${q}_{D}\le 1$$\end{document} (e.g., the NEEM record) implies that the duration of the longest gap in a chronology is of the same order as the sum of all the smaller gaps.

Supplementary information

Transparent Peer Review file Time scales and gaps, Haar Fluctuations and multifractal Geochronologies Supplementary Material