51.1. Introduction to Seismic Wave Propagation link

Seismic waves, generated by abrupt releases of energy from natural events such as earthquakes or volcanic eruptions, serve as vital probes into the Earth's interior. These waves yield indispensable data for identifying and analyzing seismic events, exploring subsurface structures in the search for resources such as oil and gas, and informing the design of structures to withstand seismic impacts. In this section, we explore the essential characteristics of seismic waves, the physics behind their propagation through the Earth's stratified layers, and several practical applications. We also present how these phenomena can be modeled and simulated using the Rust programming language, taking advantage of its robust memory safety and concurrency features.

Seismic waves can be broadly divided into two categories: body waves and surface waves.

Body waves travel through the interior of the Earth and are primarily categorized into two types. Primary waves (P-waves) are compressional in nature, making them the fastest among all seismic waves. They can traverse both solid and liquid regions, as the particles move parallel to the wave's direction, resulting in alternating regions of compression and rarefaction. Secondary waves (S-waves), on the other hand, are shear waves that travel more slowly and are restricted to solids. Their motion is perpendicular to the direction of propagation, which generates shear stresses in the medium.

Surface waves, which propagate along the Earth's exterior, are often responsible for the greatest damage during seismic events. Rayleigh waves induce a rolling motion with particles following elliptical paths in both vertical and horizontal directions, much like ocean waves. Love waves produce horizontal shearing that is perpendicular to the direction of travel, leading to pronounced lateral ground movement.

The Earth itself is structured into several distinct layers—namely, the crust, mantle, and core—each possessing unique physical properties that influence wave propagation. As seismic waves pass through these layers, their speeds, trajectories, and intensities are modified due to variations in density, elasticity, and temperature. For example, within the crust, where seismic activity is first observed, both P-waves and S-waves can be detected, though S-waves cannot travel through liquid regions. In the mantle, composed of solid rock that slowly deforms over time, the increasing density and rigidity cause the wave velocities to rise and the waves to refract. In the core, which is divided into a solid inner core and a liquid outer core, the inability of S-waves to traverse liquid results in distinct shadow zones, while P-waves are refracted as they pass through.

A particularly important aspect of seismic wave behavior is the interaction of waves with the Earth’s layered structure. For instance, when seismic waves encounter boundaries—such as the transition between the crust and mantle—a portion of the energy is reflected back while the remainder refracts into the underlying layer. These interactions, which include reflection, refraction, diffraction, and scattering, reveal critical insights about the internal composition and structure of the Earth.

Several factors govern the propagation of seismic waves. Wave velocity is determined by the physical properties of the medium, with P-waves generally traveling faster than S-waves, especially in denser and more elastic materials. Attenuation, which arises from scattering, absorption, and geometric spreading, leads to a gradual decrease in wave amplitude as the energy dissipates over distance. Additionally, the frequency of the seismic wave influences its propagation characteristics: lower-frequency waves travel farther and provide broad information about subsurface features, while higher-frequency waves, though more quickly attenuated, offer finer resolution of geological details.

In practical applications, seismic waves are invaluable in various fields. Global seismic networks utilize the arrival times of P-waves and S-waves to accurately locate earthquake epicenters and estimate their magnitudes. In the realm of subsurface exploration, artificial seismic waves are generated and their reflections off geological layers are analyzed to identify potential resource deposits. Similarly, in structural engineering, the response of buildings and infrastructure to seismic waves is studied to improve designs for better resilience against earthquakes.

Simulating seismic wave propagation, while challenging due to the complexity of the Earth's heterogeneous and anisotropic environments, can be achieved through numerical methods such as finite-difference techniques. Rust provides a powerful framework for such simulations, owing to its strong guarantees of memory safety and efficient handling of concurrent computations. In the following example, we present a simplified one-dimensional simulation of P-wave propagation. This simulation takes into account varying velocities and attenuation as waves travel through a homogeneous medium, serving as a foundation for more elaborate models that include multi-dimensional propagation and heterogeneous materials.

Below is a robust Rust implementation that illustrates the propagation of a seismic P-wave in a one-dimensional medium using finite difference methods. The code is thoroughly commented to explain each step of the simulation.

extern crate ndarray;
use ndarray::Array1;

/// Simulate the propagation of a P-wave in a 1D medium using finite difference methods.
///
/// # Arguments
/// * `time_steps` - The number of time iterations for the simulation.
/// * `grid_points` - The number of spatial grid points representing the medium.
/// * `velocity` - The constant propagation velocity of the seismic wave.
/// * `dt` - The time step size.
/// * `dx` - The spatial step size.
///
/// # Returns
/// An Array1<f64> containing the final displacement values of the wave after simulation.
fn simulate_p_wave(time_steps: usize, grid_points: usize, velocity: f64, dt: f64, dx: f64) -> Array1<f64> {
    // Initialize displacement arrays for previous, current, and next time steps.
    let mut u_prev = Array1::<f64>::zeros(grid_points); // Displacement at time step n-1
    let mut u_curr = Array1::<f64>::zeros(grid_points); // Displacement at time step n
    let mut u_next = Array1::<f64>::zeros(grid_points); // Displacement at time step n+1

    // Set an initial condition: a pulse at the center of the grid to simulate an instantaneous energy release.
    u_curr[grid_points / 2] = 1.0;

    // Iterate over the specified number of time steps.
    for _ in 0..time_steps {
        // Update the displacement at each interior grid point using the finite difference approximation of the wave equation.
        for i in 1..grid_points - 1 {
            // Calculate the second spatial derivative using central differences.
            let spatial_derivative = (u_curr[i + 1] - 2.0 * u_curr[i] + u_curr[i - 1]) / (dx * dx);
            // Update the displacement at the next time step.
            u_next[i] = 2.0 * u_curr[i] - u_prev[i] + (velocity * velocity * dt * dt) * spatial_derivative;
        }

        // Apply reflective boundary conditions by mirroring the displacement at the boundaries.
        u_next[0] = u_next[1];
        u_next[grid_points - 1] = u_next[grid_points - 2];

        // Shift the time steps: the current state becomes the previous state, and the next state becomes the current state.
        u_prev.assign(&u_curr);
        u_curr.assign(&u_next);
    }

    // Return the final displacement values after completing the simulation.
    u_curr
}

fn main() {
    // Define simulation parameters.
    let time_steps = 500;    // Total number of time iterations
    let grid_points = 100;   // Number of spatial grid points in the medium
    let velocity = 1.0;      // Propagation velocity of the seismic wave (units per time step)
    let dt = 0.01;           // Time step size (in appropriate time units)
    let dx = 0.1;            // Spatial step size (in appropriate distance units)

    // Run the P-wave simulation.
    let result = simulate_p_wave(time_steps, grid_points, velocity, dt, dx);
    // Print the final displacement values of the wave.
    println!("Final wave displacement: {:?}", result);
}

In this simulation, the spatial domain is discretized into a fixed number of grid points representing a section of the Earth's interior. A pulse is introduced at the center of the grid to simulate the initial energy release characteristic of an earthquake. The simulation then proceeds iteratively, updating the displacement at each grid point by approximating the second derivative in space and applying a finite difference scheme for time evolution. Reflective boundary conditions are enforced at the edges of the grid, ensuring that the wave energy is contained within the computational domain.

This implementation, while basic, provides a foundation for more complex seismic simulations. It can be extended to model two-dimensional or three-dimensional wave propagation, incorporate heterogeneous material properties, and simulate the effects of various types of seismic waves, including surface waves. By fine-tuning parameters such as velocity and grid spacing, one can explore a wide range of seismic phenomena and gain deeper insights into both natural seismic activity and engineered systems designed to mitigate its impact.

Understanding the propagation of seismic waves is crucial for both scientific research and practical applications. Through numerical simulations like the one demonstrated above, researchers can analyze how seismic energy travels through different layers of the Earth, predict the behavior of waves under various conditions, and ultimately enhance our ability to monitor, interpret, and respond to seismic events.

51.2. Mathematical Modeling of Seismic Waves link

Seismic wave propagation is captured by several fundamental equations that describe the behavior of waves as they travel through different media. The wave equation, the Helmholtz equation for harmonic motions, and the elastodynamic equations each provide critical insights into the mechanics of seismic phenomena. These equations enable us to analyze how seismic waves interact with both homogeneous and heterogeneous materials, deepening our understanding of the Earth's layered structure and the dynamics within.

Wave Equation\ The wave equation is the cornerstone of seismic modeling. For a scalar displacement field u, the three-dimensional wave equation is written as

$$∂²u/∂t² = v² ∇²u$$

where v represents the seismic wave velocity, $∇²$ is the Laplacian operator that encapsulates the second spatial derivatives, and t denotes time. This formulation applies to both compressional and shear waves in elastic media, with the actual speed v varying based on whether the wave is a P-wave or an S-wave.

Helmholtz Equation\ When focusing on harmonic waveforms, the time-dependent wave equation can be transformed into the Helmholtz equation by assuming a time-harmonic dependence. The Helmholtz equation is expressed as

$$∇²u + k² u = 0$$

with k denoting the wave number, which is directly related to the wave’s frequency. This equation is particularly useful in steady-state analyses, where frequency-domain methods provide detailed insight into seismic responses.

Elastodynamic Equations\ The elastodynamic equations extend the wave equation to elastic media where the displacement field is vectorial. For a vector displacement field 𝑢, the general form of the elastodynamic equation is

$$∂²𝑢/∂t² = ∇ · σ$$

where $σ$ is the stress tensor related to the strain in the material through constitutive relations, and the density ρ of the medium plays a crucial role in the dynamics. This framework encompasses both compressional and shear wave propagation and accounts for the internal forces arising within the Earth's elastic layers.

Boundary and Initial Conditions\ In solving the wave equation, the specification of initial and boundary conditions is essential. Initial conditions define the state of the system at the beginning of the simulation, such as the initial displacement or velocity field generated by a sudden event. Boundary conditions, on the other hand, determine how the wave interacts with the limits of the modeled domain. For example, Dirichlet boundary conditions fix the displacement or velocity to a predetermined value, while Neumann boundary conditions specify the gradient of the displacement, ensuring reflective boundaries that simulate a closed system.

In some seismic applications, particularly those involving fluid-filled environments, an acoustic approximation is applied. In this approximation only P-waves are considered because the shear modulus is assumed to be zero, thus simplifying the model. This approach is valuable in studies of bodies of water or fluid-saturated rock formations.

In more realistic geological scenarios, seismic wave propagation becomes increasingly complex when dealing with anisotropic, porous, or heterogeneous structures. Anisotropic materials exhibit direction-dependent wave velocities, while porosity alters the medium’s compressibility, affecting P-wave speeds. Heterogeneity introduces spatial variations in material properties like density and elasticity, which in turn lead to scattering, diffraction, and other intricate wave phenomena.

Let us now illustrate these concepts with a practical example by implementing the wave equation using the finite difference method (FDM) in Rust. In the code below, we approximate the second spatial derivative and perform time-stepping on a discretized grid to simulate the propagation of a seismic wave in a homogeneous medium.

extern crate ndarray;
use ndarray::Array1;

