Motivation #

In MSMA 2, we have level-pool routing and the Modified Puls equation, which engineers use to compute the pond water level and discharge timestep by timestep. It’s always fun for me to find out exactly how such equations follow from more general physical principles. So in this article we will do just that. We start from the continuity equation, which expresses conservation of mass and, with constant water density, becomes a balance of water volume. We combine it with the momentum and energy relationships developed in From Newton’s Second Law to Bernoulli’s Equation, which supply the outflow. We then introduce the assumptions that turn this balance into level-pool routing, derive the Modified Puls equation and its storage indicator, work through one timestep, and calculate how long the stored water takes to drain after inflow stops.

Level-pool routing gives a differential equation for the pool level. Modified Puls solves it numerically, one timestep after another, which suits a spreadsheet or computer. But what about hand calculation? Can we limit ourselves to some special, simple pond geometry and counter-check the numerical results? As a bonus, we will also derive an analytical, closed-form solution for the drawdown of simple basins after inflow stops, and use it to check the Modified Puls results.

Take note that level-pool routing is itself a simplification of more general equations. So there are situations where it may not apply, and we will remark briefly on them where they arise. In those cases, we have to use a more general simulation or equation that considers more factors.

Water entering and leaving a short slice #

Let’s first consider a short slice of an open water body, such as a drain, channel or detention basin. Take two fixed sections perpendicular to the direction of flow. Let x measure distance downstream, t denote time, and \Delta x be the distance between the sections. Let A(x,t) be its wetted cross-sectional area in square metres and Q(x,t) its discharge in cubic metres per second. We initially exclude lateral inflow, infiltration and evaporation, and take water density as constant.

For a sufficiently thin slice, volume is approximately A\Delta x . The rate at which that volume increases is inflow minus outflow:

\displaystyle \frac{\partial A}{\partial t}\Delta x=Q(x,t)-Q(x+\Delta x,t).

To first order,

\displaystyle Q(x+\Delta x,t)=Q(x,t)+\frac{\partial Q}{\partial x}\Delta x.

Substitute this expression, cancel Q(x,t) , divide by \Delta x , and take the small-segment limit:

\displaystyle \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=0.

This is the one-dimensional continuity equation for our stated conditions. If discharge decreases downstream, \partial Q/\partial x<0 , and \partial A/\partial t>0 : water accumulates and the wetted area increases. The same slice is passing water and changing its storage.

For steady flow, \partial A/\partial t=0 , so inflow equals outflow. There is still water held in the slice; its volume simply does not change with time. The absence of a storage change does not mean the physical volume is zero.

Integrate continuity over the whole facility #

We now add the water held in all the thin slices between two fixed sections. Each slice contributes its wetted cross-sectional area multiplied by its length. As the slices become thinner, their sum becomes an integral. Let S(t) denote the water volume, in cubic metres, between fixed sections x_1 and x_2 :

\displaystyle S(t)=\int_{x_1}^{x_2}A(x,t)\,dx.

Integrating continuity over the same length gives

\displaystyle \int_{x_1}^{x_2}\frac{\partial A}{\partial t}\,dx+\int_{x_1}^{x_2}\frac{\partial Q}{\partial x}\,dx=0.

The first integral is the time derivative of the stored volume. The second is discharge at the downstream boundary minus discharge at the upstream boundary:

\displaystyle \frac{dS}{dt}+Q(x_2,t)-Q(x_1,t)=0.

Water mass is density multiplied by volume. With constant density, that common factor cancels from the mass balance, leaving a balance of water volumes and volumetric flows. Define I(t)=Q(x_1,t) as inflow and O(t)=Q(x_2,t) as outflow, both positive in their named directions. Then

\displaystyle \frac{dS}{dt}=I-O.

For a pond with several inlets and outlets, we obtain the same result by summing all inflows into I and all outflows into O .

Where momentum and Bernoulli enter #

The balance dS/dt=I-O tells us how storage changes once we know the outflow, but continuity alone, whether for a slice or the whole facility, is insufficient to uniquely determine the motion of the water. We also need the relationship between motion, pressure, elevation and resistance. This is where the Newton–Bernoulli derivation comes in.

