Measuring Water Flow of Rivers
Michael Monagan
Department of Mathematics
Simon Fraser University
mmonagan@cecm.sfu.ca
Introduction
One of my first summer jobs as a university student was working for the Hawkes Bay Catchment Board in Napier, New Zealand. Each day I would go out to measure the amount of water flowing down rivers in the province. We used water flow meters like the one shown in the figure below. We would go across the river in waders and take measurements. We did about 5 to 10 depending on the width of the river. At each point we would measure the depth and the velocity at that point.
After we were done for the day one of us would enter the data into a program to determine the flow. We did this for each river periodically. It was a fantastic job. Especially in the summer. One day the boss found out that I was doing computer science at university and he asked me to write a program to calculate the water flow from the data for their new computer which they had purchased. It had an HP plotter. He also asked me to generate plots of the cross section of the river with the data on it. So I had to figure out how to calculate the flow and how to use the HP plotter language. The point of this article is to explain what I did, that I got it wrong, and how to use Calculus to fix it, and, finally, now that I've done a bit of research, how the Geological Survey does this calculation. I am using this example for my teaching. It's nice. It's got a bit of everything. I think this is a great problem for an integral calculus class using Maple, or class on numerical methods. One useful tool that I've also learned is Maple's canvas facility which I'll use for the sketches in this article.
The Problem
Suppose we take n measurements (in meters) across the river at points where L and R are the left and right edges of the river. At each point we measure the depth in meters and the velocity in meters per second. So the depth and velocity at L and R is zero.
Because water flows faster at the surface of a river than at the bottom, there are different methods for taking a velocity measurement. We used two. For depths less than 0.5 m, we measured the velocity at 60% below the surface. For greater depths we took two measurements, one 20% below the surface and one 20% above the bottom and used the average of the two. I don't know the theory behind that. Referring to the two point method, the article by Thomas J . Buchanan and William P. Somers at
http://pubs.usgs.gov/twri/twri3a8/html/toc.html
from Techniques of Water Resource Investigations of the United States Geological Survey says "This method is based on many studies of actual observation and on mathematical theory . Experience has shown that this method gives more consistent and accurate results than any of the other methods except the vertical velocity curve method. (See p. 31.)" The article suggests also a three point method where you take three velocity measurements at depths 20%, 60% and 80% below the surface and use the weighted average Interesting. So we will have some data like this which we can visualize like the figure below.
restart; with(plots): xdata := [0,5,10,18,25,32,38,43,47,51]; ddata := [0.0,0.2,0.3,0.3,0.4,0.5,0.6,0.5,0.3,0.0]; vdata := [0.0,0.1,0.2,0.3,0.3,0.35,0.4,0.3,0.25,0.0]; n := nops(xdata)-2;
riverbed := plot( [seq( [xdata[i],-ddata[i]], i=1..n+2 )], style=line, color=blue ): readings := seq( plot( [[xdata[i],0], [xdata[i],-ddata[i]]], style=line ), i=1..n+2 ): for i to n do x := xdata[i+1]; d := ddata[i+1]; v := vdata[i+1]; s := sprintf("%4.2f m/s",v); xplot[i] := textplot( [x,0,cat("x",i)], align=above ); vplot[i] := textplot( [x,-(d+0.05),s], font=[HELVETICA,14] ); od: xtickmarks := seq( xplot[i], i=1..n ): velocities := seq( vplot[i], i=1..n ): display( [readings,riverbed,velocities,xtickmarks], view=[0..51,-0.7..0] );
What do you see in this figure? What I see is some trapezoids! Well, actually, I see trapeziums not trapezoids, because I'm from New Zealand. So I am thinking, for each section of the river, let's compute the flow on that section using Flow = Area of section TIMES Average velocity of water in that section We can approximate the area of each section as a trapezoid and approximate the average velocity on each section as the average of the two velocities on either side. Seems simple. In my program, let me use to indicate the position, depth and velocity at the left of side of the section and for the right side of the section.
n := nops(xdata);
TotalFlow := 0: for i from 1 to n-1 do (xL, xR) := xdata[i], xdata[i+1]; (dL, dR) := ddata[i], ddata[i+1]; (vL, vR) := vdata[i], vdata[i+1]; Area := (xR-xL) * (dL+dR)/2.0; # trapezoid Vel := (vL+vR)/2.0; # average velocity TotalFlow := TotalFlow + Area*Vel; od: TotalFlow;
I believe that is what I programmed for the HP plotter. The only problem with this is that taking the average of the two velocities and is wrong! Let's focus on one section. I'll sketch this with the Canvas feature in Maple. And put in some numbers to make a point.
The canvas facility allows you to insert lines, arrows, rectangles, circles, ovals, text (with font control), mathematics and as I've used, the mouse to sketch and/or write. I like the yellow paper. I actually use that for my lecture notes!
So what is the average velocity of the water on that section? It's NOT because the water is deeper on the right - there is more water flowing faster there. My first thought was oh, this is just a simple weighted average. We should take Average Velocity = . But this is not right either. For suppose . Then we would get for the average velocity which is not right. Hmm.
Riemann Sums
In a first course in integral calculus, we learn how to construct definite integrals as Riemann sums. Many instructors will also apply this to volumes of revolution. Shown in the figure I've divided the section up into subsections (the vertical lines). Each subsection is a little trapezoid. What is the velocity of the water in each little trapezoid and the area of each little trapezoid? Using linear interpolation to interpolate the velocity on the section, that is, for we have
Similarly we can approximate the depth for by a straight line as shown in the figure
Now to approximate as a Riemann sum we take the depth at positions on the interval and multiply by the velocity at that position and the width of the interval. This is just the area of each subsection multiplied by a velocity. So this just the definite integral for the area, , but we are multiplying by velocity. Therefore the flow on the section is given by
Flow =
This integral we can calculate in Maple:
x,xL,xR,dL,dR,vL,vR := 'x,xL,xR,dL,dR,vL,vR';
v := vL + (vR-vL)/(xR-xL)*(x-xL);
d := dL + (dR-dL)/(xR-xL)*(x-xL);
The flow on this section is given by .
Flow := factor( int( v*d, x=xL..xR ) );
To understand this formula as a weighted average of and let's collect in those variables:
(xR-xL) * collect( Flow/(xR-xL), [vL,vR], simplify );
Now we can calculate the Average Velocity = Flow / Area:
Area := (xR-xL)*(dR+dL)/2;
AveVel := collect( simplify(Flow/Area), [vL,vR] );
Expressed this in this way, it helps me see how the terms contribute to the average. I could not figure out how to do this transformation using Maple.
'AveVel' = 1/3*v1 + 1/3*v2 + d1/(d1+d2)*v1/3 + d2/(d1+d2)*v2/3;
So let's apply it and see if how much difference it makes.
xdata := [0,5,10,18,25,32,38,43,47,51]; ddata := [0.0,0.2,0.3,0.3,0.4,0.5,0.6,0.5,0.3,0.0]; vdata := [0.0,0.1,0.2,0.3,0.3,0.35,0.4,0.3,0.25,0.0];
n := nops(xdata); TotalFlow := 0: for i from 1 to n-1 do (xL, xR) := xdata[i], xdata[i+1]; (dL, dR) := ddata[i], ddata[i+1]; (vL, vR) := vdata[i], vdata[i+1]; Area := (xR-xL) * (dL+dR)/2.0; # trapezoid AveVel := vL/3+vR/3 + (vL*dL+vR*dR)/(dL+dR)/3; # average velocity TotalFlow := TotalFlow + Area*AveVel; od: TotalFlow;
Error := (5.3367-5.28625)/5.28625;
Another useful feature that I found in Maple is the numeric formatting features. If you right click on a number the context menu shows the option "Numeric Formatting". You can select from the following.
none
fixed (2 decimal places)
currency with control on how negative numbers are displayed e.g.
percentage
scientific
engineering with control on the number of digits
So on this example, the difference in the two approximations is 0.95 %. I suppose measurement errors might be more than 1% so this is not such a big deal.
In reflecting on this, if the river has a normal shape, so deeper in the middle, and if the river flows faster in the middle, then using the average of the two velocities and will underestimate the flow. In this example it's about 1%. I wonder, what is maximum difference between the two methods? That sounds like a good calculus exercise.
The MidPoint Rule
After I had done this I was thinking about numerical integration methods and I recalled that the Midpoint Rule is more accurate than the Trapezoidal Rule. The Stewart Calculus text says the error is half the error of the Trapezoidal Rule. Maybe we should be applying the Midpoint Rule to the river measurements? Let's replot the river using the midpoint rule. Actually, any kind of drawing like this may involve quite a bit of "fiddly programming". So I'll do the sections using blue polygons for illustration. I'll use the polygon command from the plottools package.
with(plottools);
for i from 2 to n-1 do x := xdata[i]; d := ddata[i]; xL := xdata[i-1]; xR := xdata[i+1]; poly[i] := polygon( [[(xL+x)/2,0],[(xL+x)/2,-d],[(xR+x)/2,-d],[(xR+x)/2,0]], color=cyan ); od:
polys := seq( poly[i], i=2..n-1 ):
plots[display]( [readings,riverbed,polys,velocities,xtickmarks], view=[0..51,-0.7..0] );
Well, theMidpoint Rule is what the article at
says the Geological Survey uses. This is the ``other'' solution I referred to in the introduction. But it too is inaccurate. For it too, by using the velocity at the middle of the section, incorrectly estimates the average velocity of the water on the rectangles. We can see this in a sketch
Shown in heavy black is the section at using the Midpoint Rule. Is the average velocity of the water in that rectangle 0.4 ? Well no, it's lower because more water is on the left. This Midpoint Rule would be correct if the measurements are taken at equally spaced intervals across the river. But that is rarely the case.
In the Classroom
This Fall of 2013 I assigned this problem to my MACM 204 students at Simon Fraser University. MACM 204 is a second year course entitled "Computing and Calculus" using Maple. We teach the students Maple and how to use it for solving problems in differential, integral and multivariate calcuclus, that is, Calc I, II and III. I described the problem in class and the formula for the flow on the river
Flow = Cross section of River times Average Velocity of water.
and how we divide up the river in sections, go into the river and take depth and velocity measurements, and I ask them, how should we calculate the flow. There's no point in giving them the answer because then they won't think. The students do come up with using trapezoids to approximate the area in class. I then gave them some data, asked them to plot the data (to recreate the figures in Maple) and calculate the flow. The course is about using Maple to do calculus so getting them to create a plot is part of that. For the flow calculation, almost all students use the trapezoids to approximate the area of each section and the average of the two velocities to approximate the average velocity on each section. Some try to generalized Simpson's rule to improve the estimate for the area. None have ever noticed the problem with taking the average of the two velocities. That probably means they didn't have time to reflect on what they are doing. But then, I didn't notice this either when I was a student looking at the problem.
