Continuous wavelet transformation to quantify small-scale cycles of petrophysical properties; a new approach applied in a potential disposal repository of nuclear waste, SW Hungary

Continuous Wavelet Transformation (CWT) was applied to study the small-scale repetitive oscillations of porosity distribution patterns in a 5 m silty-claystone core sample of the Boda Clay-stone Formation. We handled the fluctuations in voxel porosity averages over unequal depth distributions as signals over uneven time intervals. The strength of wavelet analysis lies in the ability to study the fluctuation of a signal in detail, i.e., the wavelet transforms permit automatic localization of the cyclic attributes’ sequences both in time (the depth domain) and according to their frequency (the frequency domain). Thereupon, three main frequency branches (cycles) were discerned: small scale (5, 6.67, and 11 cm), intermediate scale (20, 30 cm), and large scale (66.67 cm). Depending on the CWT coefficients magnitude plot, we were able to detect the developments of porosity oscillation according to the depth variable. Thus, small-scale cycles were seen throughout the core sample., the intermediate-scale cycles were strong in the upper parts of the core sample and dwindled toward greater depths, and the large cycle was predominant in the lower part of the core sample. The cross-correlation of the wavelet coefficients of porosity and rock-forming components al - lows a detailed study of the inter-dependence of such parameters as their relationship changes over time. The distinct peaks at zero lag indicates that the measured wavelet coefficient series were contemporaneously correlated; their strong positive correlations suggest that both exam - ined series respond similarly and simultaneously to other exogenous factors. The results em - phasize that cyclical porosity fluctuations at all scales would concern three main factors; sedi - ment deposition, diagenetic processes, and structural deformation (i.e., convolute laminations).


INTRODUCTION
Time series analysis is the application of mathematical and statistical tests to time-varying data to quantify the variation and gain some physical understanding of the system behaviour (TEMPLETON, 2004).Many systems are monitored and evaluated for their behaviour using time signals; additional information about the properties of a time signal can be obtained by representing the time signal by a series of coefficients based on an analysis function.One example of a signal transformation is the transformation from the time domain to the frequency domain.The oldest and probably best-known method for this is the Fourier transform (MERRY, 2005).
In the frequency domain, Fourier Transform (FT) deals with infinite sines and cosines as base functions to decompose stationary time series (POTOČKI ET AL., 2017).This transform calculates the relative intensity of each frequency component of the whole signal (the global frequency content of signals), but the small segments (stationary parts), calculates the Fourier transform of a windowed part of the signal, and shifts the window over the signal.The short-time Fourier transform gives the time-frequency content of a signal with a constant frequency and time resolution due to the fixed window length.Thus, a short temporal window provides suitable time resolution, but the resulting frequency resolution is low.If a larger temporal window is selected, frequency resolution will be improved, but time resolution will be reduced (TEMPLETON, 2004;MERRY, 2005; SHOKROL-LAHI ET AL., 2013).
The Wavelet Transform is a computationally efficient technique for extracting information about transient and non-sinusoidal signals.The Wavelet Transform (WT) and more particularly the Continuous Wavelet Transform (CWT) is calculated analogously to the Fourier Transform by the convolution between the signal and analysis function.However, the trigonometric analysis functions are replaced by a wavelet function (mother wavelet).From this, the continuous wavelet transform retrieves the timefrequency content information with an improved resolution compared to the STFT (MERRY, 2005;LABAT, 2010;PANDA ET AL., 2000;POTOČKI ET AL., 2017).
Wavelet analysis is a rapidly developing area in many disciplines of science and engineering.For example, the wavelet trans-form originated in geophysics in the early 1980's for the analysis of seismic signals and is now being exploited for the analysis of several other geophysical processes such as spatio-temporal rainfall patterns, layered structures, and climate change.Moreover, the application of one-directional wavelet analysis has become wide spread in various fields of geosciences, such as gravity analysis (MARTELET ET AL., 2001, OUADFEUL ET AL., 2010), lithofacies segmentation (OUADFEUL & ALIOUANE, 2011), geophysics (CHAMOLI, 2009), seismic data analysis (OUAD-FEUL, 2007), hydrology and water resources, discharge and suspended sediment, and reservoir-data analysis, i.e., permeability data (LABAT, 2010;GURLEY & KAREEM, 1999A;PANDA ET AL., 2000;POTOČKI ET AL., 2017), and also to detect gradual and abrupt changes in the sedimentation rate, discontinuities, and superimposed periodic cycles (PROKOPH & BAR-THELMES, 1996).
In this study, the CWT was applied to track the small-scale repetitive oscillations of porosity distribution patterns in a 5 m silty-claystone core sample.The computed tomography scan images (CT) lay the foundation for obtaining voxel datasets on the millimetre scale.Fluctuations in voxel porosity averages over unequal depth distributions were handled as signals over uneven time intervals.The cross-correlation of the wavelet coefficients of the voxels dataset allows quantification of the similarities of the wavelet coefficients of porosity and rock-forming components on large, medium, and small-scale cycles, both contemporaneously and at various lagged values.
Currently, there is no general framework for modeling or tracking the repetitive oscillations of porosity distribution patterns in clastic rocks.However, applying such a method grants a unique insight into the cyclical behaviour of the studied parameter (petrophysical properties, i.e., porosity) and provides a reference basis for interpreting the fluctuation patterns based on sediment depositions, structural deformation, and diagenetic processes constraints.

