Ejemplo sencillo del método de Hubbert

Métodos potenciales de prospección, FCAG, 2024.

$\newcommand{\vec}[1]{\mathbf{#1}}$ $\newcommand{\versor}[1]{\hat{\mathbf{#1}}}$ $\newcommand{den}{\Delta\sigma}$ $\newcommand{G}{\mathcal{G}}$ $\newcommand{M}{\dfrac{4}{3}\pi\den R^3}$ $\newcommand{cte}{\dfrac{4}{3}\pi\G\den R^3}$ En esta ejercicio implementamos el método de la integral de línea de Hubbert,

$$g_z(x=x_E) = 2\G\den\oint_{\Gamma} z(\theta)\,\mathrm{d}\theta = 2\G\den\sum_{i=1}^{4} \int_{\Gamma_i} z(\theta)\,\mathrm{d}\theta,$$

para calcular la respuesta gravimétrica, $g_z$, en puntos de estación $x_E$ para un cuerpo 2D de sección rectangular. En esta expresión, la integral de línea es recorrida en sentido horario y $\Gamma = \Gamma_1 \cup \Gamma_2 \cup \Gamma_3 \cup \Gamma_4$.

Cálculo de $g_z$ según el método de Hubbert (Talwani-Worsel-Landisman)

En primer lugar resolvemos y validamos las expresiones para $g_z$ dadas por el método de Hubbert para el caso de un cuerpo 2D de sección rectangular. Utilizamos los lineamientos de Talwani-Worsel y Landisman para codificar el método. Decidimos emplear una subrutina para hacer más versátil la aplicación a otros escenarios y facilitar la colaboración.

La subrutina requiere las estaciones de medición $x_i$ dadas en el vector $\vec{x}$ (x), en unidades de metros, el contraste de densidad $\den$ en kg/m$^3$ (sigma) y las coordenadas de los vértices (vert) dadas por $\vec{vert}=(x_1,x_2,z_1,z_2)$ con $x_1<x_2$ y $z_1<z_2$, en unidades de metros. La salida de la subrutina son los valores de $g_{z_i}$ en unidades de mGal dadas en el vector $\vec{g_z}$ (gz).

def gz_hubbert(x,sigma,vert)

    G = 6.674e-11 # [m3 / kg s2]    
    # conversión de unidades
    convert = 1e5 # m/s2 to mGal
    cte     =  2.0 * G * sigma  # [1/s2]
    cte     = convert * cte     # [mGal] 

    x1, x2, z1, z2 = vert
    ...

    integral = []    

    for xe in x:
        if (xe < x1):
            ...
            Gamma1 = ...
            Gamma2 = ...
            Gamma3 = ...
            Gamma4 = ...

        elif(xe == x1):
            ...
        elif(xe > x1 and xe < x2):  
            ...
        elif(xe == x2)
            ...
        elif (xe > x2):
            ...

        integral.append( Gamma1 + Gamma2 + Gamma3 + Gamma4 ) # [m]

    gz = cte * np.array(integral) # mGal
return  gz

Aplicación

Buscamos reproducir la Figura 3.17 (a) y (b) de libro Gravity and Magnetic Exploration, W. J. Hinze et al. (página 60). Para ello definimos el modelo y procedemos a la evaluación de $g_z$. El modelo está parametrizado por el valor del contraste de densidad $\den$ y las coordenadas de los vértices $x_1,x_2,z_1,z_2$.

...
kmtom = 1e3
...
# Figura 3.17a
sigma = 200. # [kg/m3]
x1    = +95.0  * kmtom
x2    = +105.0 * kmtom
z1    = +5.0   * kmtom
z2    = +15.0  * kmtom
...

Para verificar, calculamos $g_z$ para el caso de una placa infinita con la misma densidad y de espesor dado por $z_2-z_1$. La anomalía calculada numéricamente debería converger a este valor como límite.

gz para placa infinita de espesor z2-z1: 83.87 mGal:
Máximo de g_z: 26.28 mGal

Gráficos

Representamos gráficamente los resultados y analizamos si hemos logrado reproducir la figura del libro de referencia. En color visualizamos la componente vertical de la atracción gravitatoria para el cuerpo 2D. La linea gris es el valor de anomalía para una placa infinita de espesor $z_2-z_1$.

Visualización interactiva

Graficamos de forma interactiva para distintos valores del contraste de densidad $\den$ y la ubicación de los vértices.

interactive(children=(FloatSlider(value=200.0, description='$\\den$ [kg/m$3$]', max=200.0, min=-200.0, step=10…

Eso es todo por hoy.