Take the straight, prismatic rectangular open channel shown above with hydrostatic pressure, small bed slope, negligible vertical acceleration and a section-average velocity approximation. The right-hand inset is a pressure plot: gauge pressure increases to the right, and depth below the free surface increases downwards. Its top point is the free surface; its bottom edge is the bed. Let V=Q/A be mean velocity, y water depth, z_b bed elevation, and \eta=z_b+y water-surface elevation above a common datum. For flow in the positive x direction, let S_f be the friction slope, the positive head loss per unit length.

We obtain the momentum equation by applying Newton’s second law to a short slice of this channel. Let \rho be the constant water density, g gravitational acceleration, and b the constant channel width, so A=by . We use one mean velocity for each section, taking the momentum correction factor as one.

First consider the pressure force on either end of the slice. Let \xi (xi) measure height above the local bed, from zero to the water depth y . The yellow strip in the left-hand drawing lies at height \xi , so its depth below the surface is y-\xi . Its thickness is d\xi and its area is b\,d\xi . The right-hand plot gives its hydrostatic gauge pressure p=\rho g(y-\xi) . Pressure multiplied by strip area gives the strip force p\,b\,d\xi . Adding these forces over the wetted end face gives the resultant pressure force F_p :

\displaystyle F_p=\int_0^y \rho g(y-\xi)b\,d\xi=\frac{1}{2}\rho gb y^2.

The lower drawing shows how these resultant forces act on the slice: F_p(x) on its upstream end pushes downstream, while F_p(x+\Delta x) on its downstream end pushes upstream. Their magnitudes can differ because the water depths at the two ends differ. Across a short distance \Delta x , their net force is approximately

\displaystyle -\frac{\partial F_p}{\partial x}\Delta x=-\rho gby\frac{\partial y}{\partial x}\Delta x=-\rho gA\frac{\partial y}{\partial x}\Delta x.

For our small bed slope, the downstream component of the slice’s weight is approximately -\rho gA(\partial z_b/\partial x)\Delta x . A falling bed has \partial z_b/\partial x<0 , giving a positive downstream force. The resisting force is -\rho gAS_f\Delta x : this expresses resistance through the friction slope already defined above. Adding the forces and using \eta=z_b+y gives

\displaystyle \text{net downstream force}=-\rho gA\left(\frac{\partial\eta}{\partial x}+S_f\right)\Delta x.

Newton’s law for this fixed slice equates net force to the rate of increase of momentum inside it plus net momentum carried out through its ends. Momentum per unit length is \rho AV ; momentum carried through an end per second is \rho AV^2 . Dividing the short-slice balance by \rho\Delta x and taking its thin-slice limit gives

\displaystyle \frac{\partial(AV)}{\partial t}+\frac{\partial(AV^2)}{\partial x}=-gA\left(\frac{\partial\eta}{\partial x}+S_f\right).

Expanding the left side by the product rule groups it as

\displaystyle A\left(\frac{\partial V}{\partial t}+V\frac{\partial V}{\partial x}\right)+V\left(\frac{\partial A}{\partial t}+\frac{\partial(AV)}{\partial x}\right).

The second bracket is zero by continuity, because AV=Q . Divide the remaining balance by the wetted area A>0 and move the surface-gradient term to the left:

\displaystyle \frac{\partial V}{\partial t}+V\frac{\partial V}{\partial x}+g\frac{\partial\eta}{\partial x}=-gS_f.

The first two terms are the acceleration of the water as its velocity changes with time and along the channel. The surface gradient combines the pressure and bed-elevation effects; S_f supplies resistance. The friction law used to determine S_f is a separate hydraulic input.

For steady flow, the local acceleration term disappears. Dividing by g gives

\displaystyle \frac{V}{g}\frac{dV}{dx}+\frac{d\eta}{dx}=-S_f.

Because V\,dV/dx=d(V^2/2)/dx , we can write

\displaystyle \frac{d}{dx}\left(\eta+\frac{V^2}{2g}\right)=-S_f.

Integrating between sections 1 and 2 gives

\displaystyle \eta_1+\frac{V_1^2}{2g}=\eta_2+\frac{V_2^2}{2g}+\int_{x_1}^{x_2}S_f\,dx.

