Computer-Code Programming Algorithm of A River (A River by Nature for Water Flow)

Simulating a natural river digitally relies on modeling three dynamic forces: gravity-driven gradient flow, hydraulic erosion, and sediment transport.

Core Mathematical Framework

A standard particle-based hydraulic simulation uses discrete water droplet iterations across a 2D heightmap H(x, y):

 * Gradient Vector (\nabla H): Determines direction of steepest descent for droplet position (x, y):

   

 * Erosion & Deposition:

   * Transport Capacity (C): C = K_c \cdot v \cdot s \cdot W (where v is velocity, s is slope, W is water volume, K_c is capacity constant).

   * If sediment carried S < C, erode terrain H by \Delta h = K_e (C - S).

   * If S > C, deposit sediment onto terrain H by \Delta h = K_d (S - C).

 * Inertia & Velocity Update:

   

Python Implementation (Droplet-Based Hydraulic Simulation)

import numpy as np

def simulate_river_flow(heightmap, num_droplets=10000, max_steps=64):

    """

    Simulates natural river formation on a heightmap using particle erosion.

    """

    rows, cols = heightmap.shape

    

    # Simulation Parameters

    inertia = 0.05 # Resistance to direction changes

    capacity_coeff = 4.0 # Sediment carrying capacity scaling

    erode_speed = 0.3 # Terrain erosion rate

    deposit_speed = 0.3 # Sediment deposition rate

    evaporation = 0.02 # Water volume decrease rate

    gravity = 9.81

    

    for _ in range(num_droplets):

        # Spawn droplet at random location

        px, py = np.random.uniform(1, cols - 2), np.random.uniform(1, rows - 2)

        dir_x, dir_y = 0.0, 0.0

        speed, water, sediment = 1.0, 1.0, 0.0

        

        for step in range(max_steps):

            ix, iy = int(px), int(py)

            

            # 1. Calculate local terrain gradient via central differences

            gx = heightmap[iy, ix + 1] - heightmap[iy, ix - 1]

            gy = heightmap[iy + 1, ix] - heightmap[iy - 1, ix]

            

            # 2. Update direction with momentum/inertia

            dir_x = dir_x * inertia - gx * (1.0 - inertia)

            dir_y = dir_y * inertia - gy * (1.0 - inertia)

            

            length = np.hypot(dir_x, dir_y)

            if length == 0:

                break

            dir_x /= length

            dir_y /= length

            

            # 3. Move particle to new position

            new_px, new_py = px + dir_x, py + dir_y

            if not (1 <= new_px < cols - 2 and 1 <= new_py < rows - 2):

                break

                

            # 4. Calculate elevation differential

            h_old = heightmap[iy, ix]

            h_new = heightmap[int(new_py), int(new_px)]

            delta_h = h_new - h_old

            

            # 5. Calculate carrying capacity

            # Capacity increases with slope and speed, scales with water volume

            slope = max(-delta_h, 0.0001)

            capacity = max(slope * speed * water * capacity_coeff, 0.01)

            

            # 6. Erode or Deposit sediment

            if sediment > capacity or delta_h > 0:

                # Flow uphill or over-capacity: deposit sediment

                amount_to_deposit = (sediment - capacity) * deposit_speed if delta_h <= 0 else min(sediment, delta_h)

                sediment -= amount_to_deposit

                heightmap[iy, ix] += amount_to_deposit

            else:

                # Accelerating downhill: erode terrain

                amount_to_erode = min((capacity - sediment) * erode_speed, slope)

                sediment += amount_to_erode

                heightmap[iy, ix] -= amount_to_erode

                

            # 7. Update velocity and water volume

            speed = np.sqrt(max(0.0, speed**2 + delta_h * gravity))

            water *= (1.0 - evaporation)

            px, py = new_px, new_py

    return heightmap