Multifractal characterization of the Coniacian–Santonian OAE3 in lacustrine and marine deposits based on spectral gamma ray logs

Limited to the Atlantic and its surrounding basins, the expression of the Coniacian–Santonian oceanic anoxic event (OAE3) was discovered in the non-marine Cretaceous Songliao Basin, Eastern Asia not long ago. In this study, based on spectral gamma ray logs data recorded in three basins, the self-similarity of the OAE3 was studied through the analysis of the scaling properties of thorium–potassium and thorium–uranium distributions both in marine and terrestrial environments using the multifractal detrending fluctuation analysis. The results indicate that, in both marine and terrestrial systems, the OAE3 intervals are characterized by their multifractal nature due to long-range correlation. However, the multifractal features of the studied OAE3 intervals are different in the three basins, although some common trends were observed. By comparing the degree of multifractality of the OAE3 deposits with the clay minerals and the redox conditions, it appears that the changes of the multifractal features are controlled by local changes such as clay mineralogy and redox conditions in both milieus under different sedimentation patterns. At all sites, the left side shortened spectrum of the thorium–potassium distribution suggests the presence of local fluctuations with minor amplitudes during the OAE3. Furthermore, the shortened singularity spectrum of the thorium–uranium distribution reflects the existence of small-scale fluctuations with large amplitudes at marine sites while in the non-marine Songliao Basin, the thorium–uranium distribution suggests the presence of local fluctuations with small amplitudes during the OAE3. Therefore, a more local behavior of the event is considered although the regional character is not neglected.

Scientific RepoRtS | (2020) 10:14363 | https://doi.org/10.1038/s41598-020-71327-w www.nature.com/scientificreports/ is characterized by significant burial of organic-rich shale, its exact paleo-ecological extent, spatial distribution and development are not well established 3,[17][18][19][20] . However, most of these previous studies are one-sided i.e., firstly, they are focused on marine environments although the influence of OAEs in lacustrine environments has been described 4 . Secondly, a complete terrestrial-marine correlation of the last Upper Cretaceous OAE is not documented in the literature, mostly due to the fact that the terrestrial expression of the OAE3 (TEOAE3) is scarce. For instance, focusing on the equatorial Atlantic, Ocean Drilling Program (ODP) sites alongside the Deep Ivorian Basin (DIB) and the Demerara Rise (DR) precisely ODP Leg 159 21,22 and ODP Leg 207 19,20 , which record this event. Although comparatively rare in lacustrine environments, the recent finding of the TEOAE3 in the upper Qingshankou in the non-marine Songliao Basin, northeastern China, suggests that paleo-lakes are proven suitable for testing hypotheses about the OAE3 triggering mechanisms in terrestrial environment 23 , thus opening a door for a possible terrestrial-marine OAE3 correlation. This is furthermore motivated by the work of Chamberlain et al. 24 who found a synchronous response of OAE2 in the lower Qingshankou Formation in the Songliao Basin. So far, the Songliao Basin represents one of the nearly complete Cretaceous terrestrial records in the world 25 , making this area an optimal site to infer global, local or zonal paleoenvironmental changes. Therefore, the question of whether there are similarities in the development of OAE3 in marine and terrestrial environments arises. Thus, carrying out a joint study of the Coniacian-Santonian oceanic anoxic event in both marine and lacustrine environments in different basins can help present a better understanding of the controlling mechanisms and its regional/local environmental responses in both marine and terrestrial systems. Moreover, previous studies mostly focused on organic proxies (biomarkers) and inorganic compounds including stable isotope compositions 26,27 . However, studies of the self-similarity of OAE3 based on petrophysical data have not been conducted.  56 ). Location of sites 959D (Deep Ivorian Basin) and 1261B (Demerara Rise) are from Wagreich 56 ; location of SK1S (Songliao Basin) is from Wang et al. 26 .
In the geosciences, fractal and multifractal features have been increasingly used to describe geological events, geological processes and geological objects 28,29 . The principal advantage of fractal/multifractal theory concerns its capacity to characterize irregular and complex phenomena or processes that display similarity over a broad sequence of scales (self-similarity) 30 that the traditional Euclidean geometry method (integral dimension) fails to analyze 31 . The fractal concept originates from geometric self-similarity which was first proposed by Mandelbrot 30 . Self-similar fractals evaluate the presence of an arrangement that is similar to itself on any scale. Since, the concept has been widely applied in medicine, geology, geophysics 28,31,32 . Many occurring natural processes like OAEs are complex phenomena with their development controlled by complicated processes, and thus do not follow normal distributions. Consequently, ordinary linear dataset analysis techniques are not sufficient for a comprehensive analysis of these events. The fractal method is one of the useful concepts to be applied to analyze complex systems [33][34][35][36] . Among numerous fractal methodologies developed since Mandelbrot's findings, the concept based on the detrended fluctuation analysis (DFA) which was devised by Peng 37 is the most frequently utilized. By expanding the concept of DFA, Kantelhardt et al. 38 initiated the notion of multifractal detrended fluctuation analysis (MFDFA), which works similarly to the DFA. MFDFA is useful to reveal the multifractal behavior in any non-stationary data series. This method also allows to detect the causes of multifractality by comparing the estimated shuffled and surrogate data series derived from the original data series to the original data series [39][40][41][42] . Multifractality helps to describe more carefully and comprehensively the dynamic properties of systems, and characterize their behaviors both locally and globally 43 . MFDFA has been used by some scholars to study geophysical well log data 44,45 , which are governed by complex spatio-temporal dynamics of which nonlinearity and scaling are the dominant processes 46 .
Using core samples from boreholes, previous studies have greatly improved our understanding of major Cretaceous marine and terrestrial changes 2,4-6,47 . However, integrated studies of marine and terrestrial Cretaceous deposits are rare 26 , and the multifractal correlation based on spectral gamma ray by means of well log data have not been conducted.
Compared to the core samples analysis, well logging has several advantages 48,49 . First, for executing highprecision continuous sampling (spectral gamma ray, nuclear magnetic resonance, etc.), and measurement of in-situ formation conditions. Second, for investigating a volume of sediment that is often greater than the one represented by a core or plug, and consequently than a cutting and so more representative of the mean properties of the rock, especially in heterogeneous rocks. Third, the human factors have little effects on the measurement processes. Among the well log data, the spectral gamma ray has proved particularly useful for understanding paleoenvironmental changes [50][51][52][53] .
To understand the multifractal behavior of the Coniacian-Santonian OAE, marine sedimentary records must be integrated with terrestrial deposits. Based on data from scientific drilling projects, this work can be achieved by correlating the non-marine Cretaceous Songliao Basin in eastern Asia with the marine Deep Ivorian Basin in western Africa and Demerara Rise off the coast of South America where a well preserved response of the OAE3 has been described 10,17,20,54-57 . This paper is the first work integrating the multifractal nature of a marine and terrestrial OAE. The aim is to highlight and compare the multifractal properties of a rare and complex geological event such as OAE3 recorded in marine and terrestrial systems as well as analyze the causes of the multifractality and the possible factors affecting the multifractal features in both environments, which can help to understand the regional or local nature of the last Late Cretaceous OAE.

