Assessment of the influence of astronomical cyclicity on sedimentation processes in the Eastern Paratethys based on paleomagnetic measurements using Discrete Mathematical Analysis

The main objective of this study is to apply Discrete Mathematical Analysis (DMA) to the development of the methodology of cyclostratigraphy. This aim is supported by exploring the magnetic 15 properties of rocks, the lithology of sediments and obtained geochronological reference definitions. The analysis was based on measurements of the variability of the magnetic susceptibility of rocks, which reflects climate variations. Astronomical cycles are global; this makes it possible to carry out a correlation analysis over a large area and on different facial types of sediments, considering their lithology and other sedimentary features. The introduction of modern methods of mathematical 20 processing of geological data is one of the prospective areas for investigation and development in geoscience. Astronomical cycles can be revealed from measurements of scalar magnetic parameters of rocks (magnetic susceptibility as presented by the authors). Specific software developed by the authors allows the processing of measurement data and assessment of the presence of stable oscillation cycles based on the obtained measurement base. The present study attempts to apply mathematical methods to 25 magnetic data using the existing PAST program, which allows spectral analysis of primary data with the construction of Lomb-Scargle and REDFIT periodograms. We interpret the spectral analysis data based on paleomagnetic determinations, considering the available dates for the boundaries of direct and reverse polarity chrons on a general stratigraphic scale.


Introduction 30
In recent years, the rapid development of various techniques to study Earth processes has contributed to the generation of large amounts of heterogeneous data. Therefore, Earth Science urgently requires new methods and approaches to analyse these large datasets . This is especially true for describing complex phenomena and explaining some global events, such as paleoenvironmental changes, sedimentation rates, and fluctuations in insolation. 35 The Miocene stratigraphic scale of the Black Sea region is based on the determination of mollusc fauna, along with dates of the boundaries of regional stages and their subdivisions derived from various geochronology methods. However, some uncertainties remain. This problem provokes heated disputes and impedes the performance of full-scale interregional comparisons and paleogeographic reconstructions. Current efforts to solve the problem include improving the stratigraphic dissection of 40 the studied sediments and modernising methods for dating rocks. These methods include one of the new directions in geology: cyclostratigraphy (dating of strata by astronomical cyclicity, recorded in scalar magnetic parameters), which is often used in combination with mathematical methods. The authors propose integrating geological and mathematical methods to solve some critical issues related to the cyclostratigraphy and paleogeography of Eastern Paratethys Miocene sediments. 45 The present study describes the influence of astronomical cyclicity on the sedimentation processes of the Eastern Paratethys based on long-term fluctuations of insolation, cyclostratigraphy and paleomagnetic measurements, using Discrete Mathematical Analysis (DMA). The cyclostratigraphy methods are based on assessing astronomical cyclicity and its recording in sediments and accompanied by detailed lithological and paleomagnetic studies. Global fluctuations in 50 insolation combined with traditional lithological methods allow an estimation of sedimentation rate and allow us to date the sediments. It is particularly important considering the lack of other methods to determine the absolute age of the sediments. The influence of astronomical cycles on sedimentation is global, meaning this method could be used for large-scale correlations. Determination of astronomical cyclicity using geochronological methods, including 55 paleomagnetic reconstruction and lithological analysis, is made possible by defining the magnetic susceptibility of sediments. Magnetic susceptibility depends on the amount of solar radiation reaching the Earth's surface and reflects climatic fluctuations. Further spectral analysis of scalar magnetic parameters allows the detection of long-term period oscillations of the Earth's axis, the angle of inclination of the Earth's axis to the plane of the ecliptic, and eccentricity. This method can produce 60 absolute ages of sediments with an accuracy on the order of 20,000 to 400,000 years. The methodology proposed by the authors consists of synthesising lithological, paleomagnetic and mathematical methods for solving some cyclostratigraphy issues. The method has been validated on several sections of the Miocene sediments of Paratethys. The distribution of these outcrops on the East and Central Paratethys provides wide geographical coverage for comprehensive regional comparison 65 and allows us to extend the methodology to other sections and outcrops accumulated under different conditions and regimes. Petromagnetic and paleomagnetic studies, including field measurements of the magnetic susceptibility of sediments using a field kappameter, are also significant component.