GEOLOGICAL SETTING
We investigated a silty claystone core sample from the Boda Claystone Formation (BCF).It is located within the Tisza Mega unit, Southern Transdanubia, SW, Hungary.During the Permian, this unit was lo cated north of the equator and belonged to the southern margin of the stable European plate (CSONTOS & VÖRÖS, 2004), where a thick conti nental clastic succession accumulated.It was detached from the southern margin of Variscan Europe during the Jurassic (HAAS & PÉRÓ, 2004;BALLA, 1987;HORVÁTH, 1993) and eventually settled in the southeastern half of the Pannonian Basin (Fig. 1).Thereby, in a closed, sub- siding basin under arid-semiarid climate conditions (KONRAD ET AL., 2010; SCHNEIDER ET AL., 2006), the BCF was deposited in a shallow-water lacustrine environment (playa mudflat, playa lake).
The BCF occurs in two domains in the Western Mecsek Mountains (WM Mts): The Boda block and the Gorica block.Data from boreholes and geological mappings reveal that the extension of the BCF is around 150 km 2 ; only about 15 km 2 is exposed at the Boda village region in the WM Mts.Several deep drillings reached the BCF in the Gorica block, but only borehole Ib-4 recovered the BCF sequence in a significant thickness, c.a. 200 m (Fig. 1).The total thickness of the BCF is estimated to be 700-900 m in the peri-anticlinal structure of the WM Mts (Boda block).In contrast, the overall thickness of the Gorica block does not exceed 350 m (JÁMBOR, 1964;BARABÁS & BARABÁS-STUHL, 1998;MÁTHÉ, 1998;ÁRKAI ET AL., 2000;VARGA ET AL. 2005;LÁZÁR & MÁTHÉ, 2012;NÉMETH & MÁTHÉ, 2016).
The lithology of the BCF is characterized by small grain size (clay to silt) beds.It starts with fine-grained sandstone beds at the base overlain by albitic claystone/siltstone, with successive claystone, albitic clayey siltstone, and silty claystone, and finally dolomite at the top.The green and greenish-grey siltstones and claystones occur infrequently.Desiccation cracks, dolomite concretions, convolute laminations, and crossbedding are recognized almost throughout the entire formation.(KONRÁD ET AL., 2010).The stratigraphic position of the core sample studied within the core sections is shown in Fig. 2.
The actual investigation of the Boda Claystone Formation as a potential host rock for high-level nuclear waste (HLW) disposal began in 1989 and has continued up to the present day.Since 1998 the Public Limited Company for Radioactive Waste Management (PURAM), a Hungarian government agency, has the responsibility and monetary funds for coordinating the research.In 1999, PURAM announced the BCF as a suitable host rock for a highlevel nuclear waste repository.This is due to a set of suitable properties of the BCF influencing radionuclide migration.For instance, it is a rather massive and homogeneous host rock, with a large extent and significant thickness (700-900 m), low bulk porosity (0.6-1.4%), appropriate hydraulic conductivity, and absence of organic residues.

FUNDAMENTAL PRINCIPLES OF X-RAY COMPUTED TOMOGRAPHY (CT)
X-ray computed tomography (CT) can provide unrivaled information about the internal structure of materials, non-destructively, from the metre down to the tens of nanometres scale (ABUTAHA ET AL., 2022).It exploits the penetrating power of X-rays to obtain a series of two-dimensional (2D) radiographs of the object viewed from many different directions.This process is called a CT scan.A computed reconstruction algorithm is then used to create a stack of cross-sectional slices from the object's 2D projections (radiographs) (WITHERS ET AL., 2021).
A single CT scan image is produced using a mono-energetic X-ray.As each X-ray beam passes through the sample, it attenuates varyingly, and the transmitted X-ray is received by a detector (HOUNSFIELD, 1973).The attenuation is measured at multiple angles and reconstructed in a 3D matrix.The 3D distribution of the X-ray attenuation coefficient in reservoir or host rocks depends on mineral composition variations (atomic number and density), porosity, and saturation.X-ray attenuation is physically determined mainly by photoelectric absorption and the Compton effect.Photoelectric absorption is dependent on the effective atomic number and is especially important at low energies (YANG ET AL., 2019; WITHERS ET AL., 2021).The Compton effect predominates at high energies, and the associated X-ray attenuation is mainly controlled by density (WITHERS ET AL., 2021).Just as 2D images are made up of 2D pixels, 3D images are made up of many cubic volume elements called voxels (Withers et al., 2021).The X-ray attenuation can be determined using Beer Lambert's law (Eq.1).Each rotation of the X-ray source around the sample produces a cross-sectional image, which can then be stacked to form a 3D volume. (md)  (1 where I is the intensity of the transmitted X-ray, I 0 is the initial X-ray intensity, m the linear X-ray attenuation coefficient, and d is the length of the X-ray path inside the object.When X-ray energy and intensity are kept constant, linear attenuation of X-ray occurs as a function of density, resulting in the sensitivity of CT images to density changes (HEISMANN ET AL., 2003;DUCHESNE ET AL., 2009).
Series of X-ray attenuation measurements are numerically processed (reconstructed) to show the spatial distribution of Xray attenuation coefficients within the sample; the signals at each point in the reconstructed images are expressed in Hounsfield units.The Hounsfield unit (HU) scale is a linear transformation of the original linear attenuation coefficient measurement into one in which the radiodensity of distilled water at standard pressure and temperature (STP) is defined as zero Hounsfield units (HU), while the radiodensity of air at STP is defined as -1000 HU.The corresponding HU value is therefore given by: where m is the attenuation coefficient of the measured material, and m w is the attenuation coefficient of water.The X-ray attenuation coefficient is represented as CT numbers for a medical CT, calibrated to air with the value -1000 and water with the value of 0 according to the Hounsfield scale.Grey-scale images are generally used to visualize the differences in X-ray attenuation.This process provides a digital 3D grey-scale representation (often called a tomogram).This can be quantitatively analyzed and virtually sliced in any direction, or specific constituents can be digitally colour-coded to visualize the 3D morphology.For example, bright colours (high values) have low porosity, and dark shades (low values) have high porosity in reservoir or host rocks with constant mineralogy and saturation (FÖLDES ET AL., 2004;WESOLOWSKI & LEV, 2005;FÖLDES, 2011).
Measurements with X-ray CT are subject to various errors and image artifacts, including Beam hardening, star-shaped, positioning errors, and machine errors.Techniques used to minimize them have been discussed in full by VAN GEET ET AL. (2000), KETCHAM &CARLSON (2001), andAKIN &KOVS-CEK (2003).