For hydrostatic open-channel flow, \eta=z_b+y combines bed elevation with the pressure head represented by water depth. The equation above is the extended Bernoulli, or energy, equation for this channel, including distributed friction loss. Its lossless form has the same total head at both sections. Entrance, exit and other local losses can be included as additional head-loss terms.

The extra assumptions that give level-pool routing #

One level describes the stored volume #

We now approximate the facility’s free surface as horizontal at each instant. The momentum equation tells us when this is reasonable: in a pond, velocities are small, so the acceleration terms and the friction slope are small, leaving g\,\partial\eta/\partial x\approx0 . Let h(t) be its elevation above a selected datum. Basin geometry then supplies the storage curve S(h) .

Let A_s(h) be the horizontal water-surface area at elevation h . The highlighted layer of thickness \Delta h adds volume \Delta S\approx A_s(h)\Delta h . Starting from known storage S(h_0) at reference level h_0 , we add the layers up to h . Let \zeta denote the layer elevation inside the integral. In the thin-layer limit, dS=A_s(h)\,dh , giving

\displaystyle S(h)=S(h_0)+\int_{h_0}^{h}A_s(\zeta)\,d\zeta,\qquad \frac{dS}{dh}=A_s(h).

Outlet hydraulics provides the outflow #

In the sketch, the horizontal plan area controls storage; the outlet opening area controls discharge. Pool level and opening centre use one elevation datum. For fixed downstream conditions and outlet settings, we assign one outflow O(h) to each level, assuming the outlet adjusts quickly enough to use its steady rating at the current pool level.

For a small, fully wetted orifice discharging into air, let z_c be its centre elevation above the same datum. The driving head is H=h-z_c , not automatically the pond depth. With negligible approach velocity, the familiar rating is

\displaystyle O(h)=C_dA_o\sqrt{2g(h-z_c)}.

Here A_o is opening area and C_d is the discharge coefficient. The orifice article explains the opening-size approximation, partial wetting and submerged discharge. A weir or culvert requires the appropriate weir or culvert relationship. Several outlets can contribute to the total rating, with shared approach and downstream conditions treated consistently.

Bring the relationships together #

Applying the chain rule to storage gives

\displaystyle \frac{dS}{dt}=\frac{dS}{dh}\frac{dh}{dt}=A_s(h)\frac{dh}{dt}.

Thus our level-pool routing equation is

\displaystyle A_s(h)\frac{dh}{dt}=I(t)-O(h).

The water-surface area controls how quickly a flow imbalance changes the level. For the same net inflow, a larger surface area rises more slowly. The rising level then changes the outlet head and therefore the outflow.

When one pool level loses important behaviour #

Equal storage may give different outlet flows #

The level-pool equation assumes that one pool level fixes both the storage and the outflow. Consider the two water-surface profiles above in a rectangular detention drain. Both have an average depth of 1.00 m, so they hold the same volume for the same width, length and horizontal bed. A level-pool calculation would represent both by a uniform depth of 1.00 m. But the outlet is at the right-hand end, where the lower profile is 1.20 m deep.

For that lower profile, we use the local depth of 1.20 m to obtain the head above the orifice centre. Using the 1.00 m average would underestimate the head and free-orifice discharge. The stored volume alone cannot tell us which outlet depth applies, because it contains no information about the surface slope.

Tailwater changes the outlet head and discharge #

The water level on the receiving side of the outlet is its tailwater. Let h_d be that downstream surface elevation above the same datum as the pool level h . When a small orifice is submerged on both sides and approach velocities are negligible, the driving head is h-h_d . For positive upstream-to-downstream flow with h>h_d ,

\displaystyle O=C_dA_o\sqrt{2g(h-h_d)}.

At the same upstream level, raising h_d reduces the driving head and discharge. If downstream level changes independently, write O(h,h_d(t)) . The pool mass balance remains valid, but a fixed one-variable rating O(h) , and the storage-indicator table we build from it below, no longer contains all the required information. When we discretise the balance in time, evaluate the start-of-step outflow with the tailwater at the start of the step and the end-of-step outflow with the tailwater at its end. An interacting downstream system may require a coupled hydraulic solution.

