Shear Velocity from Seismic Ambient Noise to Determine Crustal Structure and Re-Fitting Rock Physics Modeling in Dhafriyah and Kumaite Oil Fields

Abstract


Introduction
The ambient noise(ANT) is not the earthquakes or human activity but it's related to the continuous motion of the earth surface.There are two types of random wavefields signals in the earth, the seismic coda and ambient seismic noise.The ambient seismic noise that no relation with the earthquakes, but it's generated in the earth naturally.So, this kind of noise can be recorded at anytime and anywhere and the opposite of the seismic coda (Yang and Ritzwoller, 2008).The ANT cross-correlation between two seismic stations can determine the crustal structure by using wavefields records to recapture the Rayleigh waves.This strategy is more helpful than other conventional strategies because it doesn't depend on source quakes or boundaries parameters.Generally, applying the cross-correlation of recorded wavefields can give an assessment of the empirical Green's function (EGFs), that is comprises of the Rayleigh waves that movement between the two seismic stations as if they were generated at one of the stations (Buffoni et al., 2018).This Green's function will assist with knowing the information about the earth's crust between the two seismic stations.Numerous specialists (e.g.Buffoni et al., 2018;Bensen et al., 2007;Yao et al., 2006;Sabra et.al, 2005Snieder, 2004;Wapenaar, 2004;Lobkis and Weaver, 2001) utilized empirical Green's function between two stations that can be assessed from the assortment of cross-relationships of consistent records information of ground movement at these stations.They affirmed that its development requires diffuse wavefields or systematic distribution of ANT sources.Utilizing these functions, the phase and group velocities can be assessed as a period function and can be measure by utilizing traditional time-frequency analysis (Lobkis and Weaver, 2001).Sabra et al. (2005) and Shapiro et al. (2005) designed velocity model in the Southern California at short periods (7.5-15s) by using the ambient seismic noise.They found that the igneous rocks have high-velocity and the sedimentary basins have low-velocity.The study area is located near platform flank of Mesopotamian Foredeep Basin in the southern part of Iraq, extending from Wasit to Missan Governorate, specifically in the area of Kumate, Dhafriyah oil fields.where, Kumate oil field located Missan Governorate to the north west of Amara city and, Dhafriyah oil field located Wasit Governorate ,12 Km north east of Kut city.Generally, the area is located between latitude (32° -33° N) and longitude (45.5° -47° E) (Fig. 1).The current study is tried to comprehend the geological and tectonic setting by estimating the Swave velocity that leading finally to determine the cortical structure of this region.And re-interpretation the general rock physics modeling by integration S-wave velocity with 3D seismic cubes in Dhafriyah and Kumaite oil fields.

Materials and Methods
Through six months have been prepared in AMR, ANB, KAR stations that were recorded from April, May, June to November, September, December for ANT procedure.The data was recorded in GCF format by Scream package (Seismometer Configuration Real-Time Acquisition and Monitoring).Ambient noise processing starting with raw data processing, crosscorrelation, measuring group and or phase velocity, then choose the optimum measurement.Consequently, the information's that's includes (station name, longitude, latitude, elevation, etc.) add as a header info.for each hour by script called SAC. then the data placed on folder has the Julien day number.There is one file must be extracted that represents 24 hours, this file has been executed by command called merge.After merge step, the instrument response must be removed to conjugate the frequency-dependent timings recorded at different types of seismometers, also it's necessary to isolate earthquake signal from ANT.One of the important steps in preparing the data is temporal normalization.This procedure reduces the effects of earthquakes on the cross-correlation and eliminates noise source near the station (Bensen et al., 2007).The one-bit normalization is chosen in this study due to it is the most aggressive, retaining just an original signal's sign by replacing all negative amplitudes with -1 and all positive amplitudes with 1 (Schimmel and Gallart, 2007).If the geometry between two stations Ammara and Anbar presented in Fig. 2. The medium velocity (assumed constant) and inter-station distance symbolize by c and r, respectively.Then, the travel time difference between the 2 stations is Δt = r cos θ/ c.The signal crosscorrelation function at 2 stations is a sum of a cross-correlation functions for every relevant direction, assumption the sources are non-correlated.most principally, if a seismic wavefield traveling between two points X1 andX2 with apparent velocity v controlled by elastic conditions of subsurface, so the wave package u (x, t) crossing position X1to X2 with ∆t= |x1-x2|/v.Therefore, we expect the wavefield recordings u (x1, t) and u (x2, t − ∆t) to be statistically similar, or correlated as Nakata et al., (2019): Generally, the Eq. 1 represents the cross correlations from uniformly distributed plane waves.This Eq. can be expressed as for individual plane waves in the cross-correlation structure (time-domain representation):

