Conceptual model

One of the underlying principles of dynamic TOPMODEL is that the landscape can be broken up into hydrologically similar regions, or Hydrological Response Units (HRUs), where all the area within a HRU behaves in a hydrologically similar fashion. Further discussion and the connection of HRUs is outlined elsewhere.

This document outlines the conceptual structure and computational implementation of the hillslope HRU. The hillslope HRU is a representation of an area of hill slope with, for example, similar topographic, soil and upslope area characteristics.

The hillslope HRU is considered to be a area of hill slope with outflow occurring across a specified width of it boundary. It is formed of four zones representing the surface water, which passes water between HRUs and drains to the root zone. The root zone characterises the interactions between evapotranspiration and precipitation and when full spills to the unsaturated zone. This drains to the saturated zone which also interacts with other HRUs. The behaviour of the saturated zone is modelled using a kinematic wave approximation. The zones and variables used below are shown in the schematic diagram.

In the following section the governing equations of the hillslope HRU are given in a finite volume form. Approximating equations for the solution of the governing equations are then derived with associated implicit and semi-implicit numerical schemes. These are valid for a wide range of transmissivity profiles, including those present in the code and presented in the third sections of the vignette.

Two appendices provide supporting information justifying the numerical schemes and an alternative derivation of the governing equations based on the infinitesimal vertical slab used in the original derivation of Dynamic TOPMODEL.

Notation

Table 1 outlines the notation used for describing an infinitesimal slab across the hillslope HRU.

The following conventions are used:

• All properties such as slope angles and time constants are considered uniform along the hillslope.

• Precipitation and Potential Evapotranspiration occur at spatially uniform rates.

• All fluxes \(r_{\star \rightarrow \star}\) are considered to be positive when travelling in the direction of the arrow, but can, in some cases be negative.

• All fluxes \(l_{\star}\) and \(q_{\star}\) are considered to be positive flows travelling downslope.

• The superscripts \(^+\) and \(^-\) are used to denote variables for the outflow and inflow to the hillslope.

• \(\left.y\right\rvert_t\) indicates the value of the variable \(y\) at time \(t\)

• Spatial averages are denoted with a \(\bar{}\) so for example integrating over the length of the hillslope we find

\[\begin{equation} \int w y dx = A\bar{y} \end{equation}\]

• Average values over a time step are denoted with a \(\hat{}\). If the temporal average is of a spatial average the \(\bar{}\) is dropped. For example the temporal average of the spatial average \(\bar{y}\) over a time step is \(\hat{y}\).

Table 1: Outline of notation for describing a cross section of the hillslope HRU

Quantity type	Symbol	Description	unit
Storage	\(s_{sf}\)	Surface excess storage per unit width	m
	\(s_{rz}\)	Root zone storage per unit width	m
	\(s_{uz}\)	Unsaturated zone storage per unit width	m
	\(s_{sz}\)	Saturated zone storage deficit per unit width	m
Vertical fluxes	\(r_{sf \rightarrow rz}\)	Flow from the surface excess store to the root zone per unit width	m/s
	\(r_{rz \rightarrow uz}\)	Flow from the root zone to the unsaturated zone per unit width	m/s
	\(r_{uz \rightarrow sz}\)	Flow from unsaturated to saturated zone per unit width	m/s
Lateral fluxes	\(l_{sf}\)	Lateral flow in the surface zone per unit outflow width	m\(^2\)/s
	\(l_{sz}\)	Lateral flow in the saturated zone per unit outflow width	m\(^2\)/s
	\(q_{sf}\)	Lateral flow in the hillslope surface zone	m\(^3\)/s
	\(q_{sz}\)	Lateral inflow to the hillslope surface zone	m\(^3\)/s
Vertical fluxes to the HRU	\(p\)	Precipitation rate	m/s
	\(e_p\)	Potential Evapotranspiration rate	m/s
Other HRU properties	\(A\)	Plan area	m\(^3\)
	\(w\)	Width of a hill slope cross section	m
	\(\Delta x\)	Effective length of the hillslope HRU (Area divided by outflow width)	m
	\(k_{sf}\)	Maximum flow rate from surface excess to root zone per unit width	m/s
	\(c_{sf}\)	Surface excess storage celerity	m/s
	\(s_{raf}\)	Surface runoff attenuation feature storage	m
	\(t_{raf}\)	Surface runoff attenuation feature time constant	s
	\(s_{rzmax}\)	Maximum root zone storage depth	m
	\(T_d\)	Unsaturated zone time constant per m of saturated storage deficit.	s/m
	\(\Theta_{sz}\)	Further properties and parameters for the solution of the saturated zone	:-:

Finite Volume Formulation

In this section a finite volume formulation of the Dynamic TOPMODEL equations is derived.

Approximating Equations

No analytical solution yet exists for simultaneous integration of the system of ODEs outlined above. In the following an implicit scheme, where fluxes between stores are considered constant over the time step and gradients are evaluated at the final state, is presented.

The basis of the solution is that a gravity driven system will maximise the downward flow of water within each timestep of size \(\Delta t\).

Recalling that the temporal and spatial average vertical fluxes are denoted with a \(\hat{}\) the implicit approximations are stated in the following sections using the unknown vertical fluxes \(\hat{r}_{\star \rightarrow \star}\) which are then determined by maximising the downward flux.

Numerical Algorithm

An implementation of the approximating equations which is consistent with maximising the downward flux produces a fully implicit numeric scheme.

The maximum downward flux from the surface is given by \[ \check{r}_{sf \rightarrow rz} = \min\left( k_{sf}, \left. \frac{1}{\Delta t}\left(\bar{s}_{sf} \right\rvert_{0} + \frac{\Delta t}{\Delta x} \hat{l}_{sf}^{-} \right)\right) \]

which results in a maximum downward flux from the root zone of \[ \check{r}_{rz \rightarrow uz} = \max\left(0, \frac{1}{\Delta t}\left( \left. \bar{s}_{rz} \right.\rvert_{0} + \Delta t \left( p + \check{r}_{sf \rightarrow rz} - e_{p} \right) - s_{rzmax} \right)\right) \]

Consider the pseudo-function \[ f\left(z\right) = z - \left. \bar{s}_{sz}\right.\rvert_{0} - \frac{\Delta t}{\Delta x} \left(g\left(z,\Theta_{sz}\right) - l_{sz}^{-}\right) + \Delta t \min\left(\frac{1}{T_d},\frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t}\right) \]

adapted from the approximating equation for the saturated zone storage with the minima ensuring the limit on the flux from the unsaturated zone. If \(f(0) \geq 0\) the downward flux is enough to produce saturation in the subsurface and \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\). If \(f(0) < 0\) then \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) is given by a solution of \(f(z) = 0\). In Appendix A it is shown that \(f\left(z\right)\) is a monotonic increasing function of \(z > 0\); hence if \(f\left(0\right) \leq 0\) and \(f\left(D\right) \geq 0\) a unique solution for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) can be found.

A direct implementation of the above results takes the form of Algorithm 1.

---

Algorithm 1: Direct Implementation

1. Compute \(\check{r}_{sf \rightarrow rz}\) and \(\check{r}_{rz \rightarrow uz}\)
2. Test for saturation. If \(f\left(0\right) \geq 0\) then go to Step 3 else go to Step 4
3. Saturated. Set \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\) and \[\hat{r}_{rz \rightarrow uz} = \frac{1}{\Delta t}\left( \left. \bar{s}_{sz}\right.\rvert_{\Delta 0} + \frac{\Delta t}{\Delta x} \left(g\left(0,\Theta_{sz}\right) - l_{sz}^{-}\right) - \left. \bar{s}_{uz} \right.\rvert_{0} \right) \] Go to Step 5.
4. Not saturated. Find \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) such that \(f\left(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\right) = 0\) and take \[ \hat{r}_{rz \rightarrow uz} = \min\left(\check{r}_{rz \rightarrow uz},\frac{1}{\Delta t} \left( \left. \bar{s}_{sz} \right.\rvert_{\Delta t} + \frac{\Delta t}{T_d} - \left. \bar{s}_{uz} \right.\rvert_{0}\right)\right)\] Go to Step 5.
5. Compute \[ \hat{r}_{sf \rightarrow rz} = \min\left(\check{r}_{sf \rightarrow rz}, \frac{1}{\Delta t}\left( s_{rzmax} - \left. \bar{s}_{rz} \right.\rvert_{0} - \Delta t \left( p - e_p - \hat{r}_{rz \rightarrow uz} \right) \right)\right)\]
6. Solve for \(\left. \bar{s}_{uz}\right.\rvert_{\Delta t}\), \(\left. \bar{s}_{rz}\right.\rvert_{\Delta t}\) and \(\left. \bar{s}_{sf}\right.\rvert_{\Delta t}\) using the computed vertical fluxes.

---

The challenge with Algorithm 1 is that maintaining a water (mass) balance and correct limits on the state values depends upon the accuracy of the solution of \(f\left(z\right)=0\). Since finding \(z\) requires multiple evaluations of the transmissivity profile it can rapidly become the most expensive part of the code.