For example, take C_dA_o\sqrt{2g}=0.20\ \mathrm{m^{5/2}/s} , the value used in the worked timestep below, and upstream elevation h=1.20\ \mathrm{m} , downstream levels of 0.60 m and 0.90 m give O=0.20\sqrt{1.20-0.60}=0.1549\ \mathrm{m^3/s} and O=0.20\sqrt{1.20-0.90}=0.1095\ \mathrm{m^3/s} , respectively. Both sides cover the opening in this example. The same upstream storage now drains at different rates. The tailwater and backwater derivation follows how a receiving-water level affects upstream levels along a drain.

Derive the Modified Puls equation #

Integrate one timestep #

We divide the storm duration into short time intervals. During each interval, we calculate the volume entering and leaving, update the stored volume and pool level, then carry that state into the next interval. This is time discretisation. Suppose we know the initial level and storage, the outlet rating and the inflow hydrograph. Let subscripts 1 and 2 denote the beginning and end of one timestep, and let \Delta t=t_2-t_1 be its duration in seconds. In particular, S_1=S(t_1) is the stored volume at the start, and S_2=S(t_2) is the stored volume at the end, both in cubic metres. The diagram above shows these two states of the same basin.

Integrating the whole-facility balance gives the exact volume relationship

\displaystyle S_2-S_1=\int_{t_1}^{t_2}I(t)\,dt-\int_{t_1}^{t_2}O(t)\,dt.

We approximate each integral by the average of its endpoint flows multiplied by the interval. This is the trapezoidal rule:

\displaystyle S_2-S_1\approx\frac{I_1+I_2}{2}\Delta t-\frac{O_1+O_2}{2}\Delta t.

The inflow integral is exact if inflow varies linearly over the interval. The outlet response generally is nonlinear, so timestep selection still matters.

Using this approximation as our numerical routing equation, multiplying by 2/\Delta t gives

\displaystyle \frac{2S_2}{\Delta t}-\frac{2S_1}{\Delta t}=I_1+I_2-O_1-O_2.

Move the new outflow to the left and the old storage to the right:

\displaystyle \frac{2S_2}{\Delta t}+O_2=I_1+I_2+\frac{2S_1}{\Delta t}-O_1.

This is the Modified Puls, or storage-indication, routing equation. The unknown end-of-step storage and outflow are on the left. The right contains the known inflow ordinates and the previous state. HEC-HMS uses this formulation in its reservoir-routing reference.

MSMA 2 (2012), §2.5, Equations 2.16–2.17, uses the same trapezoidal balance and its rearrangement. We use seconds and flows in cubic metres per second throughout: a flow multiplied by a timestep must give a volume in cubic metres.

For the calculation in practice, see our stage-storage, stage-discharge and storage-indicator example and the Pond Routing guide, which follows the tabular procedure. The MSMA 2 Appendix 7A benchmark compares a worked detention-pond design with the software calculation.

Find the next pool level from one lookup #

The new storage S_2 and outflow O_2 are both unknown, but they belong to the same new pool level h_2 . For each trial level, the storage curve supplies S(h) and the outlet rating supplies O(h) . We can therefore tabulate the entire left side of the routing equation in advance. Call this tabulated quantity F(h) :

\displaystyle F(h)=\frac{2S(h)}{\Delta t}+O(h).

This is the storage indicator: a lookup value used to find the next water level. It adds 2S/\Delta t to outflow; both terms have units of discharge, \mathrm{m^3/s} . It does not represent an additional physical flow through the pond. From the known right-hand side, we solve

\displaystyle F(h_2)=I_1+I_2+\frac{2S(h_1)}{\Delta t}-O(h_1).

We then obtain both S_2=S(h_2) and O_2=O(h_2) . For an ordinary passive outlet with nondecreasing outflow and increasing storage, F(h) increases with level, permitting a unique lookup or numerical solution over its valid range.

In the graph above, R denotes the known right-hand side of the routing equation. We move horizontally from R to the indicator curve, then vertically down to find h_2 . The same level gives both the new storage and the new discharge. The curve is schematic; its shape and scale depend on the storage curve, outlet rating and timestep.

Some tables use S/\Delta t+O/2 instead. It is half of F , obtained by halving both sides of the same routing equation. The two conventions are equivalent, provided we use one consistently. Changing the timestep changes the indicator curve, because its storage term contains \Delta t .

