<jupyter_start><jupyter_text>ARC-seismic | demo Anisotropic Reflection Coefficient seismic modelling software.***  1. IntroductionARC-seismic uses the Zoeppritz equations to compute exact reflection coefficients for anisotropic elastic media across all angles of incidence (0$^{\circ}-$90$^{\circ}$), including post-critical angles. It consists of 3 individual programs:1. **_zoeppritz.py_** - Comprehensive code for modelling reflection and transmission coefficients for isotropic and anisotropic elastic media in three dimensions.2. **_montecarlo.py_** - A code that uses zoeppritz.py to conduct Monte Carlo simulations of reflection coefficients.3. **_kirchhoff.py_** - Kirchhoff program that uses zoeppritz.py to generate synthetic seismograms of the reflected finite frequency wave-field. The zoeppritz.py code has been validated using exact plane-wave reflection coefficient solutions, and the only condition it requires is that both media possess at least **_monoclinic symmetry_**, defined by an elastic tensor of the form:$$\mathbf{C}_{ij} = \begin{bmatrix}C_{11} & C_{12} & C_{13} &  &  & C_{16}\\ C_{12} & C_{22} & C_{23} &  &  & C_{26}\\ C_{13} & C_{23} & C_{33} &  &  & C_{36}\\  &  &  & C_{44} & C_{45} & \\  &  &  & C_{45} & C_{55} & \\ C_{16} & C_{26} & C_{36} &  &  & C_{66}\end{bmatrix}$$*** 2. Background & MotivationARC-seismic was initially developed to use the anisotropic Zoeppritz equations to investigate a seismic data set, the Buzzard Field, where acoustic Full-Waveform Inversion (FWI) fails (see Figure 1). I was initially surprised when I couldn't find a program online that accurately solved the anisotropic Zoeppritz equations for pre- and post-critical incidence angles. It turns out that developing an accurate and robust anisotropic Zoeppritz modelling program is a challenging task$^{\textbf{(2)}}$, requiring the solution to multiple eigenvalue problems and the correct handling of complex numbers. Using ARC-seismic, it was discovered that the failure of acoustic FWI on the Buzzard Field was a consequence of post-critical elastic effects originating from an unconformable top-chalk interface, which resulted in almost zero post-critical energy being reflected. It was also discovered that post-critical reflection coefficients are highly sensitive to weak/moderate anisotropy within the reflecting medium of an interface.For a detailed case study example, please refer to my thesis: *'The Limitations of Acoustic Full-Waveform Inversion in an Elastic World'*.   Figure 1. A) Starting reverse time migrated p-wave velocity model. The high chalk velocities have been verified using common-image gathers.  B) Acoustic FWI recovered p-wave velocity model. Acoustic FWI clearly fails to recover the high p-wave velocities of the chalk in the starting model.*** 3. TheoryThe Zoeppritz equations$^{\textbf{(1)}}$ govern the proportion of seismic energy that is reflected and transmitted from an interface between two homogenous elastic media (see Figure 2). They are formulated by assuming the solution to the wave equation: $$\left (  \rho\ \delta_{ij}\ \frac{\partial^2 }{\partial t^2}\ -\ \mathbf{C}_{ikjm}\ \frac{\partial^2 }{\partial x_{k}\partial x_{m}} \right )u_{j}=0,\ \ \ i=1,2,3,$$takes the form of an elastic plane wave:$$u_{i} = q_{i}\ exp\left [ i\omega\left ( t\ -\ s_{k}\ x_{k} \right ) \right ],\ \ \ i=1,2,3,$$where $\rho$ is the density, $\mathbf{\delta}_{ij}$ is the Kronecker delta function, $\mathbf{C}_{ikjm}$ is the fourth order elastic tensor, $u_{ij}$ is the displacement, $q_{ij}$ is the polarization vector, $s_{km}$ is the slowness vector and $\omega$ is the angular frequency. This yields the following equation: $$\rho\ q_{i} = \mathbf{C}_{ikjm}\ s_{k}\ s_{m}\ q_{j},\ \ \ i=1,2,3,$$$\mathbf{C}_{ikjm}$ governs the three-dimensional relationship between stress and strain in a medium. It was condensed into a more convenient second order 6$\times$6 matrix, $\mathbf{C}_{ij}$, by exploiting its symmetrical properties. The previous equation can then be rearranged and expressed as:$$\left ( \mathbf{\Gamma}-\rho \mathbf{I}\right )\mathbf{q} = 0$$where $\mathbf{\Gamma}$ is the 3$\times$3 Christoffel matrix, defined for media with monoclinic symmetry as follows:$$\begin{matrix} \Gamma_{11}=C_{11}s_{1}^{2}+C_{66}s_{2}^{2}+C_{55}s_{3}^{2}+2C_{16}s_{1}s_{2}\\ \Gamma_{22}=C_{66}s_{1}^{2}+C_{22}s_{2}^{2}+C_{44}s_{3}^{2}+2C_{26}s_{1}s_{2}\\ \Gamma_{33}=C_{55}s_{1}^{2}+C_{44}s_{2}^{2}+C_{33}s_{3}^{2}+2C_{45}s_{1}s_{2}\\ \Gamma_{12}=\Gamma_{21}=C_{16}s_{1}^{2}+C_{26}s_{2}^{2}+C_{45}s_{3}^{2}+(C_{12}+C_{66})s_{1}s_{2}\\ \Gamma_{13}=\Gamma_{31}=(C_{13}+C_{55})s_{1}s_{3}+(C_{36}+C_{45})s_{2}s_{3}\\ \Gamma_{23}=\Gamma_{32}=(C_{36}+C_{45})s_{1}s_{3}+(C_{23}+C_{44})s_{2}s_{3}.\end{matrix}$$Explicit solutions for the reflection and transmission coefficients of incident P and S-waves were then found by imposing continuity of displacement and traction across the interface$^{\textbf{(1)}}$.   Figure 2. Diagram of the system defined by the anisotropic Zoeppritz equations. $Pi$ is the incident p-wave, $Pr$ and $Pt$ are the reflected and transmitted p-waves, $qSHr$, $qSHt$, $qSVr$, $qSVt$ are the reflected and transmitted s-waves polarized in the quasi-horizontal/vertical planes, $\theta$ is the incidence angle and $\phi$ is the azimuth angle.For a comprehensive discussion of the theory and methods used, please refer to my masters thesis.*** References1. <NAME>. (1919), ‘Erdbebenwellen VIII B, Uber Reflexion and Durchgang seismischer wellen durch Unstetigkeisflachen’, *Gottinger Nachr* **1**, 66–84.2. <NAME>. & <NAME>. (1990), ‘Zoeppritz rationalized, and generalized to anisotropic media’, *The Journal of the Acoustical Society of America* **88**(S1), S46–S46.3. <NAME>. & <NAME>. (1980), *Quantitative seismology; theory and methods*, Freeman Co., San Francisco.4. <NAME>. (1986), ‘Weak elastic anisotropy’, *Geophysics* **51**(10), 1954–1966.5. <NAME>. (1943), ‘The mathematics of the physical properties of crystals’, *Bell Labs Technical Journal* **22**(1), 1–72.6. <NAME>. (2009), *Introduction to seismology*, Cambridge University Press. ***  Demo  1. Isotropic Reflection Coefficients<jupyter_code>from zoeppritz import * # import zoeppritz.py<jupyter_output><empty_output><jupyter_text>First, we need to define the isotropic physical parameters of the upper (1) and lower (2) media. E.g. $vp1$ is the p-wave velocity of the upper medium and $p2$ is the density of the lower medium. Note that all values should be given in **SI units** (e.g. $ms^{-1}$ and $kgm^{-3}$). These are the estimated parameters for the top-chalk interface from Figure 1.<jupyter_code>vp1 = 2000; vp2 = 4000 # p-wave velocities
vs1 = 400; vs2 = 2150 # s-wave velocities
p1 = 2000; p2 = 2600 # densities<jupyter_output><empty_output><jupyter_text>To generate a scattering matrix$^{\textbf{(4)}}$, call the **isotropic_zoeppritz()** function, passing in the angle of incidence (i_angle) in degrees. The first element of the first row of the scattering matrix is the complex valued p-wave reflection coefficient ($Rpp$). As shown, the magnitude and phase can be easily extracted.<jupyter_code>Rpp = isotropic_zoeppritz(vp1, vp2, vs1, vs2, p1, p2, i_angle=35)[0][0] # p-wave reflection coefficient