To limit the computational cost suppose that \(f\left(z\right)=0\) is solved approximately to give \(\tilde{z}\) such that for some numerical tolerance \(\epsilon\) we have \(\left\lvert f\left(\tilde{z}\right) \right\rvert \leq \epsilon\). Using \(\tilde{z}\) to evaluate the outward lateral flux from the saturated zone gives the pseudo-function \[ h\left(z\right) = z - \left. \bar{s}_{sz}\right.\rvert_{0} - \frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) - l_{sz}^{-}\right) + \Delta t \min\left(\frac{1}{T_d},\frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t}\right) \]

Since \(h\left(z\right)\) is a monotonic increasing function of \(z > 0\) the solution for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) is similar to that using \(f\left(z\right)\). Taking \(z \in \left[0,D\right]\) the solution for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\) is given by

• If \(h(0) \geq 0\) the inward flux is enough to produce saturation in the subsurface and \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\).
• If \(h(z_{max}) < 0\) (possible only due to the outflow flux approximation) then \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = D\) and the lateral outflow is adjusted such that \[ \frac{\Delta t}{\Delta x} l_{sz}^{+} = D - \left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} l_{sz}^{-} + \Delta t \min\left(\frac{1}{T_d},\frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} D + \Delta t}\right) \]
• If \(h(0) < 0\) and If \(h(D) \geq 0\) then with \(z_{C} = \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t\left( \check{r}_{rz \rightarrow uz} - \frac{1}{T_{d}}\right)\)
  • If \(h\left( z_{C} \right) \geq 0\) then the inflow from the unsaturated zone is \(1/T_{d}\) and \[ \left. \bar{s}_{sz}\right.\rvert_{\Delta t} = \left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) - l_{sz}^{-}\right) - \frac{\Delta t}{T_d} \]
  • else (see Appendix C) \[ z = \frac{1}{2}\left(-b + \sqrt{ b^2 - 4c }\right) \\ = \frac{1}{2} \left( z_{C} - h\left(z_{C}\right) + \sqrt{ \left(h\left(z_{C}\right)-z_{C}\right)^2 - 4 \frac{\Delta t}{T_{d}} h\left(z_{C}\right) } \right) \]

The implementation of this approach is presented in Algorithm 2. it can be seem that the exchange for decreasing the accuracy of the solution of the zero finding problem is a potential to evaluate a power and a square root.

---

Algorithm 2: Stable Implementation

1. Compute \(\check{r}_{sf \rightarrow rz}\) and \(\check{r}_{rz \rightarrow uz}\)
2. Compute the lateral outflow
   • If \(f\left(0\right) \geq 0\) then \(l_{sz}^{+} = l_{szmax}\)
   • else find \(\tilde{z}\) such that \(\left\lvert f\left(\tilde{z}\right) \right\rvert \leq \epsilon\). Set \(l_{sz}^{+} = g\left(\tilde{z},\Theta\right)\)
3. Solve for \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t}\)
   • If \(h(0) \geq 0\) then \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\).
   • If \(h(z_{max}) < 0\) then \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = D\) and \[ \frac{\Delta t}{\Delta x} l_{sz}^{+} = D - \left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} l_{sz}^{-} + \Delta t \min\left(\frac{1}{T_d},\frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} D + \Delta t}\right) \]
   • else with \(z_{C} = \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t\left( \check{r}_{rz \rightarrow uz} - \frac{1}{T_{d}}\right)\)
     • If \(h\left( z_{C} \right) \geq 0\) then \[ \left. \bar{s}_{sz}\right.\rvert_{\Delta t} = \left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) - l_{sz}^{-}\right) - \frac{\Delta t}{T_d} \]
     • else \[ \left. \bar{s}_{sz}\right.\rvert_{\Delta t} = = \frac{1}{2} \left( z_{C} - h\left(z_{C}\right) + \sqrt{ \left(h\left(z_{C}\right)-z_{C}\right)^2 - 4 \frac{\Delta t}{T_{d}} h\left(z_{C}\right) } \right) \]
4. Solve for \(\hat{r}_{rz \rightarrow uz}\)
   • If \(\left. \bar{s}_{sz}\right.\rvert_{\Delta t} = 0\) then \[\hat{r}_{rz \rightarrow uz} = \frac{1}{\Delta t}\left( \left. \bar{s}_{sz}\right.\rvert_{\Delta 0} + \frac{\Delta t}{\Delta x} \left(g\left(0,\Theta_{sz}\right) - l_{sz}^{-}\right) - \left. \bar{s}_{uz} \right.\rvert_{0} \right) \]
   • else \[ \hat{r}_{rz \rightarrow uz} = \min\left(\check{r}_{rz \rightarrow uz},\frac{1}{\Delta t} \left( \left. \bar{s}_{sz} \right.\rvert_{\Delta t} + \frac{\Delta t}{T_d} - \left. \bar{s}_{uz} \right.\rvert_{0}\right)\right)\]