A worked timestep #

The shaded areas are volumes: inflow adds 24.00 m³ and outflow releases about 12.04 m³, leaving about 11.96 m³ in storage. Each line joins the endpoint flows used in this 60 s trapezoidal calculation; the outlet line approximates the continuous response.

Consider an illustrative rectangular basin with vertical sides and horizontal surface area A_s=1000\ \mathrm{m^2} . Take its floor as the level datum, so S(h)=1000h . The small orifice centre is 0.20 m above the floor. Over the fully wetted range considered here, use the illustrative rating

\displaystyle O(h)=0.20\sqrt{h-0.20}\quad\mathrm{m^3/s}.

Levels are in metres, and the coefficient 0.20 has units \mathrm{m^{5/2}/s} . It represents the chosen C_dA_o\sqrt{2g} for this example; it is not a universal orifice coefficient. The outlet discharges freely into air and its approach velocity is negligible.

At the beginning, h_1=1.20\ \mathrm{m} , giving S_1=1200\ \mathrm{m^3} and O_1=0.200\ \mathrm{m^3/s} . Inflow rises linearly from I_1=0.200 to I_2=0.600\ \mathrm{m^3/s} over \Delta t=60\ \mathrm{s} .

The known right-hand side is

\displaystyle I_1+I_2+\frac{2S_1}{\Delta t}-O_1=0.200+0.600+\frac{2(1200)}{60}-0.200=40.600\ \mathrm{m^3/s}.

Our indicator is

\displaystyle F(h)=\frac{2000h}{60}+0.20\sqrt{h-0.20}.

Tabulating nearby levels gives:

The target 40.600 lies between the last two rows. Interpolation gives approximately h_2=1.21196\ \mathrm{m} . Solving the indicator equation directly gives S_2=1211.9642\ \mathrm{m^3} and O_2=0.2011929\ \mathrm{m^3/s} .

We can check the volume balance independently:

\displaystyle \text{inflow volume}=\frac{0.200+0.600}{2}(60)=24.0000\ \mathrm{m^3},

\displaystyle \text{outflow volume}\approx\frac{0.200+0.2011929}{2}(60)=12.0358\ \mathrm{m^3}.

The retained difference is 24.0000-12.0358=11.9642\ \mathrm{m^3} , agreeing with S_2-S_1 . Dividing by the 1000 m² surface area gives the level rise of approximately 0.011964 m.

How long does the stored water take to drain? #

We need the time to fall from the starting level to a stated recovery level. The end of rainfall does not necessarily mark the end of runoff entering the pond; first consider the period after inflow has ceased, with no other significant gains or losses.

Let h_i be the initial pool elevation, h_f the lower target elevation, and T the elapsed drawdown time in seconds. Take t=0 when the pool is at h_i . For a fixed outlet rating that remains positive throughout this interval, the routing equation reduces to

\displaystyle A_s(h)\frac{dh}{dt}=-O(h).

Separate the time and level terms, then integrate from the initial state to the target:

\displaystyle dt=-\frac{A_s(h)}{O(h)}\,dh,\qquad T=-\int_{h_i}^{h_f}\frac{A_s(h)}{O(h)}\,dh=\int_{h_f}^{h_i}\frac{A_s(h)}{O(h)}\,dh.

The reversed limits make the elapsed time positive. The integrand has units \mathrm{m^2}/(\mathrm{m^3/s})=\mathrm{s/m} ; integrating over elevation gives seconds. A graded basin uses its actual water-surface area at each level. If several outlets operate, their combined rating belongs in the denominator, split into the applicable flow regimes where necessary.

An orifice and a basin with vertical walls #

Return to a small, fully wetted orifice discharging freely into air, with fixed opening area and coefficient and negligible approach velocity. Its centre remains at elevation z_c . Define H_i=h_i-z_c and H_f=h_f-z_c as the initial and final heads above that centre. These are heads on the outlet, not depths above the basin floor.

For brevity, let k=C_dA_o\sqrt{2g} , with units \mathrm{m^{5/2}/s} . The orifice rating is O=k\sqrt{h-z_c} . Substituting it into the general drawdown integral gives

\displaystyle T=\frac{1}{k}\int_{h_f}^{h_i}\frac{A_s(h)}{\sqrt{h-z_c}}\,dh.

