{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "
\n", "
\n", "
Computational Seismology
\n", "
Seismic Wavefield of a Double-Couple Point Source
\n", "
\n", "
\n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", " \n", "\n", "\n", "

\n", "\n", "\n", "---\n", "\n", "This notebook is part of the supplementary material \n", "to [Computational Seismology: A Practical Introduction](https://global.oup.com/academic/product/computational-seismology-9780198717416?cc=de&lang=en&#), \n", "Oxford University Press, 2016.\n", "\n", "##### Authors:\n", "* David Vargas ([@dvargas](https://github.com/davofis))\n", "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic Equations\n", "\n", "The fundamental analytical solution to the problem of a double couple point source in infinite homogeneous media (Aki and Richards - 2002) is implemented in this Ipython notebook. This solution of seimic waves in an infinite homogeneous medium provide fundamentamental information used as benchmark to understand kinematic properties of seismic sources, quasi-analytical solutions to wave propagation problems, and influence of earthquakes on crustal deformation. \n", "\n", "Simulations of 3D elastic wave propagation need to be validated by the use of analytical solutions. In order to evaluate how healty a numerical solution is, one may recreate conditions for which analytical solutions exist with the aim of reproduce and compare the different results. In this sense, the fundamental solution for the double couple point source offers an alternative to achieve this quality control\n", "\n", "We which to find the displacement wavefield $\\mathbf{u}(\\mathbf{x},t)$ at some distance $\\mathbf{x}$ from a seismic moment tensor source with $M_xz = M_zx = M_0$. According to (Aki and Richards 2002), the displacement $\\mathbf{u}(\\mathbf{x},t)$ due to a double-couple point source in an infinite, homogeneous, isotropic medium is\n", "\n", "\\begin{align*}\n", "\\mathbf{u}(\\mathbf{x},t) &= \\dfrac{1}{4\\pi\\rho} \\mathbf{A}^N \\dfrac{1}{r^4} \\int_{{r}/{\\alpha}}^{{r}/{\\beta}} \\tau M_o(t-\\tau)d\\tau +\\\\\n", "&+\\dfrac{1}{4\\pi\\rho\\alpha^2}\\mathbf{A}^{IP}\\dfrac{1}{r^2} M_o(t-{r}/{\\alpha}) +\\dfrac{1}{4\\pi\\rho\\beta^2}\\mathbf{A}^{IS}\\dfrac{1}{r^2} M_o(t-{r}/{\\beta})+\\\\\n", "&+\\dfrac{1}{4\\pi\\rho\\alpha^3}\\mathbf{A}^{FP}\\dfrac{1}{r} \\dot M_o(t-{r}/{\\alpha}) +\\dfrac{1}{4\\pi\\rho\\beta^3}\\mathbf{A}^{FS}\\dfrac{1}{r} \\dot M_o(t-{r}/{\\beta})\n", "\\end{align*}\n", "\n", "where the radiation patterns $\\mathbf{A}^N$ (near-field), $\\mathbf{A}^{IP}$ (intermediate-field P wave), $\\mathbf{A}^{IS}$ (intermediate-field S wave), $\\mathbf{A}^{FP}$ (far-field P wave) and $\\mathbf{A}^{FS}$ (far-field S wave) are:\n", "\n", "\\begin{align*}\n", "\\mathbf{A}^N &= 9sin(2\\theta)cos(\\phi)\\hat{\\mathbf{r}} - 6(cos(2\\theta)cos(\\phi)\\hat{\\mathbf{\\theta}} - cos(\\theta)sin(\\phi))\\hat{\\mathbf{\\phi}}\\\\\n", "\\mathbf{A}^{IP} &= 4sin(2\\theta)cos(\\phi)\\hat{\\mathbf{r}} - 2(cos(2\\theta)cos(\\phi)\\hat{\\mathbf{\\theta}} - cos(\\theta)sin(\\phi))\\hat{\\mathbf{\\phi}}\\\\\n", "\\mathbf{A}^{IS} &= -3sin(2\\theta)cos(\\phi)\\hat{\\mathbf{r}} + 3(cos(2\\theta)cos(\\phi)\\hat{\\mathbf{\\theta}} - cos(\\theta)sin(\\phi))\\hat{\\mathbf{\\phi}}\\\\ \n", "\\mathbf{A}^{FP} &= sin(2\\theta)cos(\\phi)\\hat{\\mathbf{r}}\\\\\n", "\\mathbf{A}^{FS} &= cos(2\\theta)cos(\\phi)\\hat{\\mathbf{\\theta}} - cos(\\theta)sin(\\phi)\\hat{\\mathbf{\\phi}}\n", "\\end{align*}\n", "\n", "The parameters one have to consider include: receiver coordinates $\\mathbf{x}$, density of the medium $\\rho$, S-Wave velocity $\\beta$, p-wave velocity $\\alpha$, and the desired time dependent seismic moment function $M_o(t)$.On the other hand, integrations limits are determined by the propagation time from source to receiver for both p-waves and s-waves ie. ${r}/{\\beta}$ and ${r}/{\\alpha}$ respectively.\n", "\n", "This a solution in spherical coordinates. Since we normally measure displacements in cartesian coordinates, it is necessary to implement a change of coordinates if we want to visualize the solution in cartesian coordinates. \n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "code_folding": [ 0 ] }, "outputs": [], "source": [ "# Please run it before you start the simulation!\n", "import matplotlib.pyplot as plt\n", "from scipy.special import erf\n", "from scipy.integrate import quad\n", "from numpy import sin, cos, arccos, arctan, pi, sign, sqrt\n", "from numpy import vectorize, linspace, asarray, outer, diff, savetxt\n", "\n", "# Show the plots in the Notebook.\n", "plt.switch_backend(\"nbagg\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Coordinate transformation methods" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "code_folding": [ 0, 9 ] }, "outputs": [], "source": [ "def sph2cart(r, th, phi):\n", " '''\n", " Transform spherical coordinates to cartesian\n", " '''\n", " x = r * sin(th) * cos(phi)\n", " y = r * sin(th) * sin(phi)\n", " z = r * cos(th) \n", " return x, y, z\n", " \n", "def cart2sph(x, y, z):\n", " '''\n", " Transform cartesian coordinates to spherical\n", " '''\n", " r = sqrt(x**2 + y**2 + z**2)\n", " th = arccos(z/r)\n", " phi = arctan(y/x)\n", " return r, th, phi " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. COMPUTE AKI & RICHARDS SOLUTION" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "#%% Initialization of setup\n", "# -----------------------------------------------------------------------------\n", "x = 4000 # x receiver coordinate \n", "y = 4000 # y receiver coodinate\n", "z = 4000 # z receiver coodinate\n", "\n", "rho = 2500 # Density kg/m^3 \n", "beta = 3000 # S-wave velocity\n", "alpha = sqrt(3)*beta # p-wave velocity\n", "\n", "stf = 'gauss' # Set the desired source time function 'heaviside' , 'gauss'\n", "Trise = 0.25 # Rise time used in the source time function \n", "Mo = 4*10E16 # Scalar Moment \n", " \n", "r, th, phi = cart2sph(x, y, z) # spherical receiver coordinates \n", "\n", "tmin = r/alpha - 2*Trise # Minimum observation time \n", "tmax = r/beta + Trise + 2*Trise # Maximum observation time \n", "\n", "# SOURCE TIME FUNCTION\n", "# ----------------------------------------------------------------------------- \n", "if stf == 'heaviside':\n", " M0 = lambda t: 0.5*Mo*0.5*(sign(t) + 1)\n", "if stf == 'gauss':\n", " M0 = lambda t: Mo*(1 + erf(t/Trise))\n", "\n", "#******************************************************************************\n", "# COMPUTE AKI & RICHARDS SOLUTION\n", "#******************************************************************************\n", "# Scalar factors int the AKI & RICHARDS solution\n", "# -----------------------------------------------------------------------------\n", "CN = (1/(4 * pi * rho)) \n", "CIP = (1/(4 * pi * rho * alpha**2))\n", "CIS = (1/(4 * pi * rho * beta**2))\n", "CFP = (1/(4 * pi * rho * alpha**3))\n", "CFS = (1/(4 * pi * rho * beta**3))\n", "\n", "# Radiation patterns: near(AN), intermedia(AIP,AIS), and far(AFP,AFS) fields \n", "# -----------------------------------------------------------------------------\n", "def AN(th, phi): \n", " AN = [[9*sin(2*th)*cos(phi), -6*cos(2*th)*cos(phi), 6*cos(th)*sin(phi)]]\n", " return asarray(AN)\n", " \n", "def AIP(th, phi): \n", " AIP = [[4*sin(2*th)*cos(phi), -2*cos(2*th)*cos(phi), 2*cos(th)*sin(phi)]]\n", " return asarray(AIP)\n", " \n", "def AIS(th, phi): \n", " AIS = [-3*sin(2*th)*cos(phi), 3*cos(2*th)*cos(phi), -3*cos(th)*sin(phi)]\n", " return asarray(AIS)\n", " \n", "def AFP(th, phi): \n", " AFP = [sin(2*th)*cos(phi), 0, 0 ]\n", " return asarray(AFP)\n", " \n", "def AFS(th, phi): \n", " AFS = [0, cos(2*th)*cos(phi), -cos(th)*sin(phi)]\n", " return asarray(AFS)\n", "\n", "# Calculate integral in the right hand side of AKI & RICHARDS solution\n", "# -----------------------------------------------------------------------------\n", "integrand = lambda tau, t: tau*M0(t - tau)\n", "\n", "def integral(t):\n", " return quad(integrand, r/alpha, r/beta, args=(t))\n", "\n", "vec_integral = vectorize(integral)\n", "\n", "# Assemble the total AKI & RICHARDS solution\n", "# -----------------------------------------------------------------------------\n", "t = linspace(tmin, tmax, 1000) \n", "UN = CN * (1/r**4) * outer(AN(th, phi), vec_integral(t))\n", "UIP = CIP * (1/r**2) * outer(AIP(th, phi), M0(t - r/alpha)) \n", "UIS = CIS * (1/r**2) * outer(AIS(th, phi), M0(t - r/beta))\n", "\n", "t, dt = linspace(tmin, tmax, 1001, retstep=True) # diff() return N-1 size vector \n", "UFP = CFP * (1/r) * outer(AFP(th, phi), diff(M0(t - r/alpha))/dt)\n", "UFS = CFS * (1/r) * outer(AFS(th, phi), diff(M0(t - r/beta))/dt)\n", "t = linspace(tmin, tmax, 1000) \n", "\n", "U = UN + UIP + UIS + UFP + UFS\n", "\n", "Ur, Uth, Uphi = U[0,:], U[1,:], U[2,:] # spherical componets of the field u \n", "Ux, Uy, Uz = sph2cart(Ur, Uth, Uphi) # spherical to cartesian coordinates" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Plot displacement components " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "