Como é o ângulo de reflexão e refração das ondas P e S para diferentes ângulos de incidência? E como isso varia dependendo das velocidades de onda sísmica do meio? A Lei de Snell explica:
$$\dfrac{\sin \theta_1}{V_1} = \dfrac{\sin \theta_2}{V_2}$$Nessa prática, vocês deverão usar a Lei de Snell para prever e explicar o comportamento das ondas P e S ao se depararem com uma interface.
Utilizaremos as simulações de ondas da biblioteca Fatiando a Terra. Essas simulações utilizam o método de diferenças finitas para calcular soluções da equação da onda.
Esse documento é um Jupyter notebook, um documento interativo que mistura texto (como esse), código (como abaixo), e o resultado de executar o código (números, texto, figuras, videos, etc).
O notebook te fornecerá exemplos interativos que trabalham os temas abordados no questionário. Utilize esses exemplos para responder as perguntas.
As células com números ao lado, como In [1]:
, são código Python. Algumas dessas células não produzem resultado e servem de preparação para os exemplos interativos. Outras, produzem gráficos interativos. Você deve executar todas as células, uma de cada vez, mesmo as que não produzem gráficos.
Para executar uma célula, clique em cima dela e aperte Shift + Enter
. O foco (contorno verde ou cinza em torno da célula) deverá passar para a célula abaixo. Para rodá-la, aperte Shift + Enter
novamente e assim por diante. Você pode executar células de texto que não acontecerá nada.
Rode as células abaixo para carregar os módulos necessários para essa prática.
%matplotlib inline
from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as ipw
from fatiando.seismic import RickerWavelet, GaussianWavelet, FDElasticPSV, FDAcoustic2D
Vamos simular uma onda P incidindo na interface entre dois meios fluidos com velocidades aumentando com a profundidade.
shape = (200, 600)
spacing = 5
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]
vp1 = 1500
vp2 = 2000
increase_density = np.zeros(shape, dtype='float32') + 1000
increase_velocity = np.zeros(shape, dtype='float32') + vp1
increase_density[100:,:] = 1500
increase_velocity[100:,:] = vp2
Agora vamos criar o nosso simulador de ondas P com uma fonte explosiva na superfície do nosso modelo.
p_increase = FDAcoustic2D(increase_velocity, increase_density, spacing=spacing)
fonte = (0, shape[1]//4)
p_increase.add_point_source(fonte, RickerWavelet(1, 60))
Agora que temos nosso simulador pronto, rode a célcula abaixo para avançar a simulação no tempo.
p_increase.run(1500)
Essa simulação pode demorar um pouco para terminar.
p_increase.animate(every=20, embed=True, dpi=60, cutoff=0.4, fps=7)
Rode a próxima célula para explorar fotos de cada etapa da simulação, uma de cada vez. A figura também mostrará 3 raios (onda incidente, refletida e refratada) para um determinado ângulo de incidência que você poderá controlar. Os ângulos de refração e reflexão são calculados com a Lei de Snell.
def plot_with_p_rays_increasing(tempo, incidencia):
fig = plt.figure()
ax = plt.subplot(111)
p_increase.snapshot(frame=tempo, ax=ax, cutoff=0.2, cmap='Greys')
fig.set_size_inches(14, 5.5)
y_bottom = shape[0]*spacing
y_interface = 100*spacing
y_source = fonte[0]*spacing
x_source = fonte[1]*spacing
x_incidence = (np.tan(np.radians(incidencia))*(y_interface - y_source)
+ x_source)
x_reflect = 2*x_incidence - x_source
arg = (vp2/vp1)*np.sin(np.radians(incidencia))
if arg <= 1:
refract = np.arcsin(arg)
x_refract = (np.tan(refract)*(y_bottom - y_interface) + x_incidence)
ax.plot([x_incidence, x_refract], [y_interface, y_bottom], '-r', linewidth=2)
ax.plot([x_source, x_incidence], [y_source, y_interface], '-k', linewidth=2)
ax.plot([x_incidence, x_reflect], [y_interface, 0], '-b', linewidth=2)
ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
ipw.interactive(plot_with_p_rays_increasing,
tempo=ipw.IntSlider(value=0, min=0, max=p_increase.it, step=50),
incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5))
Agora vamos ver o que acontece quando a velocidade diminui com a profundidade.
vp1_diminui = 1500
vp2_diminui = 1000
decrease_density = np.zeros(shape, dtype='float32') + 1500
decrease_velocity = np.zeros(shape, dtype='float32') + vp1_diminui
decrease_density[100:,:] = 2000
decrease_velocity[100:,:] = vp2_diminui
Novamente, crie a simulação com uma fonte explosiva no topo e a avance no tempo.
p_decrease = FDAcoustic2D(decrease_velocity, decrease_density, spacing=spacing)
p_decrease.add_point_source(fonte, RickerWavelet(1, 60))
p_decrease.run(1500)
p_decrease.animate(every=20, embed=True, dpi=60, cutoff=0.4, fps=7)
Rode abaixo para gerar a mesma figura interativa da simulação anterior.
def plot_with_p_rays_decreasing(tempo, incidencia):
fig = plt.figure()
ax = plt.subplot(111)
p_decrease.snapshot(frame=tempo, ax=ax, cutoff=0.4, cmap='Greys')
fig.set_size_inches(14, 5.5)
y_bottom = shape[0]*spacing
y_interface = 100*spacing
y_source = fonte[0]*spacing
x_source = fonte[1]*spacing
x_incidence = (np.tan(np.radians(incidencia))*(y_interface - y_source)
+ x_source)
x_reflect = 2*x_incidence - x_source
arg = (vp2_diminui/vp1_diminui)*np.sin(np.radians(incidencia))
if arg <= 1:
refract = np.arcsin(arg)
x_refract = (np.tan(refract)*(y_bottom - y_interface) + x_incidence)
ax.plot([x_incidence, x_refract], [y_interface, y_bottom], '-r', linewidth=2)
ax.plot([x_source, x_incidence], [y_source, y_interface], '-k', linewidth=2)
ax.plot([x_incidence, x_reflect], [y_interface, 0], '-b', linewidth=2)
ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
ipw.interactive(plot_with_p_rays_decreasing,
tempo=ipw.IntSlider(value=0, min=0, max=p_decrease.it, step=50),
incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5))
Vamos fazer agora uma simulação que inclui tanto ondas P quanto ondas S propagando em um meio sólido. A onda P irá passar para um meio com velocidades maiores.
vp1_solido = 3000
vp2_solido = 4000
vs1_solido = 2000
vs2_solido = 3000
shape2 = (300, 600)
l = 150
twosolid_density = np.zeros(shape2, dtype='float32') + 1800
twosolid_density[l:, :] = 2100
twosolid_vs = np.zeros(shape2, dtype='float32') + vs1_solido
twosolid_vs[l:, :] = vs2_solido
twosolid_vp = np.zeros(shape2, dtype='float32') + vp1_solido
twosolid_vp[l:, :] = vp2_solido
Rode as células abaixo para criar e rodar a simulação utilizando uma fonte explosiva na superfície.
fonte_explosiva = (70, shape[1]//4)
ps_twosolid = FDElasticPSV(twosolid_vp, twosolid_vs, twosolid_density, spacing=spacing)
ps_twosolid.add_blast_source(fonte_explosiva, GaussianWavelet(1, 100))
ps_twosolid.run(1200)
ps_twosolid.animate(every=20, plottype=['vectors', 'wavefield'],
cutoff=1e-6, scale=1e7, every_particle=10,
dpi=60, fps=6, embed=True)
Novamente, use a célula abaixo para ver cada etapa da simulação separadamente e os raios correspondentes as reflexões e refrações.
def plot_with_ps_rays(tempo, incidencia):
fig = plt.figure()
ax = plt.subplot(111)
ps_twosolid.snapshot(frame=tempo, ax=ax, plottype=['wavefield'],
cutoff=1e-6, every_particle=10, cmap='Greys')
fig.set_size_inches(14, 6.5)
y_bottom = shape2[0]*spacing
y_interface = l*spacing
y_source = fonte_explosiva[0]*spacing
x_source = fonte_explosiva[1]*spacing
x_incidence = (np.tan(np.radians(incidencia))*(y_interface - y_source)
+ x_source)
ax.plot([x_source, x_incidence], [y_source, y_interface], '-k', linewidth=2)
# P reflection
x_preflect = (np.tan(np.radians(incidencia))*(y_interface - 0) + x_incidence)
ax.plot([x_incidence, x_preflect], [y_interface, 0], '-b', linewidth=2)
# S reflection
sreflect = np.arcsin((vs1_solido/vp1_solido)*np.sin(np.radians(incidencia)))
x_sreflect = (np.tan(sreflect)*(y_interface - 0) + x_incidence)
ax.plot([x_incidence, x_sreflect], [y_interface, 0], '--b', linewidth=2)
# P refraction
arg = (vp2_solido/vp1_solido)*np.sin(np.radians(incidencia))
if arg <= 1:
prefract = np.arcsin(arg)
x_prefract = (np.tan(prefract)*(y_bottom - y_interface) + x_incidence)
ax.plot([x_incidence, x_prefract], [y_interface, y_bottom], '-r', linewidth=2)
# S refraction
arg = (vs2_solido/vp1_solido)*np.sin(np.radians(incidencia))
if arg <= 1:
srefract = np.arcsin(arg)
x_srefract = (np.tan(srefract)*(y_bottom - y_interface) + x_incidence)
ax.plot([x_incidence, x_srefract], [y_interface, y_bottom], '--r', linewidth=2)
ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
ipw.interactive(plot_with_ps_rays,
tempo=ipw.IntSlider(value=0, min=0, max=ps_twosolid.it, step=50),
incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5))
Course website: https://github.com/leouieda/geofisica2
Note: This notebook is part of the course "Geofísica 2" of Geology program of the Universidade do Estado do Rio de Janeiro. All content can be freely used and adapted under the terms of the Creative Commons Attribution 4.0 International License.