/// Simulates the propagation of a seismic wave in a 1D homogeneous medium using finite difference methods.
///
/// # Arguments
/// * `time_steps` - The number of time iterations for the simulation.
/// * `grid_points` - The number of spatial points representing the medium.
/// * `velocity` - The constant wave propagation velocity.
/// * `dt` - The time step size.
/// * `dx` - The spatial step size.
///
/// # Returns
/// A one-dimensional array of f64 values representing the final displacement field.
fn simulate_wave_1d(time_steps: usize, grid_points: usize, velocity: f64, dt: f64, dx: f64) -> Array1<f64> {
    // Create arrays to hold the displacement values at the previous, current, and next time steps.
    let mut u_prev = Array1::<f64>::zeros(grid_points); // Displacement at time step n-1
    let mut u_curr = Array1::<f64>::zeros(grid_points); // Displacement at time step n
    let mut u_next = Array1::<f64>::zeros(grid_points); // Displacement at time step n+1

    // Set the initial condition: a pulse at the center of the grid representing an initial disturbance.
    u_curr[grid_points / 2] = 1.0;

    // Perform the time-stepping for the specified number of iterations.
    for _ in 0..time_steps {
        // Update displacement for interior grid points using a central difference approximation for the second spatial derivative.
        for i in 1..grid_points - 1 {
            let spatial_derivative = (u_curr[i + 1] - 2.0 * u_curr[i] + u_curr[i - 1]) / (dx * dx);
            u_next[i] = 2.0 * u_curr[i] - u_prev[i] + (velocity * velocity * dt * dt) * spatial_derivative;
        }

        // Apply reflective boundary conditions by mirroring the adjacent grid point values.
        u_next[0] = u_next[1];
        u_next[grid_points - 1] = u_next[grid_points - 2];

        // Shift the arrays for the next iteration: the current becomes previous and next becomes current.
        u_prev.assign(&u_curr);
        u_curr.assign(&u_next);
    }

    // Return the final displacement field after completing the simulation.
    u_curr
}

fn main() {
    // Simulation parameters for the homogeneous medium.
    let time_steps = 500;     // Number of time iterations
    let grid_points = 100;    // Number of spatial grid points
    let velocity = 1.0;       // Constant propagation velocity
    let dt = 0.01;            // Time step size
    let dx = 0.1;             // Spatial step size

    // Execute the simulation and print the final displacement.
    let result = simulate_wave_1d(time_steps, grid_points, velocity, dt, dx);
    println!("Final wave displacement: {:?}", result);
}

This implementation discretizes a one-dimensional medium into a fixed number of grid points. A pulse, introduced at the center, serves as the initial disturbance. The simulation updates the displacement values at each grid point iteratively by approximating the second derivative in space and applying the finite difference formula for time evolution. Reflective boundary conditions are enforced by mirroring adjacent values at the domain boundaries.

The complexity of seismic wave propagation increases in heterogeneous media, where the wave velocity can vary spatially. The following code extends the previous implementation by introducing a spatially varying velocity profile. In this scenario, an array of velocities is defined to represent differing material properties along the grid.

extern crate ndarray;
use ndarray::Array1;

/// Simulates seismic wave propagation in a 1D heterogeneous medium using finite difference methods.
///
/// # Arguments
/// * `time_steps` - The number of time iterations for the simulation.
/// * `grid_points` - The number of spatial points in the simulation.
/// * `velocities` - A one-dimensional array of wave velocities at each grid point.
/// * `dt` - The time step size.
/// * `dx` - The spatial step size.
///
/// # Returns
/// A one-dimensional array of f64 values representing the final displacement field.
fn simulate_wave_heterogeneous(time_steps: usize, grid_points: usize, velocities: &Array1<f64>, dt: f64, dx: f64) -> Array1<f64> {
    // Initialize arrays for displacement at different time steps.
    let mut u_prev = Array1::<f64>::zeros(grid_points); // Displacement at time step n-1
    let mut u_curr = Array1::<f64>::zeros(grid_points); // Displacement at time step n
    let mut u_next = Array1::<f64>::zeros(grid_points); // Displacement at time step n+1

    // Set the initial condition with a pulse in the center.
    u_curr[grid_points / 2] = 1.0;

    // Perform the simulation over the defined number of time steps.
    for _ in 0..time_steps {
        // Update the displacement using spatially varying velocities.
        for i in 1..grid_points - 1 {
            // Extract the local velocity from the provided velocities array.
            let local_velocity = velocities[i];
            let spatial_derivative = (u_curr[i + 1] - 2.0 * u_curr[i] + u_curr[i - 1]) / (dx * dx);
            u_next[i] = 2.0 * u_curr[i] - u_prev[i] + (local_velocity * local_velocity * dt * dt) * spatial_derivative;
        }

        // Apply reflective boundary conditions.
        u_next[0] = u_next[1];
        u_next[grid_points - 1] = u_next[grid_points - 2];

        // Update arrays for the next iteration.
        u_prev.assign(&u_curr);
        u_curr.assign(&u_next);
    }

    // Return the final displacement field.
    u_curr
}

fn main() {
    // Define simulation parameters.
    let time_steps = 500;     // Number of time iterations
    let grid_points = 100;    // Number of spatial grid points
    let dt = 0.01;            // Time step size
    let dx = 0.1;             // Spatial step size

    // Define a heterogeneous velocity profile.
    // For example, the first half of the grid has a velocity of 1.0 and the second half has a velocity of 2.0.
    let velocities = Array1::from(
        vec![1.0; grid_points / 2]
            .into_iter()
            .chain(vec![2.0; grid_points - grid_points / 2].into_iter())
            .collect::<Vec<_>>()
    );

    // Run the heterogeneous simulation and output the final displacement field.
    let result = simulate_wave_heterogeneous(time_steps, grid_points, &velocities, dt, dx);
    println!("Final wave displacement in heterogeneous media: {:?}", result);
}

In this extended implementation, the wave velocity is defined as an array that varies across the grid, allowing the simulation to capture realistic scenarios where material properties differ spatially. The finite difference scheme remains similar to the homogeneous case, but the use of a local velocity at each grid point introduces variations that lead to more complex wave behavior. Reflective boundary conditions continue to ensure that the wave remains within the computational domain.

The mathematical models and numerical implementations presented here lay the foundation for more advanced simulations. Rust’s integration with external libraries such as nalgebra and ndarray provides a robust environment for scaling these simulations to two or three dimensions, handling larger data sets, and applying more sophisticated numerical techniques such as the finite element method (FEM) or spectral element methods (SEM). This approach ultimately equips researchers and engineers with the computational tools needed to explore and predict seismic phenomena in increasingly realistic settings.

51.3. Numerical Methods for Seismic Simulation link

Seismic wave simulations depend on robust numerical methods to solve the governing partial differential equations that describe wave propagation. Among the most widely used methods are the Finite Difference Time-Domain (FDTD) method, the Finite Element Method (FEM), and the Spectral Element Method (SEM). These techniques enable the simulation of seismic wave propagation in complex geological settings. Equally important is the choice of time-stepping algorithms, which are essential to accurately update the wavefield over time, and the careful management of boundary conditions to maintain stability and minimize artifacts.

The FDTD method discretizes both time and space to solve the wave equation using finite difference approximations. It is well suited for structured grids and simple geometries, offering a straightforward and computationally efficient approach. Nevertheless, its relative inflexibility when dealing with irregular domains or heterogeneous material properties has led to the development of alternative methods.

The FEM divides the simulation domain into discrete elements, such as triangles in 2D or tetrahedra in 3D, and solves the wave equation within each element using variational techniques. This approach is particularly effective for complex geometries and heterogeneous media, though it typically comes with a higher computational cost compared to FDTD.

The SEM combines the geometric flexibility of FEM with the high accuracy of spectral methods by employing high-order polynomial basis functions within each element. While SEM delivers excellent accuracy for simulations that involve complex material properties and boundary conditions, it is computationally more demanding and involves a more intricate implementation.

Time-stepping in these numerical schemes is critical to capturing the dynamic evolution of the seismic wavefield. Explicit time-stepping methods are often favored for their simplicity and efficiency, but they require adherence to stability criteria, such as the Courant–Friedrichs–Lewy (CFL) condition, which restricts the maximum time step based on the spatial grid spacing and maximum wave velocity in the medium. In contrast, implicit methods allow for larger time steps at the cost of solving a system of equations at each iteration, making them better suited for stiff problems.

In addition to the discretization schemes and time-stepping algorithms, appropriate boundary conditions are necessary to emulate realistic wave behavior. Reflecting boundaries are straightforward to implement but can introduce artificial interference patterns. To counteract this, absorbing boundary conditions – for example, via Perfectly Matched Layers (PML) – are employed to simulate an open domain by damping outgoing waves.

Below is a Rust implementation of the FDTD method for simulating 2D seismic wave propagation with absorbing boundary conditions. The code discretizes the wave equation in both time and space and includes inline comments to explain each step.

extern crate ndarray;
use ndarray::Array2;

/// Simulate 2D seismic wave propagation using the Finite Difference Time-Domain (FDTD) method.
/// This function applies absorbing boundary conditions to dampen outgoing waves.
///
/// # Arguments
/// * `time_steps` - The total number of time iterations to perform.
/// * `nx` - Number of grid points in the x-direction.
/// * `ny` - Number of grid points in the y-direction.
/// * `velocity` - Constant wave propagation velocity.
/// * `dt` - Time step size.
/// * `dx` - Spatial step size in the x-direction.
/// * `dy` - Spatial step size in the y-direction.
/// * `absorbing_factor` - Damping factor applied at the boundaries (value between 0 and 1).
///
/// # Returns
/// A 2D array representing the final wavefield after simulation.
fn fdtd_2d(
    time_steps: usize,
    nx: usize,
    ny: usize,
    velocity: f64,
    dt: f64,
    dx: f64,
    dy: f64,
    absorbing_factor: f64,
) -> Array2<f64> {
    // Initialize arrays for previous, current, and next time steps.
    let mut u_prev = Array2::<f64>::zeros((nx, ny));
    let mut u_curr = Array2::<f64>::zeros((nx, ny));
    let mut u_next = Array2::<f64>::zeros((nx, ny));

    // Set initial condition: a pulse in the center of the grid.
    u_curr[[nx / 2, ny / 2]] = 1.0;

    // Time-stepping loop.
    for _ in 0..time_steps {
        // Iterate over interior grid points (avoiding boundaries).
        for i in 1..nx - 1 {
            for j in 1..ny - 1 {
                // Compute second spatial derivatives using central difference approximations.
                let laplacian_x = (u_curr[[i + 1, j]] - 2.0 * u_curr[[i, j]] + u_curr[[i - 1, j]]) / (dx * dx);
                let laplacian_y = (u_curr[[i, j + 1]] - 2.0 * u_curr[[i, j]] + u_curr[[i, j - 1]]) / (dy * dy);
                // Update the wavefield using the discretized wave equation.
                u_next[[i, j]] = 2.0 * u_curr[[i, j]] - u_prev[[i, j]] + (velocity * velocity * dt * dt) * (laplacian_x + laplacian_y);
            }
        }

        // Apply absorbing boundary conditions on the edges to dampen outgoing waves.
        for i in 0..nx {
            u_next[[i, 0]] *= absorbing_factor;
            u_next[[i, ny - 1]] *= absorbing_factor;
        }
        for j in 0..ny {
            u_next[[0, j]] *= absorbing_factor;
            u_next[[nx - 1, j]] *= absorbing_factor;
        }

        // Shift the wavefields: current becomes previous, and next becomes current.
        u_prev.assign(&u_curr);
        u_curr.assign(&u_next);
    }

    // Return the final wavefield.
    u_curr
}

fn main() {
    // Define simulation parameters for the 2D FDTD method.
    let time_steps = 500;
    let nx = 100;
    let ny = 100;
    let velocity = 1.0;
    let dt = 0.01;
    let dx = 0.1;
    let dy = 0.1;
    let absorbing_factor = 0.9;

    // Run the 2D FDTD simulation.
    let result = fdtd_2d(time_steps, nx, ny, velocity, dt, dx, dy, absorbing_factor);
    println!("Final 2D wavefield: {:?}", result);
}

The FEM approach is especially beneficial when the simulation domain is irregular or when the material properties vary spatially. In the example below, a simplified 1D FEM implementation is provided using Rust. The simulation divides the domain into finite elements and uses matrix operations to solve the discretized wave equation.

extern crate nalgebra as na;
use na::{DMatrix, DVector};

