Utilizing Maximum Entropy Spectral Analysis (MESA) to identify Milankovitch cycles in Lower Member of Miocene Zhujiang Formation in north slope of Baiyun Sag, Pearl River Mouth Basin, South China Sea

1 Introduction

Quasi-periodicity or periodicity is the common phenomenon for almost every object in the universe. Derived from the astronomical cyclical movement, the deposition processes and sediment sequences of the Earth are within the cyclicity and periodicity as well. The quasi-periodic variations of the precession, obliquity and eccentricity of the Earth orbital parameters generate different features and rhythmicity records in the sedimentary rocks, called Milankovitch cycles in geology [1]. However, the outcrops in fields and cores in petroleum wells just provide a small part of the formation. It is impossible to directly view and observe the cycles in the deep burial sedimentary rocks continuously. Fortunately, logs in the boreholes can continuously collect the features of the underground sediments, which offers other indirect materials to research the Milankovitch cycles [2, 3, 4]. These digital signals may consist by various stacked cyclical indexes, and one of the challenges is how to analyze the sedimentary records and to find the information that indicate the cyclicity in the stratigraphy and cyclostratigraphy.

Recent advances in signal analysis provide many approaches in different theories for this problem, such as Multitaper method (MTM) [5, 6, 7], Average Spectral Misfit (ASM) [8, 9], first-order autoregressive (AR1) [10], Evolutive Harmonic Analysis (EHA) [11], Time Scale Optimization (TimeOpt) [12, 13], Wavelet Transform (WT) [14, 15, 16, 17, 18], Fast Fourier Transform (FFT) [2, 4, 19] and so on. However, MTM, ASM, AR1, EHA and TimeOpt prefer discrete astronomical and climatological experimental data sampled in regularly-space [5, 6, 7, 8, 9, 10, 11, 12, 13], and WT and FFT, spectral analysis methods, are generally applied to continuous date in boreholes [14, 15, 16, 17, 18]. The primary reason for using spectral methods is the understanding that cyclical climatic changes are recorded in strata, which can therefore be expected to contain information in the frequency domain [19, 20, 21]. If a wireline log is viewed as a series of data points and a regularly-spaced sampling of a complex waveform, then spectral analysis methods are required to fully understand it.

However, WT methods have advantages in the original data with very spiky-large spikes which would dominate the frequency spectra to the point and lead to the obscure of the interesting frequencies. Although FFT may be the most familiar method in other applications, it attempts to decompose the data into exact sine waves, which are very rarely in the geological data. Maximum Entropy Spectral Analysis (MESA), another widely used spectral analysis approach, is much more tolerant of any imperfections in the expression of an underlying periodicity, such as Milankovitch cyclicity, by generating a model of the data rather than an exact decomposition mathematically [22]. Based on the MESA, a method of identifying the Milankovitch cycles and estimating the cyclical sedimentation rate in the Lower Member of Miocene Zhujiang Formation in north slope of Baiyun Sag, Pearl River Mouth Basin, is introduced in this paper.

2 Milankovitch cycles

The Earth’s orbit around the Sun is influenced by gravitational attractions of the Moon and the other planets in the solar system, which produces the quasi-periodic variations, termed Milankovitch theory, in the orbital parameters: precession, obliquity and eccentricity (Figure 1) [23, 24]. At present, the precession of equinoxes refers to the wobble of the Earth on its axis with major periods of ~23.00 and ~19.00 kyr. The obliquity, varies between 22.1^∘ and 24.5^∘, is oscillation of the Earth’s axial tilt with major periods of ~41.00 kyr and secondary period of ~29.00 and ~54.00 kyr. As the elliptical orbit, eccentricity varies between 0.0005 and 0.0607 with long period of ~413.00 kyr and short period of ~95.00 and ~123.00 kyr (average in 100.00 kyr) [23, 24]. At different geological age, these orbital parameters would be slightly changed (Table 1) [25].

Figure 1

Sketch of the Earth’s orbit. The precession, obliquity and eccentricity are three important orbital parameters [24].

The quasi-periodic variations of the orbital parameters directly control the all of the environment changes on the Earth [26]. Then, the environment influences all intimately linked factors eventually that directly or indirectly control sediment production, transportation and accumulation. Depending on the sensitivity of the sedimentary systems, environment forcing may exist in the stratigraphic records in continuous depositional settings where the preservation potential is high. It is clear that the strata thus indirectly reflect the orbital parameters. Although the orbital changes are not truly cyclical but quasi-periodic, it is commonly accepted by the geological community to call as “Milankovitch cycles” in the stratigraphy and cyclostratigraphy [25].

3 Geological setting

Pearl River Mouth Basin, in the continental shelf northern South China Sea, is a representative rift basin of passive continental margin. Strike in NE-SW, the basin is approximately paralleled with the shoreline of south China continent [27, 28, 29]. It is a vital offshore petroliferous basin of China. In tectonic, Pearl River Mouth Basin is divided into five major structural units from north to south: North Faultterrace zone, North Depression zone, Central Uplift zone, South Depression zone and South Fault-terrace zone (Figure 2a) [27, 28, 29]. Each unit include several sags and uplifts furtherly. As one of the significant oil-rich sags, Baiyun Sag is in North Depression zone, encircled by Panyu Low Uplift, Yunkai Low Bulge and South Fault-terrace zone (Figure 2b).