APPLIED X-RAY COMPUTED TOMOGRAPHY
A core sample of BCF (Ib-4), about five m-long, was scanned at a high-resolution X-ray CT facility at the Institute of Diagnostic Imaging and Radiation Oncology, University of Kaposvar, Hungary.The instrument operates at 120 kVp (peak kilovoltage), with a 250 mAs (milliampere-seconds) current, and sampling intervals of 1.0 s.The lateral resolution was (0.1953 x 0.1953) mm 2 with a 1.25 mm scan-slice thickness.The image reconstruction matrix was 512 x 512 pixels.The field of view (FOV) was approximately 9.99 cm (ABUTAHA ET AL., 2021A).
Scans were made using a modified dual-scanning approach (BALAZS ET AL., 2018).Usually, rock samples are dried in a vacuum oven at temperatures of 120 to 210 °F (50 to 100 °C).Drying is terminated when the samples reach a stable weight (SO-EDER, 1986).After six hours of vacuuming the sample, all pore water is removed, and CT measurements were acquired (scan of the dry core).The next phase constituted pumping water into the dried sample (saturation process).After an hour of relaxation, those slices that went under vacuum and saturated were rescanned.CT images were stored in a DICOM (Digital and Imaging Communications in Medicine) format.
A DICOM file contains the scanning parameters and the scanned object identification under different attributes in its metadata.Of these metadata, the pixel spacing, and slice thickness attributes are important for geoscientific applications as they record the dimension (in millimetres) of each voxel in the x, y, and z-directions.Each CT number can be assigned a real-world distance (or depth), allowing CT number profiles to be constructed to calculate the depth and geometrical measurements.DICOM images can easily be read by 'classical' 3D volume rendering software (ABUTAHA ET AL., 2021A, 2021B, 2022).
The laboratory guaranteed that the DICOM files were free of any artifacts and that during the second scan, the same pixels were measured as during the first one.We have calculated the porosity values for each voxel of the image slices from both dry and saturated scans (MOSS ET AL., 1990).

CONTINUOUS WAVELET ANALYSIS (CWT)
The time-frequency resolution problem caused by the Heisenberg uncertainty principle exists regardless of the analysis technique used.Simply, this principle states that one cannot know the exact time-frequency representation of a signal, i.e., one cannot know what spectral components exist at what instances of time.What one can know is the time intervals in which certain bands of frequencies exist, which is a resolution problem (POLIKAR, 1999).
By using an approach called multi-resolution analysis (MRA) it is possible to analyze a signal at different frequencies with different resolutions (WENG &LAU, 1994;WENG ET AL., 2017).Every spectral component is not resolved equally, as was the case in the STFT.The change in resolution is schematically displayed in Fig. 3A, it is assumed that low frequencies last for the entire duration of the signal, whereas high frequencies appear from time to time as short bursts.This is often the case in practical applications.
The wavelet analysis calculates the correlation between the signal under consideration and a wavelet function ψ(t), or the mother wavelet.The similarity between the signal and the analyzing wavelet function is computed separately for different time intervals and on different frequencies, resulting in a two-dimensional representation (scalogram).The continuous wavelet transform is defined as (ADDISON, 2002;POLIKAR, 1999): The transformed signal X WT (τ, s) is a function of the translation parameter τ and the scale parameter s.The mother wavelet is denoted by ψ; the * indicates that the complex conjugate is used in case of a complex wavelet.The signal energy is normalized at every scale by dividing the wavelet coefficients by 1/ s .This ensures that the wavelets have the same energy at every scale.
The mother wavelet is contracted and dilated by changing the scale parameters.The scale s is used instead of the frequency for representing the results of the wavelet analysis.The translation parameter τ specifies the wavelet's location in time (MERRY, 2005).By changing τ, the wavelet can be shifted over the signal as shown in Fig. 3B.For constant scale s and varying translation τ, the rows of the time-scale plane are filled.Varying the scale s and keeping the translation τ constant fills the columns of the time-scale plane.The elements in X WT (τ, s) are called wavelet coefficients.Each wavelet coefficient is associated with a scale (frequency) and a point in the time domain.
The scale s is inversely proportional to the frequency (SHOKROLLAHI ET AL., 2013;FOSTER, 1996).Thus, a large scale corresponds to a low frequency (giving global information about the signal), while small scales correspond to high frequencies (providing detailed signal information) (TEMPLETON, 2004).

METHODS
The workflow followed is shown in Fig. 4. The major aspects of the methodology are briefly addressed in three areas: Pre-processing, wavelet implementation, and cross-correlation.

PRE-PROCESSING
Due to the cooperation between the University of Szeged, GEO-CHEM, and the Public Limited Company for Radioactive Waste Management (PURAM), the raw data sets of a (5 m) vacuumed dry and saturated CT measurements of the BCF core sample were transferred to the University of Szeged for further analyses application.Only the texturally intact, unbroken parts of the scanned sample were used.Five CT volume bricks of the core sample thus served as the basis of the present work (Fig. 3).
A 3D-nearest neighbour algorithm was used to build the 3D volumes of the five scanned core bricks.This process resulted in two Hounsfield lattices, one for the vacuum dried and one for the saturated core volumes.The voxel porosities of the scanned slices were computed from the dual scan lattices subtraction.Their re- Although the output of computed tomography lends itself to straightforward interpretation, so-called scanning artifacts may obscure details of interest or cause the CT value of a single material to change in different parts of an image.The most encountered artifact in CT scanning is beam hardening.Various methods have been developed to reduce or remove the effects of beam hardening (e.g., VAN GEET ET AL., 2000;KETCHAM & CARLSON, 2001;AKIN & KOVSCEK, 2003).One of these is the so-called "subset" CT volumes in which the image's outer edges are removed, and only central volumes of the original three-dimensional images are used for quantitative analysis.
Determination of the HU intervals and calculating the relative percentages of the rock-forming components for the dry-vacuum scan was necessary for finding out and determining smallscale layers (CT layers).So, as soon as the compositional data of the core sample was defined by layers, the dominant rock-forming components were loaned the name of the rock type.