/// Simulate 1D seismic wave propagation using the Finite Element Method (FEM).
/// This function solves the wave equation by constructing simplified stiffness and mass matrices.
///
/// # Arguments
/// * `time_steps` - Total number of time iterations.
/// * `elements` - Number of finite elements in the domain.
/// * `velocity` - Constant wave propagation velocity (used in scaling the system).
/// * `dt` - Time step size.
///
/// # Returns
/// A vector representing the final displacement field.
fn fem_1d_wave(time_steps: usize, elements: usize, velocity: f64, dt: f64) -> DVector<f64> {
    // Initialize displacement vectors for previous, current, and next time steps.
    let mut u_prev = DVector::zeros(elements);
    let mut u_curr = DVector::zeros(elements);
    let mut u_next = DVector::zeros(elements);

    // Initial condition: a pulse at the center of the domain.
    u_curr[elements / 2] = 1.0;

    // Construct simplified stiffness and mass matrices.
    // Here, the stiffness matrix approximates the second derivative operator.
    let mut stiffness_matrix = DMatrix::<f64>::zeros(elements, elements);
    for i in 0..elements {
        stiffness_matrix[(i, i)] = -2.0;
        if i > 0 {
            stiffness_matrix[(i, i - 1)] = 1.0;
        }
        if i < elements - 1 {
            stiffness_matrix[(i, i + 1)] = 1.0;
        }
    }
    let mass_matrix = DMatrix::<f64>::identity(elements, elements);

    // Time-stepping loop using an implicit method to ensure stability.
    for _ in 0..time_steps {
        // Construct the right-hand side vector using the current and previous displacement vectors.
        let rhs = &mass_matrix * (2.0 * &u_curr - &u_prev) - &stiffness_matrix * &u_curr;
        // Solve the linear system to obtain the next time step displacement.
        // Using LU decomposition; the 'expect' call will panic if the system cannot be solved.
        u_next = mass_matrix.clone().lu().solve(&rhs).expect("Failed to solve FEM system");
        // Update displacement vectors for the next iteration.
        u_prev = u_curr.clone();
        u_curr = u_next.clone();
    }

    // Return the final displacement field.
    u_curr
}

fn main() {
    // Simulation parameters for the 1D FEM model.
    let time_steps = 500;
    let elements = 100;
    let velocity = 1.0;
    let dt = 0.01;

    // Execute the FEM simulation and output the final displacement.
    let result = fem_1d_wave(time_steps, elements, velocity, dt);
    println!("Final 1D displacement (FEM): {:?}", result);
}

The SEM further enhances accuracy by integrating high-order polynomial basis functions into the finite element framework. This approach offers superior precision with fewer elements when modeling complex material properties and boundary conditions. The following code provides a simplified SEM implementation for a 1D seismic wave simulation.

extern crate nalgebra as na;
use na::{DMatrix, DVector};

/// Simulate 1D seismic wave propagation using the Finite Element Method (FEM).
/// This function solves the 1D wave equation by constructing simplified stiffness and mass matrices,
/// and then time-stepping the solution using an explicit scheme for demonstration purposes.
/// The wave equation is given by uₜₜ = c² uₓₓ, which is discretized in time as:
///
///     uⁿ⁺¹ = 2 uⁿ − uⁿ⁻¹ + dt² * c² * uₓₓ
///
/// In this implementation the stiffness matrix approximates the second spatial derivative (uₓₓ)
/// and the mass matrix is taken as the identity matrix for simplicity.
///
/// # Arguments
/// * `time_steps` - Total number of time iterations.
/// * `elements` - Number of finite elements in the domain.
/// * `velocity` - Wave propagation velocity (c), used in scaling the system.
/// * `dt` - Time step size.
///
/// # Returns
/// A vector representing the final displacement field.
fn fem_1d_wave(time_steps: usize, elements: usize, velocity: f64, dt: f64) -> DVector<f64> {
    // Initialize displacement vectors for previous, current, and next time steps.
    let mut u_prev = DVector::zeros(elements);
    let mut u_curr = DVector::zeros(elements);
    let mut u_next = DVector::zeros(elements);

    // Initial condition: a pulse at the center of the domain.
    u_curr[elements / 2] = 1.0;

    // Construct simplified stiffness matrix approximating the second derivative.
    // For each node, uₓₓ ≈ u[i-1] - 2 u[i] + u[i+1].
    let mut stiffness_matrix = DMatrix::<f64>::zeros(elements, elements);
    for i in 0..elements {
        stiffness_matrix[(i, i)] = -2.0;
        if i > 0 {
            stiffness_matrix[(i, i - 1)] = 1.0;
        }
        if i < elements - 1 {
            stiffness_matrix[(i, i + 1)] = 1.0;
        }
    }
    // Use an identity matrix as the mass matrix.
    let mass_matrix = DMatrix::<f64>::identity(elements, elements);

    // Time-stepping loop using an explicit update scheme.
    // The update formula for the wave equation is:
    //     u_next = 2 u_curr - u_prev + dt² * velocity² * (stiffness_matrix * u_curr)
    for _ in 0..time_steps {
        // Compute the right-hand side:
        // Note: mass_matrix * (2 u_curr - u_prev) is simply (2 u_curr - u_prev) since mass_matrix is identity.
        let rhs = &mass_matrix * (2.0 * &u_curr - &u_prev)
            + (dt * dt * velocity * velocity) * (&stiffness_matrix * &u_curr);
        // Solve the linear system: mass_matrix * u_next = rhs.
        // Since mass_matrix is the identity, this simply yields u_next = rhs.
        u_next = mass_matrix.clone().lu().solve(&rhs).expect("Failed to solve FEM system");
        // Update displacement vectors for the next iteration.
        u_prev = u_curr.clone();
        u_curr = u_next.clone();
    }

    // Return the final displacement field.
    u_curr
}

fn main() {
    // Simulation parameters for the 1D FEM model.
    let time_steps = 500;
    let elements = 100;
    let velocity = 1.0;
    let dt = 0.01;

    // Execute the FEM simulation and output the final displacement.
    let result = fem_1d_wave(time_steps, elements, velocity, dt);
    println!("Final 1D displacement (FEM): {:?}", result);
}

In this section, the numerical methods underpinning seismic simulations are discussed in detail. Each method – FDTD, FEM, and SEM – provides a distinct balance between computational efficiency and accuracy. The FDTD method offers a straightforward explicit approach ideal for structured grids, whereas the FEM is highly flexible for complex and heterogeneous domains. The SEM pushes accuracy further by using high-order basis functions, making it well-suited for large-scale simulations. Practical Rust implementations provided for each method demonstrate how these techniques can be effectively utilized to simulate seismic wave propagation in both 1D and 2D, addressing complex boundary conditions and ensuring numerical stability.

51.4. Seismic Wave Attenuation and Dispersion link

Seismic wave attenuation and dispersion are critical phenomena that fundamentally influence the propagation characteristics of waves through the Earth’s subsurface. Attenuation describes the gradual reduction in wave amplitude as energy is lost during propagation, while dispersion refers to the frequency-dependent variation in wave velocity that leads to the spreading of a wave pulse over time. These effects carry vital information about the physical properties of geological layers, including absorption due to internal friction, scattering from heterogeneities in the medium, and anelastic behavior where the material does not fully recover after deformation. Understanding and accurately modeling these processes is essential for both theoretical investigations and practical applications such as resource exploration and seismic hazard assessment.

Attenuation occurs as wave energy is converted into heat or redistributed by interactions with microstructures in the rock. The quality factor Q is commonly used to quantify attenuation; a lower Q signifies stronger attenuation and vice versa. Mathematically, the attenuation coefficient α is related to the angular frequency ω, the wave velocity v, and the quality factor Q via the relationship

$$α = ω / (2 v Q)$$

This relation underscores that higher frequency components of a seismic wave are generally more affected by attenuation than lower frequency components. In contrast, dispersion arises when different frequency components of a wave travel at varying speeds due to the medium’s elastic and viscous properties. This frequency dependence not only modifies the shape of the wave packet over time but also provides insights into the material’s density, viscosity, and elastic moduli.

The following Rust implementations demonstrate how to incorporate these phenomena into seismic wave simulations. The first code block illustrates a one-dimensional simulation that introduces an attenuation term into the finite difference approximation of the wave equation using a specified Q factor. The second code block presents a simplified approach to modeling dispersion by transforming the wavefield into the frequency domain, applying a frequency-dependent velocity, and then returning to the time domain using a Fourier transform.

extern crate ndarray;
use ndarray::Array1;

/// Simulate 1D seismic wave propagation with attenuation using a finite difference method.
/// The attenuation is modeled through an attenuation coefficient derived from the wave velocity,
/// quality factor Q, and spatial resolution. This function updates the wave amplitude at each grid point,
/// gradually reducing it over time according to the specified Q factor.
///
/// # Arguments
/// * `time_steps` - The number of time iterations for the simulation.
/// * `grid_points` - The number of spatial grid points representing the medium.
/// * `velocity` - The constant propagation velocity of the seismic wave.
/// * `dt` - The time step size.
/// * `dx` - The spatial step size.
/// * `q_factor` - The quality factor representing the level of attenuation (higher Q means lower attenuation).
///
/// # Returns
/// A one-dimensional array representing the final displacement field with attenuation applied.
fn simulate_wave_with_attenuation(
    time_steps: usize,
    grid_points: usize,
    velocity: f64,
    dt: f64,
    dx: f64,
    q_factor: f64,
) -> Array1<f64> {
    // Initialize arrays for the displacement field at previous, current, and next time steps.
    let mut u_prev = Array1::<f64>::zeros(grid_points);
    let mut u_curr = Array1::<f64>::zeros(grid_points);
    let mut u_next = Array1::<f64>::zeros(grid_points);

    // Calculate an attenuation coefficient based on the wave velocity, Q factor, and spatial step size.
    // This coefficient scales the amplitude reduction per time step.
    let attenuation_coefficient = velocity / (q_factor * dx);

    // Set an initial condition: a pulse introduced at the center of the grid.
    u_curr[grid_points / 2] = 1.0;

    // Time-stepping loop to update the wavefield.
    for _ in 0..time_steps {
        for i in 1..grid_points - 1 {
            // Compute the second spatial derivative using a central difference scheme.
            let spatial_derivative = (u_curr[i + 1] - 2.0 * u_curr[i] + u_curr[i - 1]) / (dx * dx);
            // Apply an attenuation term that reduces the amplitude based on the attenuation coefficient.
            let attenuation_term = 1.0 - (attenuation_coefficient * dt);
            // Update the displacement using the finite difference formulation with attenuation.
            u_next[i] = (2.0 * u_curr[i] - u_prev[i] + (velocity * velocity * dt * dt) * spatial_derivative)
                        * attenuation_term;
        }

        // Apply reflective boundary conditions to the edges of the grid.
        u_next[0] = u_next[1];
        u_next[grid_points - 1] = u_next[grid_points - 2];

        // Shift the arrays to prepare for the next iteration.
        u_prev.assign(&u_curr);
        u_curr.assign(&u_next);
    }

    // Return the final displacement field incorporating attenuation.
    u_curr
}

fn main() {
    // Define simulation parameters.
    let time_steps = 500;
    let grid_points = 100;
    let velocity = 1.0;      // Propagation velocity of the seismic wave.
    let dt = 0.01;           // Time step size.
    let dx = 0.1;            // Spatial step size.
    let q_factor = 50.0;     // Quality factor indicating the level of attenuation.

    // Execute the simulation and print the final displacement field.
    let result = simulate_wave_with_attenuation(time_steps, grid_points, velocity, dt, dx, q_factor);
    println!("Final wave displacement with attenuation: {:?}", result);
}