Geological Setting of the drilling sites
Songliao Basin (SLB). The SLB is an elongated large non-marine Cretaceous basin located in northeastern China ( Fig. 1). It represents a Mesozoic-Cenozoic intracratonic basin with the thickness of Cretaceous strata reaching up to 6,000 m 58,59 , making the SLB one of the largest and longest-lasting Cretaceous non-marine basins on the Earth 25 . A fluvio-deltaic and lacustrine sedimentary succession ranging from the Turonian to Campanian has been recorded through two scientific drilling sites (SK1S and SK1N). The Late Cretaceous formations at site SK1S consists of the older to the younger strata, of six formations namely the Quantou ( K 2 q ), Qingshankou ( K 2 qn ), Yaojia ( K 2 y ), Nenjiang ( K 2 n ), Sifangtai ( K 2 s ), and Mingshui ( K 2 m ) formations controlled by local tectonic motions and climatic variations 60,61 . K 2 qn is subdivided into lower K 2 qn (Member 1) composed mainly of organic-rich laminated mudstone and upper K 2 qn (Members 2 and 3) essentially formed of undifferentiated grey shales 62,63 , with presences of thin layers of colored marl in the lower section 58,59 . Members 2 and 3 of K 2 qn in hole SK1S are semi-deep lacustrine sediments, and mainly late Turonian to late Coniacian in age 25 (Fig. 2). The basin experienced a general shallowing trend from deep lacustrine deposits in Member 1 of K 2 qn 62 , to fairly deep and shallow lake in Members 2 and 3 of K 2 qn 63 .
Deep Ivorian Basin (DIB). The DIB belongs to that of Togo-Benin-Nigeria, to a vast sedimentary complex whose subsidence began in the Cretaceous. It is a quasi-covered basin, and situated off equatorial West Africa (Fig. 1); whose formation is related to the expansion of the opening Equatorial Atlantic Gateway 57 (Fig. 2). Deposition of black shales started during the Turonian when the margin differentiation commenced and continuous subsidence generated a semi-enclosed sub-basin 57 , and was controlled by the climate changes in equatorial Africa.  68 . These sediments are mainly constituted of calcareous claystone with organic matter combined with clayey chalk with nannofossils and clayey limestones (Fig. 2). DR experienced a consistent deepening trend during the accumulation of the organic-rich sediments from Cenomanian to early Campanian 69 .
Results. The multifractal characteristics of the Th/U and Th/K ratios of the OAE3 and non-OAE sections at SLB, DIB and DR were studied by the MFDFA. Prior to the MFDFA application, the datasets were first segmented and integrated random walk time series into non-overlapping segments with different scaling (except DIB before the OAE3 due to the logging depth). By using Eq. (5) and for − 5 < q < 5, the fluctuation function (Fq(s) versus s) of the Th-K and Th-U distributions have been computed. The scaling behavior obtained by the log-log plots of the fluctuation function of the Th/U and Th/K ratios of the OAE3 and non-OAE units in different boreholes at SLB (SK1S), DIB (959D) and DR (1261B) are shown in Figs. 3 and 4, respectively. All the generalized Hurts exponents (Hq) of the Th-K and Th-U distributions in each OAE3 and non-OAE segments are q dependent (decrease with the increasing of q). The slope of all the fluctuation functions at different scale decrease with the increase of q, implying that the Th-U and Th-K distributions of the studied OAE3 and non-OAE intervals have a multifractal property 34,70,71 . By analyzing the slopes, it appears that the scaling properties of the Th-K and Th-U distributions at SLB, DIB and DR are closer in the OAE3 interval than the non-OAE3 intervals (Figs. 3, 4; Table 1). Moreover, the scaling behavior of the Th/U ratio in SK1S and 959D are nearly uniform for all q values which may infer some similarities in the thorium and uranium distributions in the OAE3 units in both boreholes (Fig. 3). Also considering the Th/K ratio, some similarities in the scaling properties of the Th/K ratio recorded in the studied OAE3 intervals in SK1S and 1261B can be noticed (Fig. 3).    Table 2 shows that on average, the highest slope differences of the mass exponent τ(q) between the negative (q < 0) and positive (q > 0) q-order of Th/U and Th/K ratios are found in the OAE3 interval, indicating that the Th-U and Th-K distributions in this depositional interval possess a higher degree of multifractality in their scaling properties. In addition, the slope difference of Th/U and Th/K ratios in SK1S and 959D are close within this unit, which may infer some similarities in their degree of multifractality. The singularity spectrum D(α) and the singularity exponent α describing the multifractal spectrum of the Th/U and Th/K ratios of the different formations in the OAE3 and non-OAE3 units are displayed in Figs. 5c,f and 6,e,j, respectively. The curves described by all the multifractal spectrums are asymmetric unimodal with different singularity spectrum widths, implying that under fluctuating environmental settings, the Th-K and Th-U distributions in the studied OAE3 and non-OAE3 intervals present varying multifractal characteristics. Furthermore, the singularity spectrum widths defined by Δα = α max − α min , which reflect the degree of multifractality/complexity were calculated using Eq. (11) ( Table 3). When Δα is low, the whole system is described to possess low heterogeneity in its local scaling and vice versa 73 . Thus, the results show that on average, the OAE3 interval has the larger multifractal spectrum widths of both Th/U and Th/K ratios, which may indicate that the temporal fluctuation of both Th-U and Th-K distributions in this unit display a high degree of multifractality/ complexity, with a more complex trend in SK1S for Th-U distribution and a more complex trend in 959D for Th-K distribution.