Figure 2

Tectonic features of Pearl River Mouth Basin. (a) Five major structural units in Pearl River Mouth Basin. The Baiyun Sag is in the South Depression zone. (b) Fourteen wells in research region.

Under adequate sediment supplying, steady formations are developed in Baiyun Sag from Eocene to Quaternary. From 32 myr, for the continuous expansion of the South China Sea, the depositional environment changed from lakes to shoreland in Eocene and Oligocene [30]. Then the Pearl River Mouth Basin turned into marine sedimentary stage, when shelf slope break zone maintained at the south of Panyu Low Uplift and north slope of Baiyun Sag (Figure 3a) [30]. Therefore, between 23.8 myr and 18.5 myr, shelf marginal delta was primarily developed to form Lower Member of Zhujiang Formation in Lower Miocene [30]. Constituted by various sands, silts and clays, the thickness of Lower Member of Zhujiang Formation is ranged from 150 m to 400 m (Figure 3b) [30, 31].

Figure 3

Formations in Baiyun Sag. (a) Formations from Eocene to Quaternary. The Zhujiang Formation is mainly formed by shelf marginal delta. (b) Lower Member of Zhujiang Formation in north slope of Baiyun Sag, Well BY9. Delta front is the dominated depositional environment.

4 Database and methodologies

4.2 Maximum Entropy Spectrum Analysis

In information theory, entropy is defined as a measure of signal uncertainty. Proposed by Burg [33], MESA extracts the periods of different frequency from the time series, which generates the minimum of the power spectrum error and increases the resolution. It is equivalent with the autoregressive model [34]. The prediction error of logging record x_i (i=1, 2, ..., N) could be regarded as

ej=xj−∑k=1mαkxj−k,j=m+1,m+2,…,N.

Taken 1, −α₁, . . . , −α_m as filtering factors and e_j as the output of the digital filter, a prediction error filter is just designed. In frequency domain, the filter is expressed by

Eω=Aω⋅Xω

E (ω) is the frequency spectrum of prediction error e_j, A(ω) is the frequency spectrum of prediction error filter and X (ω) is the frequency spectrum of the logging x_i. The frequency spectrum of the x_i is thus

Sω=Xω2=Eω2Aω2=2σα21−∑k=1mαke−i2πωk2

σα2is the power of the prediction error filter and S(ω) is the maximum entropy frequency spectrum of the stochasticprocess. It is outside the key points of this paper to derive all detail procedures that illustrate such spectral theory.

For example, a simulate log has been processed by MESA (Figure 4). The simulate log is synthetized by five logs with constant wavelengths of 41.3 m, 10 m, 4.1 m, 2.3 m and 1.9 m respectively. In the MESA spectral band, frequency spectrum generated by MESA, there are five peaks whose wavelengths are 41.3 m, 10 m, 4.1 m, 2.3 m and 1.9 m respectively. It suggests that MESA is effective to unmix the cycle property in a composite log.

Figure 4

Spectral band of simulate log. From second to sixth pane show synthetic logs with wavelengths of 41.3 m, 10 m, 4.1 m, 2.3 m and 1.9 m respectively. The seventh pane shows the simulate log, sum of the synthetic logs. The eighth pane shows the two-dimension frequency spectral band calculated by MESA (window=30 m). High spectral amplitudes correspond to the wavelengths of 41.3 m, 10 m, 4.1 m, 2.3 m and 1.9 m respectively.

5 Results

Without exposure and denudation, the sustaining marine sedimentary successions can often preserve more complete record of orbital cycles than lacustrine strata. From 23.8 myr, the north slope of Baiyun Sag was inundated by marine water and the slope break belt is appeared. Until 18.5 myr, Lower Member of Miocene Zhujiang Formation is developed in delta front facies. Between both times, there is another chronologic age recorded by seismic reflection, 21 myr, which divided Lower Member of Zhujiang Formation into two parts. Taking Well BY6 for example, each part comprises lower delta outer front and upper delta inner front (Figure 5). Thick sands would be deposited in delta inner front, and interbed of thick mudstone and thin sands may be developed in delta outer front.

Figure 5

MESA result of Lower Member of Zhujiang Formation in Well BY6. The sandstones are in low GR values and siltstones and mudstones are in high GR values. The two-dimension frequency spectral band is calculated by MESA (window=30 m). Six typical one-dimension frequency spectrum at different depth match the Milankovitch periods.