This is the drawdown relationship in MSMA 2 §2.4.4, Equation 2.14, with the integration direction made explicit. For vertical walls, the horizontal area A_s is constant. Set H=h-z_c , so dH=dh , and evaluate the integral:

\displaystyle T=\frac{A_s}{k}\int_{H_f}^{H_i}H^{-1/2}\,dH=\frac{A_s}{k}\left[2\sqrt{H}\right]_{H_f}^{H_i}.

Therefore,

\displaystyle T=\frac{2A_s}{C_dA_o\sqrt{2g}}\left(\sqrt{H_i}-\sqrt{H_f}\right).

This is MSMA 2 Equation 2.15. The coefficient and the wetted-opening approximation remain hydraulic inputs; conservation of mass supplies the drawdown calculation once that rating is accepted.

We can also stop the integration at any intermediate time t , when the head is H(t)=h(t)-z_c :

\displaystyle t=\frac{2A_s}{k}\left(\sqrt{H_i}-\sqrt{H(t)}\right),\qquad H(t)=\left(\sqrt{H_i}-\frac{kt}{2A_s}\right)^2.

Use this expression only while the quantity inside the parentheses is nonnegative and the fully wetted orifice model remains applicable. The square-root head decreases linearly with time. Head itself falls progressively more slowly as the discharge reduces.

Drawdown of the worked basin #

Use the earlier basin: A_s=1000\ \mathrm{m^2} , z_c=0.20\ \mathrm{m} above the floor, and k=0.20\ \mathrm{m^{5/2}/s} . Take h_i=1.20\ \mathrm{m} and h_f=0.45\ \mathrm{m} , giving H_i=1.00\ \mathrm{m} and H_f=0.25\ \mathrm{m} . For this recession, retain the rating but realise the outlet as sixteen separate, equal circular openings at the same centre elevation. With an illustrative C_d=0.62 for each, diameters of approximately 76 mm supply the combined rating, neglecting interaction between jets. Here A_o is their total opening area. Their crowns are about 0.238 m above the floor, so the openings remain fully wetted and small relative to the final centre head. This is an illustrative hydraulic model; the adopted coefficient and practical outlet arrangement require their own design checks.

\displaystyle T=\frac{2(1000)}{0.20}\left(\sqrt{1.00}-\sqrt{0.25}\right)=5000\ \mathrm{s}=83.33\ \mathrm{min}.

The released volume is 1000(1.20-0.45)=750\ \mathrm{m^3} . Discharge falls from 0.200 to 0.100\ \mathrm{m^3/s} . Dividing the volume by the initial discharge would give only 750/0.200=3750\ \mathrm{s} : 25% shorter than the calculated time, because that shortcut holds the initial, higher outlet rate throughout the recession.

The 83.33 minutes is the time to lower the water level from 1.20 m to 0.45 m. At that final level, S(h_f)=1000(0.45)=450\ \mathrm{m^3} of water remains in the basin.

Check Modified Puls against the analytical solution #

We can now check Modified Puls by hand: calculate the same drawdown analytically and with Modified Puls, then compare the pool levels at the same time. Both use the same basin, orifice and zero inflow, so any difference comes from the Modified Puls timestep. The analytical solution gives the head at any chosen time from one equation; Modified Puls must march there one timestep after another.

Such a check is possible only when:

• there is no inflow, so A_s(h)\,dh/dt=-O(h) separates and integrates;
• outflow depends on pool level alone: free discharge, steady tailwater and fixed outlet settings, with any several outlets combined into one O(h) ;
• A_s(h)/O(h) integrates in closed form, as for vertical walls or a linearly widening basin with one orifice.

These conditions hold for drawdown after inflow stops, when the highest level is simply the starting level. The peak level during a storm has no general closed form, which is why we use Modified Puls to find it.

The basin with vertical walls #

First use the earlier basin with vertical walls: A_s=1000\ \mathrm{m^2} , k=0.20\ \mathrm{m^{5/2}/s} , z_c=0.20\ \mathrm{m} and initial head H_i=1.00\ \mathrm{m} . The analytical answer comes from the head-versus-time expression derived above, H(t)=\left(\sqrt{H_i}-kt/(2A_s)\right)^2 . Here k/(2A_s)=0.20/2000=0.0001\ \mathrm{m^{1/2}/s} , so after 3600 s

