Model Sediment Erosion and Sorting in Matter.js

This guide explains how to simulate fluvial sediment dynamics, riverbed erosion, and grain-size sorting using Matter.js. By coupling Matter.js's 2D rigid-body physics engine with custom hydrodynamic forces—such as fluid drag, buoyancy, and lift—you can model how variable water discharge mobilizes finer grains while leaving heavier gravel behind to form natural armor layers.

1. Representing Sediment Grains

Matter.js lacks native fluid dynamics, so sediment particles are modeled as discrete circular rigid bodies with varying dimensions and masses representing a grain-size distribution (e.g., sand, pebbles, and cobbles).

const { Bodies, Body, Composite } = Matter;

function createSedimentGrain(x, y, radius) {
  // Density scales mass proportional to volume in 2D
  const density = 0.00265; // ~2.65 g/cm³ (quartz equivalent)
  
  return Bodies.circle(x, y, radius, {
    density: density,
    friction: 0.8,         // High friction replicates inter-particle locking
    frictionStatic: 0.9,
    restitution: 0.1,      // Low bounciness for natural gravel packing
    render: {
      fillStyle: radius > 15 ? '#5A5A5A' : radius > 8 ? '#8A8A8A' : '#C2B280'
    }
  });
}

Scatter these grains inside a confined channel formed by static boundary bodies to establish the initial unworked riverbed.

2. Defining the Flow Field Velocity Profile

Water flow velocity in an open channel is not uniform. Due to boundary layer friction, fluid velocity approaches zero near the channel bed and peaks near the water surface. A logarithmic or power-law velocity profile can be approximated based on the vertical coordinate \(y\):

function getFlowVelocity(particleY, riverSurfaceY, bedY, baseFlowRate) {
  if (particleY > bedY) return 0; // Below bed level
  
  // Normalize depth: 0 at bed, 1 at surface
  const relativeDepth = Math.max(0, Math.min(1, (bedY - particleY) / (bedY - riverSurfaceY)));
  
  // Logarithmic-style velocity profile
  return baseFlowRate * Math.sqrt(relativeDepth);
}

Modulating baseFlowRate over time allows you to simulate hydrographs, such as seasonal high flows or flash floods.

3. Calculating Hydrodynamic Forces

During every engine update step, calculate and apply the net fluid forces acting on each grain: buoyancy, drag, and lift.

Buoyancy

Buoyancy counteracts gravity. In Matter.js, subtract a portion of the downward gravitational force based on the displaced fluid density:

\[\vec{F}_{\text{buoyancy}} = -\rho_{\text{fluid}} \cdot \text{Area} \cdot \vec{g}\]

Fluid Drag

Drag is the primary driver of horizontal transport. It depends on the relative velocity between the water and the grain:

\[\vec{F}_{\text{drag}} = \frac{1}{2} \rho_{\text{fluid}} C_d A (\vec{v}_{\text{fluid}} - \vec{v}_{\text{grain}}) |\vec{v}_{\text{fluid}} - \vec{v}_{\text{grain}}|\]

Lift Force

Vertical fluid velocity gradients create upward hydrodynamic lift (Bernoulli/Saffman effect), allowing bed grains to become entrained in the flow:

\[\vec{F}_{\text{lift}} = \frac{1}{2} \rho_{\text{fluid}} C_L A v_{\text{fluid}}^2\]

4. Applying the Force Integration Loop

Hook the force calculations into the beforeUpdate event of your Matter.js engine.

const { Events, Vector } = Matter;

const WATER_DENSITY = 0.001; // Water density relative to grain density
const GRAVITY = 1;          // Match your engine's world.gravity.y
const DRAG_COEFF = 0.47;    // Spherical drag coefficient
const LIFT_COEFF = 0.2;

Events.on(engine, 'beforeUpdate', () => {
  const particles = Composite.allBodies(engine.world).filter(b => !b.isStatic);

  particles.forEach(grain => {
    const r = grain.circleRadius;
    if (!r) return;

    const area = Math.PI * r * r;
    const flowSpeed = getFlowVelocity(grain.position.y, surfaceY, bedY, currentFlowRate);
    
    // 1. Buoyant force (upward)
    const buoyantMagnitude = WATER_DENSITY * area * GRAVITY;
    const fBuoyancy = { x: 0, y: -buoyantMagnitude };

    // 2. Drag force (horizontal)
    const relativeVx = flowSpeed - grain.velocity.x;
    const dragMagnitude = 0.5 * WATER_DENSITY * DRAG_COEFF * (2 * r) * (relativeVx * Math.abs(relativeVx));
    const fDrag = { x: dragMagnitude, y: 0 };

    // 3. Lift force (upward, active under strong local shear)
    let fLift = { x: 0, y: 0 };
    if (grain.position.y > surfaceY) {
      const liftMagnitude = 0.5 * WATER_DENSITY * LIFT_COEFF * (2 * r) * (flowSpeed * flowSpeed);
      fLift = { x: 0, y: -liftMagnitude };
    }

    // Apply combined forces to the particle center
    const totalForce = Vector.add(Vector.add(fBuoyancy, fDrag), fLift);
    Body.applyForce(grain, grain.position, totalForce);
  });
});

5. Simulating Incipient Motion and Shields Criterion

Sediment only mobilizes when fluid drag overcomes static bed friction—a threshold governed physically by the Shields parameter. In Matter.js, this emerges mechanically from the combination of frictionStatic, inter-particle contact normal forces, and grain mass.

• Low Flow Rates: Hydrodynamic drag remains below the static friction threshold of all particles. The bed remains stable.
• Moderate Flow Rates (Winnowing): Drag exceeds the threshold for small, low-mass particles. Fine grains lift and saltate downstream, while heavy cobbles stay stationary.
• High Flow Rates (Pervasive Erosion): High shear stresses overcome the friction of all grain classes, causing bulk channel degradation.

6. Emergence of Sediment Sorting and Armoring

As the simulation runs:

1. Vertical Sorting (Armoring): Fine sediment is swept away from the surface layer, exposing coarse gravel beneath. The coarse gravel protects the underlying subsurface fines from further erosion, accurately simulating bed armoring.
2. Downstream Sorting: Transported grains settle in areas where local flow velocity decreases, such as behind obstacles, in channel depressions, or where the channel widens. Heavy pebbles deposit first, followed by lighter sands further down-gradient.