In order to identify probable Milankovitch cycles in this borehole, GR log is processed utilizing MESA to generate two-dimension spectral band (Figure 5). The GR are high values in sandstone intervals and low values in siltstone and mudstone intervals. In the two-dimension spectral band, the ratio of wavelengths at each depth point should be check carefully to judge whether match the Milankovitch periods approximating wavelengths ratio in 41.3: 10: 4.1: 2.3: 1.9. Finally, six typical one-dimension frequency spectrum at different depth point, who’s wavelengths ratio of amplitude peak matches the Milankovitch periods well (Figure 5, Table 2). The maximum amplitude suggests the dominating Milankovitch period. Between the 23.8 myr and 21 myr, the dominating period is mainly ~100.00 kyr, while between the 21 myr and 18.5 myr, the dominating period is mainly ~413.00 kyr.

The abscissa axis of the amplitude peak is the wavelength, implying the sedimentation thickness during the dominating period. Taking the spectral amplitude at depth of 3432.24 m for example, five amplitude peaks are found distinctly. Their wavelengths are 40.14 m, 10.4 m, 4.03 m, 2.27 m, 1.89 m respectively, closing to the ratio of 41.3: 10: 4.1: 2.3: 1.9. What’s more, the amplitude in wavelength of 40.14 m is the maximum value, which means that the major long eccentricity period ~413.00 kyr may dominate the cycles. Finally, in this interval, the sedimentation rate controlled by cycle is estimated about 40.14/413=0.0972 m/kyr., Therefore, the cyclical sedimentation rates could be estimated at these six depth points (Table 2). Vertically, the cyclical sedimentation rates of thick sands in delta inner front are slightly greater than that of the delta outer front.

6 Discussions

The north slope of Baiyun Sag is a gently dipping monoclinal slope whose tilts towards southeast with abundant sandy sediments. Profit from three-dimension seismic data prospecting, three sedimentary boundaries are interpreted as 23.8 myr, 21 myr and 18.5 myr in the Lower Member of Miocene Zhujiang Formation. From depositional zone to bathyal region, a correlation wells profile involving six wells is built to discussed the Milankovitch cycles in this formation (Figure 6). The linear distance between Well BY1 and Well BY14 is about 95 km. The depositional microfacies are extracted from conventional electrofacies [30, 31]. From Well BY1 to Well BY14, depositional environment changed from delta inner front to delta outer front, prodelta and bathyal gradually. The GR of five wells are processed by MESA, except Well BY14 for its partially interval is eruptive deposit. Using the method proposed in this paper, there are obvious seven, eight, six, eight and thirteen typical Milankovitch cycles found respectively in Well BY1, BY3, BY8, BY10 and BY12 (Table 2). The results suggest that: 1) the dominating cycles frequency in delta inner and outer front is mainly ~413.00 kyr and ~100.00 kyr, and the concluded cyclical sedimentation rate is about 0.1 m/kyr; 2) in the mudstone of prodelta and bathyal, dominating cycles frequency is occasionally ~41.00 kyr and ~23.00 kyr, and the implied cyclical sedimentation rate is about 0.085 m/kyr.

Figure 6

MESA result and microfacies of a correlation wells profile whose location can be found in Figure 2b. No MESA in Well BY14 for the eruptive deposit.

In the stratigraphy, formation developed in the delta inner and outer front are thin than that in prodelta and bathyal, the average sedimentation rate of the former is less than that of the latter. For example, the average sedimentation rate of Lower Member of Zhujiang Formation in Well BY1 and Well BY12 is 0.0321 m/kyr and 0.0661 m/kyr respectively. However, the cyclical sedimentation rate estimated from Milankovitch cycles imply the inverse consequence, which is more accord with principles of sedimentology.

7 Conclusions

A method of Milankovitch cycles identification and sedimentation rate estimation is proposed based on Maximum Entropy Spectral Analysis using GR log in the petroleum boreholes, which could be divided into four stages. Choosing the calculating windows is the first stage. Calculating the frequency and generates the two-dimension frequency spectrum in short data sections is the second stage. Finding the potential Milankovitch cycles in the one-dimension frequency spectrum at each depth point serves as the third stage. Estimating the sedimentation rate of formation in the borehole sites is as the fourth stage. In one-dimension frequency spectrum of Milankovitch cycles, the wavelength of average amplitude spectrum peak would be matched to the ratio of the Milankovitch period. The Milankovitch period of the maximum amplitude spectrum is regarded as the dominant periodicity that controls the cyclicity, and the ratio of the wavelength and period offers an approximated cyclical accumulation rate of the depth.

The proposed approach has been applied to the Lower Member of Miocene Zhujiang Formation in north slope of Baiyun Sag, Pearl River Mouth Basin. The dominating cycles frequency in delta inner and outer front is mainly ~413.00 kyr and ~100.00 kyr, and the concluded cyclical sedimentation rate is about 0.1 m/kyr. In the mudstone of prodelta and bathyal, dominating cycles frequency is occasionally ~41.00 kyr and ~23.00 kyr, and the implied cyclical sedimentation rate is about 0.085 m/kyr. This method may be feasible and effective to identify Milankovitch cycles and estimate the cyclical sedimentation rate in the marine depositional environment and be useful to other geological researches.