Hoje vamos estudar o fenômeno da refração no ângulo crítico (ângulo de incidência que gera refração a 90°). Essa onda possui propriedades distintas que possibilitam seu uso para inferir propriedades da subsuperfície (o método da sísmica de refração). Veremos quais são essas propriedades e possíveis fatores limitantes para a aplicação do método.
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.
Sua resposta deve conter no mínimo:
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, FDAcoustic2D
Vamos simular uma onda P incidindo na interface entre dois meios gerada por uma explosão na superfície.
Dessa vez, vamos ignorar as ondas S para não complicar a vida.
shape = (250, 900)
spacing = 30
extent = [0, shape[1]*spacing, shape[0]*spacing, 0]
vp1 = 3000
vp2 = 4000
print("Vp1 = {} m/s".format(vp1))
print("Vp2 = {} m/s".format(vp2))
density = np.zeros(shape, dtype='float32') + 2000
velocity = np.zeros(shape, dtype='float32') + vp1
interface = 80
density[interface:,:] = 2300
velocity[interface:,:] = vp2
print("Interface entre os dois meios a {} m de profundidade.".format(interface*spacing))
A fonte será uma explosão na superfície com as seguintes coordenadas:
fonte = (0, shape[1]//5)
print('Coordenadas da fonte:')
print(' x = {} m'.format(fonte[1]*spacing))
print(' z = {} m'.format(fonte[0]*spacing))
Agora vamos criar o nosso simulador de ondas P com uma fonte explosiva na superfície do nosso modelo.
simul = FDAcoustic2D(velocity, density, spacing=spacing)
simul.add_point_source(fonte, RickerWavelet(1, 30))
Agora que temos nosso simulador pronto, rode a célcula abaixo para avançar a simulação no tempo. Essa simulação vai demorar um pouco para terminar.
simul.run(2000)
Rode a célula abaixo para gerar a animação.
simul.animate(every=30, embed=True, dpi=60, cutoff=2, fps=6, cmap='Greys')
Note que na simulação existem ondas refletindo nas bordas da área. Ignorem essas ondas. Elas são artefatos da simulação que são difíceis de evitar.
Rode a próxima célula para gerar uma figura interativa com:
def plot_with_rays(tempo, incidencia, x, angulo):
fig = plt.figure()
ax = plt.subplot(111)
simul.snapshot(frame=tempo, ax=ax, cutoff=2, cmap='Greys')
fig.set_size_inches(14, 5.4)
ax.set_title('Tempo em segundos: {} s'.format(tempo*simul.dt))
y_bottom = shape[0]*spacing
y_interface = interface*spacing
y_source = fonte[0]*spacing
x_source = fonte[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)
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)
tmp = np.tan(np.radians(angulo))
x_marc_topo = (y_interface - 0)*tmp + x
x_marc_base = (y_interface - y_bottom)*tmp + x
ax.plot([x_marc_base, x_marc_topo], [y_bottom, 0], '--b', linewidth=2)
ax.hlines(y_interface, 0, shape[1]*spacing, colors='grey')
ipw.interactive(plot_with_rays,
tempo=ipw.IntSlider(value=0, min=0, max=simul.it, step=50),
incidencia=ipw.FloatSlider(value=45, min=0, max=90, step=0.5),
x=ipw.FloatSlider(value=10e3, min=0, max=shape[1]*spacing, step=500),
angulo=ipw.FloatSlider(value=0, 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.