5. Compute \[ \hat{r}_{sf \rightarrow rz} = \min\left(\check{r}_{sf \rightarrow rz}, \frac{1}{\Delta t}\left( s_{rzmax} - \left. \bar{s}_{rz} \right.\rvert_{0} - \Delta t \left( p - e_p - \hat{r}_{rz \rightarrow uz} \right) \right)\right)\]
6. Solve for \(\left. \bar{s}_{uz}\right.\rvert_{\Delta t}\), \(\left. \bar{s}_{rz}\right.\rvert_{\Delta t}\) and \(\left. \bar{s}_{sf}\right.\rvert_{\Delta t}\) using the computed vertical fluxes.

---

Transmissivity Profiles

The solution and numerical scheme have been presented with a general transmissivity function \(g\left(\bar{s}_{sz},\Theta_{sz}\right)\). Table 2 present the transmissivity profiles present as options within the dynatop package, the corresponding value to use for the transmissivity_profile value in the model options vector and the additional parameters or properties required. The additional parameter and properties are defined in Table 3.

Table 2: Transmissivity profiles available with the dynatop package.

Name	\(g\left(\bar{s}_{sz},\Theta_{sz}\right)\)	Celerity	\(\Theta_{sz}\)	transmissivity_profile value	Notes
Exponential	\(T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m}s_{sz}\right)\)	\(\frac{\cos\beta}{m}g\left(\bar{s}_{sz},\Theta_{sz}\right)\)	\(T_0,\beta,m\)	"exp"	Originally given in Beven & Freer 2001
Bounded Exponential	\(T_{0}\sin\left(\beta\right)\left(\exp\left(-\frac{\cos\beta}{m}s_{sz}\right) - \exp\left(-\frac{\cos\beta}{m}D\right)\right)\)	\(\frac{cos \beta}{m}T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m}s_{sz}\right)\)	\(T_0,\beta,m,D\)	"bexp"	Originally given in Beven & Freer 2001
Constant Celerity	\(c_{sz}\left(D-s_{sz}\right)\)	\(c_{sz}\)	\(c_{sz},D\)	"cnst"
Double Exponential	\(T_{0}\sin\left(\beta\right)\left( \omega\exp\left(-\frac{\cos\beta}{m}s_{sz}\right) + \left(1-\omega\right)\exp\left(\frac{\cos\beta}{m_2}s_{sz}\right)\right)\)	\(\omega\frac{cos \beta}{m}T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m}s_{sz}\right) + \left(1-\omega \right)\frac{cos \beta}{m_2}T_{0}\sin\left(\beta\right)\exp\left(-\frac{\cos\beta}{m_2}s_{sz}\right)\)	\(T_0,\beta,m,m_2,\omega\)	“dexp”

Table 3: Additional parameters used in the transmissivity profiles.

Symbol	Description	unit
\(T_0\)	Transmissivity at saturation	m\(^2\)/s
\(m\),\(m_2\)	Exponential transmissivity constants	m\(^{-1}\)
\(\beta\)	Angle of hill slope	rad
\(c_{sz}\)	Saturated zone celerity	m/s
\(D\)	Depth at which zero lateral flow occurs	m
\(\omega\)	weighting parameter	-

Appendix A - Monotonicity of \(f\left(z\right)\)

By definition \(\left. \bar{s}_{uz} \right.\rvert_{0}\), \(\Delta t\), \(\Delta X\), \(T_d\), \(\check{r}_{rz \rightarrow uz}\) and \(z\) are all greater or equal to 0. Rewrite the expression for \(f\left(z\right)\) as \[ f\left(z\right) = \left\{ \begin{array}{cl} z - \left. \bar{s}_{sz}\right.\rvert_{\Delta 0} - \frac{\Delta t}{\Delta x} \left(g\left(z,\Theta_{sz}\right) - l_{sz}^{-}\right) + \frac{\Delta t}{T_d} & \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} \geq \frac{1}{T_d} \\ z - \left. \bar{s}_{sz}\right.\rvert_{\Delta 0} - \frac{\Delta t}{\Delta x} \left(g\left(z,\Theta_{sz}\right) - l_{sz}^{-}\right) + \Delta t \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} & \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} < \frac{1}{T_d} \end{array}\right. \]