In order to model dispersion, we introduce frequency-dependent velocity variations into the wave equation. Dispersion is simulated by transforming the wavefield to the frequency domain, modifying each frequency component with a velocity that increases with frequency, and then transforming the wavefield back to the time domain. The following code utilizes the rustfft crate along with ndarray and num_complex to perform Fourier and inverse Fourier transforms, thereby incorporating dispersive effects into the seismic wave simulation.

extern crate ndarray;
extern crate num_complex;
extern crate rustfft;

use ndarray::Array1;
use num_complex::Complex;
use rustfft::FftPlanner;

/// Simulate 1D seismic wave dispersion by transforming the wavefield into the frequency domain,
/// applying frequency-dependent velocity modifications, and then transforming back to the time domain.
/// The frequency-dependent velocity is modeled as an adjustment to the base velocity that increases linearly
/// with frequency, simulating dispersive behavior.
///
/// # Arguments
/// * `time_steps` - The number of time iterations for the simulation.
/// * `grid_points` - The number of spatial grid points representing the medium.
/// * `base_velocity` - The base wave propagation velocity.
/// * `dt` - The time step size.
/// * `dx` - The spatial step size.
///
/// # Returns
/// A one-dimensional array representing the final wavefield after applying dispersion.
fn simulate_dispersion(
    time_steps: usize,
    grid_points: usize,
    base_velocity: f64,
    dt: f64,
    dx: f64,
) -> Array1<f64> {
    // Initialize the time-domain wavefield with zeros and set a pulse at the center.
    let mut u_time = Array1::<Complex<f64>>::from_elem(grid_points, Complex::new(0.0, 0.0));
    u_time[grid_points / 2] = Complex::new(1.0, 0.0);

    // Prepare an array for the frequency-domain representation.
    let mut u_freq = u_time.clone();

    // Initialize FFT planners for forward and inverse transforms.
    let mut planner = FftPlanner::<f64>::new();
    let fft = planner.plan_fft_forward(grid_points);
    let ifft = planner.plan_fft_inverse(grid_points);

    // Time-stepping loop to simulate dispersion.
    for _ in 0..time_steps {
        // Perform forward Fourier transform to convert the wavefield to the frequency domain.
        fft.process(u_time.as_slice_mut().unwrap(), u_freq.as_slice_mut().unwrap());

        // Apply frequency-dependent modifications to each frequency component.
        // The modification simulates dispersion by increasing the effective velocity for higher frequencies.
        for (i, freq_component) in u_freq.iter_mut().enumerate() {
            // Compute the frequency corresponding to this index.
            let freq = (i as f64) / (grid_points as f64 * dx);
            // Define a dispersive velocity that increases linearly with frequency.
            let dispersive_velocity = base_velocity * (1.0 + 0.1 * freq);
            // Scale the frequency component by the dispersive velocity.
            *freq_component = *freq_component * dispersive_velocity;
        }

        // Perform inverse Fourier transform to revert the wavefield back to the time domain.
        ifft.process(u_freq.as_slice_mut().unwrap(), u_time.as_slice_mut().unwrap());
    }

    // Extract the real part of the final wavefield to represent physical displacement.
    let result: Array1<f64> = u_time.mapv(|c| c.re);
    result
}

fn main() {
    // Define simulation parameters.
    let time_steps = 500;
    let grid_points = 100;
    let base_velocity = 1.0; // Base wave propagation velocity.
    let dt = 0.01;           // Time step size.
    let dx = 0.1;            // Spatial step size.

    // Run the dispersion simulation and print the final wavefield.
    let result = simulate_dispersion(time_steps, grid_points, base_velocity, dt, dx);
    println!("Final wave displacement with dispersion: {:?}", result);
}

The models and simulations presented above illustrate how seismic wave attenuation and dispersion can be incorporated into numerical simulations. By integrating physical parameters such as the quality factor and frequency-dependent velocities into the finite difference and Fourier-based frameworks, these examples provide a pathway to simulate complex seismic phenomena. Adjusting parameters like the Q factor or the dispersive scaling factor enables the modeling of various geological conditions and enhances our understanding of wave behavior in the Earth's subsurface.

51.5. Seismic Wave Reflection and Refraction link

Seismic waves encounter interfaces between geological layers, and at these boundaries, they experience reflection and refraction that are central to seismic imaging and exploration techniques. When a wave reaches a boundary such as a transition from sediment to bedrock or a subsurface fault, the contrast in material properties causes a portion of the wave energy to reflect back into the original medium while the remainder is transmitted into the new medium. The behavior of the wave at such an interface is governed by Snell’s law, which relates the angles of incidence, reflection, and refraction to the velocities in the two media. In mathematical terms, Snell’s law is expressed as

$$sin(θ₁) / v₁ = sin(θ₂) / v₂$$

where θ₁ is the angle of incidence, θ₂ is the angle of refraction, and v₁ and v₂ denote the wave velocities in the respective media. The acoustic impedance of each medium, defined by Z = ρv with ρ representing density, plays a crucial role in determining the amount of energy reflected or transmitted at the boundary.

The reflection coefficient R and the transmission coefficient T quantify these proportions. They are given by the expressions

$$R = (Z₂ - Z₁) / (Z₂ + Z₁)$$

$$T = 2Z₁ / (Z₂ + Z₁)$$

with Z₁ and Z₂ being the impedances of the first and second media, respectively. These relationships indicate that larger impedance contrasts, such as those between soft sediments and solid bedrock, result in stronger reflections that are key to delineating subsurface features. In addition, at certain angles known as critical angles, total internal reflection may occur, and mode conversion between P-waves and S-waves can take place at the interface.

The following Rust implementations provide practical examples of how to simulate seismic wave reflection and refraction at a geological boundary. In the first example, a function is defined to compute the reflection and transmission coefficients based on input velocities, densities, and the incident angle. In the second example, a finite difference method is used to simulate the propagation of a wavefield through a one-dimensional medium containing an interface between two layers. At the interface, the computed coefficients are applied to adjust the wavefield in a manner that mimics the reflection and transmission of the incident wave.

extern crate ndarray;
use ndarray::Array1;
use std::f64::consts::PI;

/// Compute the reflection and transmission coefficients based on acoustic impedance contrast.
/// The function uses Snell's law to compute the refracted angle given an incident angle.
/// 
/// # Arguments
/// * `v1` - Wave velocity in the first medium.
/// * `v2` - Wave velocity in the second medium.
/// * `rho1` - Density of the first medium.
/// * `rho2` - Density of the second medium.
/// * `incident_angle` - Incident angle in radians.
///
/// # Returns
/// A tuple containing the reflection coefficient and the transmission coefficient.
fn compute_reflection_transmission(
    v1: f64,
    v2: f64,
    rho1: f64,
    rho2: f64,
    incident_angle: f64,
) -> (f64, f64) {
    // Compute acoustic impedances for both media.
    let z1 = rho1 * v1;
    let z2 = rho2 * v2;
    
    // Use Snell's law to compute the refracted angle.
    // Ensure that the value inside asin does not exceed 1.0.
    let sin_refracted = (v2 / v1 * incident_angle.sin()).min(1.0);
    let _refracted_angle = sin_refracted.asin();
    
    // Calculate the reflection and transmission coefficients.
    let reflection_coefficient = (z2 - z1) / (z2 + z1);
    let transmission_coefficient = 2.0 * z1 / (z2 + z1);
    
    (reflection_coefficient, transmission_coefficient)
}

/// Simulate the propagation of a seismic wave in a 1D medium with an interface between two layers.
/// The wavefield is updated using a finite difference approximation of the wave equation.
/// The simulation applies reflection and transmission at the interface based on the calculated coefficients.
///
/// # Arguments
/// * `time_steps` - Number of time iterations for the simulation.
/// * `grid_points` - Total number of spatial grid points.
/// * `velocity1` - Wave velocity in the first layer.
/// * `velocity2` - Wave velocity in the second layer.
/// * `dx` - Spatial step size.
/// * `interface` - Index representing the interface between the two layers.
///
/// # Returns
/// A one-dimensional array representing the final wavefield after simulation.
fn simulate_wave_propagation_with_interface(
    time_steps: usize,
    grid_points: usize,
    velocity1: f64,
    velocity2: f64,
    dx: f64,
    interface: usize,
) -> Array1<f64> {
    // Initialize the wavefield arrays for the previous, current, and next time steps.
    let mut wavefield_prev = Array1::<f64>::zeros(grid_points);
    let mut wavefield = Array1::<f64>::zeros(grid_points);
    let mut wavefield_next = Array1::<f64>::zeros(grid_points);
    
    // Set initial condition: a pulse at the start of the grid.
    wavefield[0] = 1.0;
    
    // Simulation loop over time steps.
    for _ in 0..time_steps {
        // Update interior grid points using finite difference approximation.
        for i in 1..grid_points - 1 {
            // Determine the local wave velocity based on the position relative to the interface.
            let local_velocity = if i < interface { velocity1 } else { velocity2 };
            
            // Calculate the second derivative with respect to space.
            let second_derivative = (wavefield[i + 1] - 2.0 * wavefield[i] + wavefield[i - 1]) / (dx * dx);
            
            // Update the wavefield using the finite difference time-stepping scheme.
            wavefield_next[i] = 2.0 * wavefield[i] - wavefield_prev[i] + (local_velocity * local_velocity) * second_derivative;
        }
        
        // Compute reflection and transmission coefficients at the interface.
        // Example densities for the two layers are chosen; these can be adjusted for different scenarios.
        let (reflection_coefficient, transmission_coefficient) =
            compute_reflection_transmission(velocity1, velocity2, 2.5, 3.0, 0.0);
        
        // At the interface, combine the incoming wave from the left and the adjacent value on the right.
        wavefield_next[interface] = reflection_coefficient * wavefield[interface - 1] + transmission_coefficient * wavefield[interface + 1];
        
        // Update the wavefields for the next time iteration.
        wavefield_prev.assign(&wavefield);
        wavefield.assign(&wavefield_next);
    }
    
    // Return the final wavefield.
    wavefield
}

fn main() {
    // Define simulation parameters.
    let time_steps = 500;
    let grid_points = 200;
    let velocity1 = 1.0;  // Wave velocity in the first layer.
    let velocity2 = 2.0;  // Wave velocity in the second layer.
    let dx = 0.1;         // Spatial step size.
    let interface = 100;  // Index of the interface between the two layers.
    
    // Run the simulation of wave propagation with an interface.
    let result = simulate_wave_propagation_with_interface(time_steps, grid_points, velocity1, velocity2, dx, interface);
    println!("Final wavefield: {:?}", result);
}

Seismic wave reflection and refraction at geological interfaces are fundamental to imaging the Earth’s subsurface. As a seismic wave propagates, changes in acoustic impedance due to variations in density and velocity lead to the partitioning of energy between reflected and transmitted waves. In seismic surveys, these phenomena enable the reconstruction of subsurface structures by analyzing travel times and amplitudes of the returning waves. The Rust code provided above demonstrates how to calculate reflection and transmission coefficients using Snell’s law and impedance contrasts, and how to integrate these calculations into a finite difference simulation of wave propagation. By adjusting velocities, densities, and the position of the interface, various subsurface scenarios can be explored, highlighting the sensitivity of seismic responses to changes in material properties. Advanced simulations may include additional complexities such as absorbing boundaries, mode conversion, and multiple interfaces, leveraging Rust’s performance and concurrency features to model large-scale seismic phenomena accurately.

51.6. Seismic Tomography and Inversion Techniques link

Seismic tomography and inversion techniques offer powerful approaches for imaging the Earth's interior and inferring subsurface structures by analyzing seismic wave travel times and waveform data. These methods are analogous to medical CT scans, where waves are used to construct three-dimensional images of internal structures. In seismic tomography, the travel times of seismic waves recorded at various stations are used to reconstruct the spatial variations in seismic velocity. This information, in turn, helps to delineate different geological formations and detect anomalies in the subsurface.

