The discretized sediment transport, bed change, and bed material sorting equations are solved in a coupled manner but decoupled from the flow equations during each time step. This is referred to semi-coupling. The coupled equations are derived by treating the erosion implicitly. The erosion rate is written as E_{tk}^{*n+1} = f_{1k}^{n+1} E_{tk}^{n+1}.  Inserting this equation into the bed change and sorting equations and then substituting the bed sorting equation into the bed change equation. The sediment transport, bed change, and bed sorting equations are computed iteratively. Each iteration is computed as follows:

First the total bed change is computed as

1) \Delta z_{bi}^{m+1} = \frac{ \sum_k \frac{A_{ki}}{B_{ki}} }{1 - \sum_k \frac{G_{ki}}{B_{ki}} }

where

n = time step counter

m = sediment outer-loop iteration counter

A_{ki} = D_{tki}^{n+1} - (m_{1ki}^n \delta_{1i}^n - m_{*ki}^n \Delta \delta_{1i} ) \frac{E_{tki}^{*n+1} }{\rho_{d1i}^{m} \delta_1^{n+1}}+S_{bki}^{n+1

B_{ki} = (1 - \phi_{bi} ) \rho_{sk} \left( \frac{1}{f_M \Delta t} + \frac{E_{tki}^{*n+1} }{\rho_{d1i}^{m} \delta_{1i}^{n+1}} \right)

G_{ki} = m_{*ki}^n \frac{E_{tki}^{*n+1} }{\rho_{d1i}^{m} \delta_{1i}^{n+1}}

\Delta \delta_{1i} = \delta_{1i}^{n+1} - \delta_{1i}^{n}

The fractional bed change is then computed as

2) \Delta z_{bki}^{m+1} = \frac{ A_{ki} + \Delta z_{bi}^{m+1} G_{ki} }{B_{ki} }

The fractional mass bed change is computed as 

3) \Delta M_{bki}^{m+1} = \Delta z_{bki}^{m+1} ( 1 - \phi_{bi} ) \rho_{sk}

The active layer fractional mass concentrations are computed as

4) m_{1ki}^{m+1} = \frac{\Delta M_{bki}^{m+1} + m_{1ki}^{n}\delta_{1i}^{n} - m_{*ki}^{n} \Delta \delta_{2i}}{\delta_{1i}^{n+1}}

The active layer bulk dry density is computed as

5) \rho_{d1i}^{m+1} = \sum_k m_{1ki}^{m+1}

The active layer grain class fractions are then computed as

6) f_{1ki}^{m+1} = \frac{ m_{1ki}^{m+1} } { \rho_{d1i}^{m+1} }

The convergence of the algorithm is monitored using the grain class fractions. At convergence the m+1 variables become n+1