Re-Fitting Multiple Filter Technique (MFT)
From noise-correlation traces, this technique was utilized to compute the phase and group velocity dispersion curves as a computer programs seismology part (Herrmann, 2013).the computing analysis is done by multiplying the time-domain cross-correlation traces to the Gaussian filters, and when compared to a function of Green's between 2 stations at different frequencies.This process' time-domain representation is as: Where: k = wavenumber r = distance t = time φ= source phase The filtered signal is produced, after applying the narrow-band Gaussian filter around the central frequency ωo and some algebra, by: Where: α = the filter width Uo = the group velocity The last term is the signal's envelope, which peaks at the time of the group's arrival tg = r/Uo.At a time of the group's arrival (tg), the phase will be: The phase term must be π/4 if a signal is linked to a function of Green's between 2 stations.As a result, based on an analysis of numerous filters, the phase velocity may be derived from: Extracted Rayleigh wave dispersion curves using stacked cross-correlations of phase and group velocity are shown in Figs. 3 and 4. MFT may be used to separating the single dispersion and removing other path effects and arrival overlap.The MFT calculates the velocity of a group or phase as a frequency function (Dziewonski et al., 1969).In current research just group velocity is utilized to conduct the results because phase velocity gives us higher values doesn't commensurate with the magnitude of the problem as shown in Fig. 4. Fig. 5 depicts the Rayleigh-wave dispersion curves generated using the surf 96 inversion program's velocity of group and phase values.The contrast between the observed dispersion data and the predicted dispersion curves in solid red is shown in these figures.The average Root Mean Square (RMS) between an observed and predicted dispersions is a same.This means that an S-wave velocity model we used to produce a projected dispersion curve is closer to the crustal velocity structures along the path.

So
• if field is the homogeneous source • Ignoring the scalar amplitude coefficient • Then, in a frequency domain, cross-correlation is calculated by integrating all plane waves over an angle:  Where (ω) is angular frequency.This can be solved this integral using the zeroth-order Bessel function of a first kind, for a homogeneous medium, the cross-correlation of uniformly dispersed plane waves is proportionate to a part of imaginary of a Two-dimension Green's function (Sanchez-Sesma et al., 2006) as shown in the fallowing Eq: For high frequency and far area limited (ωr/c≫1), the previous Eq.may be approximated by: Cross correlations between uniformly distributed plane waves can be calculated using this equation (Abramowitz and Stegun, 1970;Boschii et al., 2012;Boschii and Weemstra, 2015).