\displaystyle H=\left(1-0.0001\times3600\right)^2=0.4096\ \mathrm{m},\qquad h=z_c+H=0.609600\ \mathrm{m}.

For Modified Puls, the zero-inflow trapezoidal balance with S=1000h and O=0.20\sqrt H becomes 1000H_2+90\sqrt{H_2}=1000H_1-90\sqrt{H_1} for 900 s steps. Four steps from H_1=1.00\ \mathrm{m} give the same level, 0.609600 m, apart from rounding. In this zero-inflow case, discharge decreases linearly with time, so the trapezoidal rule integrates it exactly. This agreement confirms the calculation, but it cannot show how the timestep affects the result. For that, we need a case where discharge varies nonlinearly with time, so we keep the same orifice and change to the widening basin below.

A widening basin #

Take a basin 40 m long with vertical end walls and the transverse section shown above. At the outlet-centre elevation its width is 20 m; above that elevation each side widens by 5 m per metre of rise. Its width at head H=h-z_c is therefore b_0+10H , where b_0=20\ \mathrm{m} . Let L=40\ \mathrm{m} be its length. Multiplying length by width gives

\displaystyle A_s(h)=L(b_0+10H)=A_0+BH,\qquad A_0=800\ \mathrm{m^2},\quad B=400\ \mathrm{m}.

Here A_0=Lb_0 is surface area at the centre elevation and B=10L is the increase in area per metre of rise. Integrating the area above that level gives

\displaystyle S(h)-S(z_c)=A_0H+\frac{B}{2}H^2.

We retain O=k\sqrt{H} with k=0.20\ \mathrm{m^{5/2}/s} , no inflow, and initial head H_i=1.00\ \mathrm{m} . Our question for both methods is: what is the pool level after 3600 s?

The analytical answer #

Let \zeta denote head above z_c within the following integral. The drawdown integral gives the elapsed time to reach a lower head H :

\displaystyle t=\frac{1}{k}\int_H^{H_i}\frac{A_0+B\zeta}{\sqrt{\zeta}}\,d\zeta=\frac{2}{k}\left[A_0(\sqrt{H_i}-\sqrt H)+\frac{B}{3}(H_i^{3/2}-H^{3/2})\right].

The two terms follow by integrating A_0\zeta^{-1/2} and B\zeta^{1/2} . Solving the expression for t=3600\ \mathrm{s} gives H=0.445118\ \mathrm{m} , or pool elevation h=0.645118\ \mathrm{m} . This is our reference answer, obtained directly at 3600 s with no intermediate steps. The same expression gives the analytical head at any other time, which we will list beside the Modified Puls heads below.

The Modified Puls answer, step by step #

We take the Modified Puls equation derived earlier and run it on the same basin, so that we have a numerical answer to compare with the analytical one.

Modified Puls marches forward one timestep at a time, as in the diagram above: h_1 and h_2 are pool elevations at the beginning and end of one timestep, and their outlet heads are H_1=h_1-z_c and H_2=h_2-z_c , measured above the fixed outlet centre z_c=0.20\ \mathrm{m} . These subscripts label the two ends of whichever step we are calculating; H_i remains the initial head at the start of the complete recession.

With no inflow, I_1=I_2=0 , and the Modified Puls equation becomes

\displaystyle \frac{2S_2}{\Delta t}+O_2=\frac{2S_1}{\Delta t}-O_1.

For this basin, S(h)=S(z_c)+800H+200H^2 and O=0.20\sqrt H . Substituting them, multiplying by \Delta t/2 and cancelling the common S(z_c) gives

\displaystyle 800H_2+200H_2^2+0.10\Delta t\sqrt{H_2}=800H_1+200H_1^2-0.10\Delta t\sqrt{H_1}.

For each step, we know H_1 , solve for the lower head H_2 , and carry it into the next step as its new H_1 .

We first use a large 900 s timestep, so four steps reach 3600 s and the hand calculation stays short. The outflow changes slowly during this drawdown, so even this step lands within about 0.2 mm of the analytical level, as the comparison below shows.

