Flow Routing and Mass Balance in QUAL2K River Segments

As described in the EPA QUAL2K technical documentation, the model divides the river into discrete reaches (computational elements) and marches downstream from the headwater boundary. At each reach, it performs a conservative mass balance to mix upstream flow with point sources and diffuse sources before applying kinetic reactions.

Downstream Flow Routing Algorithm

The simulation proceeds sequentially from the headwater (reach 0) to the terminus (reach n). At each reach $i$ , three steps are performed in order:

Figure 1: Three-step downstream march for each reach

Step 1: Mass Balance by Conservative Mixing

The mass entering reach $i$ from all upstream sources is summed and divided by the total flow to produce mixed concentrations:

C_{mix,k} = \frac{Q_{up} \cdot C_{up,k} + \sum Q_{ps} \cdot C_{ps,k} + \sum Q_{ds} \cdot C_{ds,k}}{Q_{total}}

where:

$Q_{up}$ — upstream flow from reach $i-1$
$C_{up,k}$ — upstream concentration of constituent $k$
$Q_{ps}, C_{ps}$ — point source flow and concentration
$Q_{ds}, C_{ds}$ — diffuse source flow and concentration
$Q_{total} = Q_{up} + \Sigma Q_{ps} + \Sigma Q_{ds} - \Sigma Q_{abs}$

Point Source Injection

Point sources represent discrete discharges such as wastewater treatment plant effluent, tributary inflow, or industrial discharges. Each point source has a location $x_{ps}$ (km) along the river and is assigned to the reach where it falls:

x_{dn,i-1} < x_{ps} \leq x_{dn,i}

Figure 2: Point source and diffuse source injection into river reaches

Diffuse Source Distribution

Diffuse (nonpoint) sources are distributed over a range of river kilometers. The contribution to each reach is proportional to the overlap length between the diffuse source range and the reach:

Q_{ds,i} = Q_{ds} \cdot \frac{L_{overlap}}{L_{ds}}

where $L_{overlap} = \min(x_{ds,max}, x_{r,max}) - \max(x_{ds,min}, x_{r,min})$ and $L_{ds}$ is the total diffuse source length.

Step 2: Kinetic Reactions

After mixing, first-order kinetic reactions are applied over the travel time $\Delta t$ of the reach. See the BOD Decay and Nitrogen Cycle pages for detailed kinetic equations.

Step 3: Dissolved Oxygen Balance

The final DO concentration is computed by balancing sources (reaeration, photosynthesis) against sinks (BOD oxidation, nitrification, SOD). See the DO Balance page.

Travel Time and Spatial Discretization

The travel time through each reach determines how long kinetic reactions act on the water:

\Delta t = \frac{A_c \cdot L}{Q \cdot 86400} \quad \text{(days)}

Longer reaches have longer travel times and more kinetic decay. The choice of reach discretization affects the spatial resolution of the model output but does not change the total mass balance (which is conservative).

Python Implementation

pythonsimulation.py — Downstream march loop

for i in range(1, nr + 1):
    reach = reaches[i]
    q_up = reaches[i - 1].get('q', 0)

    # Step 1: Conservative mixing
    mass_in = [c * q_up for c in conc[i - 1]]
    q_total = q_up

    for ps in point_sources:
        if _in_reach(ps['xptt'], x_up, x_dn):
            q_ps = ps['Qptt']
            for k in range(num_vars):
                mass_in[k] += q_ps * ps['cpttMean'][k]
            q_total += q_ps

    c_mix = [m / q_total for m in mass_in]

    # Step 2: Kinetics
    c_kin = apply_kinetics(c_mix, rates, reach, temp, trav_time)

    # Step 3: DO balance
    c_kin[2] = calculate_do_balance(c_kin[2], ...)