Materials and Methods
The stages of development of mathematical geology show a clear trend towards solving 70 geological problems using increasingly complex mathematical models. Multidimensional methods are replacing simple statistics. The concept of stochastic processes changed our understanding of geological history and led to the development of geostatistics. Nonlinear models are replacing linear models, and the introduction of fractal sizes has led to the notion of chaotic behaviour. The methods used in cyclostratigraphic studies utilise achievements in all of these areas. It should be noted that cyclic 75 sequences are predictable [Schwarzacher, 1993]. Over the past 100 years, mathematical methods have evolved along with the available hardware, from simple computations in the early twentieth century to laborious mainframe algorithms in the midtwentieth century and the microcomputer revolution in the late twentieth century. The use of rapid algorithms on modern computers has made even the most complex multidimensional methods available 80 to all researchers. The most frequently used mathematical methods in geology can be roughly divided into three groups: • Time series analysis, including spectral analysis (Fourier, spectrogram, and wavelet analysis), autocorrelation, cross-correlation, smoothing, filtering, and extremum search. 85 • Multivariate data analysis, including multivariate distributions and cluster analysis.
• Statistical methods, such as statistical distributions, correlation, regression, and chi-square tests.
Cyclostratigraphy is a new scientific domain in stratigraphy and sedimentology that deals with identifying, describing, correlating, and interpreting cyclic variations in the stratigraphic sequence. In particular, cyclostratigraphy involves applying this knowledge to geochronology by increasing the 90 accuracy and resolution of chronostratigraphic units. It uses the astronomical cycles of the currently known periodicity and an interpretation of the sedimentation conditions. Orbital cycles are the most important cycles. These cycles translate into climatic, oceanographic, sedimentary and biological changes recorded in sedimentary deposits over geological time. Numerous case studies have shown that detailed analysis of sedimentary rocks allows these cycles to be identified with a high degree of 95 confidence. An unprecedented high temporal resolution is available once the relationship between the sedimentary record and orbital forcing is established. This relationship provides a basis for timing the processes occurring in the Earth system [Strasser et al., 2006]. The main task of cyclostratigraphy is to analyse he structure of sedimentary deposits and identify stable signals reflecting the orbital influence. An essential aspect of time series analysis in 100 cyclostratigraphy is transforming raw data into the time domain using calibration points. Once the time series is established, mathematical methods can be used to detect orbital periodicities. Of the groups of mathematical methods listed above, the group 'Time series analysis', especially spectral analysis, is invaluable for solving cyclostratigraphy problems. A brief description of the spectral methods most commonly used in cyclostratigraphy is provided below. 105 Any time series can be analysed in terms of its description in the frequency domain. The classical way to detect frequency components in time series is Fourier spectral analysis. The importance of each frequency component in the time series is established using paired sine and cosine waves. Cyclostratigraphic data are usually discrete. For this reason, a discrete Fourier transform is applied. In a discrete Fourier transform, a time series is multiplied, point by point, by a cosine wave of a specific 110 frequency. The results are summed and multiplied by a constant (2/N, where N is the number of points in the data series). This calculation gives the average amplitude of the cosine frequency component. The calculations continue assuming half the spectrum exists as a mirror image of the actual spectrum at 'negative frequencies'. Since negative frequencies have no physical meaning for time series observations, the average amplitude must be doubled using a constant. The time series is then multiplied 115 by a sine wave of the same frequency. The results are again summed and multiplied by a constant. For each frequency component investigated, the relative size of the average amplitude of the sine and cosine waves determines the average phase of the time series oscillations. Fourier transform can be thought of as reorganising time series data to a different location based on frequency.
To summarise, the Fourier transform is a mathematical operation that takes any waveform and 120 breaks it down into separate sinusoidal components with different frequencies and amplitudes. The components are then presented as peaks in the frequency spectrum. Spectral analysis of data series can be performed using the periodogram method. The periodogram is the square of the modulus of the amplitude of the Fourier spectrum. The Lomb periodogram is a frequency analysis technique for non-uniform data series; it is more suitable for 125 cyclostratigraphic data than the Fourier transform. The Lomb periodogram is invariant concerning the time scale shift, and its statistical properties for non-uniform samples are equivalent to the Fourier transform properties for uniform samples. Frequency analysis by the Lomb method solves the problem of detecting oscillatory processes in data series. However, to study the evolution of the observed phenomena requires a time-frequency analysis. In this analysis, a subset of the sample is selected by a 130 sliding observation window. When using an observation window, the data series is multiplied by a specific weight function. Time series spectra are often characterised by a continuous decrease in spectral amplitude with increasing frequency (red noise). Additionally, time series are often unevenly distributed in time, making it difficult to estimate their red noise spectra. The REDFIT mathematical method [Schulz and 135 Mudelsee, 2002], an advanced version of the simple periodogram, solves this problem by directly fitting the first-order autoregressive process to unevenly distributed time series. In this way, the method avoids interpolation in the time domain and its inevitable shift. REDFIT is used to test if peaks in the spectrum of a time series are significant against a background of red noise from a first-order autoregressive process. 140 As mentioned above, knowledge of the periods can provide important information about the function and its phenomenon. Traditionally, spectral analysis was used to determine the period of a given function. At the Geophysical Center of the Russian Academy of Sciences (GC RAS), a novel technique has been developed, allowing assessment of the consistency of any positive number as a period of the original function. The most accurate answers, if they exist, will be the required periods. 145 The proposed technique is based on DMA [Agayan et al., 2018]. The new method is an original technique for analysing discrete data, developed at the GC RAS. DMA is a series of algorithms united by a common formal basis: fuzzy comparisons of numbers, a measure of proximity in discrete spaces, and a discrete limit. DMA was developed to create discrete equivalents of the concepts of classical mathematical analysis: for example, limit, continuity, smoothness, connectivity, monotonicity, and 150 extremum. DMA methods and algorithms have proven to be useful in numerous studies related to the processing and analysis of various geological , geophysical [Gvishiani et al., 2008a], geomagnetic , seismological [Gvishiani et al., 2016] and other data [Gvishiani et al., 2008b]. In the following section, we provide a strict formulation of the problem and a description of the 155 methodology. Suppose that a function f is given on the segment [a, b], and we want to understand what periods it has. The period Т, ideal for f, is defined as follows: T is the period for f if ∀ t, t+T ∈ [a, b] is true for f(t) = f(t+T). In reality, the function f may not have any ideal period T, but it may have approximate periods. A description of approximate periods and how to determine them is provided below. 160 It is clear that the applicant for the period Т must lie in the interval (0, (b-a)/2). For such a T and t ∈ [a, a+T), we denote by {f, T, t} the finite sequence {f(t), f(t+T), f(t+2T),…}. Firstly, it is necessary to understand to what extent the sequence {f, T, t} can be considered constant, i.e., to determine the exponent of its constancy C{f, T, t} ≥ 0. The equality C{f, T, t} = 0 must correspond to the constancy of the sequence {f, T, t}. 165 Here we present three different constructions of the exponent C{f, T, t}. Firstly, the variance D of the sequence {f, T, t} relative to the usual or average median: C{f, T, t} = D{f, T, t} or C{f, T, t} = Dm{f, T, t}. (1) The second and third constructions presuppose the determination of the modulus of the difference sequence |∆{f, T, t}| for consistency {f, T, t}: The question of the constancy of {f, T, t} is reduced to the question of the triviality of 170 |∆{f, T, t}|. The second construction is the Kolmogorov mean Mp with a positive exponent p of the sequence |∆{f, T, t}|: C{f, T, t} = Mp|∆{f, T, t}|. ( As mentioned previously, using the Kolmogorovsky mean Mp, C indicates the proximity of the sequence |∆{f, T, t}| to zero. 175 Another solution is given by the third construction, using the distribution function F|∆{f, T, t}|: C{f, T, t} = α(|∆{f, T, t}|), where α(|∆{f, T, t}|) is α-quantile of the distribution F|∆{f, T, t}|. The estimate C{f, T, t} characterizes T as a period on the sequence {t, t+T, t+2T,…}. The next stage is the development of a unified estimate T as a period independent of t. Such an estimate C(f, T) should be an indicator of the smallness of the collection (C(f, T, t), t ∈ [a, a+T)). For 180 this, we again use the Kolmogorov mean Mp: The final stage is the search for strong minima of the estimate C(f, T). Provided that they exist, we can obtain the necessary periods for the function f using the apparatus of minimality measures [Agayan et al., 2018].

185
In the following section, we present the results of applying the above spectral methods to an example time series taken from [Lisiecki and Raymo, 2005]. The time series represents a series of oxygen isotopes of calcite in the shells of Pliocene and Pleistocene benthic foraminifera collected from 57 core samples worldwide. This series contains the last 400,000 years of the record in 1,000-year increments (Figure 1). This series was chosen for demonstrating the application of spectral methods as 190 it has been comprehensively studied.  the Fourier and Lomba periodograms [Carbonell et al., 1992]. Secondly, the differences in periodograms are caused by the presence of minor noise in the original data. The Lomb method [Lomb, 1976] better distinguishes periods of 106,380 and 65,360 years from this noise. The REDFIT spectrum (Figure 2b) shows five peaks at the frequencies (in descending order of spectral power): 0.0094, 0.0250, 0.0163, 0.0425, and 0.0325. Note that the two main peaks in the 210 REDFIT spectrum are located at the same frequencies as the main peaks in the Fourier and Lomb periodograms. In this case, two REDFIT peaks are shifted in the frequency domain, and one peak combines two Fourier and Lomb peaks. The discrepancies may be because REDFIT considers the spectrum of red noise in the time series, which is not accounted for in the Fourier and Lomb transforms. The results shown in Figure 2 demonstrate that spectral methods for studying time series are 215 valuable tools in mathematical geology, particularly in cyclostratigraphy for detecting orbital periodicities. Academy of Sciences is demonstrated using the example of the time series shown in Figure 1. Figure 3 shows the calculated constancy exponent of the series in blue and its smoothed version in red. As mentioned previously, the required periods of the time series are the points of the minimum constancy exponent C(f, T). Figure 3 shows that the strong periods are 115,000, 44,000 and 85,800 years. Weakly expressed periods are 58,600 and 67,600 years. Note that the identified periods are relatively close to 225 the periods identified by spectral methods. Four out of five periods in Figure 3 closely match the periods in Figure 2. A possible explanation for the appearance of the 85,800-year period may be that it is a multiple of the period of 44,000 years. The results shown in Figure 3 indicate that the proposed method for determining the periods of time series can be used simultaneously with spectral methods in practice. 230 We apply the above spectral methods, and the new approach developed at GC RAS to search for periods in magnetic susceptibility data of the Pontian (the total thickness of the sediments is 41.2 m; Figure 4a) and Lower Maeotian deposits (the total thickness of the deposits is 30 m; Figure 4b) of 235 Zhelezny Rog Cape. The section is located on the Black Sea coast of the Taman Peninsula (Russia) and is a reference section for the Pontian regional stage of the European part of Russia. Measurements of the magnetic susceptibility of rocks were carried out using a field kappameter across the strike of the layers at regular intervals of 20 cm. Three measurements were taken at each point to ensure accuracy. For each point, the average susceptibility was calculated from the three measured values. The averaged 240 magnetic susceptibility values are shown in Figure 4.  Figure 5 shows the results of applying spectral methods to the susceptibility data for the Pontian deposits. Of the three prominent peaks on the Lomb periodogram, two peaks fall within the 95% 245 confidence interval (Figure 5a). These peaks are located at frequencies of 0.0425 and 0.2519, corresponding to periods of 23.53 m and 3.97 m, respectively. This result is confirmed by the Fourier periodogram (Figure 5a), where the main peaks are located at frequencies of 0.0437 and 0.2519periods of 22.88 m and 3.97 m. In contrast to the example described above (Figure 2a), the Fourier and Lomb periodograms in Figure 5a differ slightly. These differences can be explained by white noise in 250 the original data ( Figure 4a). Figure 5b shows the spectrum obtained using the REDFIT algorithm. Two peaks fall within the 95% confidence interval. They occur at frequencies of 0.0444 and 0.2525, which correspond to periods of 22.52 m and 3.96 m. These periods coincide with the periods identified by the Fourier and Lomb periodograms. 260 Figure 6 shows the Fourier and Lomb periodograms and the REDFIT spectrum for the magnetic susceptibility data from the Lower Maeotian deposits. Only one peak of the Lomb periodogram falls within the 95% confidence interval (Figure 6a). This peak is located at a frequency of 0.1408, which corresponds to a period of 7.1 m. This result is partially confirmed by the Fourier periodogram, where the main peak is at a frequency of 0.1516 (period 6.6 m). We are inclined to believe that the differences 265 in the periodograms can be explained by noise in the data. The REDFIT spectrum is shown in Figure  6b. In this spectrum, just one peak is located in the 95% confidence interval. It occurs at a frequency of 0.1497, corresponding to a period of 6.68 m. The similarity in the periods obtained from the Fourier periodogram and the REDFIT spectrum may indicate more white noise than red noise in the initial data. We apply the new algorithm created in the Geophysical Center to selected periods. In Figure 7, the calculated constancy values for the magnetic susceptibility data of the Pontian and Lower Meotian 275 deposits are shown in blue. As discussed previously, periods are considered to be the points of 'global' minima of the constancy exponent. For this reason, the calculated constancy exponents were 'strongly' smoothed to obtain their global trends. The smoothed constancy exponents are shown as red lines in Figure 7. Figure 7a shows that the prominent period in the magnetic susceptibility data series of the Pontian deposits is 3.6 m. The minimum of the constancy exponent at 15.2 m can be considered a 280 'weak' period. The absence of the ~23 m period in Figure 7a is due to the algorithm considering values in half of the studied segment to be contenders for the period, which in this case is 41.2 m. For the Lower Meotian deposits, only one period of 7 m can be distinguished (Figure 7b). Strong periods in Figure 7 emerge as very close to the periods obtained using spectral methods of analysing data series. In the case of a 'softer' smoothing of the constancy exponent (Figure 7), weaker periods can be obtained. 285 In Figures 5-6, these were not valid above the 95% confidence intervals or could not be obtained due to the low spatial resolution of the initial data.

Conclusions
This study presents new experimental data of the Pontian and Lower Maeotian of the Black Sea Basin (Paratethys) obtained by Time series analysis of magnetic susceptibility data from relatively deep-water sediments exposed in the Zheleznyi Rog Cape section (Taman Peninsula, Russia).
The DMA-based algorithm developed by the authors allows us to carry out the main tasks of 295 cyclostratigraphy effectively. Periodograms were constructed to further identify repetitive cycles for carrying out geological reconstructions and correlating sediments of the same age over a large area. In the studied interval, a ∼71-m-long sedimentary sequence, spectral analysis revealed statistically significant signals with some highland peaks. These signals correspond to the precession and obliquity cycles. The 3.6 m peak corresponds to the precession periodicity (19-24 thousand years). The 7 m peak 300 corresponds to the periods of changes in the angle of inclination of the Earth's axis (41,000 years). The 15.2 m peak corresponds to 100,000-year cycles; however, its validity is questionable due to the length of the interval (15×3 = 45 m). The 23 m peak is not valid, as the sampling interval is 41.2 m (cycle lengths are valid when they are three times the thickness of the interval: 23×3 = 69 m). The analyzed parameters of the magnetic susceptibility of the Zhelezny Rog section are consistent with the data from 305 the study of these strata using the Past and AnalySeries programs, the selection of astronomical cycles was carried out [Rybkina and Rostovtseva, 2014;Rostovtseva and Rybkina, 2017]. This study correlates the main steps of Messinian Salinity Crisis (MSC) of the Mediterranean to the Black Sea Pontian record based on astronomical tuning of the study sequence and evaluation of integrated biostratigraphic, paleomagnetic and sedimentological data [Hsü et al., 1973]. 310 It is necessary to say in conclusion a few words, about the mathematical apparatus of DMA [Agayan et al., 2018]. It is being developed at the Geophysical Center of the Russian Academy of Sciences and forms the basis of the methodology for searching periods in the series of cyclostratigraphic data presented in the article. It is important to emphasize that DMA is a direction of modern applied systems analysis [Zgurovsky and Pankratova, 2007]. 315 DMA has all the necessary tools to generate mining algorithms for geological and geophysical data, including searching for hidden periods / cycles. Based on fuzzy sets and fuzzy logic, DMA has the ability to convey expert ideas about the structure, morphology, monotony, and other of studied data series. Thus, DMA enables a systematic approach to the analysis of complex data series of Earth sciences. 320 The Geophysical Center of the Russian Academy of Sciences has a vast experience of application DMA-structures for solving various geophysical problems. Among them are recognition of low-amplitude geomagnetic pulsations of various frequency ranges from 1 to 30 MHz (Pc3, Pi2, Pc5) on records of geomagnetic field variations; recognition of anomalies on magnetic records of the INTERMAGNET network; recognition of areas prone to strong earthquakes; recognition of anomalies 325 in the records of residual gravity field variations of the global network of superconducting gravimeters GGP Network; recognition of anomalies on the records of the intrinsic electric potential of the Piton de la Fournaise volcano (Reunion Island, France); identifying areas of stability in radar interferometry images for monitoring the activity of Mount Etna (Sicily); analysis of a complex interference magnetic field in the Gulf of Saint-Malo (Brittany, France); direction assessment of the magnetization vector for 330 geological bodies in the area of the Ahaggar massif (Algeria) and others Gvishiani et al., 2008a;Gvishiani et al., 2002;Zelinskiy et al., 2014;Soloviev et al., 2005;Zelinskiy et al., 2014;Soloviev et al., 2015;Gvishiani et al., 2003;Gvishiani et al., 2004;Kulchinsky et al., 2010;Mikhailov et al., 2003;Zlotnicki et al., 2005]. Thus, the results of this work show that DMA is also applicable for solving cyclostratigraphy problems. 335

Code availability
The code presented in this study is available on request from the corresponding author. The code is temporarily not publicly available due to research policies and implementing research programmes.

Data availability 340
The data presented in this study are available on request from the corresponding author. The data are temporarily not publicly available due to research policies and implementing research programmes.