Step 1, from 0 to 900 s. At the start, H_1=1.00\ \mathrm{m} . The right side is 800(1)+200(1)^2-0.10(900)\sqrt{1}=910 , so we solve

\displaystyle 800H_2+200H_2^2+90\sqrt{H_2}=910.

A numerical root calculation gives H_2=0.852122513\ \mathrm{m} . For example, in a spreadsheet, set the left side minus 910 to zero by varying the cell containing H_2 , keeping 0<H_2<H_1 . The left side increases with head, so there is only one root in that interval.

Step 2, from 900 to 1800 s. The head at the end of step 1 becomes the new starting head, H_1=0.852122513\ \mathrm{m} . The right side is now 800(0.852122513)+200(0.852122513)^2-90\sqrt{0.852122513}=743.841 , so we solve

\displaystyle 800H_2+200H_2^2+90\sqrt{H_2}=743.841,

giving H_2=0.709284064\ \mathrm{m} . Steps 3 and 4 repeat the same calculation, reaching 0.573049637 m at 2700 s and 0.445329371 m at 3600 s.

To see how the timestep affects the answer, we repeat the whole calculation with 450 s steps (eight steps, with 0.10\Delta t=45 in place of 90) and with 225 s steps (sixteen steps, using 22.5). The table lists the head at the end of every step of each Modified Puls run, beside the analytical head at the same time. A dash means that the run has no step ending at that time. Calculations retain full precision; the displayed heads are rounded to nine decimal places.

In every row, the Modified Puls heads sit slightly above the analytical head, and the gap narrows as the timestep shrinks.

Compare the two answers #

Adding the outlet-centre elevation z_c=0.20\ \mathrm{m} to the last-row Modified Puls heads gives the pool elevations at 3600 s:

\displaystyle \begin{aligned} \Delta t=900\ \mathrm{s}:\quad h&=0.20+0.445329371=0.645329371\ \mathrm{m},\\ \Delta t=450\ \mathrm{s}:\quad h&=0.20+0.445170929=0.645170929\ \mathrm{m},\\ \Delta t=225\ \mathrm{s}:\quad h&=0.20+0.445131243=0.645131243\ \mathrm{m}. \end{aligned}

The table compares them with the analytical elevation; errors are calculated from the unrounded values.

Halving the timestep reduces the error by approximately four. All three runs balance their discrete inflow, outflow and storage volumes, yet their predicted levels differ. Volume balance alone therefore cannot establish timestep accuracy.

For an actual pond, reduce the timestep and compare maximum level, outflow peak and timing, and time to the selected recovery level. Resolve changes in rainfall, tailwater, outlet regime and controls, and continue through the recession while retaining any ongoing inflow.

Conclusion #

We set out to see how the level-pool routing and Modified Puls equations in MSMA 2 follow from more general physical principles. Continuity, integrated over the whole facility, gives dS/dt=I-O , whether the internal water surface is horizontal or sloping. The momentum and Bernoulli relationships supply the outlet rating. Level-pool routing adds one approximation: a single horizontal pool level describes both the stored volume S(h) and the outlet head that gives O(h) . Trapezoidal integration over each timestep then gives the Modified Puls equation, and tabulating 2S/\Delta t+O against level gives its storage indicator.

For a hand check, the drawdown after inflow stops has a closed-form solution for simple basins. For vertical walls and a free orifice, it reproduces MSMA 2 Equations 2.14–2.15, and Modified Puls matches it apart from rounding, because the discharge then falls linearly with time. For the widening basin, the Modified Puls error falls by about four each time we halve the timestep, even though every run balances its volumes. The same drawdown calculation also tells us whether the pond recovers in time for the next storm, which a low peak level alone does not show.

Finally, level-pool routing is no longer applicable when one pool level stops fixing the outflow. This happens when:

• the water surface slopes materially, as along a long detention drain, so the depth at the outlet differs from the average depth that the stored volume implies;
• the tailwater rises or falls during the event, changing the driving head for the same pool level;
• the flow accelerates enough that the momentum terms we neglected in assuming a horizontal surface become significant.

The mass balance dS/dt=I-O still holds in these cases, but the outflow can no longer come from a single O(h) curve. We then need a more general hydraulic model that resolves the water-surface profile and the downstream level.