In the process of seismic tomography, forward modeling plays a crucial role. Forward modeling involves generating synthetic travel times or waveforms from a given velocity model by solving the wave equation. These synthetic data are then compared to observed seismic measurements. The discrepancies, or residuals, between the predicted and observed data drive the inversion process, which iteratively updates the subsurface model to reduce these differences. This inversion process is challenging due to the non-linear relationship between the observed data and the model parameters, and the inherent ill-posedness of the problem where many models can explain the data. Regularization techniques such as Tikhonov regularization or smoothing constraints are often employed to stabilize the inversion and yield physically meaningful solutions.

Travel-time tomography focuses on reconstructing the velocity structure of the subsurface using the travel times of seismic waves between sources and receivers. The basic procedure involves two steps. First, forward modeling computes synthetic travel times based on an initial velocity model; then, inversion techniques update the model by correcting the differences between observed and predicted travel times. The following Rust implementation demonstrates these basic principles. In this example, a simple forward model computes travel times using straight-line distances from a source to multiple receiver locations. The inversion procedure then iteratively updates the velocity model by applying a correction proportional to the travel-time residuals.

extern crate ndarray;
use ndarray::{Array1, Array2};

/// Compute synthetic travel times for seismic waves propagating from a source to multiple receivers
/// using a simple straight-line distance approximation. The travel time is computed as the distance
/// divided by the local velocity at the source.
///
/// # Arguments
/// * `velocity` - A 2D array representing the subsurface velocity model.
/// * `source` - A tuple representing the grid indices of the seismic source.
/// * `receivers` - A vector of tuples, each representing the grid indices of a receiver location.
/// * `dx` - The spatial grid spacing.
///
/// # Returns
/// An array of travel times corresponding to each receiver.
fn travel_time_forward_model(
    velocity: &Array2<f64>,
    source: (usize, usize),
    receivers: &Vec<(usize, usize)>,
    dx: f64,
) -> Array1<f64> {
    let mut travel_times = Array1::<f64>::zeros(receivers.len());

    for (i, &(rx, ry)) in receivers.iter().enumerate() {
        // Compute the Euclidean distance from the source to the receiver and scale by dx.
        let distance = (((rx as f64 - source.0 as f64).powi(2)
            + (ry as f64 - source.1 as f64).powi(2))
            .sqrt()) * dx;
        // Use the velocity at the source location for a basic approximation.
        let velocity_at_source = velocity[[source.0, source.1]];
        travel_times[i] = distance / velocity_at_source;
    }

    travel_times
}

/// Perform a simple travel-time inversion to update the velocity model. The function compares observed
/// travel times with those predicted by the current model and applies a correction to the velocity
/// at the receiver locations based on the residuals.
///
/// # Arguments
/// * `velocity` - A mutable 2D array representing the subsurface velocity model.
/// * `observed_times` - An array of observed travel times for each receiver.
/// * `source` - A tuple representing the grid indices of the seismic source.
/// * `receivers` - A vector of tuples representing the receiver locations.
/// * `dx` - The spatial grid spacing.
/// * `iterations` - Number of inversion iterations to perform.
fn travel_time_inversion(
    velocity: &mut Array2<f64>,
    observed_times: &Array1<f64>,
    source: (usize, usize),
    receivers: &Vec<(usize, usize)>,
    dx: f64,
    iterations: usize,
) {
    for _ in 0..iterations {
        // Compute predicted travel times based on the current velocity model.
        let predicted_times = travel_time_forward_model(velocity, source, receivers, dx);
        // Calculate residuals between observed and predicted travel times.
        let residuals = observed_times - &predicted_times;

        // Update the velocity model at each receiver location using a small correction factor.
        for (i, &(rx, ry)) in receivers.iter().enumerate() {
            let correction = residuals[i] * 0.1; // Correction factor for stability.
            velocity[[rx, ry]] += correction;
        }
    }
}

fn main() {
    let nx = 100;
    let ny = 100;
    let dx = 1.0; // Spatial step size.
    // Initialize a simple velocity model (e.g., 3.0 km/s across the domain).
    let mut velocity = Array2::<f64>::from_elem((nx, ny), 3.0);
    let source = (50, 50);
    let receivers = vec![(20, 20), (80, 20), (50, 80)]; // Receiver locations.
    // Mock observed travel times for demonstration purposes.
    let observed_times = Array1::<f64>::from_vec(vec![12.5, 10.0, 15.0]);

    // Run the inversion process over a set number of iterations to update the velocity model.
    travel_time_inversion(&mut velocity, &observed_times, source, &receivers, dx, 100);

    println!("Updated velocity model: {:?}", velocity);
}

For more advanced imaging techniques, Full-Waveform Inversion (FWI) uses the complete waveform, including both amplitudes and phases, to update the velocity model. FWI is computationally demanding as it requires iterative forward modeling and adjustment of the model based on the residuals between observed and synthetic waveforms. A simplified example of FWI is provided below. In this demonstration, a basic forward wave simulation function is defined (for illustrative purposes), and then the inversion updates the velocity model using the computed residuals from the observed waveforms.

extern crate ndarray;
use ndarray::{Array1, Array2};

/// A simplified forward simulation function for the wave equation that generates a synthetic wavefield
/// from a given velocity model and source location. This implementation is a basic placeholder to illustrate
/// the forward modeling process for FWI.
///
/// # Arguments
/// * `velocity` - A 2D array representing the subsurface velocity model.
/// * `source` - A tuple representing the grid indices of the seismic source.
/// * `time_steps` - The number of time iterations for the simulation.
/// * `dx` - The spatial grid spacing.
/// * `dt` - The time step size.
///
/// # Returns
/// A 2D array representing the synthetic wavefield after the simulation.
fn forward_wave_simulation(
    velocity: &Array2<f64>,
    source: (usize, usize),
    time_steps: usize,
    dx: f64,
    dt: f64,
) -> Array2<f64> {
    // Create an array for the wavefield with the same dimensions as the velocity model.
    let mut wavefield = Array2::<f64>::zeros(velocity.raw_dim());
    // Set an initial pulse at the source location.
    wavefield[[source.0, source.1]] = 1.0;

    // A simple and highly approximate wave propagation loop for demonstration purposes.
    for _ in 0..time_steps {
        // In a full implementation, update the wavefield based on the wave equation.
        // This placeholder does not perform an actual finite difference update.
        // The focus is on demonstrating the inversion update process.
    }

    wavefield
}

/// Perform a basic Full-Waveform Inversion (FWI) to update the velocity model by minimizing the difference
/// between observed and synthetic waveforms. The inversion loop updates the velocity model by applying a
/// correction proportional to the residuals at each grid point.
///
/// # Arguments
/// * `velocity` - A mutable 2D array representing the initial subsurface velocity model.
/// * `observed_waveforms` - A 2D array containing the observed waveform data.
/// * `source` - A tuple representing the grid indices of the seismic source.
/// * `time_steps` - The number of time iterations for the forward simulation.
/// * `dx` - The spatial grid spacing.
/// * `dt` - The time step size.
/// * `iterations` - The number of inversion iterations to perform.
fn full_waveform_inversion(
    velocity: &mut Array2<f64>,
    observed_waveforms: &Array2<f64>,
    source: (usize, usize),
    time_steps: usize,
    dx: f64,
    dt: f64,
    iterations: usize,
) {
    for _ in 0..iterations {
        // Generate synthetic waveforms using the current velocity model.
        let predicted_waveforms = forward_wave_simulation(velocity, source, time_steps, dx, dt);
        // Compute residuals between observed and predicted waveforms.
        let residuals = observed_waveforms - &predicted_waveforms;

        // Update the velocity model by applying a small correction based on the residuals.
        for i in 0..velocity.shape()[0] {
            for j in 0..velocity.shape()[1] {
                velocity[[i, j]] += residuals[[i, j]] * 0.1; // Correction factor for stability.
            }
        }
    }
}

fn main() {
    let nx = 100;
    let ny = 100;
    let time_steps = 500;
    let dx = 1.0;
    let dt = 0.01;
    // Initialize the velocity model with a uniform value (e.g., 3.0 km/s).
    let mut velocity = Array2::<f64>::from_elem((nx, ny), 3.0);
    let source = (50, 50);

    // For demonstration purposes, create a mock observed waveform as a zero matrix.
    // In practice, this would be replaced with real seismic waveform data.
    let observed_waveforms = Array2::<f64>::zeros((nx, ny));

    // Run the Full-Waveform Inversion (FWI) process over a number of iterations.
    full_waveform_inversion(&mut velocity, &observed_waveforms, source, time_steps, dx, dt, 100);

    println!("Updated velocity model after FWI: {:?}", velocity);
}

Seismic tomography and inversion techniques, including travel-time tomography and full-waveform inversion, are fundamental for reconstructing the subsurface velocity structure from seismic data. The forward modeling steps generate synthetic data based on an initial model, and the inversion process iteratively refines this model by minimizing the residuals between observed and synthetic data. This iterative procedure, while computationally demanding, can reveal detailed images of the Earth's interior and provide insights into geological features and anomalies. Rust's performance and concurrency features make it a robust platform for implementing these advanced numerical techniques, and integration with visualization libraries further aids in the interpretation of the reconstructed subsurface models.

51.7. Earthquake Source Modeling link

Earthquake source modeling plays a critical role in understanding how seismic waves are generated during fault ruptures. By constructing models of the earthquake source, one can simulate the generation of seismic waves, analyze how they propagate through the Earth’s crust, and study the effects of rupture dynamics on wave characteristics. Different models, such as point sources, finite fault models, and kinematic rupture models, are used depending on the scale and complexity of the seismic event. While point source models are suitable for small earthquakes, larger events require finite fault models that divide the fault plane into sub-fault segments, each contributing to the overall energy release. Kinematic rupture models go further by prescribing the slip distribution and rupture velocity along the fault, capturing the evolution of the rupture and the resulting seismic energy release.

The generation of seismic waves is influenced by key source parameters including stress drop, fault slip, rupture velocity, and scaling relationships that link fault dimensions and energy release. A higher stress drop leads to stronger seismic waves, while greater fault slip increases the amplitude of the waves. The rupture velocity, typically a fraction of the shear wave velocity, affects both the frequency content and duration of the seismic waves. These parameters combine to determine the energy radiated during an earthquake and the characteristics of the seismic waves that are ultimately recorded.

To model earthquake sources in Rust, we simulate fault rupture by dividing the fault plane into a grid, where each grid point represents a sub-fault. As the rupture propagates across the fault at a prescribed velocity, each sub-fault releases a certain amount of slip, generating seismic waves. In the following Rust implementation, a simplified finite fault model is presented. In this model, the fault plane is represented by a two-dimensional array, and a corresponding slip distribution is applied to simulate the energy release during rupture. The wavefield is updated over time as the rupture front reaches each sub-fault, with the generated wave amplitude being proportional to the slip and the elapsed time since rupture at that location.

extern crate ndarray;
use ndarray::{Array2, Array1};