m = abs(Rpp) # magnitude
p = np.degrees(cm.phase(Rpp)) # phase

print(f'magnitude = {m}, phase = {p}')<jupyter_output>magnitude = 0.21887645438940628, phase = -87.20829659064667<jupyter_text>In order to make a plot of magnitude and phase vs angle of incidence, call the **isotropic_plot()** function. The solid line represents the magnitude and the dotted line represents the phase shift.<jupyter_code>isotropic_plot(vp1, vp2, vs1, vs2, p1, p2) # angle plot<jupyter_output><empty_output><jupyter_text>Assuming we know the depth to the interface, we can easily convert incidence angle to offset by specifying the depth (d) in metres. The maximum offset can be set by using the maxoff parameter.<jupyter_code>isotropic_plot(vp1, vp2, vs1, vs2, p1, p2, d=1000, maxoff=5000) # offset plot<jupyter_output><empty_output><jupyter_text>The above plot is an elastic model of reflection coefficients vs offset. If we want to invetigate whether acoustic FWI will be effective, as discussed in section 2, we can easily generate an acoustic model by setting the upper and lower s-wave velocities to zero (note that we have to use very small numbers instead of zero to ensure that matrix P is non-singular and invertible).<jupyter_code>isotropic_plot(vp1, vp2, 1e-10, 1e-10, p1, p2, d=1000, maxoff=5000) # acoustic offset plot<jupyter_output><empty_output><jupyter_text>As demonstrated above, there is a huge disparity between the elastic and acoustic post-critical reflection coefficients. This means that acoustic FWI will likely fail. Elastic FWI, or methods for mitigating elastic effects, are required in order to proceed with FWI.  ***  2. Anisotropic Reflection Coefficients To model media with vertical transverse isotropy (VTI), we can specify Thomsen's anisotropy parameters ($\epsilon$, $\delta$, $\gamma$)$^{\textbf{(4)}}$ and generate elastic tensors for the upper and lower media using the **thomsen_c()** function.<jupyter_code>e1 = 0.1; e2 = 0.2; d1 = 0.05; d2 = 0.1; g1 = 1e-10; g2 = 1e-10 # define some random but realistic Thomsen parameters
C1 = thomsen_c(vp1, vs1, p1, e1, d1, g1) # generate elastic tensors
C2 = thomsen_c(vp2, vs2, p2, e2, d2, g2)<jupyter_output><empty_output><jupyter_text>Similarly to before, we can calculate an individual reflection coefficient using **anisotropic_zoeppritz()**, with the addition of the azimuth angle parameter (a_angle) in degrees. The first element of the first row of the first matrix that is returned is the complex valued p-wave reflection coefficient. To generate a plot of magnitude and phase vs incidence angle, call the **anisotropic_plot()** function. Just like before, to convert incidence angle to offset, simply specify the depth (d) in metres. If the results are unstable, try increasing the pre-whitening factor (p_white) to stabilise the solution.<jupyter_code>Rpp = anisotropic_zoeppritz(C1, C2, p1, p2, i_angle=10, a_angle=0)[0][0][0] # p-wave reflection coefficient

m = abs(Rpp) # magnitude
p = np.degrees(cm.phase(Rpp)) # phase

print(f'magnitude = {m}, phase = {p}')
anisotropic_plot(C1, C2, p1, p2, a_angle=0, p_white=1e-7) # angle plot<jupyter_output><empty_output><jupyter_text>To model horizontal transverse isotropy (HTI), without losing the convenience of Thomsen's anisotropy parameter, we can apply a 90$^{\circ}$ bond transformation to the elastic tensors$^{\textbf{(5)}}$. This rotates the tensors from VTI to HTI, whilst retaining the same magnitude of anisotropy. Note that for HTI media, the reflection coefficients depend on the azimuth angle (a_angle).<jupyter_code>anisotropic_plot(bond_transformation(C1,90), bond_transformation(C2,90), p1, p2, 
                                                        a_angle=70, p_white=1e-7) # hti plot<jupyter_output><empty_output><jupyter_text>***  3. Monte Carlo Simulations Monte Carlo simulations can be used to investigate the full range of possible reflection coefficients from a geological interface where the precise lithology and physical parameters are unknown or poorly constrained. Additionally, if you have high quality seismic data available, Monte Carlo simulations can provide constraints on the range of physical parameters that could produce reflections with the amplitude behaviour that is observed in the data. They are also a very effective tool for sensitivity analyses. For instance, to determine the sensitivity of reflection coefficients to small changes in the density, velocity, or anisotropy of the system.<jupyter_code>from montecarlo import * # import montecarlo.py<jupyter_output><empty_output><jupyter_text>3.1 Isotropic Simulation To demonstrate an isotropic Monte Carlo simulation, I am going to use the same example as before; investigating the top-chalk interface from Figure 1. The estimated parameters that I used to generate the isotropic reflection coefficient profiles were very poorly constrained. In this case, it makes more sense to estimate the minimum and maximum possible values for each parameter, and then fully explore how the reflection coefficients behave in this possible parameter space using Monte Carlo simulations. In order to conduct an isotropic Monte Carlo simulation, we must specify the **minimum** and **maximum** values for each parameter, and the **number of samples**.<jupyter_code>vp1_min = 2000; vp1_max = 2000 # upper and lower p-wave velocity range 