WEIGHTED WAVELET Z-TRANSFORM (WWZ)
The wavelet transform shows great promise as a method for period analysis in time series, particularly for detecting the time evolution of the parameters describing periodic signals.However, when the wavelet transform is applied to unevenly sampled time series, the response of the wavelet transformation is often more dependent on irregularities in the number and spacing of available data than on actual changes in the parameters of the signal.In 1996 Foster brought forward the idea of vector projection concerning the wavelet transformation process of non-equal-interval data.He pointed out that the analysis result could be greatly improved if wavelet transformation was taken as a vector projection.Not only could a desired period be more accurately calculated, but the period's stability could be disclosed.He adopted a Morlete mother wavelet which fluctuates according to a term of the form e iz .The wavelet function used in Fosters' algorithm includes both a periodic, sinusoidal test function of the form exponential (e iw(t-t) ) and a Gaussian window function of the form (e -cw 2 (t-t) 2 ) as a wavelet generating function for transformation ( f(z), see equation 4.
1 4 (4) The constant c determines how rapidly the analyzing wavelet decays; for example, a popular choice for c is 1/8p 2 .However, c could be treated as a parameter, where the wavelet transform would be characterized by three parameters instead of two.Foster (1996) viewed the wavelet transform not just as a projection onto a set of trial functions, but also as a weighted projection (the weighted wavelet transform (WWT)) in equation 5.
N eff the effective number of data points.Vx is the weighted variation of the data, and Vy is the weighted variation of the model function.A better estimate of the frequency of a significant peak would be obtained using a test statistic that was less sensitive to the effective number of data.Fortunately, there is such a statistic for projections; the weighted wavelet Z transform (WWZ).(Equation The ratio of the WWZ to the WWT has been given in the equation ( 7).
the S/N is the estimated signal-to-noise ratio.As S/N is much less than 1, there is almost no difference between the WWZ and the WWT, even in the limit of large data size N. Thus, at low S/N, the WWZ shares the drawbacks of the WWT.Otherwise, on a moderate S/N ratio (1), the WWZ would significantly be improved than WWT.
As the wavelet function fits data (signal), it weighs the data points by applying the sliding window function to the data.One point of the time-scale plane is computed for every scale and time (interval).Thereby, the computations at one scale construct the rows of the time-scale plane, and the computations at different scales construct the columns of the time-scale plane.(FOSTER, 1996;TEMPLETON, 2004).
The (WWZ) algorithm applied in this study is fixed at c = 1/72, corresponding to the wavenumber of 6 that is usually used in the equal (even) Wavelet Transform module.As the PAST 3.2 software package was used, we were able to handle the fluctuations of voxel porosity averages of CT layers over unequal depth distributions as signals over uneven time intervals.

CROSS-CORRELATION
We have calculated the cross-correlations of the continuous wavelet transform spectrum at six representative frequencies, using IBM SPSS Statistical Software, to quantify the similarities of the wavelet coefficients of porosity and rock-forming components on large, medium, and small-scale cycles contemporaneously and at various lagged values.The scales in both variables are defined at uniform frequencies, have a temporal dependency, and the sampling interval in one succession is the same in the other (have the same depth).Two information items may emerge from such a comparison: the strength of the relationship between the two series and the lag or offset in time or distance between them at their position of maximum equivalence.
The cross-correlation is simply the dot product of vectors x and y (Davis, 1986).The number of data points the signal shifts is called the lag.The version of the cross-correlation for the matching position is given in Equation 8.
n* is the number of overlapped positions between the X and Y series.The summations here are understood to extend only over the segments of the two sequences overlapped at the match position.r m , represents the cross-correlation for the matched position.the calculated cross-correlation coefficients series were delineated by the 5% significance level (The 95.0% confidence level) The zero-lag is often set where the origins of the two series are aligned; negative lags represent an arbitrary choice of the sense of the movement of one sequence past another.Since the two series are not identical, the lags in which series A leads series B differ from those in which B leads A (Fig. 5).In our case, wavelet coefficients of voxel-porosity lead the rock-forming components series on positive lags.
The linear trend was removed as a pre-processing step of cross-correlation by subtracting a linear regression line from the transformed wavelet frequencies of the rock-forming components.

QUANTITATIVE ANALYSES OF ROCK TYPES
Based on the most frequent rock-forming constituent, we classified the stacked slices of CT volumes into 80 thin consecutive layers (abbreviated as L).The entire estimation's grid number of the 80 CT layers is more than 40 million voxels.The statistical properties as well as the averaged compositional rock components of each defined CT layer, were calculated and are shown in Table 1.Those 80 layers were used as input data to implement the wavelet transformation.
The layer-averaged compositional frequency percentages have proposed two major rock types; clayey siltstone and fine siltstone (Table 1).The widest variability of the average frequency percentages is noticed in the fine siltstone.Its range varies from 39.24% to 64.84%.In comparison, the claystone variation range is significantly smaller (from 23.76% to 41.18 %).Next, carbonate with an average frequency range extends from 7.15% to 23.41%.Then, the smallest average percentage of the detrital fragments is 0.02%, and the largest is 2.22%.Lastly, the albite average does not exceed 1% as a max, and its minimum is zero.
The total averaged percentages of rock-forming components of each rock type are found to be as follows: in the fine siltstone rock type, for the detrital fragments, it is 0.44, about 55.4 for fine siltstone, 31.5 for claystone, 12.54 for carbonate, and 0.11 for albite.However, the percentages of the rock-forming components in the clayey-siltstone rock type are: 0.22 for the detrital fragments, 45.74 for the fine siltstone, 35.63 for the claystone, 18.16 for carbonate, and 0.2 for the albite.That is to say, the coarser compositional constituents' percentages, such as detrital fragments and siltstone, dominate the fine siltstone rock type.In contrast, the finer components of claystone, carbonate, and albite show an occurrence in the clayey-siltstone rock type.

CONTINUOUS WAVELET TRANSFORMATION (CWT) RESULTS
In Fig. 6A, the wavelet power spectrum graph (scalogram) plots the absolute value of wavelet coefficients.Those X and Y axes numbers show that the CWT was computed at 400 translations (in 300 cm) and 100 scale locations on the translation-scale plane.
The parameter scale in the wavelet analysis is similar to the scale used in maps.As in the case of maps, high scales correspond to a non-detailed global view (of the signal), and low scales correspond to a detailed view.Similarly, in terms of frequency, low frequencies (high scales) correspond to global information about a signal (that usually spans the entire signal), whereas high fre-  These effects in the scalogram arise from areas where the stretched wavelets extend beyond the edges of the observation interval.As such, within the inner region delineated by the black line, you are sure that the information provided by the scalogram is an accurate time-frequency representation of the data.Outside the black line, information in the scalogram should be treated as suspect due to the potential for edge effects.The third dimension on the scalogram indicating the amplitude (power) of a particular frequency at a particular time, is represented as a colourmap (Fig. 6A).The colourmap of the scalograms ranges from blue to red, and passes through the colours cyan, yellow, and orange.The bar on the bottom of the scalogram plot (Fig. 6A) indicates the percentage of energy for each wavelet coefficient i.e., the darker red, the higher amplitude and more considerable value of wavelet coefficient.whereas the corresponding amplitude and energy of the great depression are represented by the blue colour.
In Fig. 6B, we detected the magnitude oscillation of the porosity scale cycles by plotting the CWT coefficients against the depth series changes.Although the small-scale fluctuations are seen all through the core sample (Fig. 6B), the oscillation has peaked at the lower parts of the core sample (deeper depths), i.e., the cycle of 6.67 cm.Conversely, the intermediate oscillations are pretty strong in the upper parts of the core sample and reduce toward greater depths.This corresponds to the gradual colour changes in Fig. 6A on the scalogram profile; for instance, the intense colour of the 30 cm scale has gradually changed from dark red at the shallower depths to almost orange at deeper depth.However, the large-scale cycle in Fig. 6B predominates in the lower part of the core sample.So, in Fig. 6A, its redness intensity has intensified (gotten darker) at the lower part of the core sample.To further validate the observations above, statistical characteristics of wavelet coefficients of the defined scale cycles were calculated and are summarized in Table 2.
In Table 2, the 6.67 cm small-scale cycle shows the highest maximum value of wavelet coefficients (5.24) and the largest coefficient of variation (128.24%).The two peaked spikes at 202 and 250 cm depth in Fig. 6A are almost the cause behind this exceptionally high coefficient of variation, as the wavelet coefficients are dispersed over a larger range.In other words, the wider the range (the difference between the min and max), the more the standard deviation is affected, and the greater the coefficient of variation.In comparison, on the 30 cm scale of the intermediate cycle, the coefficient of variation was the lowest (52.74%).This is relevant to the narrower range standard deviation spreads around (notice the max and min values) and the relatively high mean value (1.78).On the 20 cm scale cycle, it is evident that the standard deviation and mean are almost equal (1.277 and 1.24, respectively), so the variation of the coefficient was close to 100% (102.30%).
Back to Fig. 6A, as was previously mentioned, the continuous wavelet transform analysis is based on integrating the multiplication results of a signal and a complex function (mother wavelet) over time.The product is non-zero when the signal falls in the region of support of the wavelet and zero elsewhere.This raises the question of why the transformation results of the continuous wavelet at the time interval (depth) of 148 -176 (Fig. 6A) is around zero at all scales.
To answer the above question, we need to refer to Table 1.It can be seen that the layers from L39 to L42 cover the exact depth interval of 148-176 cm.The coarser compositional constituents (detrital fragments and siltstone) were highest in L42 and L41, unlike the finer rock-forming components (claystone and carbonate) that were dominant in L40 and L39.The sedimentary signature of these layers (L39-L42) takes the form of normal grading.So, the coarser sediments occur at the base (L42 and L41) and grade upward into progressively finer ones (L39 and L40).The influential occurrence of the claystone (authigenic and allogenic) in L40 and L39 (from the depth of 148 to 176 cm) might enhance a full or partial plugging of the pore throats (ABUTAHA ET AL., 2022).This considerably reduces the average porosity (~2%) and decreases the pore volume variations (low Std.Dev < 3.5).Therefore, the low average porosity in this portion of the core sample, delineated in L39 and L40, is basically the cause behind the weak porosity fluctuations (weak frequency component of the signal) at a depth of 148-176 cm.
On the scale-translation plane, when the window width changes with increasing scale (decreasing frequency), the transform starts picking up the lower frequency components of the porosity signal.As soon as the wavelet window starts covering a depth of 148-176 cm, where the frequency is not a major signal component, the transformation product (wavelet coefficient) gives a (relatively) small value.This means that the frequency component in the signal has a small amplitude (See Fig. 6A, B).However, when the wavelet window picks the higher frequencies (smaller scales), the transformation at the time interval of 148-176 cm would always be zero due to the minority or non-existent frequency component of the porosity signal.This explanation supports the zeros of the minimum wavelet coefficients for the smaller scale cycles (11, 6.67, and 5) and the non-zero minimum values for the larger ones (66.68, 30, and cm cycles) in Table 2. Thereupon, the wavelet coefficients' variation is larger at the smaller scales and lower at the higher scales.

THE CROSS-CORRELATION OF CWT
In the large cycle (66.76cm scale) in Fig. 7, a positive relationship arises between wavelet coefficients of porosity and carbonate and porosity and detrital fragments series at zero lags.The correlation at these match positions is 0.4.However, in the intermediate cycle (30 cm), the highest correlation of the wavelet coefficients of porosity and rock-forming components occurs in two different lags; the first one is at the lag of 100 for porosity and siltstone, and porosity and carbonate, with r = 0.9 and r = 0.7 respectively.The second position is at lag 50 for wavelet coefficients of porosity with albite, porosity with claystone, and porosity with detrital fragments (r ≃ 0.9).Furthermore, at the 20 cm scale, the porosity wavelet coefficient series correlates positively with the wavelet  coefficients of claystone, albite, and detrital fragments at zero lags (r = 0.8).
In the small-scale cycle (11 cm) in Fig. 8, the strongest conformity between the wavelet coefficients series happens in porosity and claystone and porosity and detrital fragments series at zero lag (r = 0.7).The second highest correlation peak, r = 0.5, appears between the wavelet coefficients of porosity and siltstone series.In the 6.67 cm scale cycle at zero lag, the wavelet coefficients correlation exceeds 0.8(r ≃ 0.85) of the porosity and the albite series and decreases to 0.6 in the porosity and detrital fragments.Otherwise, at the 5 cm scale, the porosity wavelet coefficients significantly correlate (r > 0.5) with the wavelet coefficients of all rock-forming components.

DISCUSSION
The distinct peaks of cross-correlation at zero lags indicate that the measured wavelet coefficient series were contemporaneously correlated.Their strong positive correlations suppose that both examined series respond similarly and simultaneously to other exogenous factors.For better understanding of those factors, we need to comprehend the linkage between the geological and petrophysical heterogeneity in each scale cycle.
Heterogeneity, which is the complexity and/or variability of the system property of interest in a space (FITCH ET AL., 2015), does not necessarily refer to the overall system or individual rock unit, but instead it may be dealt separately for a particular parameter, such as porosity (FRAZER ET AL., 2005).In general, the heterogeneity, which reflects the variance of the pore distribution throughout rock types, is caused by sediment deposition, structural deformation, and diagenetic processes.By carefully comparing porosity and the standard deviation in Table 1, it can be noticed that the porosity variation is highly consistent with the standard deviation values.Since standard deviation is a tool to measure the data dispersion around the mean value, we could expect that the higher porosity average, the more heterogeneous pores volume, and the larger the standard deviation.Then, the porosity wavelet's oscillatory power increases.To give ground for this, we ought to consider the porosity classification of this core sample in particular.
In a former study, the CT-layers were classified into three genetic groups based on data-mining techniques (ABUTAHA ET AL., 2021A; 2022); macro-porosity, micro-porosity, and no-porosity clusters.The mean average of the matrix cluster porosity was 3.39%, 10.77% for the macro-porosity, and zero for the noporosity.Ratios of voxel-porosity of the three clusters were as follows: 30.37% for the matrix cluster, 14.65% for the macropores cluster, and 55% for the no-porosity cluster.Comparing the claystone and carbonate ratios across the micro and macro-po- rosity clusters illustrates that the matrix cluster has the higher proportion of claystone and carbonate.Diagenetic clay was also highly prospected to be found in the matrix cluster's claystone component.Authigenic minerals were influential in decreasing and increasing the pore volume ratio.For instance, smectite might swell and become free to move.This suggestion could explain the full or partial plugging of the pore throats.Beyond that, collapse of smectite could lead to the adsorption of less water during the saturation phase, which together with particles released from the detrital deposit (i.e., albite and carbonate cement) structure, could create additional pore space.Alternatively, the macro-porosity cluster has the lowest claystone content, albite is absent, and showed the highest detrital fragments occurrence.Moreover, the pore spaces (macro-pores) might be created within the portion of the detrital fragment because of the larger grain size (ABUTAHA ET AL., 2021B).Eventually, the higher micro and macro-porosity rate, the greater the standard deviations.
For a deeper discussion of the positive harmonizing oscillations of porosity, carbonate, and detrital fragments on the large scale-cycle, we've nominated two cases representing the highest CWT coefficients on the magnitude oscillation plot (Fig. 6B) at depths of 218 cm and 233 cm.They are symbolized as L55 and L58, respectively in Table 1.
Examination of the rock-forming components of L58, where the porosity average (2.86%) has demonstrated the high dispersion value of the standard deviation (4.56) shows that this clayey siltstone layer is composed of 0.3% detrital fragments (consisting of coarse siltstone or/and fine sandstone), 44% siltstone, 33% claystone, 22% carbonate, and approximately 0.5% albite.By way of explaining our example, we could propose three potential scenarios.
At first, as the basement rocks in the source area of the BCF weathered and transported away, dissolved material was carried in solution and concentrated in the upper Permian depositional basin by evaporation to provide a sodium-rich brine.Added Na ions into Boda sediments were responsible for the partial to complete albitization of the detrital plagioclase.The Ca ions released from plagioclase albitization were available to be removed in solution and react with HCO 3-and Mg 2+ (MILLIKEN & LAND,  1991, 1993).The albitization of detrital plagioclase can therefore result primarily in the formation of calcite.some of the microporosity within the albite pseudomorphs is caused by volume reduction during albitization of the detrital plagioclase (MORAD ET AL., 1989).However, dissolving calcite (calcite cement) during periodic ephemeral inundations of the fluvial influx of the playa mudflat sediments of the BCF could be the primary reason for the cyclical pattern of the spatial porosity distribution.Another possible explanation is the intergranular porosity associated with detrital grains, in which microporosity (matrix porosity) is reinforced to exist between different grains such as carbonate and sand/siltstone grains.One more assumption might be related to the so-called fracture porosity, in which void space is associated with unfilled extension fractures developed in grains and/or cement (ADEYILOLA ET AL., 2022).
In the L55 example, the averages of the rock-forming components, where the maximum porosity is 3%, illustrate an exceptionally high standard deviation (4.64), and are: 0.81% for the detrital fragments, 54% for the fine siltstone, 30% for the claystone, 15% for the carbonate and 0.15% for the albite.It is worth noting the discernibly high value of the detrital fragments, which represents the maximum value throughout the core sample.If we compare the detrital fragments and the standard deviation values, it is noticeable that the higher the detrital fragments, the higher the standard deviation (the increasing of the detrital fragments leads to larger standard deviations).To put it differently, the presence of the detrital fragments might increase the variance of pore volume sizes, particularly in the matrix porosity medium.As an illustration, the pore opening between the detrital fragments and matrix grains caused by the dissolution of calcite cement in the initial stage of weathering would gradually connect and widen; thus, creating a larger pore volume and causing a higher standard deviation.Conversely, the absence of detrital fragments increases the proportion of the claystone-siltstone components.In such a manner, microporosity predominantly occurs (low standard deviation).Briefly, albitization of the detrital plagioclase, the dissolution of the calcite cement, and the detrital fragments abundance are seemingly the fundamental reasons behind the heightened porosity oscillation in the large scale-cycle.
According to the 20 cm scale cycle, where the porosity oscillation shows good consistency with clay, albite, and detrital fragments fluctuations (Fig. 7), L15 is illustrative as an example of the highest CWT coefficients on the magnitude oscillation plot (Fig. 6B).The rock-forming components averages (Table 1) in this layer are as follows: 45% of the fine siltstone, 34% of claystone, and 20% carbonate.The detrital fragments and the albite have almost an equivalent value; they are (0.24% and 0. 25%).
To enhance our understanding of the concurrent occurrence of these variables, it is crucial to walk briefly through the (coupling albitization and Illitization formation) process.That is to say, since the detrital fragments -support matrix grains, the rocks can retain more pore space (macro-porosity) (e.g., SAIGAL ET AL., 1988).However, when sodium was diagenetically added from brines to the Boda sediments (e.g., VARGA ET AL., 2005;MATHE, 2015), the predominance of the clay mineral content resulted in the relatively great extent of albitization of the K-feldspar.The reaction of albitization (i.e., K-feldspar) is described below (CHOWDHURY & NOBLE, 1993;SAIGAL ET AL., 1988;NORBERG, ET AL., 2011): The K-feldspar loses K + , gains Na + , and transforms into albite.The different molar volumes of the K-feldspar and albite (the molar volume of K-feldspar is 109.1 cm 3 /mol, and that of albite is 100.2 cm 3 /mol), can generate an approximately 8% volume reduction after albitization of the K-feldspar occurred (NORBERG ET AL., 2011; WILKIN ET AL., 1996).This expresses the feldspar albitization as a volume-reducing reaction (NORBERG ET AL., 2011), by which irregular intergranular pores often form between the replacement albite crystals.Further, the allogenic claystone occurrence (clay minerals formed by depositional process), is expected here.Pores hosted by clay minerals would be held rigidly by silt-sized dolomite and quartz crystals preventing the clays from compaction and ultimately preserve pores (ADEY-ILOLA ET AL., 2022).In conclusion, as detrital fragments might support the macro-pores occurrence in the 20 cm cycle scale, the microporosity (intergranular pores) whether generated within the albitized feldspar grains, or preserved in the allogeneic claystone could promote the porosity oscillatory power and increase the heterogeneity of the pores volume distribution (high std.Dev).
The K-feldspar albitization is often coupled with the alteration of smectite to illite (diagenetic clay) as follows (BOLES &FRANKS, 1979): The equations above demonstrate that the release Na + , Ca 2+ , Fe 3+ , Mg 2+ , and Si 4+ .The Si 4+ would precipitate in the form of quartz.The alteration of kaolinite to illite would consume K+ (BJØRLYKKE ET AL., 1986).The Ca 2+ and Mg 2+ would be involved in the carbonate cement evolution.However, the Na + could be evolved to precipitate diagenetic albite (albite cement) (e.g., VARGA ET AL., 2005;MÁTHÉ., 1998;ÁRKAI ET AL., 2000).Diagenetic albite cement could be precipitated as albite pore filling or preferentially along grain margins as overgrowths (MÁTHÉ, 2015).Once the grain margin is completely albitized, the margin of the grain would be composed of stable albite (MIN ET AL., 2019).This suggests that these types of cement (carbonate and albite) would be able to fully fill the textural composition's pores; consequently, porosity is expected to be almost zero.Notwithstanding, at the small cycle (6.67cmscale) the porosity oscillatory power heightens when albite fluctuation increases.This leads us to suppose two scenarios.The first is dissolution of the albite cement during late diagenetic processes.Or, during the saturation phase, the difference between the pore pressure and the injected water pressure may force some small movable particles to get out from some semi-filled pore space.In this way, the pore volume increases in the exhausted pores (SIMON & AN-DERSON 1990;ZHOU ET AL., 1995;HAYATDAVOUDI & GHALAMBOR 1996;AL-YASERI ET AL., 2015).Along these lines, releasing albite particles from the pore volume could generate an additive micro pore space resulting from particle migration.Those additive micro pores might cause a higher porosity measurement, resulting in a higher overall porosity oscillation (ABUTAHA ET AL., 2022).
At the 11 cm scale cycle, the porosity oscillation correlates positively with claystone, siltstone, and detrital fragments fluctuations.L35 is an instance of the highest CWT coefficients on the magnitude plot (Fig. 6B) at a depth of 124 cm.Here it is observed that this layer selected in the fine siltstone rock type (L35), comprises of 53% of fine siltstone, 30% of claystone, 16% of carbonate, 0.3% of albite, and around 0.6% of detrital fragments.It shows the highest porosity average (3.6%) and the greater standard deviation (5.26).Due to the low percentages of both carbonate and clay, and the dominant occurrence of siltstone, as well as the relatively high detrital fragments ratio, one could expect a low diagenetic process (if any) and low movement of particles during the saturation process.This suggests that this inherent porosity ought to be accurate.So, the porosity fluctuation might be related to the primary textural pore heterogeneity rather than the diagenetic processes (secondary porosity) (ABUTAHA ET AL., 2022).Simply put, porosity fluctuation in such a homogeneous rock compositional layer (see sedimentary feature column in Table 1), could be attributable to the variability of the relevant rock-forming components.That is, as the siltstone component in this core sample was the significant parameter to getting the micro-pores to thrive, the detrital fragments, again, were fundamental in developing/controlling the macro-pores.Thereupon, the standard deviation of porosity has shown the greatest value.
Intriguingly, in the small scale-cycle the porosity wavelet coefficients were positively correlated with all the rock-forming components.However, the multiple short spikes on the magnitude coefficients plot in Fig. 6B are evoked by subsequent short, individual, events that happened at specific depths for short durations.For further explanation, we consider three major spikes as three different cases selected from the magnitude oscillation plot (Fig. 6B); they are: L9, L49, and (L59-L60).
Going through the compositional components of the L9, it is apparent that this layer is similar to the example mentioned above (L35), where Low carbonate and clay are aligned with the dominant occurrence of siltstone and detrital fragments (Table 1).This is mostly spiking the porosity fluctuation to be 2.93% at a depth of 20 cm.In comparison with the L49, it is noticeable that the L9 has higher claystone, carbonate and albite (32%, 19% ,0.4%) and lower fine siltstone and detrital fragments (45%, 0 .5%).In addition, the presence of the convolute lamination structure suggests well mixed, poorly sorted grains so likely to have a low porosity sediment developed.Nonetheless, it has the maximum porosity average (3%) and one of the highest standard deviation values (4.7).In fact, the presence of sedimentary struc tures in this core sample is more commonly interpreted as a type of heterogeneity.Seemingly, internal void space in the layers' sedimenta ry features might be the actual reason for developing the high average porosity of this CT layer.More precisely, the tight siltstone-claystone core samples preparation required drying out and heating in an oven before analyses to remove trapped moisture within the pores i.e., connate water, (the water that remains trapped in pores as the sediments compact and bind together).The drying process is terminated when the samples reach a stable weight, which indicates that all pore water has been removed (HINAI ET AL., 2014).Thereby, the vacuum-dry process altogether with macro-porosity resulted from the detrital fragments and the siltstone/claystone matrix porosity, could affect the available pore space and cause higher-porosity measurements, especially where connate water is expected to occur primarily.
Otherwise, L59-L60 encompass claystone, (~37%).The albite and the detrital fragments are almost zero, carbonate is around 16% and the siltstone average is around 45%.The abnormally high average of the claystone is predominantly attributed to the allogenic and diagenetic claystone.Basically, clay minerals formed by depositional processes tend to have relatively higher porosity and lower permeability initially which is often lessened or eliminated by compaction.On the other hand, clay minerals formed by diagenetic processes reduce pore space and block pore-throats.However, in these samples, there is a positive trend with increasing clay minerals and porosity (ADEYILOLA ET AL., 2022).So, this is either related to the pores hosted by clay minerals, or the sample preparation technique.Drying in the oven at relatively high temperature (about 60 ℃) has been observed to induce potential microcracks and shrinkage in rocks containing clay minerals (HINAI ET AL., 2014).Since our technique involves saturating the samples fully with fluids that are not in chemical equilibrium with the clays, it is possible that the saturating fluids will fill the induced microcracks as well as causing the clay minerals to swell (i.e., smectite is collapsed).The swelling clay effect and microcracks artificially introduces some new pores which ultimately increase the porosity average and provide the spiking peak at a depth of 236-237 cm.

SUMMARY AND CONCLUSION
The general purpose of this research is to investigate the smallscale repetitive oscillations of porosity distribution patterns in a 5 m silty-claystone core sample in the BCF using CT scan images.The methodology applied had three categories: Pre-processing, wavelet implementation, and cross-correlation.
In the pre-processing stage, we classified the stacked slices of CT volumes into 80 thin consecutive layers depending on the most frequent rock-form ing constituent.The statistical properties as well as the averaged compositional rock components of each defined CT layer were calculated.The porosity averages of the 80 CT layers were later used as input data to implement the continuous wavelet transformation.To locate the repetitive oscillations of porosity temporally and spatially, the fluctuations of porosity over unequal depth were handled as signals over uneven time intervals based on the WWZ algorithm.Thereafter, we were able to quantify the similarities of the wavelet coefficients of porosity and rock-forming components on all different scale cycles contemporaneously and at various lagged values, using cross correlation methods.Lastly, the distinct peaks of cross-correlation at zero lags were interpreted by defining the highest magnitude oscillation of the CWT coefficients in terms of geological and petrophysical (i.e., porosity) heterogeneity.
Based on the procedures above, the following conclusions can be drawn: • By localizing the porosity high-energy spectrum within the cone of influence, three scale cycles were defined: small (5, 6.67, and 11 cm), intermediate (20, 30 cm), and large (66.67 cm).Small-scale fluctuations were almost seen throughout the core sample.In comparison, the intermediate oscillations were pretty strong in the upper parts of the core sample (at a depth interval of 0-100 cm) and less strong toward greater depths.However, the large cycle predominated the lower part of the core sample (at a depth interval of 200-300 cm).• The fluctuation magnitude of the CWT coefficients at each scale cycle could be attributed to three main factors; sediment depositions, structural deformation, and diagenetic processes.For illustration: • At the large scale-cycle, the albitization of the detrital plagioclase, the dissolution of the calcite cement, and the detrital fragments abundance were the fundamental reasons behind the heightened porosity oscillation.• On the intermediate cycle (20 cm), the detrital fragments supported the macro-pores occurrence, however, those micro-porosity (intergranular pores) whether generated within the albitized feldspar grains, or preserved in the allogeneic claystone would promote the porosity oscillatory power.
• At the 11 cm scale cycle: the porosity fluctuation might be related to the primary textural pore heterogeneity rather than the diagenetic processes.Thus, the siltstone component was significant to evolution of the micropores, and the detrital fragments were fundamental in developing/controlling the macro-pores.• At the small cycle (6.67cm scale), releasing albite particles from pore volumes during saturation process or even dissolving albite cement in late diagenesis could generate an additive micro pore space, which may increase the overall porosity oscillation.• At the 5 cm scale cycle: spiking porosity peaks were related to short, limited events such as sedimentary structures i.e., convolute lamination, or may even have been produced due to the data records glitches, including microcracks and shrinkage that occurred during the saturation process etc. • The diagenetic products must be figured out as a rule rather than an exception in studying the cyclical porosity behaviour, because the variability of the diagenetic process could modify or even obscure the original depositional porosity and reorganize the spatial pore volume distribution (pore system).In this study, the diagenetic processes i.e., albitizations, had an important influence to enhance the microporosity occurrence, however dissolution processes were, more often than not, associated with macro-porosity.• The detrital fragments were the decisive rock component, improving the porosity oscillatory power relationships at all scales.Thus, when the detrital fragments support matrix grains, the rocks could almost retain more pore space because of the larger grain size (coarse silt to fine sand particles).• Moreover, the presence of the detrital fragments might increase the variance of pore volume sizes, particularly in the matrix porosity medium; so the higher the detrital fragments, the higher the standard deviation.Conversely, the absence of detrital fragments increases the proportion of the claystone-siltstone components.In such a manner, microporosity predominantly occurs (low standard deviation).Therefore, the cyclical porosity fluctuations at all scales would somehow concern the detrital fragments component (intraclasts); it provides a reference basis for the porosity cyclical behaviour.• Sedimentary structures i.e., convolutions would enhance the micro-porosity abundance in a minor way.• In the fine siltstone rock type, where the coarser rock-forming components dominated the composition, i.e., detrital fragments and siltstone components are >50%, the textural heterogeneity would play an influential role to elevate the porosity oscillation magnitude.
A future proposed endeavor is to analyze the spatial distribution of porosity and permeability at each scale cycle.This can add an essential insight into the prediction and understanding of the fluid flow behaviour at different scales of heterogeneity.

Fig. 1 :
Fig. 1: Geological map of the Boda Claystone Formation (modified after Konrad et al., 2010); the star in the upper left corner shows the location of the studied core sample (Ib-4).

Fig. 5 :
Fig. 5: Cross-correlogram of two dissimilar series (modified after Devis, 1986).Shaded intervals indicate portions of A and B that are compared at two lag positions.Zero lag positions indicate both series are aligned at their origin.

Fig. 7 :
Fig. 7: Cross-correlogram resulting from the wavelet coefficients series of porosity and rock-forming components at large and medium scale cycles.

Fig 8 :
Fig 8: Cross-correlogram resulting from the wavelet coefficients series of porosity and rock-forming components at small scale cycles.

Table 1 :
Summary for general characters of CT layers defined, including averaged compositional frequency percentages of rock-forming components, average porosity mean, and standard deviations of porosity.They typically appear from time to time as short bursts or spikes.Therefore, three main frequency branches (cycles) are discerned in bold dashed lines in Fig.6A: small scale (5, 6.67, and 11 cm), intermediate scale (20, 30 cm), and large scale (66.67 cm).The solid black line marks what is known as the cone of influence.The cone of influence shows areas in the scalogram potentially affected by edge-effect artifacts.

Table 2 :
Summary statistics of wavelet coefficients in each scale cycle defined.