Discussion
Sources of multifractality. In general, the multifractality observed in time or spatial series have two main causes 38,42,43,73 : (1) the multifractality originating from a wideness of the probability density function (PDF), and (2) the multifractality caused by long-range correlation for small and large fluctuations in the pattern arrangement. To highlight the dominant one between these two types of multifractality in the Th-K and Th-U distributions in the studied OAE3 and non-OAE3 intervals, we analyzed the corresponding randomly shuffled and surrogate datasets generated using the original Th-K and Th-U datasets based on Kimiagar et al. 41 , Movahed et al. 39 and Niu et al. 40 . Indeed, the shuffled process removes the multifractality due to long-range correlation while retaining the multifractality caused by the broad PDF signal. Thus, the shuffled dataset will exhibit nonmultifractal behavior, with h(2) = 0.5 if the multifractality properties are solely controlled by the long-range correlation. In case the multifractality of the analyzed Th-K and Th-U datasets originates from both types of multifractality, the shuffled series will show weaker multifractality than the original series 39 . By considering the surrogate method, the generalized Hurst exponents h(q) generated have the capability to be q-independent if the multifractality observed in the Th-K and Th-U distributions is mainly from the wideness of PDF. However, if the multifractality properties are triggered by both sources, then weaker multifractality would be found in both shuffled and surrogate data of the Th-K and Th-U datasets 39,42 .
Hence, to explore the origins of multifractality in boreholes SK1S, 959D and 1261B, firstly, the shuffled and surrogate procedure were applied to the original Th-K and Th-U datasets of the OAE3 interval, and the results are shown in Fig. 7. Afterwards, the same procedure was applied to the non-OAE3 units (Table 3).
From Fig. 7, it is observed that the generalized Hurst exponents h(q) of the shuffled data set are slightly q dependent for all Th/U and Th/K ratios with monotonically decreasing curves. Also, the analysis of the generalized Hurst exponents indicates that for q = 2, h(q) slightly converge to 0.5. This trend shows that the multifractal Table 1. Slope difference of the fluctuation functions of the thorium-uranium and thorium-potassium distributions between localities before, within and after the OAE3 intervals. SK1S (Songliao Basin), 959D (Deep Ivorian Basin), and 1261B (Demerara Rise).
Scientific RepoRtS | (2020) 10:14363 | https://doi.org/10.1038/s41598-020-71327-w www.nature.com/scientificreports/ behavior described in the Th/U and Th/K ratios is due to the broad PDF and the long-range correlation. However, by analyzing the surrogate curves, we observe that the generalized Hurst exponents h(q) of the surrogate series are all quasi q-dependent, suggesting that the multifractality due to long-range correlation is dominant in the Th-K and Th-U distributions in the OAE3 intervals. Tanna and Pathak 74 , and Wu et al. 42 described a similar trend where both the broad PDF and the long-range correlation were the cause of the multifractality when studying an ionospheric scintillation time series and hydrological data, respectively. To identify the influence of the broadness of the PDF and the long-range correlation in their dataset, the dimension of the multifractal spectrum width of the original time series was correlated to the estimated dimensions of the multifractal spectrum width of the shuffled and surrogate time series. Their results indicated a narrow multifractal spectrum width of the shuffled time series compared to that of the surrogate time series, which also had a narrow multifractal spectrum width compared to that of the initial dataset. Subsequently, they concluded that the contribution of the broad PDF on the multifractal behavior observed in their dataset was weaker than the long-range correlation. www.nature.com/scientificreports/ Following the same procedure, the multifractal spectrum widths of the original, shuffled and surrogate series of the Th-U and Th-K distributions were examined for the OAE3 interval as well as the non-OAE3 units (Table 3). Through Table 3, it is obvious that the multifractality triggered by the long-range correlation is dominant in the OAE3 intervals in all the boreholes, while it is not the case in the non-OAE3 units. In SK1S, the �α Surrogate of Th-U distribution is stronger than �α Original before and after the OAE3 interval, indicating that the PDF dominates the multifractality 75 . Th-K distribution exhibits similar behavior before the OAE3 interval at www.nature.com/scientificreports/ site 1261B. Moreover the multifractal behavior of the Th-U and Th-K distributions after the OAE3 intervals at sites 1261B and 959D, respectively is dominated by the PDF since �α Shuffled is stronger than �α Surrogate 42 . Thus, the main cause of the multifractal behavior of the Th-U and Th-K distributions in the non-OAE3 intervals is non-homogeneous, however, it is homogeneous in the OAE3 interval (Table 3). Consequently, is there a common factor controlling the multifractality in the OAE3 intervals in the three basins?
Possible factors affecting the multifractal properties of Th-U and Th-K distributions in the OAE3 interval. The Th/U and Th/K ratios in sediments fluctuate under a number of different controls such as clay mineral content, paleo-redox conditions and paleoclimate 53,76 . Th/K ratio is for the most part a function of clay mineral content in shale formations 77 . Previous studies have shown that the presence of shale and variations in the sub-surface sedimentation pattern largely influence the multifractal behavior of the gamma ray log response 44,45 .
Based on the thorium and potassium distributions, we analyzed the likely relationship between clay minerals and studied OAE3 intervals multifractality (complexity) in each borehole through a Th-K cross-plot (Fig. 8g) and singularity spectrum width correlation (Table 2). However, since thorium and potassium can also be associated with non-clay minerals, we first evaluated the correlation between Th and K, and the estimated weight percent of clay (WT% clay) based on the empirical expression developed by Bhuyan and Passey 78 (Fig. 8a-f). WT% clay represents the total clay obtained from the total GR. Meanwhile, the dispersions between the total GR and Th, K and U elements are useful to investigate the likely sources of these elements 79 . Therefore following Ito et al. 79 , a strong correlation between WT% clay, and Th and K suggests that clay is the main contributor of Th and K responses, while weak correlation suggests that the main source of Th and K may not be the clay minerals. From Fig. 8a-f, a good positive correlation appears between the weight percent of clay, and Th and K variations at sites 959D and 1261B, while in SK1S, the correlation is moderate for K. The analysis of the correlation coefficients supposes that clay minerals are the predominant contributors of Th and K responses, and therefore led to clay Table 2. Slope values of mass exponent τ(q) for thorium-uranium and thorium-potassium distributions before, within and the OAE3 intervals at Songliao Basin (SK1S), Deep Ivorian Basin (959D) and Demerara Rise (1261B).  www.nature.com/scientificreports/ type evaluation by a Th-K cross-plot (Fig. 8g). The results show that in the three boreholes, the nature of dominating clay minerals are different just as it differs from the singularity spectrum widths Δα. In 1261B, the OAE3 segment is dominated by a single type of clay mineral mainly montmorillonite (smectite) and the estimated dimension of the singularity spectrum width is the lowest (Δα = 1.37). In borehole 959D, three types of clay minerals (montmorillonite, chlorite heavy thorium bearing minerals) dominate the OAE3 unit and the estimated dimension of the singularity spectrum width is the highest (Δα = 1.82) ( Table 3). At SK1S, although the TEOAE3 interval is dominated by three categories of clay minerals i.e. mixed layer clays (illite-montmorillonite mainly), illite and micas (Fig. 8g), the estimated dimension of the spectrum width is less than that of 959D (Δα = 1.51) ( Table 3). This partially falls in line with the findings of Gao et al. 80 , where similar clay distribution results with illite, illite-smectite and chlorite were found by studying SK1S clay minerals distribution. These results suggest that the multifractal behavior of the Th-K distribution in the studied OAE3 intervals is influenced by the nature and diversity of clay minerals, and changes in the deposit environment. This result is consistent with the work of Dashtian 44 . Looking at the above differences of thorium and potassium distributions in these three OAE3 intervals, it suggests that the main factors controlling their scaling behavior during the development of the Coniacian-Santonian OAE are different. Therefore, the long-term persistence of the Th-K distribution might be related to the persistence in the deposition of clay minerals and/or Th-K rich minerals under different sedimentation patterns. Furthermore, the change in the shape of the singularity spectrum for Th/K and Th/U ratios may indicate some variations in the thorium-potassium and thorium-uranium distributions. When the multifractality is sensitive to the local fluctuations with small amplitudes, the singularity spectrum will be found with left side shortness 74 . In all boreholes, the left side shortened spectrum of the Th/K ratio of the studied OAE3 intervals (Fig. 5f) suggests the presence of local changes with small amplitudes in  www.nature.com/scientificreports/ theorized that periodical shift of the intertropical convergence zone controlled by minor orbital-driven changes triggered an increasing rainfall and weathering leading to the genesis of smectite-rich clay 84 . This climatic control is supported by the works of Chamley 81 , who correlated the kaolinite abundance found in the equatorial Eastern Atlantic Ocean to a strong climatic control regulated by the intensity of the continental hydrolysis. Furthermore, the left-side truncated spectrum of Th/K ratio (Fig. 5f) indicates the existence of local intermittency in the thorium and potassium deposition, which may suggest the existence of local fluctuations in the climate evolution. Previous studies based on climate tracers at site 959 showed the presence of highly fluctuating climate modulated at different timescales during the Upper Cretaceous 54 . Thus, since the Th/K ratio depends on the nature of the clay mineral assemblage, which are controlled by climatic fluctuations, we correlate the change in the multifractal features of the Th-K distribution to the change in the climate dynamics during the Coniacian-Santonian OAE3 in the DIB, since several scholars have shown that the multifractality properties of geological data change with the changes in the climatic processes 73,[85][86][87] . At DR (borehole 1261B), the Late Cretaceous clastic sediments derive from weathered granitoid basement rocks of the Guyana Shield 67,88 . These sediments afterwards were accumulated under seasonal and long-term (mainly eccentricity and obliquity) bands 89 . The existence of smectite (montmorillonite) as the predominant type of clay mineral in the Coniacian-Santonian OAE3 interval (Fig. 8g) suggests an increase in chemical weathering 81 under small seasonal climate fluctuation during its depositional time, because smectite is known as a product derived from a seasonal climate, minor changes in continental climate 84 . Based on carbonate analysis recorded alongside DR, Nederbragt et al. 90 also reported a tropical seasonality on DR throughout the Turonian to Santonian. The local intermittency in the thorium and potassium deposition revealed by the multifractal spectrum width of the Th/K ratio may also be due to this seasonal change. Therefore, the abundance of one type of clay mineral implies that the seasonal fluctuation is far from the main factor controlling the thorium-potassium ratio. Consequently, the multifractal nature of the thorium-potassium distribution observed in the OAE3 segment in borehole 1261B may be due to persistence in Al-rich weathered granitoid basement rocks from the Guyana Shield resulting from minor seasonal climate fluctuations chiefly controlled by eccentricity and obliquity. Figure 8g shows that illite, mixed layer clays (illite-montmorillonite) and micas dominate the TEOAE3 unit in SK1S without smectite. This trend was described by Gao et al. 80 as a consequence of burial diagenesis due to the fact that illite persists during burial diagenesis, while smectite is transformed in illitic layers-associated with geochemical changes involving K and Al dissolved from proximate feldspars and mica 91 . Since SLB is enclosed by crystalline igneous outcrops 92 , their weathering would have released Al and K which get integrated in the smectite layers to form illitic clays. This is consistent with the medium correlation between the weight percent of clay and K in SK1S (Fig. 8d) indicating that a substantial part of K may derive from non-clay minerals. Furthermore, based on climatically sensitive proxies, Wang et al. 26 found that SLB did not undergo large climatic variations throughout the Late Cretaceous but experienced abundant rainfall inferring that local fluctuations were common phenomena in the paleo-basin. This seems consistent with the left-side truncated spectrum of Th/K ratio (Fig. 5f) indicating the existence of local intermittency in the thorium and potassium deposition, which may suggest the existence of local fluctuations during the OAE3. Since, the burial diagenesis is the main process controlling the clay mineralogical change below 1,100 m in SK1S 80 , and the basin did not experience large climatic shifts 26 , and as the TEOAE3 interval in SLB ranges from 1,380 to 1,440 m (Fig. 9a), we therefore correlate the multifractality of the lacustrine OAE3 Th/K ratio to a long-term persistence of the local burial diagenesis during the basin subsidence under the control of local tectonic and climatic variations 80,93 , which played an important role during the basin structuration.
Similar to the Th/K ratio, the probable factors affecting the multifractal properties of the thorium-uranium distribution in the studied OAE3 intervals were analyzed in all boreholes by correlating the thorium-uranium cross-plot (Fig. 10) and the singularity spectrum results. Th/U ratio is a useful proxy to track the paleo-redox conditions of the original sedimentary environment and/or subsequent diagenetic processes 94 since it is often strongly connected with the changes in the depositional environment 95 . The depositional environment is probably described as reducing when Th/U < 2 (commonly marine), or oxidizing when Th/U > 7 (possibly terrestrial) 95,96 . Based on Fig. 10, the three OAE3 intervals were deposited in different paleo-redox conditions. In boreholes 959D and 1261B, the Th/U ratio shows that the OAE3 unit is mainly characterized by reducing environment while in SK1S it is characterized by a period ranging from oxidation to reduction sedimentation.
According to San José Martínez et al. 97 and Liu et al. 34 , the shape features of the multifractal spectrum defined by ΔD = D(α max ) − D(α min ) can estimate the probability of the dominate subset in time series. When ΔD < 0, a small probability subset dominates, while a large probability subset dominates when ΔD > 0. Therefore, a period of strong and long-term reducing environment dominates in the OAE3 interval at sites 959D and 1261B and correlates well with ΔD > 0 (Table 4). This is supported by the computed average of Th/U ratio in both boreholes ((Th/U) mean = 0.66 and (Th/U) mean = 1.3, respectively). Previous studies based on geochemical data analysis found a similar result at both sites 17,57,67,89 . Furthermore, the analysis of the singularity spectrum of the Th/U ratio indicates a right-side shortness of the multifractal spectrum (Fig. 5c), reflecting the existence of small-scale intermittency in the thorium-uranium distribution at both sites 74 , which may reveal some shifts in the redox conditions during the deposition of the OAE3 interval. This result is consistent with the works of März et al. 17 and Hofmann and Wagner 67 who described periodic changes between anoxic and euxinic bottom water condition under periodic shifts of the Intertropical Convergence Zone during the deposition of the Coniacian-Santonian OAE3 intervals at DR and DIB. Hence, the multifractal nature of the thorium-uranium distribution in the OAE3 intervals may be due to long-term persistence of reducing conditions at sites 959D 98 and 1261B 55 , which contributed to the deposition of black shale at DR and DIB.
In SK1S, ΔD < 0 (Table 4) suggesting that there is a small probability that neither the oxidizing nor reducing environment dominates during the OAE3. This result is consistent with the computed mean of Th/U ratio ((Th/U) mean = 4.04) in the TEOAE3 interval, indicating that the OAE3 occurred during a transitional period Scientific RepoRtS | (2020) 10:14363 | https://doi.org/10.1038/s41598-020-71327-w www.nature.com/scientificreports/ between reducing and oxidizing sedimentation as also revealed by Jones et al. 23 . As mentioned before, the multifractal singularity spectrum displays left-side shortness when the dataset has a multifractality that is sensitive to the local fluctuations with minor amplitudes. In the case of SK1S, the thorium-uranium distribution suggests  (Fig. 5c). The multifractal behavior of the thorium-uranium distribution may be a consequence of long-term persistence of local changes like increased seasonality, freshening of bottom-water, and/or lake depth reducing which likely contributed to the reoxygenation of the water column and consequently, the demise of the reducing condition 23 .
The computed multifractal spectrum widths of Th/U ratio differ from OAE3 segment to OAE3 segment indicating distinct redox-potential strength ( Table 3). The expression of the OAE3 unit recorded in SK1S displays the highest singularity spectrum width of Th/U ratio implying that this interval was deposited in more complex paleo-redox environment compared to those in boreholes 959D and 1261B, as also revealed by the thorium-uranium cross-plot (Fig. 10) and by the fractal dimensions (Table 5) estimated using the commonly used box dimension (D B ) 99 and Higuchi dimension (D H ) 100 defined by: where N is the number of thorium-uranium data series and δ the scale; L(k) the normalized length of the data series.
We hypothesize that the changes in the paleo-redox conditions influenced the multifractal properties of the thorium-uranium distribution in the studied OAE3 intervals.
Furthermore, the similarities of the paleo-redox conditions during the OAE3 were explored through the root-mean-square (RMS) variation which helps to study the average variation of time series. Indeed, the RMS can be used to segregate the fluctuation function between the amplitudes of the local changes 71 based on the q-order variation. The scaling properties of RMS quantifies the variation of the fluctuation function over the interval with large and minor changes. In marine environment (959D and 1261B), the fluctuation function of Rise based on thorium-uranium cross-plot. Th/U < 2 represents a reducing environment; Th/U > 7 represents an oxidizing environment; 2 < Th/U < 7 represents a transitional period between oxidizing and reducing environment. Table 4. The multifractal spectrum shape features of the thorium-uranium distribution in the studied OAE3 intervals.

