Model Seismic Waves in Layered Soil Using Matter.js

This article explains how to simulate seismic wave propagation through stratified soil layers using a discrete particle-lattice approach in the 2D physics engine Matter.js. By representing soil strata as an interconnected network of mass particles and viscoelastic constraints, you can model how seismic body and shear waves refract, attenuate, and amplify across varying geological media. The following guide covers domain setup, layer parameterization, boundary conditions, and seismic excitation.

Representing Soil as a Discrete Element Network

Matter.js is a rigid-body physics engine, so continuous soil mechanics must be discretized using a bonded-particle or lattice-spring model. In this approach, soil mass is concentrated at nodal particles (Matter.Bodies.circle), while inter-particle forces (elasticity and damping) are handled via distance constraints (Matter.Constraint.create).

To model different soil strata—such as bedrock, dense sand, and soft clay—you vary three fundamental physical parameters across vertical depth:

• Mass/Density: Assigned to individual particles via density.
• Stiffness: Controlled by the stiffness property of the connecting constraints, proportional to the soil layer's shear modulus (\(G\)) and elastic modulus (\(E\)).
• Damping: Controlled by constraint damping, representing material damping and energy dissipation.

Constructing Layered Strata

Generate the domain on a regular triangular or hexagonal grid to ensure isotropic wave propagation. Divide your domain vertically into distinct strata:

1. Bedrock Layer (Bottom): High particle mass, maximum constraint stiffness (e.g., stiffness: 0.9 to 1.0), and low damping.
2. Intermediate Layer (e.g., Dense Sand/Gravel): Moderate mass, moderate constraint stiffness (e.g., stiffness: 0.4 to 0.6).
3. Near-Surface Layer (e.g., Soft Clay/Alluvium): Lower mass, low constraint stiffness (e.g., stiffness: 0.05 to 0.2), and higher damping.

Each particle is connected to its nearest neighbors using cross-braced constraints to support both compressional (P) waves and shear (S) waves. Setting collisionFilter.group = -1 on all soil particles disables default contact collisions between bonded neighbors, allowing the constraints alone to dictate the linear elastic response without erratic contact chatter.

Applying Non-Reflecting Boundary Conditions

A common issue in computational geomechanics is the reflection of waves off domain edges. To mitigate artificial boundary reflections:

• Lateral Boundaries (Left and Right): Attach edge particles to fixed anchor bodies using viscous dampers (constraints with high damping and zero rest-length deviation) to approximate absorbing boundaries.
• Top Surface: Leave completely unconstrained to form a free boundary, allowing Rayleigh wave generation and natural surface amplification.
• Bottom Boundary: Kept mobile horizontally or vertically depending on the type of wave being induced.

Simulating Seismic Input

Seismic excitation is introduced by driving the base layer particles. Instead of fixing the bottom particles in place, manipulate their kinematics directly in each simulation tick of the beforeUpdate event:

• Harmonic Shear Waves (S-Waves): Modulate the horizontal velocity of the bedrock particles using a sine function:

  const velocityX = amplitude * Math.sin(2 * Math.PI * frequency * engine.timing.timestamp);
  Matter.Body.setVelocity(baseParticle, { x: velocityX, y: 0 });

• Compressional Waves (P-Waves): Apply the harmonic signal to the vertical velocity (y component) instead.
• Earthquake Time-Histories: Read an external accelerogram file (e.g., acceleration data sampled at 0.01-second intervals), integrate to velocity, and set the base particle velocities accordingly.

Tuning the Engine for Wave Fidelity

Standard physics engines prioritize speed over numerical precision. To accurately resolve dynamic wave propagation without numerical instability or artificial dispersion:

• Sub-stepping: Increase solver fidelity by increasing engine.positionIterations and engine.velocityIterations to at least 10 or 15.
• Timestep Stability: Keep the simulation delta fixed and small (e.g., \(\Delta t \le 16.6\text{ ms}\)). If stiffness values are near 1.0, reduce the timestep to prevent numerical explosion.
• Wavelength Resolution: Ensure that the shortest physical wavelength expected in your simulation spans at least 8 to 10 particles in the grid to minimize lattice dispersion.