Métodos potenciales de prospección, FCAG, 2024.
Implementamos un algoritmo sencillo para determinar el perfil en profundidad de una cuenca a partir de observaciones de gravedad.
Bott (1960) describió un método para estimar la forma transversal de una cuenca sedimentaria. La cuenca se asume extendida infinitamente en una dirección y con un contraste de densidad $\Delta\rho$ uniforme con respecto a las rocas vecinas. Geométricamente, la cuenca se divide en $N$ bloques rectangulares infinitamente extendidos paralelos a la cuenca, de espesor (thickness) $t_j$, $j=1,\ldots,N$.
Se consideran $N$ puntos del campo observado $g_j$, $j=1,\ldots,N$, a lo largo de un perfil perpendicular a la cuenca. Se asume que cada medición $g_j$ en la estación de coordenada $x_j$ se encuentra centrada sobre un bloque de la cuenca a determinar.
Inicialmente, se supone cada bloque $j$ como una placa plana infinita en la dimensión horizontal. La siguiente ecuación provee el espesor de una placa infinita basada en una medición de gravedad $g_j$,
con $j=1,2,\ldots,N$.
Sea $\mathbf{x}$ el vector columna con las estaciones de medición $x_j$, $\mathbf{T}$ el vector columna con los espesores de cuenca $t_j$, $\mathbf{g}^o$ el vector columna con los valores de gravedad observada $g_j$ en $x_j$ y $\mathbf{g}^c$ el vector columna con los valores de gravedad dada en cada $x_j$ por un modelo directo que utiliza el perfil $\mathbf{T}$ y un contraste de densidad $\Delta\rho$ dado por el intérprete: $\mathbf{g}^c(\mathbf{x},\mathbf{T},\Delta\rho)$.
Para decidir la convergencia, se observan el ajuste entre las curvas de los datos observados y calculados, el perfil de la cuenca y la norma del residuo.
Para obtener $\mathbf{g}^c(\mathbf{x},\mathbf{T},\Delta\rho)$ en esta práctica se calcula la componente vertical de la anomalía de la gravedad según el método de Talwani-Worzel-Landisman analizado en la notebook del método de Hubbert.
Escribimos el código del método iterativo de Bott.
def bott(x,go,drho,maxiter=5):
G = 6.674e-11 # [m3 / kg s2]
convert = 1e5 # [m/s2] to [mGal]
cte = (2*np.pi*G*drho) * convert # [mGal]
T = np.copy(go) / cte
for niter in np.arange(maxiter):
gc = model(x,T,drho) # modelo directo utilizando el T obtenido
residue = go - g_c
T = T + residue / cte # recursión
return T, np.sum(np.square(residue))
Debemos además utilizar un código para el modelo directo con el fin de calcular las anomalias verticales de gravedad para los espesores dados.
def model(x,T,drho):
...
N = len(x)
gz = np.zeros(N)
# dx = (dx[0],dx[0], ..., d[N], d[N])
dx = np.diff(x)/2
dx = np.append(dx[0],dx) # el primer prisma va de x[0]-dx[-1] a x[0]+dx[1]. Pedimos dx[-1] = dx[0]
dx = np.append(dx,dx[-1]) # el último prisma va de x[N]-dx[N] a x[N]+dx[N+1]. Pedimos dx[N+1] = dx[N]
for i in np.arange(N): # aporte de cada prisma
x1, x2, z1, z2 = x[i]-dx[i], x[i]+dx[i+1], 0., T[i] # el prisma se extiende desde x1 a x2 y de z=0 a z=T[i]
gz = gz + gz_hubbert(x,drho,[x1,x2,z1,z2]) # para cada x_i en x, sumamos el aporte de todas la capas.
return gz
Utilizamos el modelo de Hubbert (Talwani-Worsel-Landisman) gz_hubbert() para calcular la atracción en cada estación $x_i$ debida a cada prisma y luego sumamos para obtener la superposición.
Para aplicar el método utilizamos datos que provienen de la cuenca New Red Sandstone en el sur de Escocia (Bott, 1959). Los datos proveen las posiciones de las estaciones de medición, $\mathbf{x}$, y la anomalía observada, $\mathbf{g}^o(\mathbf{x})$.
data = np.loadtxt('datos-bott.txt') # las lineas con '#' no son leidas
x , g_o = data[:,0], data[:,1] # estaciones [m] | gz observado [mGal]
print("Número de observaciones:",len(x))
Número de observaciones: 17
Graficamos $\mathbf{g}^o$ vs. $\mathbf{x}$. Notamos que los puntos de estación no están equiespaciados.
Utilizamos los siguientes valores de contraste de densidad drho y número máximo de iteraciones maxiter:
drho = -400. # [kg/m3]
maxiter = 3
Aplicamos el método:
T, error = bott(x, g_o, drho=drho, maxiter=maxiter) # espesores obtenidos y norma del error de ajuste
Y calculamos los valores de gravedad en los espesores obtenidos para visualizar:
gc = model(x,T,drho) # modelo directo utilizando los espesores en "T".
Visualizamos los resultados del método iterativo:
Si se aumenta el número de iteraciones, el ajuste con el dato observado es cada vez menor (overfitting). Sin embargo, el perfil de la cuenca obtenido puede no ser aceptable.
Podríamos justificar el perfil de la cuenca obtenida con un número menor de iteraciones en base al principio de parsinomia, ya que el perfil resulta ser más simple en esa situación.
Eso es todo por hoy.
Geophysical Journal International, Volume 3, Issue 1,Pages 63–67, doi:10.1111/j.1365-246X.1960.tb00065.x
for Two-Dimensional Bodies with Application to the Mendocino Submarine Fracture Zone, J. Geophys. Res., 64(1), 49-59, doi:10.1029/JZ064i001p00049.
Con esta forma de calcular los límites de los prismas no se genera superposición de bordes entre prismas consecutivos. Sin embargo, no necesariamente todas la estaciones caen en el centro de un prisma (como pide la hipótesis del método).
x = data[:,0]
dx = np.diff(x)/2
dx = np.append(dx[0],dx) # dx[-1] = dx[0]
dx = np.append(dx,dx[-1]) # dx[N+1] = dx[N]
Si x tiene longitud N, dx tiene longitud N+2.
x[i]-dx[i] a x[i]+dx[i+1], entonces no hay superposición de bordes entre prismas. Pero no necesariamente las estaciones caen en el centro de cada prisma.x[i]-dx[i] a x[i]+dx[i], entonces las estaciones caen siempre en el centro de un prisma pero hay superposición de bordes entre prismas.
Una forma de hacer un gráfico del perfil de la cuenca es la siguiente:
# Gráfico simple de espesores vs. estaciones
dx = np.diff(x)/2
dx = np.append(dx[0],dx) # dx[-1] = dx[0]
dx = np.append(dx,dx[-1]) # dx[N+1] = dx[N]
X, dX, TT = x * mtokm, dx*mtokm, T * mtokm
plt.figure(figsize=(8,5))
for (i,xi) in enumerate(X):
plt.scatter(xi,0,color="gray",marker="v",s=100) # estaciones
plt.scatter([xi-dX[i],xi+dX[i+1]],[0,0],color="limegreen",marker="|",s=100,alpha=0.5) # bordes entre prismas
plt.fill_between([xi-dX[i],xi+dX[i+1]],0,-TT[i],color="goldenrod",alpha=0.5,zorder=0) # relleno
plt.ylim([-1.2*np.max(TT),np.mean(TT)])
plt.xlabel("x [km]")
plt.ylabel("z [km]")
plt.show()