vp2_min = 4000; vp2_max = 4000
vs1_min = 200; vs1_max = 600 # upper and lower s-wave velocity range
vs2_min = 1800; vs2_max = 2500
p1_min = 1800; p1_max = 2200 # upper and lower density range
p2_min = 2400; p2_max = 2800

num_samples = 300 # number of samples used for the simulation<jupyter_output><empty_output><jupyter_text>Then, call the **isotropic_monte_carlo()** function. This creates a random uniform distribution of samples for each parameter, which are then used to generate a reflection coefficient profile for each sample. The **isotropic_monte_carlo()** function returns an array of p-wave reflection coefficient profiles of length num_samples and the mean profile (generated using the mean value of each parameter).<jupyter_code>SIM_isotropic, mean = isotropic_monte_carlo(vp1_min, vp1_max, vp2_min, vp2_max, \
                                            vs1_min, vs1_max, vs2_min, vs2_max, \
                                            p1_min, p1_max, p2_min, p2_max, num_samples)<jupyter_output>100%|██████████| 300/300 [00:05<00:00, 53.51it/s]<jupyter_text>To visualise the result, use the **monte_carlo_plot()** function. The blue lines are the simulated reflection coefficient profiles and the red line is the mean profile. Probability density functions are calculated and plotted at incidence angles of 15$^{\circ}$, 45$^{\circ}$ and 75$^{\circ}$.<jupyter_code>monte_carlo_plot(SIM_isotropic, mean, pdf=True)<jupyter_output><empty_output><jupyter_text>3.2 Anisotropic Simulation I am going to demonstrate how to implement an anisotropic Monte Carlo simulation by showing you how to conduct a sensitivity analysis of p-wave reflection coefficients to weak/moderate amounts of VTI anisotropy. A VTI Monte Carlo simulation is conducted using the same general approach. The only difference is that we need to define the azimuth angle, the minimum and maximum Thomsen parameters, and use the **anisotropic_monte_carlo()** function. Note that this will take longer than the isotropic simulation due to the added complexities and computational cost of solving the anisotropic Zoeppritz equations!<jupyter_code>vp1_min = 2000; vp1_max = 2000 # hold isotropic parameters constant