Shear Velocity Model and Final Results
The next phase in the preparation of data is cross-correlation and stacking.Collecting tens of crosscorrelations is a typical research method when ambient seismic noise data is processed across wide geographical ranges.The highest distance between the stations in this study was 420km.Three paths have been considered here, AMR-ANB, AMR-KAR, and KAR-ANB, to cover study area.Correlations between -300-300sec.were saved in this phase, after which the E and N horizontal components were rotated into radial (R) and transverse (T) directions respectively for every pair of the station to generate one cross-correlation (T-T and R-R) (Lin et al., 2008).From 3 routes, was obtained 9 components.The 3 paths represent AMR-ANB, AMR-KAR, and KAR-ANB, the interstation was 20km, 56km, and 420km respectively.All components of paired stations were band-pass filtered daily between 0.01Hz and 8Hz to consider a signal-to-noise ratio (SNR) (Fig. 6).Applied the scale of a signal-to-noise ratio (SNR) to waveforms to aid in separating the needed signal level from background noise, where the presence of this noise causes data corruption.Every day, the cross-correlations has been calculated between station pairs and then added or stacked them for 4 months to generate two-time series.The stacks for the 6 months are April, May, June, November, September, and December, in that order.After obtained a stack file, the Rayleigh waves in a period range of interest were chosen.The group and phase velocities were measured using the multiple filter technique (MFT) once a Rayleigh wave was chosen (Dziewonski et al., 1969).After processing the waveform record with only R and Z components, utilized an MFT to calculate the Rayleigh waves group and phase velocity dispersion curves.Because love waves analysis was outside the scope of this study, the T component was not employed.There are four steps to solve inverse problem (Zhao, 2015).First, model parametrization: by divide the spatial slowness perturbation into cells of constant slowness.So, if the travel time along a ray is a path integral of slowness (one over velocity) along the ray path: Where u is slowness (one over velocity), s is the length of the ray path as it incremental that linked receivers to sources along medium with velocity (v).that is means the Eq. 5 is nonlinear because of the ray path itself determined by structure velocity; this can usually linearize by truncated Taylor expansion: Where ray0 is the reference ray path.There are two assumptions to linearize reference model, one if area have no heterogeneous velocity structures, then path will be determined by reference model with straight circles.The other one if there is lateral heterogeneity then the problem can be solved iteratively, until reaching the optimum model.Second, forward computation: for this purpose, we can use wavefront tracking or Ray tracing is used to calculate data from model parameters and is based on the eikonal equation that is the high-frequency assumption: where (x) is the position vector and (Tt) is the wavefront travel-time.then the pervious Eq. can be totally recall in another form called the ray-path equation (Lay & Wallace, 1995): Third, inverse calculation: the gradient-based methods are utilized to solve the inversion step, which inverts observable data for model parameters, minimizing the mismatch between data predicted from a data of observed and second.Four, uncertainty and resolution estimation model: the quality condition of gained model parameters is evaluated by two approaches; one is to compute the covariance and resolution matrices based on a linearization of an inversion's final iteration.Another one is to utilizing different synthetic tests, for instance spike tests or checkerboard tests (Crosson, 2007).
After reaching 2D velocity maps, the depth inversion process can be initialized, following (Herrmann, 2013) for depth inversion and the production of depth sensitivity kernels, the variational principle for Rayleigh waves is applied.The following is a solution to the Rayleigh wave's eigenvalue problem: Where I0, I1, I2, and I3 are: Where: μ is shear modulus, λ is Lame's parameter, and ρ is density.R1 and r2 denote the horizontal and vertical Rayleigh-wave displacement eigenfunctions.Eq.14 yields the general form of an eigenvalue equation in terms of k and eigenfunctions r1 and r2 for a given frequency.The linearized relationship between phase velocity disturbances and material property perturbations μ, λ, and ρ is determined to invert for depth based on Eq.14.An inversion is nonlinear that may be solved iteratively given the starting model.Strata with constant density and velocity are split into dm-thickness layers in the medium.The partial derivatives of shear-wave velocity (β), P-wave velocity (α), and a phase velocity concerning ρ are found in layer m are calculated: The group velocity U in the previous relations is calculated by: In each progression of the iterative inversion, the distinction between the noticed and anticipated dispersion curves is calculated.The process is continued till it converges to an elastic model (density compressional and shear velocities) which satisfies data (the local phase-velocity dispersion curve).The information of Rayleigh-wave dispersion at seismic station can be utilized and inverted to the shear velocity model along the path between the source and the receiver.In this research, the initial model was designed and then iterated (usually thirty times) until the best fit was accomplished between the predicted Rayleigh-wave dispersion and observed Rayleigh-wave dispersion.The inversion results show that the crustal structure has three discontinuities; these are a discontinuity within the sedimentary thickness 10 km, the basement rocks discontinuity at 12km.The first path AMR-ANB, Fig. 8A show sediments cover 10 km, lower crust 28 km, upper crust 14 km, The Conrad discontinuity 20 km.The second path AMR-KAR Fig. 8B reveal that the sedimentary cover at 4 km to 10 km, lower crust 27 km, upper crust 14 km, The Conrad discontinuity 19 km.while in third path ANB-KAR fig.(8, C), the sediments cover end at 10 km, lower crust 22 km, upper crust 12 km, Conrad discontinuity 20 km.The final velocity model for whole area is represented by Fig. 8, D where the sediments cover starts from 1-4 km and 4-10 km, lower crust 27 km, upper crust 14 km, Conrad discontinuity 20 km. the output of final velocity model is represent by Table 1.The results for S-velocity model in current research are closely related to (Ramthan et al., 2020).In fact, the difference in resolution between seismic and ANT must be taking in our consideration, for this reason there's some ambiguity in small intervals measurements of ANT, but the average picture of rock physics modeling by ANT is closely related to results from seismic technique.The intervals of rock physics modeling in Kumaite field, Figs.9&10 showed decreasing in ROHB, Vp,Vs and Vp/Vs within Mishrif(-1900 sec.) and Nuhr-Umr(-2150 sec.)formations zone , where density value ranges from 2.2 GM/CC to 2.4 GM/CC , Vs value is within 2400 m/s to 2800 m/s , acoustic impedance value approximate to 9500 Kpa.s/m-10000 Kpa.s/m .The Vp/Vs ratio also decrease within productive layers reaching out (1.5 -1-6).The coverage of ANT will give us a lot of details deeper than seismic well logs reached out 8000 m, within Kumaite seismic cube there are two zones could be recognized below well log coverage, the first one began with 3000 sec.to 3500 sec.characterized with decreasing in RHOB, AI, Vs and Vp/Vs.that refer to probable fluid layers.
The second one began within time 3500 sec.to 4500 sec.characterized with higher RHOB, AI, Vs and Vp/Vs.In Dhafriyah oilfield the situation is different Figs.11&12 where the productive layers are within the Zubair Fn. (-2300 sec.).Where RHOB decreasing from 2.2 to 2.3 GM/CC, Vs (3000 m/s-3400 m/s), AI (9500 Kpa.s/m -10000 Kpa.s/m), Vp/Vs is higher than Kumaite range from (1.5 to 1.7).There are two zones in Dhafriyah oil field have been recognized below well logs coverage, it covered by ANT.The first one start with 2500 sec.to 3600 sec.characterized by general decreasing in RHOB, AI, Vs and Vp/Vs. it was located within Jurassic period and it was probably containing major layers of fluid.The other zone is located within 4000 sec.to 5000 sec.and characterized by higher RHOB, AI, Vs and Vp/Vs.That was probably refer to higher compaction circumstances and lack of fluid presence.