Differentiating with respect to \(z\) gives \[ \frac{d}{dz} f\left(z\right) = \left\{ \begin{array}{cl} 1 - \frac{\Delta t}{\Delta x} \frac{d}{dz}g\left(z,\Theta_{sz}\right) & \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} \geq \frac{1}{T_d} \\ 1 - \frac{\Delta t}{\Delta x} \frac{d}{dz}g\left(z,\Theta_{sz}\right) - T_d \Delta t \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { \left(T_{d} z + \Delta t\right)^2} & \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} < \frac{1}{T_d} \end{array}\right. \]

The definition of the function \(g\) states that \[ -\frac{d}{dz}g\left(z,\Theta_{sz}\right) \geq 0 \] Using this relationship and the expression for the range of \(z\) gives \[ \frac{d}{dz} f\left(z\right) \geq \left\{ \begin{array}{cl} 1 & \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} \geq \frac{1}{T_d} \\ 1 - \frac{\Delta t}{T_{d} z + \Delta t} & \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} < \frac{1}{T_d} \end{array}\right. \]

From this it can be seen that \(\frac{d}{dz} f\left(z\right) \geq 0\) with equality possible when \(z=0\).

Appendix B - Comparison to the Infinitesimal slab derivation

Appendix C Solution to \(h\left(z\right)\) when \(h(0) < 0\) and \(h(D) \geq 0\)

In this case we see that \(0 < z < D\). Taking \(z_{C} = \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t\left( \check{r}_{rz \rightarrow uz} - \frac{1}{T_{d}}\right)\) the function \(h\) can be rewritten as \[ h\left(z\right) = \left\{ \begin{array}{cl} z - \left. \bar{s}_{sz}\right.\rvert_{\Delta 0} - \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) + \frac{\Delta t}{T_d} & z_{C} \geq z \\ z - \left. \bar{s}_{sz}\right.\rvert_{0} - \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) + \Delta t \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} & z_{C} < z \end{array} \right. \]

Using the monotonicity of \(h\left(z\right)\) we solve each part separately. Firstly if \(h\left( z_{C} \right) \geq 0\) then the inflow from the unsaturated zone is \(1/T_{d}\) and \[ \left. \bar{s}_{sz}\right.\rvert_{\Delta t} = \left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} \left(g\left(\tilde{z},\Theta_{sz}\right) - l_{sz}^{-}\right) - \frac{\Delta t}{T_d} \]

if \(h\left( z_{C} \right) < 0\) the value of \(z\) satisfies \[ z - \left. \bar{s}_{sz}\right.\rvert_{0} - \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) + \Delta t \frac{ \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} } { T_{d} z + \Delta t} = 0 \] Rearrangement to a quadratic gives \[ 0 = T_{d} z^2 - T_{d} z \left(\left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) \right) + \Delta t z - \Delta t \left(\left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) \right) + \Delta t \left( \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} \right) \] \[ 0 = z^2 - z \left(\left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) - \frac{\Delta t}{T_{d}}\right) + \frac{\Delta t}{T_{d}} \left( \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} - \left. \bar{s}_{sz}\right.\rvert_{0} - \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) \right) \]

The commonly quoted solution for the quadratic \[ 0 = z^2 + b z +c \] is \[ z = \frac{-b \pm\sqrt{ b^2 - 4c }}{2} \] Comparing terms we see that \[ -b = \left. \bar{s}_{sz}\right.\rvert_{0} + \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) - \frac{\Delta t}{T_{d}} \]

and \[ c = \frac{\Delta t}{T_{d}} \left( \left. \bar{s}_{uz} \right.\rvert_{0} + \Delta t \check{r}_{rz \rightarrow uz} - \left. \bar{s}_{sz}\right.\rvert_{0} - \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) \right) \]

Since \(-b = z_{C} - h\left(z_{C}\right)\) gives \(-b > 0\). Since \(z > z_{C}\) the subsurface (\(s_{uz}\) and \(s_{sz}\)) does not fill so \[ \left. \bar{s}_{sz}\right.\rvert_{0} - \left. \bar{s}_{uz} \right.\rvert_{0} \gt \Delta t \check{r}_{rz \rightarrow uz} - \frac{\Delta t}{\Delta x} \left(l_{sz}^{+} - l_{sz}^{-}\right) \] and \(c < 0\)

Combining these conditions results in \(\sqrt{b^2 - 4c} > \left\lvert b \right\rvert\). Since \(z>0\) only one root of the quadratic is valid (as always positive) so \[ z = \frac{1}{2}\left(-b + \sqrt{ b^2 - 4c }\right) \]