Multifractal spectrum shape
Th/U
Since the OAE3 is known as associated with a period of long-term significant change in the Earth's climate dynamism with a major cooling trend spanning from the mid-Cretaceous mega-greenhouse to the normal greenhouse state 10 , this transition likely triggered changes in the depositional systems. Furthermore, astronomical forcing seems to control the black shale deposition in the OAE3 interval at DIB 21,54 , DR 101 as well as the depositional processes in SK1S 25 . Thus, in light of the aforementioned detailed discussion, we surmise that the climatic and environmental disturbances during the OAE3 (probably due to an orbital forcing) induced changes in both the weathering mechanism and the deep-water sediments supply which in turn influenced the clay mineralogy and the bottom water redox conditions in marine and continental basins-and consequently the scaling behavior of Th-K and Th-U distributions.
conclusions. Based on spectral gamma ray log data recorded in the non-marine Cretaceous Songliao Basin and the marines Deep Ivorian Basin and Demerara Rise, for the first time, the multifractality of the Coniacian-Santonian oceanic anoxic event was explored through analysis of the multifractal behavior of thorium-potassium and thorium-uranium distributions both in marine and terrestrial environments using the MFDFA. The following conclusions can be drawn: 1. The q-dependence of the generalized Hurst exponents h(q), classical scaling exponents and multispectral spectrum indicate that the thorium-potassium and thorium-uranium distribution in the OAE3 interval from both marine and terrestrial records exhibit a multifractal nature. The results of the estimated shuffled and surrogate data series show that the multifractality due to long-range correlation is dominant in the Coniacian-Santonian OAE interval in both marine and terrestrial environments. 2. The multifractal behavior of the studied OAE3 intervals is associated with the presence of clay minerals and change in the paleo-redox conditions. At all sites the left side shortened spectrum of the Th/K ratio estimated based on the spectral gamma ray indicate the presence of local changes with small amplitudes in the deposition of clay minerals during the OAE3. Moreover, the shortened singularity spectrum of thorium-uranium distribution reflects the existence of small-scale fluctuation with large amplitudes in the redox conditions at the marine sites (DIB and DR). Whereas in the Songliao Basin, the singularity spectrum indicates the presence of local fluctuations with small amplitudes during the OAE3. 3. The multifractal features, the fractal dimensions and the RMS of the described Coniacian-Santonian oceanic anoxic events based on the spectral gamma ray differ from each other and correlate with the changes in the sedimentation pattern under different paleoenvironmental conditions in both marine and terrestrial environments, suggesting a more local behavior of the event even though the regional character is not neglected.