/// Simulate seismic wave generation from a finite fault rupture using a simplified model.
/// This function updates the fault plane and generates a corresponding wavefield based on a prescribed
/// slip distribution and rupture velocity. The rupture propagates along the fault plane, and each sub-fault
/// generates a seismic signal as slip is applied.
///
/// # Arguments
/// * `fault_plane` - A mutable 2D array representing the cumulative slip on the fault plane.
/// * `slip_distribution` - A 2D array representing the prescribed slip (in meters) at each sub-fault.
/// * `rupture_velocity` - The speed at which the rupture propagates along the fault (in grid units per time unit).
/// * `time_steps` - The total number of time iterations for the simulation.
/// * `dt` - The time step size.
///
/// # Returns
/// A 2D array representing the simulated seismic wavefield generated by the fault rupture.
fn finite_fault_model(
    fault_plane: &mut Array2<f64>,
    slip_distribution: &Array2<f64>,
    rupture_velocity: f64,
    time_steps: usize,
    dt: f64,
) -> Array2<f64> {
    // Initialize the wavefield with the same dimensions as the fault plane.
    let mut wavefield = Array2::<f64>::zeros(fault_plane.raw_dim());
    
    // Loop over each time step to simulate the rupture process.
    for t in 0..time_steps {
        let current_time = t as f64 * dt;
        
        // Iterate over each sub-fault (grid point) on the fault plane.
        for ((i, j), fault_value) in fault_plane.indexed_iter_mut() {
            // Compute the rupture arrival time at this sub-fault based on its position.
            // Here, we assume that the rupture starts at the top of the fault (i = 0) and propagates downward.
            let rupture_time = (i as f64) / rupture_velocity;
            if current_time >= rupture_time {
                // Apply the prescribed slip to this sub-fault.
                let slip = slip_distribution[[i, j]];
                *fault_value += slip;
                
                // Generate a seismic signal at this location.
                // The amplitude is modeled as proportional to the slip and the elapsed time since rupture.
                wavefield[[i, j]] = slip * (current_time - rupture_time);
            }
        }
    }
    
    wavefield
}

fn main() {
    // Define grid dimensions for the fault plane.
    let nx = 100;
    let ny = 100;
    let time_steps = 500;
    let dt = 0.01;
    let rupture_velocity = 2.0; // Rupture propagation speed in grid units per time unit.
    
    // Initialize the fault plane with zero slip.
    let mut fault_plane = Array2::<f64>::zeros((nx, ny));
    // Define a uniform slip distribution where each sub-fault is prescribed a slip of 1.0 meter.
    let slip_distribution = Array2::<f64>::from_elem((nx, ny), 1.0);
    
    // Simulate the finite fault model to generate the seismic wavefield.
    let wavefield = finite_fault_model(&mut fault_plane, &slip_distribution, rupture_velocity, time_steps, dt);
    
    println!("Final seismic wavefield: {:?}", wavefield);
}

In this implementation, the finite_fault_model function simulates the progression of a fault rupture over time. The fault plane is represented as a 2D grid, and the rupture is assumed to start at the top and propagate downward. As the rupture front reaches each sub-fault, the prescribed slip is applied and a corresponding seismic signal is generated. The wavefield reflects the cumulative effect of slip across the fault, providing a basic approximation of the seismic energy release.

For scenarios that mirror historical earthquake events, the slip distribution may be non-uniform, with higher slip values concentrated in regions where maximum energy was released. The following example demonstrates how to simulate a historical earthquake scenario using a non-uniform slip distribution.

fn historical_earthquake_simulation() {
    let nx = 200;
    let ny = 100;
    let time_steps = 1000;
    let dt = 0.01;
    let rupture_velocity = 2.5; // Rupture velocity in grid units per time unit.
    
    // Initialize the fault plane with no initial slip.
    let mut fault_plane = Array2::<f64>::zeros((nx, ny));
    // Create a non-uniform slip distribution to mimic a historical earthquake.
    let mut slip_distribution = Array2::<f64>::zeros((nx, ny));
    // Assign higher slip values to a central region of the fault.
    for i in 50..150 {
        for j in 30..70 {
            slip_distribution[[i, j]] = 2.0; // A higher slip of 2.0 meters in the center.
        }
    }
    
    // Simulate seismic wave generation using the finite fault model.
    let wavefield = finite_fault_model(&mut fault_plane, &slip_distribution, rupture_velocity, time_steps, dt);
    
    println!("Simulated wavefield for historical earthquake: {:?}", wavefield);
}

fn main() {
    historical_earthquake_simulation();
}

In the historical earthquake simulation, the fault plane is modeled with a non-uniform slip distribution where the central region experiences greater slip, representing the area of maximum energy release during the rupture. The simulation updates the fault plane over time and computes the resulting wavefield accordingly. This approach illustrates how variations in rupture dynamics—such as changes in rupture velocity, slip magnitude, and stress drop—can influence the characteristics of the seismic waves generated during an earthquake.

Rust’s high performance and efficient concurrency make it well suited for large-scale earthquake source modeling. By integrating more advanced numerical methods for wave propagation, such as finite difference or spectral element methods, one can simulate realistic seismic scenarios that account for complex fault geometries and dynamic rupture processes. These models form the basis for studying earthquake mechanics, assessing seismic hazards, and developing improved early warning systems.

51.8. Seismic Hazard Assessment and Mitigation link

Seismic hazard assessment is a critical discipline that seeks to quantify both the likelihood of earthquake occurrence and the severity of ground shaking at specific locations. This field employs both probabilistic and deterministic approaches to evaluate seismic risks. In probabilistic seismic hazard analysis, one estimates the probability of various levels of ground motion occurring over a given period by incorporating earthquake recurrence rates, fault activity, and historical seismicity. Deterministic seismic hazard analysis, on the other hand, focuses on a particular earthquake scenario—often a worst-case event—and computes the expected ground motion based on that scenario. Ground motion prediction models, which estimate parameters such as peak ground acceleration and velocity, play an integral role in both approaches by taking into account earthquake magnitude, distance from the fault, and local soil conditions.

The influence of local site conditions is particularly significant in hazard assessment. Variations in soil type, rock stiffness, and other geological characteristics can amplify or attenuate seismic waves, leading to different levels of ground shaking even over short distances. For example, soft soils generally amplify seismic waves, whereas bedrock tends to diminish their amplitude. Such site-specific responses are assessed using computational simulations, which are essential for designing structures that meet modern building codes and for implementing effective mitigation strategies.

The following Rust implementations illustrate key components of seismic hazard assessment. In the first code example, a simplified ground motion prediction model estimates ground shaking based on earthquake magnitude, distance from the fault, and a site factor that represents local amplification effects. The second example simulates site response by adjusting ground motion values according to the soil type at various locations. Finally, a third code example demonstrates how mitigation strategies, such as base isolation systems, can be modeled to reduce the effective ground motion experienced by structures.

extern crate ndarray;
use ndarray::Array1;

/// A simplified ground motion prediction model that estimates the level of ground shaking
/// based on earthquake magnitude, distance from the fault, and a site factor representing
/// local amplification or attenuation. The model computes a base motion using a logarithmic
/// scaling of magnitude and then adjusts the value based on distance and site conditions.
///
/// # Arguments
/// * `magnitude` - The earthquake magnitude.
/// * `distance` - The distance from the fault in kilometers.
/// * `site_factor` - A multiplicative factor representing local site conditions.
///
/// # Returns
/// The estimated ground motion value.
fn ground_motion_model(magnitude: f64, distance: f64, site_factor: f64) -> f64 {
    // Calculate base motion using a simplified scaling law.
    let base_motion = 10f64.powf((magnitude - 6.0) / 2.0) / (distance + 1.0);
    base_motion * site_factor
}

/// Simulate ground motion values for an array of distances and corresponding site factors.
/// This function applies the ground motion model to each pair of distance and site factor,
/// returning an array of predicted ground motion values.
///
/// # Arguments
/// * `magnitude` - The earthquake magnitude.
/// * `distances` - An array of distances from the fault.
/// * `site_factors` - An array of site factors corresponding to each distance.
///
/// # Returns
/// An array of predicted ground motion values.
fn simulate_ground_motion(magnitude: f64, distances: &Array1<f64>, site_factors: &Array1<f64>) -> Array1<f64> {
    let mut ground_motions = Array1::<f64>::zeros(distances.len());
    for (i, &distance) in distances.iter().enumerate() {
        let site_factor = site_factors[i];
        ground_motions[i] = ground_motion_model(magnitude, distance, site_factor);
    }
    ground_motions
}

fn main() {
    let magnitude = 7.0; // Example earthquake magnitude.
    let distances = Array1::<f64>::from(vec![5.0, 10.0, 20.0, 50.0]); // Distances from the fault in kilometers.
    let site_factors = Array1::<f64>::from(vec![1.2, 1.0, 0.8, 0.6]); // Site factors reflecting local conditions.
    
    // Simulate ground motions based on the provided parameters.
    let ground_motions = simulate_ground_motion(magnitude, &distances, &site_factors);
    println!("Simulated ground motions: {:?}", ground_motions);
}

Local site conditions can dramatically alter the amplitude of seismic waves. The following code simulates site response by applying an amplification or attenuation factor based on the type of soil or rock present at each location. For instance, soft soils can significantly amplify ground motion, whereas rock sites typically lead to reduced shaking.

extern crate ndarray;
use ndarray::Array1;

/// Adjust the ground motion value based on the local soil type. This function applies
/// a multiplicative factor to the ground motion depending on whether the site is classified
/// as rock, soft soil, or stiff soil.
///
/// # Arguments
/// * `ground_motion` - The original ground motion value.
/// * `soil_type` - A string slice representing the soil type.
///
/// # Returns
/// The adjusted ground motion value.
fn site_amplification(ground_motion: f64, soil_type: &str) -> f64 {
    match soil_type {
        "rock" => ground_motion * 0.8,     // Rock tends to attenuate ground motion.
        "soft soil" => ground_motion * 1.5,  // Soft soil amplifies seismic waves.
        "stiff soil" => ground_motion * 1.2, // Stiff soil provides moderate amplification.
        _ => ground_motion,                // Default: no adjustment.
    }
}

/// Simulate site response by applying site amplification factors to an array of ground motion values.
/// Each ground motion value is adjusted according to the corresponding soil type provided.
///
/// # Arguments
/// * `ground_motions` - An array of original ground motion values.
/// * `soil_types` - A vector of string slices representing the soil type at each location.
///
/// # Returns
/// An array of ground motion values after applying site-specific amplification.
fn simulate_site_response(ground_motions: &Array1<f64>, soil_types: &Vec<&str>) -> Array1<f64> {
    let mut amplified_motions = Array1::<f64>::zeros(ground_motions.len());
    for (i, &ground_motion) in ground_motions.iter().enumerate() {
        let soil_type = soil_types[i];
        amplified_motions[i] = site_amplification(ground_motion, soil_type);
    }
    amplified_motions
}

fn main() {
    let ground_motions = Array1::<f64>::from(vec![0.2, 0.4, 0.6, 1.0]); // Example ground motions.
    let soil_types = vec!["rock", "stiff soil", "soft soil", "rock"]; // Site classifications for each location.
    
    // Simulate the effect of local site conditions on ground motion.
    let amplified_motions = simulate_site_response(&ground_motions, &soil_types);
    println!("Amplified ground motions: {:?}", amplified_motions);
}

Seismic hazard mitigation involves applying engineering solutions to reduce the impact of seismic shaking on structures and populations. One common approach is the implementation of base isolation systems that decouple a building from ground motion, thereby reducing the seismic forces transmitted to the structure. The following code provides a basic simulation of the effectiveness of a base isolation system by reducing the ground motion according to an isolation factor.

extern crate ndarray;
use ndarray::Array1;

/// Simulate the reduction in ground motion due to a base isolation system. The function reduces
/// the original ground motion by a factor determined by the isolation efficiency.
///
/// # Arguments
/// * `ground_motion` - The original ground motion value.
/// * `isolation_factor` - A factor representing the effectiveness of the base isolation system,
///                        where a higher value indicates greater reduction.
///
/// # Returns
/// The ground motion value after applying the base isolation reduction.
fn base_isolation_reduction(ground_motion: f64, isolation_factor: f64) -> f64 {
    ground_motion * (1.0 - isolation_factor)
}