vp2_min = 4000; vp2_max = 4000
vs1_min = 600; vs1_max = 600
vs2_min = 1800; vs2_max = 1800
p1_min = 2200; p1_max = 2200
p2_min = 2400; p2_max = 2400
e1_min = -0.02; e1_max = 0.2 # upper and lower epsilon ranges
e2_min = -0.02; e2_max = 0.2 
d1_min = -0.1; d1_max = 0.2 # upper and lower delta ranges
d2_min = -0.1; d2_max = 0.2 
g1_min = 1e-10; g1_max = 1e-10 # upper and lower gamma ranges
g2_min = 1e-10; g2_max = 1e-10 

a_angle = 0 # azimuth angle - note that this does not affect the results for VTI anisotropy

num_samples = 300 # number of samples used for the simulation
SIM_vti, mean = anisotropic_monte_carlo(vp1_min, vp1_max, vp2_min, vp2_max, \
                                        vs1_min, vs1_max, vs2_min, vs2_max, \
                                        p1_min, p1_max, p2_min, p2_max, \
                                        e1_min, e1_max, e2_min, e2_max, \
                                        d1_min, d1_max, d2_min, d2_max, \
                                        g1_min, g1_max, g2_min, g2_max, \
                                        a_angle, num_samples)
monte_carlo_plot(SIM_vti, mean, pdf=True)<jupyter_output><empty_output><jupyter_text>***  4. Kirchhoff Synthetic Seismograms Zoeppritz equations assume seismic energy travels in rays with infinite frequency$^{\textbf{(1)}}$. However, seismic waves have finite frequencies and accounting for this will effectively average out the reflection coefficient profile over the Fresnel zone$^{\textbf{(3)}}$. The kirchhoff.py program was developed in order to account for finite frequencies and geometrical spreading by modelling the reflected wave-field. The reflected wave-field for a specific receiver, $\phi_{R}$, is calculated using a Kirchhoff integral of the form: $$\phi_{R}=\frac{1}{4\pi c}\int_{S}\delta\ (t-\frac{r+r_{0}}{c})\ \frac{R_{pp}(\theta_{0})}{r r_{0}}\ (cos\theta_{0}+cos\theta)\ dS\ *\ \frac{\partial w(t,f)}{\partial t}$$where c is the p-wave velocity of the incident medium, $\delta$ is the Dirac delta function, $R_{pp}(\theta_{0})$ is the plane-wave p-wave reflection coefficient at incidence angle $\theta_{0}$, $r_{0}$ and $r$ are the source to surface grid-point and surface grid-point to receiver distances respectively, $\theta$ is the angle between the scattered ray and the surface normal, $w$ is the source time function (Ricker wavelet) and the integral is conducted over surface $S\ ^{\textbf{(6)}}$. For post-critical incidence angles, $\phi_{R}$ becomes complex, requiring an additional phase shift. Therefore, the final Kirchhoff seismogram is given by:$$\phi_{R}^{final}=Re(\phi_{R})\ +\ \textbf{H}\{Im(\phi_{R})\}$$where $\textbf{H}$ is the Hilbert transform, and the real and imaginary component are denoted by $Re$ and $Im$ respectively$^{\textbf{(6)}}$.<jupyter_code>from kirchhoff import * # import kirchhoff.py<jupyter_output><empty_output><jupyter_text>In order to model the reflected wave-field, we first need to specify the geometry of the system we wish to model.<jupyter_code>rec_min = 50 # minimum reciever distance from source
rec_max = 5000 # maximum reciever distance from source
drec = 100 # reciever spacing (every 100 meters in this case)
d = 1000 # interface depth<jupyter_output><empty_output><jupyter_text>Next, we choose some modelling parameters. Note that if the time interval spacing (dt) and grid-point spacing (ds) are too large, the synthetic seismograms will be incorrect. To run this a bit faster you can use a larger grid-point spacing (ds) value. It is recommended that you start small (using the default value or lower), and increase until the resultant seismograms are no longer the same. At this point you have reached the maximum spacing allowed for the system you have defined.<jupyter_code>w = 500 # width of 3D Kirchhoff integral
ds = 50 # grid-point spacing over the interface
dt = 5e-4 # time interval spacing
f = 6 # frequency in hz<jupyter_output><empty_output><jupyter_text>4.1 Isotropic Synthetic Seismograms To generate an isotropic synthetic seismic trace for each reciever in the specified geometry, call **isotropic_synthetic()**, which will return an array of traces (one for each reciever) and a time array. I am using the same isotropic physical parameters (defined in section 3) for modelling the top-chalk interface in Figure 1.<jupyter_code>traces_elastic, time = isotropic_synthetic(rec_min, rec_max, drec, 
                             vp1, vp2, vs1, vs2, p1, p2, d, w, f, ds=ds, dt=dt) # isotropic elastic traces<jupyter_output>100%|██████████| 50/50 [00:07<00:00,  6.95it/s]<jupyter_text>We can plot these traces by passing them into the **plot_synthetic()** function along with the time array. For better visualisation, you can amplify the result by setting the scale_fac parameter to a value greater than 1. The ymin and ymax parameters control the time axis limits.<jupyter_code>scale_fac = 1