Data and analysis method
Data. The log gamma records the total natural gamma radiation while the spectral gamma ray log not only indicates the whole intensity of radioactivity, but also estimates quantitively the content of the main radioactive constituents (Th, U and K) in the formation. Because the concentration of natural radioactive elements is mainly a function of the depositional environment and diagenesis 102 , the spectral gamma ray log have been widely used for paleoenvironmental research [51][52][53] . The data used in this work are chiefly the spectral gamma ray logs data ( Fig. 9) originating from downhole logging using the gamma ray spectrometry tool during the first stage of the Continental Scientific Drilling Project of Cretaceous Songliao Basin (SK1), Ocean Drilling Program (ODP) Leg 159 and ODP leg 207. The Th/K and Th/U ratios were subsequently computed prior to the multifractal analysis due to the fact that these two proxies are sensitive to paleoenvironmental changes. Analysis method. The approach applied for the analysis of the spectral gamma ray data in the present paper is the MFDFA. MFDFA has been widely employed to explore the self-similarity nature of nonstationary time series 42,45,73,103,104 . The analysis procedure can be found as follows 38 : • Determine the fluctuation profile Y(t) of the Th-K and Th-U distributions considering x t the data of Th-K and Th-U ratios of length N where x is the mean of the data series x t . • Segregate the N-length Y(t) of the Th-K and Th-U distributions obtained from Eq. (2) into uniform nonoverlapping intervals N s = int(N/s) of length s. Mostly, because N is not a common multiple of s, a short part of the fluctuation profile may be left at the end of the division. To have a complete sampling of the whole segment, the segmentation is repeated backwards i.e. starting from the opposite end of the data series. This process provides 2N s non-overlapping segments of Th-K and Th-U distributions at the end. • Estimate the variance of the 2N s segments of the Th-K and Th-U distributions by statistical approach: where generally, q ∈ R\{0} and s ≥ m + 2. For q = 0, F q (s) is formulated by: When q = 2, the process downgrades to DFA. Generally, the variation of F q (s) helps to understand the scaling behavior of the data series. Therefore, for different values of q, the segregation into 2N s intervals, F 2 (s, v) estimation and F q (s) determination steps are executed many times for several scale s as we want to highlight the scaling behavior of the fluctuation functions Fs of the Th-K and Th-U distributions in the studied OAE3 and non-OAE3 intervals at any scale s for different values of q . • Analyze the scaling properties of the fluctuation functions of the Th-K and Th-U distributions through the log-log plots of F q (s) versus s evaluation. If the fluctuation functions F q (s) and the scaling s are positively correlated as a power law defined by Eq. (7), the data series x t of Th-K and Th-U ratios are described as long-range power law correlated.
In Eq. (7), h(q) defined as the slope of the log(Fs) versus log(s) plot represents the generalized Hurst exponent 38 . Analyzing h(q) helps determine the fractality in the Th-K and Th-U distributions in the studied OAE3 and non-OAE3 intervals. The non-dependence of h(q) on q implies a monofractal nature of the data series while in the case of multifractal dataset, h(q) is mostly q dependence. Furthermore, the h(q) characterizes the scaling features of the intervals with large variations for positive q order (q > 0). On the other hand (for negative q order, q < 0), h(q) describes the scaling features of the intervals with minor variations. Besides, the generalized Hurst exponent and the mass exponent τ(q) define a first-order equation expressed by: To characterize the multifractality in the studied OAE3 and non-OAE3 intervals, the singularity spectrum or multifractal spectrum D(α) of the Th-K and Th-U distributions can be determined via the first-order Legendre transformation where α is the multifractal singularity exponent. By using Eq. (9), a relationship exists between the multifractal spectrum D(α) and the multifractal singularity exponent α: The degree of multifractality or complexity (also known as the width of the multifractal spectrum) of the data series is related to the dependence of h(q) on q 42 and is expressed by Eq. (11). Therefore, the degree of multifractality of the Th-K and Th-U distributions in the studied OAE3 and non-OAE3 intervals was explored based on Eq. (11).
The width of the spectra helps to quantify the strength of the multifractality in the Th-K and Th-U distributions. A narrow spectra width implies weak multifractality in the Th-K and Th-U distributions and vice versa.

Data availability
The data used in the present article are available in the Supplementary Information files.