/// Simulate the effect of mitigation measures on an array of ground motion values. Each value is
/// reduced according to the corresponding isolation factor provided in the input array.
///
/// # Arguments
/// * `ground_motions` - An array of original ground motion values.
/// * `isolation_factors` - An array of isolation factors for each location.
///
/// # Returns
/// An array of ground motion values after applying mitigation measures.
fn simulate_mitigation(ground_motions: &Array1<f64>, isolation_factors: &Array1<f64>) -> Array1<f64> {
    let mut mitigated_motions = Array1::<f64>::zeros(ground_motions.len());
    for (i, &ground_motion) in ground_motions.iter().enumerate() {
        let isolation_factor = isolation_factors[i];
        mitigated_motions[i] = base_isolation_reduction(ground_motion, isolation_factor);
    }
    mitigated_motions
}

fn main() {
    let ground_motions = Array1::<f64>::from(vec![0.5, 0.7, 1.0, 1.2]); // Example ground motion values.
    let isolation_factors = Array1::<f64>::from(vec![0.3, 0.5, 0.2, 0.4]); // Corresponding isolation factors.
    
    // Simulate the effect of base isolation systems on reducing seismic ground motion.
    let mitigated_motions = simulate_mitigation(&ground_motions, &isolation_factors);
    println!("Mitigated ground motions: {:?}", mitigated_motions);
}

Seismic hazard assessment and mitigation encompass the analysis of earthquake ground motions, the evaluation of site-specific amplification effects, and the application of mitigation strategies such as base isolation. The implementations provided above, written in Rust, demonstrate how computational models can be used to simulate earthquake scenarios, predict ground shaking intensities, assess site responses, and evaluate the benefits of seismic mitigation measures. These tools are essential for informing building codes and designing structures that can better withstand the dynamic forces generated by earthquakes.

51.9. Case Studies and Applications link

Seismic wave propagation modeling is indispensable for many real-world applications ranging from earthquake engineering to geophysical exploration and environmental monitoring. These models enable scientists and engineers to interpret seismic data, image subsurface structures, monitor volcanic activity, and even detect landslides. In earthquake engineering, detailed simulations help in understanding the impact of seismic waves on structures and in designing buildings that can withstand the forces generated by earthquakes. In the realm of geophysical exploration, seismic models are used to image subsurface layers by analyzing the reflections of waves generated at the surface, thereby identifying potential oil and gas reservoirs. Environmental monitoring also benefits from seismic modeling through the tracking of volcanic activity and landslide detection. The following examples demonstrate how seismic models can be implemented in Rust to address these diverse applications, with emphasis on performance optimization, data analysis, and interpretation.

Seismic imaging of subsurface structures is a key technique in geophysical exploration. Seismic waves generated at the surface travel through the Earth's layers and are reflected back from interfaces where there are changes in material properties. By analyzing these reflections, one can map subsurface features. The following code implements a basic simulation of seismic wave propagation in a layered medium. The wavefield is updated using a simplified update rule based on the wave velocity and time step, and the model is applicable to the exploration of oil and gas reservoirs.

extern crate ndarray;
use ndarray::Array2;

/// Simulate seismic wave propagation in a layered medium.
/// This function iteratively updates a wavefield based on a simplified wave equation.
/// The wavefield is influenced by the local wave velocity from the provided velocity model.
/// 
/// # Arguments
/// * `layers` - A 2D array representing the geological layers (can be used for additional properties).
/// * `wave_velocity` - A 2D array representing the wave velocity at each grid point (in km/s).
/// * `time_steps` - The number of time iterations for the simulation.
/// * `dx` - The spatial step size (in km).
/// * `dt` - The time step size (in seconds).
/// 
/// # Returns
/// A 2D array representing the simulated seismic wavefield after propagation.
fn seismic_wave_propagation(
    layers: &Array2<f64>,
    wave_velocity: &Array2<f64>,
    time_steps: usize,
    dx: f64,
    dt: f64,
) -> Array2<f64> {
    // Initialize the wavefield with the same dimensions as the layers.
    let mut wavefield = Array2::<f64>::zeros(layers.raw_dim());

    // Loop over each time step to simulate wave propagation.
    for _ in 0..time_steps {
        // Iterate over each grid point.
        for ((i, j), value) in wavefield.indexed_iter_mut() {
            // Retrieve the local wave velocity.
            let velocity = wave_velocity[[i, j]];
            // Update the wavefield based on a simplified propagation step.
            // This is a basic approximation where the wave amplitude is increased by a term proportional to velocity, dt, and dx.
            *value += velocity * dt * dx;
        }
    }
    wavefield
}

fn main() {
    let nx = 100;
    let ny = 100;
    let time_steps = 500;
    let dx = 1.0;   // Spatial resolution of 1 km
    let dt = 0.01;  // Time step of 0.01 s

    // Define a uniform subsurface layer.
    let layers = Array2::<f64>::from_elem((nx, ny), 1.0);
    // Define a uniform wave velocity model (for example, 2 km/s throughout the domain).
    let wave_velocity = Array2::<f64>::from_elem((nx, ny), 2.0);

    // Run the seismic wave propagation simulation.
    let wavefield = seismic_wave_propagation(&layers, &wave_velocity, time_steps, dx, dt);
    println!("Final wavefield after propagation: {:?}", wavefield);
}

For large-scale seismic simulations, computational performance is crucial. Rust’s concurrency features, such as those provided by the rayon crate, enable parallel execution of computation-intensive tasks. The following example demonstrates how to parallelize the seismic wave propagation simulation using rayon. The code updates the wavefield concurrently, thus reducing overall computation time for large grids or long time simulations.

extern crate rayon;
extern crate ndarray;
extern crate rand;

use ndarray::Array2;
use rand::Rng;
use rayon::prelude::*;

/// Parallel seismic wave propagation in a layered medium.
/// This function leverages Rayon to update the wavefield concurrently across spatial grid points for each time step.
/// It applies a simplified update rule where each grid point's amplitude is incremented based on the time step and spatial step size.
///
/// # Arguments
/// * `layers` - A 2D array representing the subsurface layers.
/// * `wave_velocity` - A 2D array representing the wave velocity at each grid point.
/// * `time_steps` - The number of time iterations for the simulation.
/// * `dx` - The spatial step size.
/// * `dt` - The time step size.
///
/// # Returns
/// A 2D array representing the final wavefield after parallel propagation.
fn parallel_wave_propagation(
    layers: &Array2<f64>,
    wave_velocity: &Array2<f64>,
    time_steps: usize,
    dx: f64,
    dt: f64,
) -> Array2<f64> {
    // Initialize the wavefield with zeros matching the dimensions of the layers.
    let mut wavefield = Array2::<f64>::zeros(layers.raw_dim());

    // Sequentially iterate over each time step.
    for _ in 0..time_steps {
        // Update each grid point in parallel by obtaining a mutable slice and using Rayon.
        if let Some(slice) = wavefield.as_slice_mut() {
            slice.par_iter_mut().for_each(|x| {
                *x += dt * dx; // Simplified update: increment by dt * dx (for demonstration purposes).
            });
        }
    }

    wavefield
}

fn main() {
    let nx = 100;
    let ny = 100;
    let time_steps = 500;
    let dx = 1.0;
    let dt = 0.01;

    // Define a uniform subsurface layer.
    let layers = Array2::<f64>::from_elem((nx, ny), 1.0);
    // Define a uniform velocity model, for example, 2 km/s.
    let wave_velocity = Array2::<f64>::from_elem((nx, ny), 2.0);

    // Execute the parallel wave propagation simulation.
    let wavefield = parallel_wave_propagation(&layers, &wave_velocity, time_steps, dx, dt);
    println!("Final wavefield after parallel propagation: {:?}", wavefield);
}

Seismic simulations are also instrumental in monitoring volcanic activity. By simulating seismic wave propagation in volcanic regions, one can track magma movement and identify precursory signals of volcanic eruptions. The following code simulates wave propagation in a volcanic region where the wave velocity is higher due to active magma movement. This model can be used as a basis for developing monitoring systems that provide early warnings of eruptive activity.

extern crate ndarray;
use ndarray::Array2;

/// Simulate seismic wave propagation in a volcanic region.
/// This function updates the wavefield in a region where both the subsurface properties and the wave velocity
/// may be affected by volcanic activity. The update rule used here is simplified for demonstration purposes.
/// 
/// # Arguments
/// * `volcano_region` - A 2D array representing the geological properties of the volcanic region.
/// * `wave_velocity` - A 2D array representing the local wave velocity, which may be higher in active volcanic zones.
/// * `time_steps` - The number of time iterations for the simulation.
/// * `dx` - The spatial step size.
/// * `dt` - The time step size.
/// 
/// # Returns
/// A 2D array representing the simulated seismic wavefield in the volcanic region.
fn volcanic_wave_propagation(
    volcano_region: &Array2<f64>,
    wave_velocity: &Array2<f64>,
    time_steps: usize,
    dx: f64,
    dt: f64,
) -> Array2<f64> {
    let mut wavefield = Array2::<f64>::zeros(volcano_region.raw_dim());

    // Simulate wave propagation over the given time steps.
    for _ in 0..time_steps {
        for ((i, j), value) in wavefield.indexed_iter_mut() {
            let velocity = wave_velocity[[i, j]];
            *value += velocity * dt * dx; // Simplified propagation update.
        }
    }
    wavefield
}

fn main() {
    let nx = 100;
    let ny = 100;
    let time_steps = 500;
    let dx = 1.0;
    let dt = 0.01;

    // Define a volcanic region where the geological property may indicate higher activity.
    let volcano_region = Array2::<f64>::from_elem((nx, ny), 1.5);
    // Define a wave velocity model that represents faster propagation in the volcanic area (e.g., 3 km/s).
    let wave_velocity = Array2::<f64>::from_elem((nx, ny), 3.0);

    // Run the volcanic wave propagation simulation.
    let wavefield = volcanic_wave_propagation(&volcano_region, &wave_velocity, time_steps, dx, dt);
    println!("Final wavefield in volcanic region: {:?}", wavefield);
}

In these case studies, seismic wave propagation models are applied to diverse real-world scenarios including subsurface imaging for geophysical exploration, performance-optimized large-scale simulations through parallel computing, and monitoring of volcanic regions to track magma movement. These examples underscore the versatility and power of computational physics in addressing complex seismic problems. Rust's efficiency, memory safety, and robust concurrency capabilities make it particularly well suited for developing high-performance seismic models that are critical for both industry applications and advanced research.

51.10. Conclusion link

Chapter 51 of CPVR equips readers with the knowledge and tools to explore seismic wave propagation using Rust. By integrating mathematical models, numerical simulations, and seismic inversion techniques, this chapter provides a robust framework for understanding the complexities of seismic wave phenomena. Through hands-on examples and real-world applications, readers are encouraged to leverage computational tools to advance geophysical research and contribute to seismic hazard mitigation efforts.

51.10.1. Further Learning with GenAI link

The following prompts are designed to help learners explore the complexities of seismic wave propagation using Rust. These prompts focus on fundamental concepts, mathematical models, computational techniques, and practical applications related to seismology.

• Discuss the significance of seismic wave propagation in understanding Earth’s interior. How do different types of seismic waves, including primary (P-waves), secondary (S-waves), and surface waves (Rayleigh and Love waves), provide detailed and complementary information about the Earth's subsurface structure? Include a discussion on how wave velocities, attenuation, and phase shifts reveal the composition and behavior of Earth's layers, such as the crust, mantle, and core.

• Explain the role of mathematical models in simulating seismic wave propagation. How do wave equations, elastodynamic principles, and the application of boundary and initial conditions contribute to the accuracy and realism of seismic simulations? Explore the challenges in solving partial differential equations (PDEs) in both homogeneous and heterogeneous media and their impact on capturing realistic seismic events in simulations.

• Analyze the importance of numerical methods in solving seismic wave propagation problems. Compare and contrast the finite difference, finite element, and spectral element methods in terms of their accuracy, stability, convergence, and computational efficiency. How do these methods address the trade-offs between resolution and computational cost in large-scale seismic simulations of complex geological environments?