plot_synthetic(traces_elastic, time, scale_fac, ymin=0.4, ymax=1.4) # plot isotropic elastic traces<jupyter_output><empty_output><jupyter_text>Just like in section 3, we can also generate an acoustic version of the reflected wave-field by setting the s-wave velocities to (practically) zero.<jupyter_code>traces_acoustic, time = isotropic_synthetic(rec_min, rec_max, drec, 
                        vp1, vp2, 1e-10, 1e-10, p1, p2, d, w, f, ds=ds, dt=dt) # isotropic acoustic traces
plot_synthetic(traces_acoustic, time, scale_fac, ymin=0.4, ymax=1.4) # plot isotropic acoustic traces<jupyter_output><empty_output><jupyter_text>We can also plot the difference between the acoustic and elastic traces to visualise the disparity between the elastic and acoustic reflected wave-fields.<jupyter_code>plot_synthetic(traces_acoustic-traces_elastic, time, scale_fac, ymin=0.4, ymax=1.4) # acoustic - elastic plot<jupyter_output><empty_output><jupyter_text>4.2 Anisotropic Synthetic Seismograms To generate anisotropic synthetic traces, we can use the same system geometry and modelling parameters. The only difference is that we need to call the **anisotropic_synthetic()** function which requres elastic tensors, only the upper p-wave velocity, and the azimuth angle. I am using the same anisotropic physical parameters that were defined in section 4.<jupyter_code>traces_anisotropic, time = anisotropic_synthetic(rec_min, rec_max, drec, 
                C1, C2, p1, p2, d, w, f, vp1, a_angle=0, ds=ds, dt=dt, p_white=1e-7) # anisotropic traces
plot_synthetic(traces_anisotropic, time, scale_fac, ymin=0.4, ymax=1.4) # plot anisotropic traces<jupyter_output><empty_output><jupyter_text>In order to visualise the minor differences between the anisotropic and elastic traces, we can increase the scale factor.<jupyter_code>scale_fac = 4  # increase scale factor to better visualise minor differences
plot_synthetic(traces_anisotropic-traces_elastic, time, scale_fac, ymin=0.4, ymax=1.4) # anisotropic - isotropic elastic<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09980187070911856, 'fraction_numerical': 0.03603188499285813, 'mean_word_length': 4.897826086956521, 'pattern_counts': {'":': 0, '<': 100, '<?xml version=': 0, '>': 98, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '408644', 'license': 'MIT', 'n_tokens_mistral': 7051, 'n_tokens_neox': 6249, 'n_words': 2555, 'path': 'ARC-seismic demo.ipynb', 'repo_name': 'GeorgeWilliamStrong/ARC-Seismic', 'size': '21741'}