Conclusions
Understanding more about the sedimentary thickness that most occurs of earthquakes is crucial.Thus, the technique of ambient seismic noise was used.The ambient seismic noise method has become increasingly popular in recent years due to the massive amounts of data recorded by seismic stations.The dispersion curves of surface waves may be extracted from the group or phase velocity, which both were utilized in this study, but group velocity gave a better result than the phase velocity.This is because of the interstation distances which is long enough.Thus, the dispersion curve was extracted from the phase velocity group velocity.The inversion of the phase alone gives a gradient of velocity with depth and a sense of the existence of several discontinuities.The results gave the final velocity model that indicates the sedimentary thickness is 12km and the continental crust thickness is 34km.The sedimentary cover has a sharp discontinuity at 4km.Based on the stratigraphic column of the Mesopotamia plain by Al-Ameri et al. ( 2011), the second discontinuity at 4.0km may represent the unconformity between the Nahr-Umr and Shuaiba formations.The ambient seismic noise results of shear velocity have been integrated with seismic cube and rock physics modeling has been utilized to reach more depths within the seismic cube.The results show a good match between seismic and ANT analysis, where Zubair and Mishrif formations reveal low physical property values in Kumaite and Dhafriyah oil fields respectively.Also, there are two potential zones in each field could contain a productive layer.

Fig. 1 .
Fig.1.The broadband seismic station location of the MP network that are used in this study.With study area location in red rectangle.

Fig. 2 .
Fig. 2.Problem geometry for the incoming plane wave at θ angle to an inter-station axis.The c is velocity, and r the distance between the seismic stations.The difference in the path from the source between stations 1 and 2 is r cos θ.

Fig. 3 .
Fig. 3. Group velocity where A, B and C shows how to use the MFT programmed from the CPS package.The color blue denotes a low amplitude, whereas the color red denotes a high amplitude.The white dots on the curve represent the best-measured group velocities at station pairings.The period span is 1-6sec.There is a label on each figure that represents the type and path component.

Fig. 4 .
Fig. 4. Rayleigh wave phase velocity dispersion curve where A, B and C show curve derived from interstation recordings.The white dots on the curve represent the best-measured group velocities at station pairings.The period span is 1-6sec.There is a label on each figure that represents the type and path component.

Fig. 5 .
Fig. 5. Comparison of predicted group velocity dispersion curves with the observed curves for three seismic stations routes (A, B, C).

Fig. 7
illustrates monthly stacked inter-station cross-correlation components functions (Z-Z, R-R and T-T) for a station pairings AMR-ANB, AMR-KAR, and KAR-ANB in the same frequency band.The stacking cross-correlations for station pairings are visible in the trace.

Fig. 6 .
Fig. 6.Stacking of seismic ambient noise cross-correlations collected at seismic broadband stations.For Rayleigh wave velocities, the red dotted lines indicate all components T, R, and Z group arrival times

Fig. 7 .
Fig. 7. A: Illustrates a distance between inter-stations.The traces fit a stacked cross-correlation and density of ANT sources in the pathways between station pairs quite well.B: frequency range utilized in cross-correlation as a function of all components Z, R, and T for all station pairings.C: Stacking of seismic ambient noise cross-correlations recorded at seismic broadband stations.

Fig. 8 .
Fig. 8. Model of the crustal structure of the study area obtained from the Rayleigh wave group velocity by using surf96 program.The red line represents the initial velocity model and the blue line represents the final model.On each figure, the station-station path is labeled

Table 1 .
Final outputs of ANT analysis