• Explore the application of seismic wave attenuation and dispersion models. How do intrinsic material properties, such as density, viscosity, and elastic moduli, along with geological structures like faults and anisotropies, influence the attenuation and dispersion of seismic waves? Discuss the modeling of anelasticity and frequency-dependent attenuation in simulating realistic seismic scenarios in Rust.

• Discuss the principles of seismic wave reflection and refraction at material interfaces. How do impedance contrasts, critical angles, and Snell’s law govern the behavior of reflected and refracted seismic waves at geological boundaries? Explain the significance of mode conversions (P-to-S and S-to-P) and total internal reflection in subsurface imaging and how Rust can model these phenomena.

• Investigate the use of seismic tomography and inversion techniques in imaging subsurface structures. How do forward modeling and inversion algorithms, such as full-waveform inversion (FWI) and travel-time tomography, contribute to reconstructing subsurface properties? Discuss the challenges of ill-posed inverse problems and how regularization techniques enhance the stability and resolution of the reconstructed models.

• Explain the significance of earthquake source modeling in understanding seismic wave generation. How do source parameters, such as stress drop, fault slip, rupture velocity, and the geometry of the fault, influence the characteristics and propagation of seismic waves? Explore the role of kinematic and dynamic source models in simulating the generation of seismic waves from real-world earthquakes.

• Discuss the role of seismic hazard assessment in earthquake engineering. How do computational models predict ground motion and assess seismic risks in urban environments? Examine the integration of probabilistic seismic hazard analysis (PSHA) with site-specific response analysis to inform the design of earthquake-resistant structures, and how Rust can enhance these simulations.

• Analyze the challenges of simulating seismic wave propagation in heterogeneous media. How does geological complexity, including anisotropy, fault zones, and heterogeneity in material properties, affect the accuracy of seismic wave simulations? Discuss the computational approaches to handling these complexities in 2D and 3D domains, and how Rust’s features contribute to optimizing these simulations.

• Explore the use of Rust in implementing seismic wave propagation models. How can Rust’s performance optimization, memory safety, and concurrency support be leveraged to implement and accelerate large-scale seismic wave propagation models? Discuss examples of Rust-based parallel computing frameworks or libraries that can enhance both simulation efficiency and scalability.

• Discuss the importance of seismic wave dispersion in identifying subsurface anomalies. How do frequency-dependent dispersion patterns of seismic waves provide insights into the composition, fluid saturation, and porosity of geological structures? Analyze how computational models simulate dispersive behavior to detect subsurface anomalies, such as oil reservoirs or geothermal sources.

• Investigate the application of numerical methods in simulating seismic wave propagation in 2D and 3D domains. How do grid generation techniques, time-stepping algorithms, and boundary condition implementations (e.g., absorbing, free-surface) impact the accuracy and efficiency of seismic simulations? Explore how Rust can be used to optimize these methods for large-scale seismic simulations.

• Explain the principles of seismic wave attenuation and how it affects seismic signal analysis. How do absorption, scattering, and anelasticity contribute to the attenuation of seismic waves as they propagate through different materials? Discuss how modeling attenuation in Rust can improve the interpretation of seismic signals, particularly for earthquake detection and subsurface exploration.

• Discuss the challenges of seismic wave inversion and the trade-offs between resolution and stability. How do regularization techniques, iterative solvers, and forward modeling strategies address the inherent trade-offs between spatial resolution and stability in seismic inversion? Analyze how inversion methods are implemented in Rust and their role in reconstructing high-fidelity subsurface models.

• Analyze the impact of seismic wave reflection and refraction on subsurface imaging. How do computational models simulate wave behavior at geological boundaries to improve the accuracy of subsurface exploration techniques, such as reflection seismology and refraction tomography? Explore how Rust can be used to enhance these models and their application in oil, gas, and mineral exploration.

• Explore the role of earthquake source modeling in predicting ground shaking intensity. How do computational models simulate the propagation of seismic waves from fault ruptures to assess ground shaking intensity and earthquake hazards? Examine the importance of source modeling in hazard prediction and mitigation, and how Rust can contribute to implementing scalable earthquake simulation frameworks.

• Discuss the application of seismic tomography in monitoring volcanic activity. How do seismic waves generated by volcanic activity reveal information about magma movement, pressure buildup, and potential eruptions? Analyze how real-time seismic data, combined with tomography techniques, can monitor volcanic processes and how Rust-based tools can facilitate these simulations.

• Investigate the use of Rust-based tools in automating seismic wave propagation simulations. How can workflow automation in Rust improve the efficiency, scalability, and reproducibility of seismic simulations? Explore specific Rust libraries or tools that support parallel processing, data visualization, and real-time simulation in seismology.

• Explain the significance of seismic hazard assessment in urban planning. How do computational models help identify high-risk areas and inform the design of earthquake-resistant infrastructure in cities? Discuss the role of site-specific hazard analysis, ground motion prediction, and mitigation strategies in reducing earthquake risks, and how Rust can optimize these workflows.

• Reflect on the future trends in seismic wave propagation research and the potential developments in computational techniques. How might advancements in Rust’s capabilities (such as support for parallelism, GPU computing, or machine learning integration) address emerging challenges in seismic wave propagation? Explore potential breakthroughs in seismic simulation techniques and their applications in geophysics, earthquake engineering, and hazard mitigation.

By engaging with these topics, you are building a strong foundation in seismology and equipping yourself with the tools to contribute to cutting-edge research and innovation. Embrace the challenges, stay curious, and let your exploration of seismic wave propagation inspire you to push the boundaries of what is possible in this dynamic field.

51.10.2. Assignments for Practice link

These self-exercises are designed to provide you with practical experience in simulating seismic wave propagation using Rust. By engaging with these exercises, you will develop a deep understanding of both the theoretical concepts and the computational methods necessary to model, analyze, and interpret seismic wave phenomena.

Exercise 51.1: Implementing Finite Difference Methods for Seismic Wave Propagation link

• Objective: Develop a Rust program to simulate seismic wave propagation using finite difference methods (FDM).

• Steps:

• Begin by researching the finite difference method (FDM) and its application in solving the wave equation for seismic wave propagation. Write a brief summary explaining the significance of FDM in computational seismology.

• Implement a Rust program that simulates seismic wave propagation in a 2D or 3D domain using FDM, including the setup of grid points, time-stepping schemes, and boundary conditions.

• Analyze the simulation results to identify wave patterns, reflection, refraction, and attenuation. Visualize the wavefield evolution over time and discuss the implications for understanding seismic wave behavior.

• Experiment with different grid resolutions, time steps, and material properties to explore their impact on the simulation accuracy and stability. Write a report summarizing your findings and discussing the challenges in simulating seismic waves using FDM.

• GenAI Support: Use GenAI to optimize the implementation of finite difference methods, troubleshoot issues in simulating seismic wave propagation, and interpret the results in the context of seismic simulations.

Exercise 51.2: Simulating Seismic Wave Reflection and Refraction at Geological Boundaries link

• Objective: Implement a Rust-based simulation to model seismic wave reflection and refraction at material interfaces, focusing on the impact of impedance contrast and critical angles.

• Steps:

• Study the principles of seismic wave reflection and refraction and their relevance in subsurface exploration. Write a brief explanation of how impedance contrast and critical angles influence wave behavior at geological boundaries.

• Implement a Rust program that simulates seismic wave propagation through multiple layers with varying material properties, including the calculation of reflection and transmission coefficients at each interface.

• Analyze the simulation results to identify reflection and refraction patterns, critical angles, and mode conversions. Visualize the wave paths and discuss the implications for subsurface imaging and exploration.

• Experiment with different material properties, layer thicknesses, and wave frequencies to explore their impact on reflection and refraction patterns. Write a report detailing your findings and discussing strategies for optimizing subsurface imaging using seismic waves.

• GenAI Support: Use GenAI to guide the implementation of reflection and refraction models, optimize the simulation of wave behavior at geological boundaries, and interpret the results in the context of subsurface exploration.

Exercise 51.3: Modeling Seismic Wave Attenuation and Dispersion in Heterogeneous Media link

• Objective: Use Rust to implement models of seismic wave attenuation and dispersion, focusing on the effects of material heterogeneity on wave propagation.

• Steps:

• Research the principles of seismic wave attenuation and dispersion, including the physical mechanisms that contribute to these phenomena. Write a brief summary explaining the significance of attenuation and dispersion in seismic signal analysis.

• Implement a Rust-based simulation that models seismic wave attenuation and dispersion in a heterogeneous medium, including the simulation of wave scattering, absorption, and anelasticity.

• Analyze the simulation results to identify the impact of heterogeneity on wave amplitude, phase, and velocity. Visualize the attenuation and dispersion effects and discuss the implications for seismic data interpretation and subsurface characterization.

• Experiment with different material heterogeneities, wave frequencies, and attenuation models to explore their impact on the simulation results. Write a report summarizing your findings and discussing strategies for improving seismic data analysis using attenuation and dispersion models.

• GenAI Support: Use GenAI to optimize the implementation of attenuation and dispersion models, troubleshoot issues in simulating wave propagation in heterogeneous media, and interpret the results in the context of seismic signal analysis.

Exercise 51.4: Implementing Seismic Tomography for Subsurface Imaging link

• Objective: Develop a Rust-based seismic tomography program to image subsurface structures by analyzing the travel times of seismic waves.

• Steps:

• Study the principles of seismic tomography and its application in imaging the Earth's interior. Write a brief explanation of how travel-time tomography reconstructs subsurface velocity structures.

• Implement a Rust program that performs seismic tomography, including the forward modeling of seismic wave travel times, inversion of travel-time data, and reconstruction of subsurface velocity models.

• Analyze the tomography results to identify subsurface anomalies, velocity contrasts, and geological structures. Visualize the reconstructed velocity model and discuss the implications for subsurface exploration and seismic hazard assessment.

• Experiment with different inversion algorithms, regularization techniques, and data sets to explore their impact on the resolution and stability of the tomography results. Write a report detailing your findings and discussing strategies for optimizing seismic tomography.

• GenAI Support: Use GenAI to guide the implementation of seismic tomography algorithms, optimize the inversion process, and interpret the results in the context of subsurface imaging.

Exercise 51.5: Case Study - Simulating Seismic Waves from an Earthquake Source link

• Objective: Apply computational methods to simulate seismic wave propagation from an earthquake source, focusing on modeling fault rupture and ground motion.

• Steps:

• Select an earthquake scenario and research the principles of earthquake source modeling, including fault mechanics and the generation of seismic waves. Write a summary explaining the importance of source modeling in seismic hazard assessment.

• Implement a Rust-based simulation to model seismic wave propagation from a fault rupture, including the calculation of stress drop, fault slip, rupture velocity, and wave radiation patterns.

• Analyze the simulation results to evaluate ground motion intensity, wave propagation patterns, and the impact of source parameters on seismic wave characteristics. Visualize the ground shaking distribution and discuss the implications for earthquake risk assessment.

• Experiment with different source parameters, fault geometries, and rupture velocities to explore their impact on ground motion and seismic hazard. Write a detailed report summarizing your approach, the simulation results, and the implications for seismic hazard mitigation.

• GenAI Support: Use GenAI to guide the design and implementation of earthquake source models, optimize seismic wave propagation simulations, and help interpret the results in the context of seismic hazard assessment.

Each exercise offers an opportunity to explore the complexities of seismology, experiment with advanced simulations, and contribute to the development of new insights and technologies in earthquake science and hazard mitigation. Embrace the challenges, push the boundaries of your knowledge, and let your passion for computational geophysics drive you toward mastering the art of seismic wave modeling. Your efforts today will lead to breakthroughs that shape the future of seismology and earthquake engineering.