.U u . 1......1 .. .QCo: . a, .J...-..., .... .7...‘ .. . w 591' - mmuuv-mm-m 12767576 |l:fiflfl:lfl:ljfi[llfl l,,,. ‘ unm‘ - Michlg. W This is to certify that the dissertation entitled FINITE ELEMENT ANALYSIS OF HYDROLOGIC RESPONSE AREAS USING GEOGRAPHIC INFORMATION SYSTEMS presented by BAX'I'ER E. VIEUX has been accepted towards fulfillment of the requirements for Ph.D. degree in Agricultural Engineering ajor ofessor Date July 21, 1988 MSUi-n-‘fi‘ .- A ‘ r1 .nr’ . , i . 042171 RETURNING MATERIALS: MSU Place in book drop to remove this ‘ checkout from your record. FINES ”BRAKES will be charged if book is returned , =- after the date stamwd below. .1m;;.. 3 _’ «if I __ W; m; o c we: C(L\ (3.9r)rl"l‘ .2 ,, J: 1*; M,- (074902 r—a ‘7 _r ABSTRACT FINITE ELEMENT ANALYSIS OF HYDROLOGIC RESPONSE AREAS USING GEOGRAPHIC INFORMATION SYSTEMS BY Baxter Ernest Vieux The methodology developed in this research utilizes a Geographic Information System and the finite element Galerkin formulation to solve the kinematic wave equation for overland flow in a watershed. The watershed studied was number 4H, located in Webster County, Nebraska, and was operated by the USDA-Agricultural Research Service. The one- and two-dimensional forms of the equation were studied and the resulting outflow hydrograph was compared to an actual storm event, May 4, 1959. Rainfall excess was calculated using the Green and Ampt infiltration equation for an unsteady rainfall. Hydrologic response areas were formed based on slope with the aid of the ARC-INFO Geographic Information System developed by ESRI, Redlands, Ca. A finite element grid representing streamlines and equipotential lines was formed such that the direction of slope forms the streamlines of flow and the elevational contours form the equipotential lines. This results in nodal slope values perpendicular and parallel to the sides of the elements. Kinematic shock was avoided due to the use of nodal slope values. This formulation allowed solution of the overland flow equations FINITE ELEMENT ANALYSIS OF HYDROLOGIC RESPONSE AREAS USING GEOGRAPHIC INFORMATION SYSTEMS BY Baxter Ernest Vieux A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Agricultural Engineering 1988 (074C902 _) r, ABSTRACT FINITE ELEMENT ANALYSIS OF HYDROLOGIC RESPONSE AREAS USING GEOGRAPHIC INFORMATION SYSTEMS BY Baxter Ernest Vieux The methodology developed in this research utilizes a Geographic Information System and the finite element Galerkin formulation to solve the kinematic wave quuation for overland flow in a watershed. The watershed studied was number 4H, located in Webster County, Nebraska, and was operated by the USDA-Agricultural Research Service. The one- and two-dimensional forms of the equation were studied and the resulting outflow hydrograph was compared to an actual storm event, May 4, 1959. Rainfall excess was calculated using the Green and Ampt infiltration equation for an unsteady rainfall. Hydrologic response areas were formed based on slope with the aid of the ARC-INFO Geographic Information System developed by ESRI, Redlands, Ca. A finite element grid representing streamlines and equipotential lines was formed such that the direction of slope forms the streamlines of flow and the elevational contours form the equipotential lines. This results in nodal slope values perpendicular and parallel to the sides of the elements. Kinematic shock was avoided due to the use of nodal slope values. This formulation allowed solution of the overland flow equations for a watershed as a continuum rather than as a series of independent cascades. The method developed through this research provides a more accurate description of the hydrologic processes in a watershed. Through more accurate description of hydrologic processes insight is provided into transport phenomena of agricultural pollution such as pesticides and nutrients in surface and subsurface water as affected by overland flow and infiltration for an agricultu COPYR I GET Baxter E. Vieux, June 17, 1988 DED I CAT I ON I wish to dedicate this work to my wife Jean and children William, Ellen, Laura, Anne, and Kimberly who helped in their own way, and to my mother and father inspired my love of science at an early age. our all who ACKNOWLEDGMENT To my advisor, Dr. Vincent Bralts, I owe a debt of gratitude for the inspiration to complete this work. To the members of my guidance committee I must thank them as a team that as a whole provided more than guidance. To the other members of the guidance committee, I individually wish to thank Dr. Larry Segerlind for his insight into the finite element method and for his more than finite patience and encouragement: Dr. Donald Edwards for his oversight and greatly appreciated support; Dr. Jon Bartholic for his efforts to provide direction to and’ synthesis of this work, and Dr. Roger Wallace for his unrelenting insight into the physical problem formulation. I also wish to thank Dr. Walter Rawls, of the USDA-ARS Hydrology Laboratory for his valued assistance in the infiltration component of the Water Erosion Prediction Project (WEPP) used in this research. vi TABLE OF CONTENTS LIST OF TABLES. O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O OViii LIST OF FIGURES. I O O O O O O O I O O O O O O O O O O O O O O O O O O O O O O O O O O O O O O I 0 ix LIST OF SYMBOLS. O I O O O O O O O O I O O O O O O O O O I O O O O O O O O O O O O O O O O O O 0 .x1 I. INTRODUCTIONOOOOOOOOOOOOO...OOOOOOOOOOOOO0.00.00.00.0001 scope and Objectives.OOOOOIOOOOOOO0.00000......0.0.00.3 II. REVIEW OF THEORY AND LITERATURE ......................8 A. Hydrologic Modeling...............................8- B. Numerical Solution of the Hydrodynamic Equations........................................28 C. Finite Element Hydrologic Models.................32 D. Geographical Information Systems.................39 E Synopsis.........................................43 III. METHODOLOGY.........................................45 A. Research Objective and Approach..................45 B. Theoretical Development..........................47 C. Finite Element Model Formulation.................61 IV. RESULTS AND ANALYSIS.................................83 A. Finite Element Formulation.......................83 B. Geographical Information System Analysis.........87 C. Hydrograph Response.............................102 D. Discussion......................................126 V O CONCLUS I ONS AND RECOWENDAT I CNS 0 O O O O O O O O O O O O O O O O O O O I O 1 4 5 A. conCluSionSOOI 0.... O... O O... O... 000...... 00.... .145 B. Recomendations. I I O O O O O O O O I O O O O O O O O O O O O O O O O O O O O O 146 C. Future Application. 0 O O O O O O O O O O O O O O O O O O O O O O O O O O O O 147 VI. REFERENCESOOOOOOOO0.0.000...OOOOOOOOOOOOOOOOOOO0.0.0150 APPENDIX MOdeling DataOOOOO0.00.0......0.0.000000000000155 vii Table Table Table Table Table LIST OF TABLES sail PropertieSOOOOIOOOOOOOOOOOOO0.0.0.....0000088 Arbitrary Grid Finite Element Input Data........95 Hydrologic Response Area Nodal Coordinates and SlopeSOOOOOOOOOOOOOOOOOOOOOOOOOOO0.00.00.00.000100 Arbitrary Finite Element Grid Outflow..........104 Isotropic Finite Element Grid Outflow. ......... 123' viii Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure \IGU'IOFOJN CD 10 ll 12 13 14 15 16 17 18 LIST OF FIGURES Geographic Information System Definition of Areas With Similar Attributes..................5 Coordinate System Definition...... .....48 Control Volume Definition.....................50 Irregular Cross-section of a Channel..........57 Element and Node Numbering Convention.........67 Five Element Region...........................74 Isoparametric Elements and Integration POintSOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO ....80 Watershed 4H Landuse..........................89 Watershed 4H Soils Map........................9O Watershed 4H Elevation Map....................92 ARC-INFO TIN Slope Map........................93 Rainfall Excess for Three Soils ............... 96 Finite Element Grid of Arbitrary Spatial FormOCOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO0.0...0.0.97 Finite Element Grid Representing Hydrologic Response Areas Based on Slope................lOl Hydrograph From Arbitrary Finite Element GridOOOOOOOO0.0.0COOCOOOOOCOCCOOOIOOOQOO0.0.0105 Comparison of Outflows for the Arbitrary Finite Element GridOOOOOOOOOOOOOOOO0......0.0106 Four Node Quadrilateral Finite Element OUtflOVOOOOO0.0....IOOOOOOOOOOOOIOOOOO0.0.0.0111 Two Element System Outflow...................ll4 ix Figure Figure Figure Figure Figure Figure Figure Figure 19 20 21 22 23 24 A1 A2 Three Node Triangle Isotropic Qutflow........1l6 Three Node Anisotropic Triangle..............119 Isotropic Finite Element Grid................121 Hydrograph From Isotropic Finite Element GridOOOOOOO0.0.0.0000...0.0.0.000000000000000122 Comparison of Outflows for the Isotropic Finite Element GridIOOOOOOOOOOOOOOOO0......0.124 Spatial Distribution of Peak Runoff Flow RateSOOI0.000000000000000000000.0.00.00000000148 Comparison of the Finite Element and Analytic Solution for Flow Depth on a Single PlaneOOOOOOOIOOOOOOOOOOOOOO0.00.0.000161 Comparison of the Finite Element and Analytic Solution for Flow Rate from a Single Plane...OOOOOOOOOOOOOOOOOOOOO00....0.162 LIST OF SYMBOLS Cross-sectional area of flow. Nodal flow area vector for a channel of irregular cross section. T Jam (mm ) man . Transpose matrix of the derivatives of the shape functions with respect to x. Transpose matrix of the derivatives of the shape functions with respect to y. Speed of a gravity wave or (5/3)V. Capacitance matrix Differential element of the control surface, cs. Differential element of control volume, cv. Force vector. Pressure in elemental control volume = pghz/Z . Gravitational force. Froude number. Frictional resistance = pthfo . Gravitational acceleration. Water depth. n = 2 Ni(x,y)hi piecewise continuous shape function 1 approximation of the flow depth h. Nodal values of the flow depth h. Vector of time derivatives of the flow depth h. Rainfall impact overpressure or induced head. Infiltration rate. Kinematic wave number. Stiffness matrix. Shape function at the ith node. Shape function at the jth node. Transpose matrix of the shape functions. Fluid density. Nodal flow rate vector. Rainfall intensity. Rainfall excess = (r-i) . Elemental residual. Angle between the horizontal axis and the x- axis. Local coordinate system in the x-direction. xi Friction slope defined by either of the Chezy or Manning equations. Friction slope in x-direction. Friction slope in y-direction. Slope of the channel bottom. Slope in the x-direction. Slope in the y-direction. ulx,y,z,t), velocity in the x-direction. Depth-averaged velocities in the x-direction. n 8 Ni(x,y)ui piecewise continuous shape function approximation of the velocity u. Nodal values of the velocity u. n v(x,y,z,t), velocity in the y-direction. Depth-averaged velocities in the y-direction. n X Ni(x,y)vi piecewise continuous shape function 1 approximation of the velocity v. Nodal values of the velocity v. Velocity vector. Dot product of the velocity vector to the normal unit vector of the control surface. w(x,y,z,t), velocity in the z-direction or flow width of the element. Time step. . Distance increment or element length. Coordinate in the primary flow direction. Coordinate perpendicular to the x-direction. Vertical coordinate or direction perpendicular to the channel bottom. Momentum correction factors for main and lateral flows. Momentum correction factor for raindrop terminal velocity. Eta natural coordinate corresponding to the global y-direction. Ksi natural coordinate corresponding to the global x-direction. Eigenvalue. Elemental control volume or mean terminal velocity of fall of raindrops. Two-dimensional element domain. Angle with the respective direction. Angle of inclination between the terminal velocity vector and the vertical axis. xii I . I NTRODUCT I ON Agricultural pollution by pesticides and nutrients threatens both surface and subsurface waters due to hydrologic transport of these contaminants. The locations within a watershed that produce similar contributions to surface and subsurface waters may be termed hydrologic response units or areas. To successfully reduce subsurface water contamination, reduce sedimentation, and' control erosion, land treatment measures must be focused on those geographic areas that will yield the greatest mitigation. Hydrologic modeling has recently been the subject of more and more attention in addressing watershed management. Modeling may be mathematical, if described by a mathematical equation, or physical if a scale model is built to represent dimensional similitude to the actual watershed. In either case the model is a conceptualization of the actual watershed. Mathematical modeling generally seeks to define the mathematical relation between a set of independent variables and a response or dependent variable. The twentieth century has witnessed a rapid acceleration in the quantitative modeling of physical processes. The mathematical description of natural phenomena in the hydrologic cycle is not, however, a child of this century. The history of quantitative hydrology has 2 been marked by important milestones. Achievements of each scientist have served as the foundation for advances by the next scientist. 0f the classical period the most notable treatise (in 66 AD) that makes mention of the hydrologic cycle was Vitruvius' ”Ten Books of Architecture" (Morgan, 1960). Vitruvius' description in Book VIII on how to find water, we read The valleys among the mountains receive the rains most abundantly, and on account of the thick woods the snow is kept in them longer by the shade of the trees and mountains. Afterwards, on melting, it filters through the fissures in the ground, and thus reaches the very foot of the mountains, from which gushing springs come belching‘ out . Though not without some misconceptions, Vitruvius essentially understood the origin of springs and groundwaters. Not until the seventeenth century did ideas of quantitative hydrology emerge. Biswas (1968) presented the following matter on the beginnings of quantitative hydrology. Pierre Perrault anonymously published the book 23 I'Origine des Fountaines in 1674 (Biswas, 1968). In this work he calculates the quantity of water that would accumulate from the rainfall in the catchment of the Seine River, France. He found that a sixth part of the rain and snow water is necessary to make the river run continually throughout the year. For the first time experimental evidence proved the pluvial origin of rivers. The greatest contribution of Perrault and other contemporaries-- including Edmond Halley, who calculated the volume of water t required to supply the rivers in the evaporation- precipitation cycle, and Edme Mariotte who expanded on the pluvial origins of rivers--was that they proved their hypotheses through quantitative methods. Watershed hydrology can be treated as either a lumped or distributed parameter model as well as by a stochastic or a deterministic method. A lumped model tends to utilize the average of a set of independent variables that represent a sub-basin or an entire watershed. A distributed model utilizes the spatial location of the‘ independent variables and computes the dependent variable directly at the spatial location of each independent variable. Distributed models represent spatial distribution as a set of grids or a series of finite elements, each with its own physical properties. The degree to which physical properties--e.g., soil infiltration parameters, surface roughness, or slope-—are averaged, represents the degree to which the model is lumped. The direct computation utilizing these parameters is a deterministic, method as opposed to a stochastic method. Scope and Objectives Early modeling efforts neglected the spatial variability of the independent variables. Lumping of these parameters does not allow visualization of the spatially distributed runoff and infiltration processes. 4 Distributed parameter models, which can readily handle complex geometry and spatial distributions of the independent variables, have typically used approximations of the actual watershed by using arbitrary grids, planes or elements. The observation of the spatial and temporal distribution of infiltration and runoff is one of the benefits claimed by proponents of distributed models. In actual practice, however, the volume of input data prevented accurate representation of the spatial character of the input data. Proper location of land treatment measures required knowing the spatial and temporal' distribution of the infiltration and runoff processes. In order to achieve an accurate characterization of the spatially distributed input parameters, an improved method of defining hydrologic response areas is needed. The value of processing a watershed by a Geographic Information System (GIS) is evident within the context of the finite element method as shown in Figure l. The GIS serves as a spatial data management system. Soils maps, landuse and other geographic data are represented in a digital media using a common coordinate system. The first step in the finite element method is to define elements over which the differential equations of overland flow may be . .mmuznwuuu4 LmHMEMm suws mmou< mo :owuwcwuon Emummm cowumau0mcH oLchLmOmo .H ousmwm mmam 4:5 Aazo_.2dm.l,iii- , 131 iazm;§2am mmwfimflsm m5¢¢cwfiotc£ HmowmoHoafimg 7 flwzmhmaaz. magnmtlil mwzficmm Mmmom mmém 4:5 mmhdjmz 024 ml: 20 ZQBUDDOmm m5 integrated. The efficiency of the Geographic Information System is realized when areas of like attribute values, such as infiltration parameters, slope, surface roughness, and so on, are aggregated within a boundary. This boundary forms a polygon consisting of vectors or arc segments. These polygons then become the set of finite elements used to model the runoff process. The purpose of this research was to develop a method that more accurately predicts the outflow hydrograph resulting from runoff during a rainstorm event for a watershed with spatially non-uniform parameters. This) method utilized the finite element method to compute the runoff rates and a Geographic Information System to more accurately represent the spatially distributed parameters for use in the computation of the runoff. The specific objectives for reaching this goal are as follows: 1) Define hydrologic response areas that exhibit similar soil infiltration parameters, surface roughness, and slope. 2) Apply the finite element method to the specific hydrologic response areas to compute and route the overland flow to the outlet. 3) Compare the accuracy of the outflow hydrographs to the actual outflow hydrograph for a given rainstorm event for the following two cases: i) A finite element grid that is of an arbitrary spatial form. ii) A finite element grid formed from hydrologic response areas defined by the Geographic Information System. It was hypothesized that the combination of the Geographic Information System and the finite element method would result in a mathematical model that predicted the. outflow hydrograph from a finite element grid formed of hydrologic response areas more accurately than from a finite element grid that was of an arbitrary spatial form. Validation of this hypothesis was accomplished if the method developed through this research more accurately predicted the actual outflow hydrograph from a watershed of non-uniform, spatially distributed parameters such as infiltration parameters, surface roughness, and slope. i) A finite element grid that is of an arbitrary spatial form. ii) A finite element grid formed from hydrologic response areas defined by the Geographic Information System. It was hypothesized that the combination of the Geographic Information System and the finite element method would result in a mathematical model that predicted the_ outflow hydrograph from a finite element grid formed of hydrologic response areas more accurately than from a finite element grid that was of an arbitrary spatial form. Validation of this hypothesis was accomplished if the method developed through this research more accurately predicted the actual outflow hydrograph from a watershed of non-uniform, spatially» distributed parameters such as infiltration parameters, surface roughness, and slope. II. REVIEW OF THEORY AND LITERATURE Since the seventeenth century, the science of hydrology evolved together with mechanical and industrial developments. The twentieth century and recent decades in particular, have seen great strides in the science of hydrology. The industrial age, which :ffected environmental degradation also necessitated more and greater strides in hydrologic methods. One such method, the mathematical model, through recent technical elaborations has allowed the direct modeling of spatially distributed hydrologic processes. Hydrologic processes, when described by physically based equations, are termed deterministic. This review expounds the theory in the literature on deterministic hydrologic models, particularly those utilizing finite difference and finite element methods to describe the spatially distributed hydrologic process of a rainfall storm event over a watershed. A. Hydrologic Modeling The physically correct representation of the surface runoff and infiltration processes in a watershed, field, or plot depends on many factors. To categorize these factors, several distinctions should be made in general as to the 9 l modeling process that seeks to represent the physical process. Abstraction is the mind's attempt to disengage the essence of something from that which is non-essential or only incidental to its make-up (V.E. Smith, 1950). The order of abstraction required in mathematical modeling is to consider the mathematical relation between cause and effect for the abstracted, conceptual model. The conceptual model strips the processes that are considered incidental or nonessential. The mathematical model then describes those essential processes contained in the conceptual model. Considering the complexity of the real“ world, it is necessary to use a conceptual model in order to successfully apply the mathematical model. This however is not without drawbacks, considering the interdependence of the many hydrologic processes. Modeling a particular process in the absence of another that is affecting the modeled process may result in a physically invalid model and would represent a poor choice of a conceptual model. These drawbacks notwithstanding, we will examine first some mathematical models that seek to model deterministically the physical process of the rainfall event. This process may include infiltration, overland flow of the rainfall excess, and channel routing of the lateral inflow from. the overland flow portion of the watershed. Smith and Woolhiser (1971) developed a mathematical model that modeled a coupled system of two complex, natural processes of an elemental watershed. The 10 conceptual model included only infiltration and overland flow. Channel flow was not considered because the conceptual model was limited to the upland portion of the watershed where channel flow is not well established. The infiltration model provides insight into the process by which rainfall becomes either runoff or subsurface water. The kinematic equations provide insight into the depth and velocity of the runoff as it accelerated down the watershed. The infiltration and the kinematic equation were coupled mathematically such that the rainfall excess as defined by the infiltration model was the boundary. value for solution of the kinematic equation. A distributed, deterministic system results when the inherent spatial nature of the processes are preserved in the solution method. A model such as Smith and Woolhiser's provides the opportunity to model the outflow of the watershed and, more importantly, the spatial and temporal distribution of the runoff-infiltration process within that watershed. 1. Distributed Parameter Models Huggins and Burney (1982) observed that hydrologic modeling is most differentiated by the manner in which parameters or input values are handled. Lumping or averaging certain parameters yields a lumped parameter model. Distributed parameter watershed models treat the individual input parameters directly without lumping. Such models avoid the errors caused by averaging of nonlinear 11 i variables or threshold values (Huggins and Burney, 1982). Distributed parameter models are relatively new and are the subject of intense research. The principal advantage of distributed models is that the geographical variation of data within the watershed is preserved. Furthermore, with measured parameters, ungaged watersheds may be investigated for the effect of landuse changes. Finally, water quality simulations on a distributed basis identify the source areas within a watershed of hydrodynamic transport of pollutants. Distributed parameter watershed models are more complex, require more computing time and increased input data. The value of the knowledge derived from distributed models must be considered vis-a-vis the increased time and costs of developing and using this class of model. High speed digital computers, however, obviate the restrictions on increased computing time and complexity allowing more and more complex, distributed models to be used. Huggins and Monke (1966) developed the ANSWERS model, a distributed parameter model that uses a grid as a method of preserving the geographic heterogeneity of the input parameters. Each cell of the grid represents the hydrologic unit over which the model equations are solved using the finite difference method. The ANSWERS model uses a fundamental method of computing distributed processes. The continuum of the watershed is resolved into discrete elements and 12 preserves the spatial heterogeneity of the input parameters. 2. Hydrodynamic Model Approach The hydrodynamic approach proposed by C.L. Chen and Ven Te Chow (1971) considered watershed hydrology as a distributed continuum, where the hydrodynamic principles of fluid flow apply. The solution of the hydrodynamic equations yields a deterministic model capable of defining distributed flow velocities and depths over the watershed. The hydrodynamic equations have been derived and) solved by various methods. There are two distinct categories of flow in a watershed--overland flow, characterized by shallow flow, and channel flow, characterized by well-defined channel geometry. The boundary between these two flows changes with time and distance and therefore are modeled with difficulty. Chen and Chow formulated a comprehensive watershed- flow model. They classified watershed hydrology by a molecular approach, a microscopic hydrodynamic approach, and a macroscopic hydrodynamic approach. The microscopic and macroscopic hydrodynamic approaches both derive the Navier-Stokes equation of motion for fluid flow with suitable boundary conditions. The difference between the micro- and the macro- approaches is that the latter utilizes averaging of variables in certain flow directions in order to simplify the Navier-Stokes equation. l3 ‘ The approach of Chen and Chow formulated the equations of continuity and momentum for flow of a Newtonian fluid in a three dimensional space. A Cartesian coordinate system is used whereby the average velocity is taken parallel to the ground surface with the x- and y- directions along the bottom of the flow. Temperature variations are not considered, so the fluid is assumed to be of a homogeneous viscosity. Their derivation is summarized as follows. a. Conservation of Mass The conservation of mass is derived by integrating‘ in the vertical direction such that the differential variation in flow is considered only in the x- and y- directions with vertical variation as defined by the depth of flow h. ah + 3(uh) + 3(vh) = r - i [1] TE TIE— ST— where :3’ I the depth of flow measured in the vertical axis 2' u = depth-averaged velocity in the x-direction, v a depth-averaged velocity in the y-direction, r a rainfall, and i a infiltration, positive when moving out of channel, negative when moving into channel. The coordinate system used assumes that the primary flow directions are parallel to the ground surface. This requires that the z-direction makes an angle 6 inclination 14 with the force of gravity. Equation [1] is in terms of vertically oriented variables h, r, and i. The x- and y- directions are not orthogonal with the vertical z-axis. b. Conservation of Momentum The conservation of momentum is derived by balancing momentum and forces for an elemental control volume. The resulting equation is the Navier-Stokes equation and is a fundamental equation of fluid mechanics. Chen and Chow (1971) begin with the Navier-Stokes equation and simplify it according to appropriate boundary conditions and‘ assumptions. This rather convoluted approach yields the general momentum equation for watershed flow in the x- direction. 3(AV) + 3(8AV2) - BrrATcosvx - BLVqL 5t t = gAsinGx - 9%[A(hcoszez + h*)] - gASfx [2] x where ex'y'z - angle with the respective direction, h a depth of the centroid of the cross- sectional area, h a rainfall impact overpressure or induced head, A - cross-sectional area of flow, Sfx a friction slope in x-direction, 8,8L - momentum correction factors for main and lateral flows, at - momentum correction factor for raindrop 15 terminal velocity, A a mean terminal velocity of fall of raindrops, and 0x = angle of inclination between the terminal velocity vector and the vertical axis. Because they are applicable to both overland and channel flow, equations [1] and [2] for the conservation of continuity and momentum, form the description of the watershed problem. They are the complete dynamic form of the shallow water equations. The continuity or conservation of mass and the one-dimensional, conservative form of the conservation of momentum are commonly known as- the St. Venant equation derived by St. Venant in 1871. The general form, ignoring the momentum of the rainfall impact and overpressure, is expressed by [3] as a set of equations aav + V3A + 3A = q 3x 5x 5t and vav + av + a('A) = g(S - s ) - Vq [3] 5x St % 5x 0 f A where cross-sectional flow area, > I depth of flow, mean water velocity, lateral inflow per unit length of channel, *< J: < «H u a distance from the water surface of the centroid, So = channel slope, and 16 Sf = friction slope. These are the conservation of mass and momentum in the x- direction. Abdel-Razaq (1967) applied the finite difference solution to the above St. Venant equations. Brutsaert (1971) experimentally verified the solution, and Yen (1973) recapitulated the open channel flow equations emphasizing the origin of the St. Venant equations as the special form of the conservation of mass and momentum. Previous researchers had simplified these equations before formulating a mathematical model. The significant- contribution is that these researchers applied these equations to a watershed flow domain. The solution of the partial differential equations requires initial and boundary conditions. If the flow domain is considered to be the watershed as a whole, then only one boundary condition on the watershed divide and one at the outlet are needed if the flow is subcritical. If supercritical, then only one boundary condition is required at the watershed divide to give a unique solution. If the flow domain is divided artificially into overland flow and channel flow, then there are two subdomains separated by an internal boundary. This boundary is not known prior to the solution yet is required in the solution. This difficulty may be removed by estimating a likely location of such an internal divide and treating the system of channel and overland flow as l7 uncoupled. This approach allows the computed values of the overland flow to act as lateral inflow for the channel flow equation. 3. Kinematic Wave Equation A comprehensive treatment of the kinematic wave equation was examined by Woolhiser (1975) as a means of computing hydrographs with the assumption of both the Chezy and Manning equations as friction relations. Grace and Eagleson (1966) developed the full dynamic equations using the control volume technique for conservation of momentum and mass. Normalized equations were then established and an order of magnitude analysis produced. This approach allows simplification of the governing equations by discarding small order terms maintaining similarity between model and prototype. Brutsaert (1968) obtained an analytical solution to the shallow water or St. Venant equations. This was done within a small solution domain bounded by the forward and backward characteristics and the x-axis from 0 to L of the plane. The series solution provides good initial values for numeric analysis of overland flow in the initial stages of the hydrograph development. Of importance is that for large slopes, for a large roughness coefficient of the plane, or for very small constant lateral inflow, the series solution reduces to the kinematic wave equation, which is an approximation to the full dynamic equation. Following is an illustration of the simplifications that 18 i .can be made to the full dynamic equation, ignoring lateral inflow: vav + av + a( A) 8 g (S - Sf) [4] 5:. at irz—x ° or s . s - 3(yA) — v av - av [51 LIV—J % X '5'; FE \KINEMATIC J \r L DIFFUSION FULL DYNAMIC The terms in equation [5] represent the possible‘ analogies for the modeling of the momentum equation. Through simplification or elimination of low order-of- magnitude terms, the significant terms are used in computing the fluid's dynamic behavior. The elimination of these terms in the above analogies introduces an error in the solution. The magnitude of the error dictates the acceptability of the analogy used under the physical conditions of the problem. Woolhiser and Liggett (1967) examined in depth the errors introduced by the kinematic wave analogy applied to the full dynamic equation for the rising hydrograph. They discovered that the kinematic analogy can be applied within a certain range of input parameters. To quantify this range, the equations of continuity and the momentum equation are first normalized, constituting a nondimensional form of the full dynamic or shallow water equations. Then, 19 . the kinematic equation or continuity plus the Chezy relation is normalized, yielding the equation of motion without the momentum terms. Kibler and Woolhiser (1970) described the kinematic cascade as a hydrologic model. The concept here is to reduce not only the full dynamic equation momentum to a simplified form but also the geometry of the watershed to simple geometric cascades, through which the overland flow is routed to the outlet of the watershed. This approach results in a distributed parameter watershed model. The kinematic equation used is the Chezy equation together) with the equation of continuity. u a ch [6] and 3h + auh = q [7] If I?" N = 3/2 for wide channels, a = CJS, where C is the Chezy roughness coefficient and S is the slope, and h = the depth of overland flow. In order to quantify the range of input parameters for which the kinematic analogy is applicable, equations [6] and [7] are normalized or made dimensionless. Dimensionless Form Equations [6] and [7] are normalized to yield 20 ‘ dimensionless equations dividing each term by a characteristic length or time. The Chezy equation is substituted into the equation for continuity and then solved using the method of characteristics. The dimensionless parametric equation is N-l 3h* + Bha 3h* a q* [8] 5E* 31* where q* a normalized lateral inflow q, h* = normalized depth, n B = N Lk / Z 1i i=1 11 = length of plane i in feet, Lk [a normalizing depth for plane k, N a Chezy parameter defined as above. The dimensionless characteristic equations derived from the above normalized, nondimensional equation is N-l dX* = Bh* [9] 3E. and db: = q* [10] 3E* The earliest solution method was the method of characteristics, which had its origins in the nineteenth century. This method is also related to a separation of variables technique. In both instances, the partial 21 ‘ differential equation is reduced to ordinary derivatives. These are then solved graphically, numerically, or by a combination of both in the x-t domain (Henderson, 1966). The value of this method is that it presents information on the partial differential equation solution. Looking at the ordinary differential equations for the characteristics, we obtain a velocity dx/dt at which information can be transmitted through the system. This velocity leads to the Courant condition, which states the limit at which a disturbance can move. This velocity is the celerity of a gravity wave. Kibler and Woolhiser (1970) made a thorough analysis of the kinematic cascade as a distributed parameter mathematical watershed model. One difficulty they encountered was the numerical phenomenon of the kinematic shock. When a disturbance occurs in an open channel system, its propagation may be described by the method of characteristics. Simple wave theory explains that when there is a change in slope between planes in the cascade, a shock or wave front is propagated within the system. The shock represents a numeric difficulty in the computation of the hydrograph. The shock parameter used to predict occurrence is defined as P5 = Yk-l “k-l > 1 [11] wk Gk w a width of the k and k-l 22 _ planes, a = Chezy coefficient or C S. Observing the shock parameter inequality will ensure that shocks will not be propagated along characteristic lines emanating from the upstream boundary and the line x = 0. Morris and Woolhiser (1980) re-examined the validity of the kinematic assumption under partial equilibrium. They found that the full dynamic or diffusion equations should be used for flat grassy slopes. The criterion Folk 3 5 should be observed when using the kinematic analogy. The physical significance of Folk or SOLO/Ho represents the ratio between the difference in elevation between the top and bottom of the plane (SOLO) and the normal flow depth at the downstream boundary(Ho). This criterion provides a convenient method of deciding the validity of the kinematic analogy. Earlier work (Woolhiser and Liggett, 1967) suggested that at full equilibrium the criterion k 3 SOLO/HO > 10 should be observed when using the kinematic analogy. If the kinematic number k is less than 10, then the full dynamic equations should be used. This condition does not normally occur in agricultural watersheds but may occur in urban areas where short, smooth watersheds with low lateral inflows prevalent. 23 4. Infiltration When modeling overland flow in a pervious watershed, it is necessary to calculate the rainfall excess that results when the rainfall intensity exceeds the infiltration capacity of the soil and the surface storage. The rainfall excess is treated as lateral inflow in the overland flow equations. The right-hand side of equation [7] is the lateral inflow and can be viewed as the forcing function. In order to characterize an infiltrating natural watershed, it is necessary that l. The conceptual model accounts for all processes of interest, e.g., unsteady rain, snow melt, etc. 2. The mathematical model adequately describes the conceptual model. 3. Soil properties within the watershed are taken ' into account. 4. Input parameters can be obtained for the domain of interest and successfully applied towards the solution. The USDA-Soil Conservation Service developed a procedure to estimate direct runoff from ungaged watersheds. Rallison (1980) gives a detailed synopsis of the development of this procedure from its inception to its final form and application to ungaged watersheds. Work on this procedure began in the mid 19505 in response to the passage of the Watershed Protection and Flood Prevention Act (P.L. 83-566). Due to the work authorized by this act, SCS anticipated the need for a simplified method of hydrologic computation. Based on extensive analyses of 24 ‘ gaged, experimental watersheds and infiltrometer studies, a relation between rainfall and runoff was developed (Andrews, 1954, and Mockus, 1949). The basic relation was derived by plotting the accumulated natural runoff versus the accumulated rainfall. It was observed that the relation is asymptotic to a line at a 45 degree slope. This shows that the runoff rate approaches rainfall rate as the accumulation of both continues. Also, the difference between rainfall and runoff, the maximum retention, approaches a constant value. Rainfall intensity and the surface sealing effects of rainfall were not considered in- the analysis. The basic hypothesis is 5': Q [12] 3 15,, where F = actual retention of precipitation during a storm, S potential maximum retention, Q a direct runoff, Pa - rainfall after initial abstraction. Curve numbers are related to S by CN = 1000 . [l3] STIU The range of curve numbers for a watershed are due primarily to variations in storm duration and intensity. In the original analysis the CN chosen was an average of a 25 ‘ range of values for the particular soil-cover complex. For a particular storm, the CN may possibly fall outside the range originally established. Infiltration patterns are not accounted for within a storm period since no time variation was incorporated into the procedure. Consequently, the curve number method is not applicable to modeling infiltration under a variable intensity rainfall. Infiltration equations that respond to variations of rainfall intensity have since been developed. Mein and Larson (1971) presented an extensive analysis of infiltration modeling as it relates to watershed modeling.- They classified models as being empirical, theoretically derived algebraic, or soil moisture flow models. The Green and Ampt equation falls into the category of theoretically derived algebraic equations. The approach of Mein and Larson was to predict the time between the inception of rainfall and the inception of runoff. Their methodology was to modify the decay function of the infiltration capacity as defined by the Green and Ampt equation. This modification is necessary in order to account for the infiltration that occurs prior to surface ponding. The result of this effort was a simple model that relates infiltration to a constant rainfall intensity, homogeneous soil properties, and uniform initial soil moisture. Brakensiek and Onstad (1977) presented a parameter estimation for the Green and Ampt infiltration equation. They observed that if watershed runoff is to be accurately 26 ‘ _predicted, an accurate estimation of the infiltration capacity is needed. The sensitivity of runoff rate and volume was very sensitive to the fillable porosity and the hydraulic conductivity. A runoff model using the kinematic wave equation for direct runoff from a plane predicted the effect of varying the Green and Ampt parameters. The volume of runoff, for the conditions modeled, was most sensitive to fillable porosity. This is termed by some researchers to be the initial soil moisture deficit. The sensitivity expressed as a ratio of the dependent to the independent variable is as follows: Qv = 5.79 C gp - 3.47 where Op peak runoff rate, Qv = volume of runoff, C fillable porosity. For each 1% error in fillable porosity there is a corresponding 5.79% error in the volume of runoff. The sensitivity of hydraulic conductivity is gv a 4.41 x = 2.68 5%? 27 where Op 8 peak runoff rate, Qv a volume of runoff, K a hydraulic conductivity. The sensitivity of these parameters indicates which parameters must be estimated with the greatest accuracy. The sensitivity indicates which parameters most easily bring the model into agreement with the observed event. The least sensitive parameter is the wetting front suction head. This would indicate that in a parameter optimization scheme, convergence may not be as rapid for the wetting. front suction head parameter as for the others. Chu (1978) studied infiltration under an unsteady rain. His study extended the application of the Green and Ampt equation to predicting the infiltration under an unsteady rainfall intensity. During an unsteady rainfall event, the intensity may recurrently shift from falling below to exceeding the infiltration capacity. The purpose of Chu's study was to extend the Green and Ampt equation to account for variable periods of rainfall intensity. The accomplished purpose is a transformed time scale that allows the computation of the infiltration with time under an unsteady rainfall. Agricultural management affects the infiltration process. Rawls, Brakensiek, and Soni (1983) present guidelines for predicting the effects of tillage on the Green and Ampt parameters. Using the soil texture data, a 28 ‘ regression analysis was made to predict the increase in porosity due to tillage of a given soil. With time, the increased porosity decreases due to consolidation. The estimated change in porosity is a basis for estimating the change in the Green and Ampt parameters of capillary pressure of the wetting front and the hydraulic conductivity parameter, which is a fraction of the saturated hydraulic conductivity. The effect of tillage is accounted for by this procedure. Brakensiek and Rawls (1983) presented the effects of surface sealing or crusting on the Green and Ampt- parameters. A two layered, hydraulic conductivity is assumed to represent crusting. The wetting front capillary head is assumed to be that of the pre-crusted soil. The crusting thickness is assumed to be 0.5 cm. The effective hydraulic conductivity is calculated for pre-ponded and post-ponded periods during a rainfall event. The model predicts the infiltration to within an order of the effects of a crust formation under a rainfall simulator. B. Numerical Solution of the Hydrodynamic Equations The solution of the full dynamic equations poses significant problems due primarily to the nonlinearity of such terms as uau/ax. Such problems give motivation to identify the domain in which reasonably accurate solutions to the simplified equations may be obtained. Further motivation to simplify the full dynamic equations is the 29 ‘ difficulty to provide the appropriate boundary _conditions and incorporate them into the solution method. The finite difference method has achieved extensive use in the computer solution of the full dynamic and simplified equations of fluid flow. Abdel-Razaq (1967) provided a finite difference solution to the surface runoff problem defined by the conservation of mass and the one- dimensional conservative form of the conservation of momentum equations. The method of characteristics is a semi-graphical solution of the full dynamic and the simplified equations - of fluid flow. Henderson (1966) gives an in-depth procedure for the solution of the full dynamic equation and the kinematic wave equations for channel flow. The method of characteristics also yields information on grid spacing in the finite differencing domain. This is related to the partial differential equation theory. The goal of these methods is to reduce the partial differential equations to ordinary and the ordinary to a linear system of equations amenable to solution. 1. Finite Difference Method The finite difference method seeks to replace a continuum with discrete points between which the differentials are approximated. By replacing the partial differential terms with a finite difference approximation, the continuous domain is replaced by a network of isolated, discrete points. This procedure reduces a continuous 30 problem to an approximating eigenvalue matrix amenable to solution (Crandall, 1956). The solution of the finite element formulation in time is most commonly done by the finite difference method. For example, the equation [C]{A} + [K]{Q} = {F} [14] requires that a temporal solution of the time derivatives {A}, aA/at be computed. This may be accomplished by writing the finite difference form of the time derivative. as aA(£) A(a) - A(b) " [15] 3t At The mean value theorem indicates that 5 must lie within the interval of a Z 5 I b Not knowing where 5 lies within this interval we must define the parameter (E - a) At and use it as A(£) = (1-0)A(a) + GA(b) [17] Substituting [16] into [17] yields the relation for g = t as A = (l-e)Aa + GAb [18] 31 5 Similarly for the right-hand force vector {F} F = (l-O)Fa + GFb _ [19] where the b values are new values at the next time step. The four common values adopted for e are 6:0, §=a, the forward difference method. Gal/2, §=At/2, the central difference method. 0:2/3, §=2At/3, the Galerkin's method. 881, 5-b, the backward difference method. Equations [14] through [19] form the basis for the finite difference solution in time commonly used in conjunction with the finite element method 2. Finite Element Method The solution of the hydrodynamic equations for fluid flow has encompassed a wide variety of disciplines, including 1. surface water equations for tidal estuaries, 2. boundary layer equations, 3. Navier-Stokes and St. Venant equations for closed and free surface fluid systems, 4. meteorological dynamics, and 5. groundwater flow. The application of the finite element method is used extensively in fluid mechanics, as evidenced by the large volume of literature dealing with the solution of the Navier-Stokes equation especially in mechanical engineering 32 \ applications. The application of the finite element method to watershed or catchment hydrology is the subject at hand and is discussed herewith. C. Finite Element Hydrologic Models The finite element method was first used by Guymon (1972) in the solution of the hydrodynamic equations for free surface water flow. He solved the equations assuming a constant depth over a region using the variational principle. He concluded that the finite element method was a suitably efficient solution technique for surface water) problems. Researchers have typically approached the watershed or catchment problem as a model composed of two distinct parts--overland and channel flow. Judah (1973) applied the finite element method to this two-part model. He used the kinematic simplification of the momentum equation and the continuity equation and the Manning equation as the friction relation. Judah's application of the Galerkin principle utilized linear shape and weighting functions to approximate depth and velocity. The element was one dimensional, having an average slope in the direction of flow. Rainfall excesses were not modeled because the model was tested for storms for which the rainfall excess was already known. Several watersheds, both experimental and natural, were modeled. Close agreement was generally found in the simulations of 33 ‘ the outflow hydrographs. It should be emphasized that rectangular strips (one-dimensional elements with an average top width) were used to represent the watersheds. The Rocky Run Branch Watershed in Brunswick County, Virginia, was subdivided into nine finite elements representing a total drainage of 555 acres. Taylor (1974), using a Navier-Stokes formulation for the momentum equation, derived the two-dimensional form for watershed flow with the kinematic wave assumption. This application of the Galerkin method resulted in a coupled set of equations of the form [Mltqlt = {Flt [20] where the vector {q} = [21] :J‘ 0 [23] where the vertical coordinate, N II 3‘ II the water depth, fluid density, ‘0 II IO II gravitational acceleration. 48 Y 12 +3qu pig—i ‘7": E7 /v- EVA ‘5y‘§ Ii Figure 2 Coordinate System Definition 49 The coordinate system is defined as follows: x = horizontal direction of primary flow, parallel to the bottom of the flow surface. y = horizontal direction perpendicular to primary flow, parallel to bottom of the flow surface. 2 = vertical direction. The time averaged local flow velocities in the elemental volume are: u t Bu Ax velocity in x-direction on up- and 5x 2— downstream face. v 1 av %y = velocity in y-direction on up- and 5y downstream face. The hydrostatic equation is the basic assumption of the first-order shallow water theory and is prevalent in engineering applications (Liggett, 1975). The hydrostatic equation implies that the flow lines possess no curvature. Other assumptions implied or not listed include: 1. There is no Coriollis acceleration. 2. Streamlines are not curved such that the pressure variation with depth is linear. 3. Turbulent velocities are time-averaged. Turbulent fluctuation velocities are not considered since on time average the‘ net effect is present as represented in the conservation of mass and momentum equations by the Reynolds stress 'UTV' (Potter and. Foss, 1982). In the one-dimensional case, velocities are depth— averaged, thus suppressing the vertical dimension in the conservation equations. U is the primary flow velocity and 50 is a function of x and t. The x- direction is such that cose a 1. The control volume and vector convention in Figure 3, is used in the Reynolds Transport equations. control surfa ce control v olu me Figure 3 Control Volume Definition The Reynolds Transport theorem provides a method of describing, among other thermodynamic quantities, the transport of mass and momentum into and out of a control volume (Potter and Foss, 1982). 51 1. Conservation of Mass‘ The Reynolds Transport equation for conservation of mass requires that a I an + I (V°n)dA = o [24] 5t cv cs where V = velocity vector, dA = differential element of control volume, cv, v-n = dot product of the velocity vector to the normal unit vector of the control surface, dA = differential element of the control surface, cs. The mass balance for a differential volume of height h(x,t) is (u - au Ax)(h -3h Ax) 5? -2 5i —2 - (u + Bu Ax)(h + 3h Ax) 8 8h Ax [25] 6i ‘2 3i .2 3? Expanding, summing like terms, and disregarding higher order terms such as an ah 5? 3i yields the discrete form - 3h Ax — hau Ax = ah Ax , [26] 3? 5? 3? which is useful in finite difference computational methods. 52 The differential form of the conservation of mass equation is 8h + 3(uh) = 0 [27] 5? fix for one dimensional, depth-averaged, time-averaged, flow. 2. Conservation of Momentum The conservation of momentum equation is derived in a similar manner to the conservation of mass. The Reynolds Transport equation for momentum requires that A a I (ov) dA + I (pV)(V-n)dA = 2F [28] 5t cv cs where V = velocity vector, dA = control volume, 0 = mass density of water, A v-n = dot product of the velocity vector to the normal unit vector of the control surface, 2F a summation of forces acting on control volume. Substitution of' the velocity components for the depth- averaged elemental control volume results in 3(pUh)Ax - [U(uh) - 3(U(uh)) 9;] + (U(uh) + 5t 3x 3(U(uh)) Ax ) = ZFx [29] 5x —_2 Upper case variables are vectors and lower case variables 53 l are scalar values. The term XFx is the vector summation of forces in the x-direction. These forces are l. Gravitational, F = pghAxsinex = pgthAx where ex = angle getween the horizontal axis and the x-axis. Sinex = tanex for small slope; Sx = tanex = sine. h h 2. Pressure, Fp = I pdz = pg f (h-z)dz = pghz/Z 0 0 3. Frictional resistance, Fs = pthfo where Sf is the frictional slope defined by either of Manning equations, u = 1.486R2/3Sf1/2 (english units) , [30] n u = l R2/3Sfl/2 (SI units) , [31] n or the Chezy equation, u = C(RSf)l/2 (SI or english units). [32] Combining these equations using the Manning equation results in p(U(uh) - 30(uh)Ax] - p(U(uh) + 3U(uh)Ax) 5x _2 _ 5x .2 + pghS Ax + 1/2pg(h2 - 3(hz)Ax - (h2 + 8(h2)Ax)] x ax ‘2 ax 7 + pgthAx = %(pUh)Ax [33] t Rearranging and disregarding higher order terms and taking the limit Ax + 0 yields 54 t hau + uah + ua(uh) + uhau + ghah = gh(Sx- Sf) [34] 8? 5? 5x 5} 5? By substitution of the differential form of the conservation of mass, 8h + 3(uh) = 0 [27] 8? 5x into the above equation [34] the conservative form of the one-dimensional, depth-averaged, time-averaged, conservation of momentum equation is obtained a0 + U8u + gah = g(Sx - Sf) [35] 5? I? 5? The two-dimensional case is derived similarly except that momentum in the x-direction may now be transported into and out of the differential control volume on both the x- and y-faces. The general form of the Reynolds Transport equation in the case of conservation of mass is A a I dA + I (v-n)dA = o [24] 5t cv cs where V = velocity vector, dA = differential element of control volume, cv, Von = dot product of the velocity vector to the normal unit vector of the control surface, dA = differential element of the control surface, cs. Substitution of the values for the two dimensional control 55 volume results in, a(ph) - [Uh - 8(uh)Ax] + [Uh + 8(uh)Ax] 5t 3x ‘2 5x _2 - (Vh- 3(vh)A J + (Vh + 3(vh)A ) = 0 [36] 5y ’2 5y .2 Combining terms yields 3h + 3(uh) + 8(vh) = 0 [37] 8? 8x 5y The Reynolds transport equation applied to the conservation of momentum yields A a I (pv) dA + I (pV)(V-n)dA = 2F [28] 5t cv cs where V = velocity vector, dA = control volume, p A = mass density of water, v-n = dot product of the velocity vector to the normal unit vector of the control surface, IF '= summation of forces acting on control volume. Substitution of the velocity components for the depth- averaged elemental control volume produces 8(Vh) - [v - av A )(uh) + (v + av A )(uh) at 37.12, 15722 - (U(uh) - 8(U(uh))Ax] + [U(uh) + 3(U(uh))Ax) 5x ‘2 5x -2 [ 1 38 56 i Simplifying, the x-momentum equation results in 3(uh) + 8(uzh) + 3(uvh) + a(h2) = gh(S -sf) [39] 5t 5x 3y §§x x and the y-momentum equation becomes gIvh) + 3(uvh) + 3(v2h) + thZ) = gh(Sy - Sf) [40] t 5x 5y Y Incorporating the continuity equation, au +‘av = 0 , [41] 3337 yields the x-momentum conservative form: an + uau + v8u+ gah = g(Sx — Sfx) . [42] 5? 5? a; 3? The y-momentum conservative form is derived similarly as 8v + uav + vav + 98h = 9(S - Sf ) . [43] a? a7 a7 a? Y y Assumptions implied in the above derivation include. 1. The above derivations do not include lateral inflow or other inflow/outflows such as rainfall or infiltration. 2. The channel is prismatic and of a rectangular cross-section. 3. Irregular Cross-section In deriving the above conservation of mass and momentum equations, the cross-section of the channel was assumed to be rectangular. However, in the application to 57 1 natural channels, a more general form is needed' in which the depth h is assumed to vary across the y-direction. The irregular cross-section is shown in Figure 4. Chmme Wafer Surface Bottom \:\/e\:/ "I; Figure 4 Irregular Cross-section of a Channel. In the derivation of the conservation of mass and momentum the depth h times the unit or differential width of the control volume was used as the dA in the Reynolds Transport equation. The Reynolds Transport equation for conservation of mass requires that A a I fill + I (V0n)dA=0 [24] 5t cv cs where V = velocity vector, dx = differential element of control volume cv, v-n - dot product of the velocity vector to the normal unit vector of the control surface, dA - differential element of the control surface cs. If instead of the depth h times the unit or differential width in the y-direction, the area A, which is a function of f(y,z), is replaced, we obtain 58 3A + a(ua) = o ' [44] 5? 5x and 3(UA) + 3(02) + F = gA(S - s ) [45] 5t 3x p x f The pressure term Fp is the force on the face of the control volume h Fp = g g(h-z)£(z)dz [46] in which C(z) is the channel width at the height 2 above the bottom of the channel. The net force in the down- stream direction is h rp - (Fp + -%:pr) = -a g pg(h - z) €(z)dz Ax 3x [47] By simplifying the right of [47] using the Leibnitz rule we have, b - a I og(h - z) £(z)dz = x 0 h -pg I a [ (h - z)€(z)]dz 0 x h h = —pg[ah I £(z)dz + I (h - 2) 35(2) dz] 5x 0 0 5x 59 The first integral on the right defines the cross- sectional area. The second term a§(z)/8x is assumed to be zero, which is the prismatic channel assumption. If the channel narrows or widens, then an additional force is exerted on the channel walls. Therefore, the pressure term in the conservation of momentum equation is interpreted as F = 9A an [48] p a7 and the conservation of momentum equation 3(UA) + a(uz) + F = gA(S - s ) [49] 5t 5x p x f becomes a(ua) + 3(02) + gA ah = gA(Sx - Sf) [50] 6t 3x 3? Using this conservation of mass equation for an irregular cross-section [50] the conservation of momentum equation becomes an + uau + gah = g(sx - Sf) , [51] 5? 5? I? which is the same as the conservation of momentum for one dimension. It was assumed that the channel width variation is negligible so that 3€(z)/ax is zero. 4. Lateral Inflow Lateral inflow must be incorporated into the conservation equations, if for example, distributed inflow 60 i from groundwater, infiltration, tributary inflow, rainfall, or overland flow occurs. In this derivation inflow is regarded as positive and outflow as negative. Let the symbol q represent lateral inflow/outflow with dimensions of [L’]/[TLz] when used in the one-dimensional equations or [L’]/[TL] when used in the two-dimensional equations. The conservation of mass equation incorporating lateral inflow would be ah + 3(Uh) = q [52] 6? 5x for the one-dimensional case and similarly for the two- dimensional case for rectangular or irregular channel cross-sections. The conservation of momentum equation must account for the momentum entering and leaving the control volume transported by the lateral inflow or outflow. The additional momentum entering is pgqux where uq is the downstream component of velocity of the lateral inflow. Upon leaving the control volume the velocity of the lateral inflow is assumed to have the same velocity as the downstream velocity of the fluid, so that the momentum of the lateral inflow leaving is quAx: haU + Uah + U3(uh) + uhgu + ghah = 5? 5? 5x 3x 6? gh(Sx - Sf) + unq - U) [53] for the two dimensional conservation of momentum 61 equation and 3(UA) + 3(02) + gA ab = ngx - 5f) + 5t 3x 5? for the two dimensional conservation of momentum equation for irregular cross-section. C. Finite Element Model Formulation The form of the partial differential equations governing. direct runoff have been derived by use of the Reynolds Transport theorem. These partial differential equations are used in the finite element method, Galerkin formulation of the surface water equations. The equation of continuity for an incompressible fluid is written as 8n + 3v + 3w = 0 [55] 5? a? 32 where u = u(x,y,z,t), velocity in the x-direction, v = v(x,y,z,t), velocity in the y-direction, w = v(x,y,z,t), velocity in the z—direction. On integration in the z—direction and using appropriate boundary conditions, we obtain: an + aIEh) + a(Vh) = r - i [56] 3? 5x 5y 62 where depth-averaged velocities, C < ll :3" II vertical depth of flow, rainfall intensity and infiltration rate respectively. '1 He II The equation of momentum, the Navier-Stokes equation for two-dimensional flow with appropriate boundary conditions, is ail + uaii' + vafi + gah = g(S - Sf ) - U(r-i) [57] a? a; a; a7 °" " a and g; + “g; + v%% + 9%; = g(Soy - Sfy) - g(r-i) [58] where Sox and Soy are the slopes of the element in the x- and y-direction respectively. Sfx and Sfy are the frictional slopes in the x- and y-direction respectively (Taylor, 1974). The continuity and momentum equations form the governing equations for watershed surface hydrology. Woolhiser and Liggett (1967) and other researchers have documented the adequacy of using the kinematic wave simplifying assumption, where all terms on the left—hand side of the momentum equation are assumed negligible. The general form then becomes V = m h“ [59] 63 where m and a are constants such as in the Chezy or the Manning equation, and v is velocity in the direction of flow. If the above assumption is made, then only the continuity and kinematic wave formulas need to be solved over each element. This in effect reduces a nonlinear set of simultaneous, partial differential equations to a linear set. The validity of this simplifying assumption is considered only if Lso > 10 [60] Fozho where L is the length of domain and F0 and ho are [the Froude number and depth of flow at the downstream location end under steady-state conditions. This restriction is dependent on the finite element grid and hydraulic parameters used to describe the watershed. The finite element formulation is applied by writing approximating functions of the form n 9e = 2 Ni(x.y)01 [61] i=1 where 9i are known values of the function 4 at each of the n nodal points. Ni(x,y) is a shape function, which approximates the function ¢e based on the n nodal values (Segerlind, 1984). 64 t The watershed surface flow variables are approximated as follows: n ue a Z Ni(x,y)fii [62] i=1 n — ve = 2 Ni(x,y)vi [63] i=1 n he = 2 Ni(x,y)hi [64] 181 where ue, Va and he are the approximate values of the velocities in the x- and y-directions, and the depth, respectively, within the finite element domain. The assembly of the elemental equations into a global matrix form for the entire domain or watershed constitutes the system of equations that models overland flow runoff. Iterative solution procedures are required (such as the Newton-Raphson technique) to obtain the solution of these equations in the time domain. 1. Element Equations . The application of the Galerkin method to the kinematic equation for overland flow is performed as follows: The convention of {] representing a vector quantity and [] representing a matrix quantity will be used. For two dimensions: 65 IIIan( £3 + a(Bh) + a(6h) - (r - i) )dxdy = o [65] t 5x 5y For one dimension: ItanI ah + 3(6h) - (r - i) )dx = o [66] 5? §§_—_ Applying the kinematic wave assumption, the momentum equation is reduced to So a Sf [67] Utilizing the Manning equation that relates the depth of flow to discharge for turbulent flow, we have n For a wide channel or overland flow, the hydraulic radius is R = A/P 2 A, so that the discharge relation becomes n Recognizing that the cross-sectional area A is, in the case of overland flow, equivalent to the flow depth defined above as h, the wide channel assumption of R = h, and the vector Q in [69] is resolved into its respective directional components in the x and y direction, with 6 as the angle with the x-axis results in Uh = Qx = Q cose [70] 66 S - I [NJT dxir-i} [72] and for two-dimensions, {R}(9) = I [N]T[N]{A]dx + I[N]T[3Ni/3x 8Nj/3X][Qx} + I[N]T[3Ni/ay aNj/aylioyl- I [n1T dx{r-i} . [73] It should be observed that the discharge values 1/2 QBX = 1.486 R2/3 SB AB [74] “8 and 1/2 QBY = 1.486 32/3 58 AB [75] “B are the B nodal values. The concept of the nodal discharge values based on nodal values of slope in the x- and y- directions and the nodal values of roughness n3 is of paramount importance if the effects of kinematic shocks are to be avoided between elements where abrupt changes in slope or roughness would otherwise occur. By causing the nodal values of discharge to be computed with nodal values of slope and roughness, a linear variation over the element of these values effects a linear variation in Q without abrupt discontinuities at the inter-element nodes. The numeric difficulty encountered by other researchers is thus avoided (personal correspondence with D.A. Woolhiser and G.A. Blandford). This is due to the elimination of the diffusion term ah/ax along with the other terms in the 67 over the element of these values effects a linear variation in Q without abrupt discontinuities at the inter-element nodes. The numeric difficulty encountered by other researchers is thus avoided (personal correspondence with D.A. Woolhiser and G.A. Blandford). This is due to the elimination of the diffusion term ahlax along with the other terms in the full dynamic equation because of the kinematic assumption. Because of this, there is no possibility other than a discontinuity in flow depth at the inter-element nodes. To avoid this difficulty, the values that relate h to Q--i.e., the slope and roughness-- must be considered as nodal values. The value of Q is a vector quantity in the direction of nodal slope subsequently resolved into x- and y-direction components according to the orientation of the slope with the global coordinate system. NODAL VALUES 01 02 Q3 Q4 81 s2 S3 84 n1 n2 n3 n4 (1) (2) (3) 2 4 ELEMENT NUMBERS = (1), (2), etc. NODE NUMBERS = 1, 2, etc. Figure 5 Element and Node Numbering Convention The convention illustrated in Figure 5 is the nodal 68 representation. by nodal values of the independent variables. With this methodology, the discontinuities in function values are avoided at the inter-element nodes. 2. Shape Functions The shape functions provide the basis of writing the linear variation across the element of the approximated functions. The local coordinate system allows easier integration over an element. For this reason the following system is defined 5 - x-direction with limits of 0 i s 5 L for a one-dimensional linear element of length L. The shape functions are N1 = 1 - s/L [76] Nj a s/L [77] or, in matrix form [N]T = [l-s/L] [78] s/L The partial derivative is computed by taking the partial derivative of the shape functions, since the approximating function is ¢(x) = N101 + NjOj [79] where the nodal values 0 are constants with respect to the x- or s-space dimension. The derivative is 69 i x- or s-space dimension. The derivative is a¢(x) = au-a- + an 4- [80] x 1(5.321 1 5i] 1 The derivatives of the shape functions with respect to x or s--which are equivalent locally within the element--are 3N1 = :1 [81] I? L aNj = l [82] I? L In matrix form this is r n ._ 3N1 F—l [beT = a} = l [83] aNj L l 3? - 1 - J Writing the individual integrals from the residual {a}(e) = I [N]T[N]{A}dx + I[N]T[3Ni/3x 3Nj/3X]{Q} - I [N]T dx{r-i} [84] and integrating over the element length L we have I [ulTIN]{A}dx = I [1 -s/L] [l-s/L s/L]{A} dx s/L = L 2 1 {A} :[1 2] = [cliil [86] 70 8 The discharge term in the residual is L IINlTIbMQ} =I [1 4/1.] [_-_1 I] {Q} 0 s/L L L = 1 {-1 1] {Q} 2 -1 1 = [bx]{Q} [87] The lateral inflow term is L L I [an dx(r-i) = I [1 -s/L] dx(r-i) 0 0 s/L = (r-i)L l [88] 2 1 Assembling the results of [86], [87], and [88] into the residual expression, we have {a}(e) = g [i i] {A} + % [:i i] {Q} - (r-i)L 1 [89] The residual {R}(E) in [89] is minimized over the system of elements only when assembled in the global form. 3. Global Matrix The global matrix form of the assembled elemental residuals may be formed by the direct stiffness procedure when local node numbering schemes are used. An expanded 71 8 form will also be examined, which uses a global nodal numbering system. The direct stiffness procedure results in a banded matrix as follows. Given a three-element, four node system as shown in Figure 5, it can be seen that node 2 receives contributions from the elements (1) and (2); node 3 from elements (2) and (3); and nodes 1 and 4 from elements (1) and (4) respectively. Assemblage by the direct stiffness procedure is as follows for elements of equal length: 2 l .1 A1 -1 l I Q1 L 1 (2+2) 1 A2 + 1 ‘1 (1‘1) 1 Q2 6 1 (2+2) 1 A3 2 '1 (1-1) 1 Q3 1 2 A“ -1 1 Qt. L - - - l rL 2 [90] _2 2 1 The nodal values of Q are related to A by the Chezy or Manning equations, resulting in a nonlinear set of differential equations with respect to time. Alternately, with [C], [K], and {F} as defined above in [86-88] without the (e) elemental designation we have {chi} + {who} = m [911 72 This matrix equation represents the system of equations to be solved. Note that .the left—hand side contains the time derivatives of the cross-sectional area {A}. A solution in time is needed such as a finite difference scheme. Note also that the discharge values of Q contain the cross-sectional area A in the Manning equation. Therefore, at each time step, when A is solved, a new set of Q's must be computed. Depending on the time weighting coefficient, an implicit or explicit solution for A and Q results. The rL values are considered constant over each element. This results in the form of a forcing function vector on the right hand side. The finite difference solution utilizing the form of equations [14-19] for the time dependent matrix equation [91] is [CllAlnew = [C]{A}01d - At[bx]((l-O){Q]01d + e{Q}new] + At[(1‘9){F}old + e{1"}new] [92] The time dependent finite difference form of equation [91] can be recast into the nonlinear form [cllalnew = {1“} [931 where {F*} is the combination of the right-hand-side terms in [93]. The right-hand-side contains terms such as {anew that are functions of the left-hand-side term {AInew- This forms a set of nonlinear equations requiring special solution techniques such as a simple iteration scheme or 73 The iterative schemes are preferred for large systems because, unlike the Newton-Raphson Method, no matrix inverses are required. At each time step the value of the rainfall excess intensities are placed in {F*}, the old and new values of {A} are assumed equal. At the first time step the old values are the initial values. Then the system of equations are solved by standard methods. The new values thus solved for become the new values for the next recursion until the solution converges to within a set tolerance. This recursion is repeated at the next time step with the new values from before becoming the old values for the present time step. 4. Expanded Form The expanded form of the global system of equations is as follows for the system of elements depicted in Figure 5. The approximating shape functions are written for each element and assembled into the expanded elemental equations: ¢(1) 3 N1(1)01 + N2(1)02 + 003 + 004 4(2) . 041 + N2(2)02 + N3(2)03 + 004 4‘3) . 001 + 042 + N3(3)O3 + N3(3)04 or, in matrix notation, “NH vvv 0 N2 N3 0 02 [94] [N1 N2 0 0 01 0 0 N3 N4 G3 MIMI-Ow 00-O- news-w 74 It is now possible to combine the elemental equations into a single region equation (Segerlind, 1976): a a . z ¢(9) [95] e-l 'rhis results in the equation ¢ = [N1(l)]01 + [N2(1) + N2(2)]¢2 + [N3(2) + N3(3)]O3 + [N3(3)]04 [96] The importance of this region equation is seen when an element region in two dimensions is formed from simpler finite elements. For example, a region as follows may form an element region, ' D2,X2,Y2 (1) ( ) Dlrxerl 2 ---11-_‘~ D3,X3,Y3 .; (5) “" D4,X4,Y4 (3) (4) ' ‘ Ds.X5:Y5 D6IX6IY6 Figure 6 Five-Element Region Writing the set of element equations, we have 75 \ {¢}(2; an N2 N3 0 0 0 01 {¢](3) 0 N2 N3 N4 0 0 T2 [4}(4) = 0 0 N3 N4 N5 0 0 43 [98] {¢}(5) 0 0 N3 0 N5 N5 04 {4} N1 0 N3 0 0 N5 05 l. - 06 This regional equation could be considered as an expanded element defined by triangular elements and area coordinate shape functions. The utility of such an approach is in the modeling of irregular patches as defined by the hydrologically homogeneous areas. By building an expanded region the irregular patch can be handled. Substructuring removes the internal node if no nodal value is desired at that point. This is done after the global matrix is constructed. By solving the matrix for removal of a nodal value and associated constants, the resulting set reflects its contribution but is not present in the set of equations. 6. Isoparametric Elements Another technique for representing finite element regions is one using isoparametric elements. This class of elements utilizes a coordinate transformation technique that maps the element in the global coordinate system into a natural coordinate system. The natural coordinate system for a linear, one-dimensional element is the i system cor- responding to the global x—coordinate system. The 5 system varies from -1 to +1 with the origin at the center of the element. The element matrices such as [C], [B], and {F} are integrated numerically in the 5 coordinate system. 76 . The coordinate transformation is written in terms of the same element shape functions as used to represent the state variables, hence, the name isoparametric. Super- and sub- parametric elements are those that use higher or lower order shape functions to perform the transformation. The transformation is done by writing the global coordinate system in terms of the nodal coordinates and the shape functions: x = N1(€)Xl + N2(€)Xz [99] where 1/2(l-€) and 2 II N 1/2(1+§) . 2 The change of variable in any integral is accomplished by writing X‘ +1 I i(x)dx = I 9(a) I d(x(§)) )dt . [100] -1 '—_af' xi The Jacobian of the transformation is d(x(€)) = — xi + xj = L [101] “3€"“ 2 2 2 Jacobian transformations are done similarly for two- dimensional elements. The four-node quadrilateral has four Shape functions in g and n. The coordinate transformation is accomplished by writing the global coordinates in terms 77 of 5 and n such that, and where The J The N ll *< II N I. acobian [J] Jacobian [J] N1(£.n)xl + Nz(5,n)xz + N3(£.n)x3 + N~(£.n)x~ [102] N1(€.n)rl . N2(€.0)Y2 + N3(€:U)Y3 + N~(5.n)Y~ = 1/4(1—5)(1-n) = 1/4(1+5)(l-n) 1/4(1+5)(l+n) 1/4(1-£)(1+n) of the transformation is 3N 3N TE’ 35‘ it: as, an an L for the linear triangle iii if if if. L z 2- 35 ag‘ 35 an3 (x - x ) (Y - Y ) I 3 I [103] aN IX Y ' g” 1 1 x Y 8N 2 2 n“ X Y .1 3 3 x Y - to 4‘ [104] finite element is 3 (X - X ) (Y - Y ) 2 3 2 3 [105] 78 where , the area coordinate and shape function 1, A A L = A , the area coordinate and shape A function 2, 1 L = - L - L , the area coordinate and shape 3 1 2 function 3, A, A1 , A2 = total area, and partial areas of local coordinate. The isoparametric, linear and quadratic quadrilateral are useful for representing odd shaped and curved boundaries. The integration of the element matrices are performed on the transformed element in the natural coordinate system. The integration method most commonly used for integrating functions is the Gauss—Legendre quadrature. This method replaces the integral with a summation of the function at mxn integration points multiplied by mxn weighting coefficients. This represented by the following equation [106] and is the method used to integrate the shape functions as used in [86-88]. 1 l m n' I I 9(€.n)d€dn = 2 Z Win g(£i.nj) [106] -1 -1 i=1 j=l The number of integration points depends on the highest power of the function to be integrated. This means that a polynomial of power (Zn-l) may be integrated exactly with n sampling points at which the element is evaluated 79 ‘ (Segerlind, 1983). Figure 7 depicts the isoparametric linear, one-dimensional and two-dimensional triangle and quadrilateral elements together with the integration points and weighting coefficients for the highest polynomial power of 52 and n2 . 5 IF -1 +1 L=l=]/3 w=|/2 o '1 O f f: 17: 2 0.577350 w t 1.0 O 0 Figure 7 Isoparametric Elements and Integration Points. 81 7. Consistent Stress Approach Nodal values of discharge Q, lepe S, roughness n, and rainfall r can be considered either as constants over an element or as nodal values. It is necessary to maintain continuity at the inter-elemental nodes. The difficulty arises whenever the state variable is either a) known at the nodes and it is required to be used as an elemental constant, or b) known as an elemental constant and it is required to be used as a nodal value. The consistent stress approach provides the correct method of computing these values. The form is [C]{¢] = r{F} [107] The vector {4} contains the nodal values, whereas r represents the constant values over the element (Segerlind, 1976). In summary, the governing equations have been developed for use in modeling the overland flow from an infiltrating watershed using a Galerkin formulation. The governing equations may be solved for the steady-state and the time-dependent cases. Computation may be done by hand for a small number of elements or by a computer program for larger systems. Before initiation of the computation, a check should be made on the applicability of the kinematic 82 wave approximation to the specific set of variables or problem. This may be done by using equation [52] to compute the kinematic number k, which must be larger than 10 for the maximum intensity. In addition, the Froude number should be computed to determine where the boundary condition should be applied--upstream or downstream--for sub- and supercritical flow, respectively. The time solution requires the selection of a time step At. This time step must not be larger than the time during which a gravity wave can travel over the length of the element. This is the Courant condition and should not be violated. Provided these conditions are met, computation of the time-dependent solution of the overland flow equations may be performed using the Galerkin finite element formulation. The method is applicable to both one-dimensional and two— dimensional finite elements using the assembled form of equation [91] and solved by any of the standard finite difference, time-domain solution techniques represented by equations [14] through [19] as shown in [92] and [93]. IV. RESULTS AND ANALYSIS The results obtained and analysis performed are presented in the following order: The finite element formulation, Geographic Information System analysis, input parameters, hydrologic response areas. The watershed modeled was Watershed 4H near Hastings, Nebraska. The rainstorm event selected was 0.33 cm on May 4, 1959. A. Finite Element Formulation The Galerkin finite element formulation was used to solve the kinematic wave equation. The Green-Ampt infiltration equation was used to compute the rainfall excess. The excess rainfall intensity was calculated using the Green and Ampt runoff model from the USDA-ARS Water Erosion Prediction Program (WEPP) project (under development). The Galerkin formulation of the kinematic wave equation has been presented in the theoretical development and will not be repeated here. The elements used for the arbitrary finite element grid were linear, one-dimensional elements. The formulation used variable width to represent trapezoidal areas. The elements used in the hydrologic response area finite element grid were the linear, two- dimensional, isoparametric four-node quadrilateral elements 83 84 4 and three-node triangles. The system of equations solved by time integration for the one—dimensional, arbitrary grid was [a]{h} + [b]{Q} = [aliiel [108] where [a] T InlN} ([N][w])[N]dn , [b] T I {N}d{N]/dx do , a {ti} i = rainfall excess, {dh/dt], Q flow rates, h = flow depth, w a flow width of the element. The system of equations solved for the two-dimensional hydrologic response area grid was [a]{h} + [bx]{Q} + [by]{Q} = i{F} [109] where T [a] = IniNliNldQ . and 85 S T [bx] = I {NldiNl/dx d9 . n T [by] = IniNldiNI/dy d9 , T {F} = I {N} do , n {a} = {dh/dt}. rainfall excess, Ho ll Q = flow rates, flow depth. Since [109] is a two-dimensional formulation and the equations are integrated over the two-dimensional element domain a , no variable width w is needed in the system. The nodal flow rates are for a unit width. The system of equations in [108] and [109] may be solved with any of the common time integration schemes-- Forward Difference Central Difference Galerkins Backward Difference depending on the time weighting coefficient chosen. Stability and accuracy of the solution will dictate the proper time weighting coefficient. The Central Difference scheme was used in this analysis. This corresponds to a value of 0 = 1/2. No numeric oscillations were observed, which may result from choices of time integration schemes 86 i that are inconsistent with the finite element grid. A criterion that must be met is the kinematic number after Woolhiser and Liggett (1967). This comes from the kinematic approximation to the full momentum equation. The approximation is accurate to within 10 % if k = LSog/V2 > 10. The steps used to check this criterion are to calculate the equilibrium outflow, which is the maximum rainfall excess, solve for the corresponding flow depth using the Manning equation, solve for velocity, and check the Froude number and the kinematic number. The actual time step used in the time integration scheme must not be longer than that time during which a gravity wave front may propagate through the system. This is known as the Courant condition: At < Ax/c . where At = time step, Ax = distance increment or element length, c = speed of a gravity wave or (5/3)V. If this condition is violated the partial differential equation theory is violated. This condition arises from the theory of the Method of Characteristics. A time step <3f one minute was selected as shown in the Appendix. 87 S B. Geographic Information System Analysis The spatial analysis of the ARS Watershed 4H near Hastings, Nebraska was performed using the ARC-INFO GIS, version 3.2, developed by Environmental Systems Research Institute, Inc., Redlands, California. Several layers were digitized in order to perform an overlay analysis of the landuse, soils, and slope data. This overlay effects a delineation of areas containing homogeneous parameters. In this case, the hydrologic response areas of homogeneous parameters consisted of only the slope interior to the watershed boundary. The landuse was digitized for the watershed according to the landuse at the time of the May 4, 1959 rainstorm event. This watershed, like many of the ARS research watersheds, was under a single crop and tillage practice. The landuse for this watershed was fallow under good residue cover at the time of the rainstorm event. Figure 8 depicts the landuse at the time of the event on May 4, 1959. The soils in this watershed are the Hastings silty loam, Hastings silty clay loam, and a Colby silt loam. These soil names are the 1939 soil survey names originally mapped for this watershed by the USDA-Soil Conservation Service. Figure 9 shows the three soil delineations for the watershed. The Hastings silty clay loam is an eroded profile of the Hastings soil series. The B-horizon properties were used in the Green and Ampt modeling for 88 this mapping unit. The Green and Ampt equation was solved for this storm and these soils using the soil properties from soil samples sampled April 21, 1965 by the USDA-Soil Conservation Service. The location of the Hastings soil sample was in Webster County, 0.15 miles west and 180 feet south of northeast corner of Sec. 1, T3N, R10W. The soil properties were taken from the Hastings Soil No. S65NE—9l-l data, Sample Nos: 20449-20456. The Colby soil prOperties were taken from the Soils-5 data sheet for the series. This information was contained in the soils data base at the USDA-Soil Conservation Service, National Soil Survey Laboratory, Lincoln, Nebraska. The soil properties are contained in Table 1. TABLE 1 Soil PrOperties 1939 SOIL S > S > S S BULK PERCENT SURVEY ' DENSITY ORGANIC CEC/ NAME 3 in 010 SAND CLAY GH/CC MATTER CLAY COLBY 81 L 0 0 10 21 1.40 1.50 0.65 HASTINGS S! L 0 0 10 24 1.27 2.50 0.79 HASTINGS 81 C L 0 0 8 34 1.26 1.77 0.67 89' LAND USE How Figure 8 Watershed 4H Landuse. 90 30: .sz mflfiom me conmuuumz m musmwm :3 an huhh hfikhfih‘fih as hflhh hfikhfih‘hfi . as huhm. Nan-.0 - MQNOW 91 The slope and topography was obtained from two-foot contour maps provided by the USDA-Soil Conservation Service in a 1942 plane table survey. Figure 10 shows the digitized elevation map derived from the 'plane table survey. The ARC-INFO GIS, utilizing the Triangular Irregular Network (TIN) subroutine, calculates the slope of each triangular area from the elevations of the vertices. The TIN slope map is shown in Figure 11. The slopes required for the finite element modeling are at the nodes of the elements. To achieve this, the finite element grid was overlaid on to the slope map an the slope at the node. tabulated in the INFO Database for each node. As the finite elements decrease in size, the elemental slopes tend toward nodal slopes in the limit. The slope map in Figure 11 shows increased complexity primarily due to the delineation of not only slope but aspect. The aspect of the slope is the direction measured in degrees between the principal slope and the north direction. The aspect and principal slope are resolved into x- and y-direction slopes for use in the finite element model. In the kinematic wave equation the friction relation, Manning or Chezy, contains the square root of the slope. When modeling the two- dimensional flow equations, the flow is caluclated in the direction of and using the principal slope. This flow rate is then resolved into x- and y-direction flow rates for use in calculating the 8(uh)/3x and 3(vh)/8y terms. To relate the orientation of all vector quantities to the watershed 92 ‘ map, such as velocity or flow rate, the aspect anlge of the principal slope was used. ELEUHT ION "‘ 197B R \ ‘\ ‘ ‘ ~ I“ \\ \ ‘ \ \\‘ \ ‘\ \‘\\\ \_\ \ \\ \ \\ \\ i \ \\ \ \ \ \ ‘. I ) ’) ‘ \ \u §T 3 \\ \ \ 'l W 2 /‘ I / I) K \ ‘/ /// fl / / / I / (V /" 1’, // / / // / / r” / r ( \ ,» L I L __/ thaw Figure 10 Watershed 4H Elevation Map. 93 Interval) X (1 SLOPE SLOPE % 7kfinvfirl 11188 _ _ _ _ _ PO71QZSQU 411.118 mmmmm 183456 811.1111: 1ACQYQEXb7AEOYI_ __ _ __ __ __ __ _ ._ _Q7I?:dagb 6183456789111111 III-EEEEEEEIl-II Figure 11 ARC-INFO TIN Slope Map. 94 Figure 12 depicts the rainfall excesses for the three soils as defined by the Green and Ampt equation--no significant difference was observed. For the purposes of computing runoff and infiltration, these soils may be treated as one for the entire watershed. While spatial nonhomogeneity may exist it is not detected by the soil mapping units interior to the watershed when considering infiltration. Figure 13 illustrates the first case analyzed, a finite element grid of arbitrary spatial form. That is, a set of linear elements of variable width. This is the same finite element grid, slopes and Manning n's as used by Peters, Blandford and Meadows (1983). Table 2 contains the arbitrary finite element grid input data. 95 Table 2 Arbitrary Grid Finite Element Input Data. ELEMENT NODE X- NODAL MANNING WIDTH NUMBER NUMBER COORD SLOPE n (PBS?) (1) 1 0 0.0406 0.035 353 2 171 0.0406 0.035 (2) 2 171 0.0406 0.035 333 3 342 0.0669 0.035 (3) 4 342 0.0669 0.035 206 4 507 0.0562 0.035 GREEN and AMPT RAINFALL EXCESS HAY 4. 1909 96 \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\‘ E/ //////// ////// \\\\\\\\\\\\\\\\\\\\\\Lfi \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ (III/III T 14.3? fill! (HOURS) HASTINGS 81 C L 10 00 12.00 - 11.00 - 9.00 -— 0.00 - 7.00 - 0.00 - 0.00 - 4.00 - 0.00 ~ 2.00 - 1.00 -« 0 00 10.00 -« (an/no) ssaaxa mam IIIIIIIII - 14.12 14.43 m 14.40 14.30 14.30 COLBY ZZ manna: s1 1. Figure 12 Rainfall Excess for Three Soils. 97 12 "GAGING STATION flovv O 50 #100 200 feet Figure 13 Finite Element Grid of Arbitrary Spatial Form. 98 1 The topography, in this particular watershed, is the one parameter of known spatial nonhomogeneity. Slopes defined by the 10-foot contour intervals, were used to form the finite element grid. Nodal points were located at the intersections of the 10-foot contour lines and the watershed boundary and the stream channel. The slope on a two-dimensional finite element grid possesses both x- and y-direction slopes. Even though two nodes lie on a contour, the element has an x- and/or y-direction slope due to the skewed spatial form with respect to the coordinate system. If the slope is spatially nonhomogeneous within the element then it should be averaged over the element to obtain the slope for the element. The scale of the variation defines the scale of the finite element grid. The 10-foot contour intervals were used in order to limit the number of elements representing the watershed. This was done in order that the -comparison of the arbitrary-grid to the hydrologic- response-area grid would not be obscured by the increased _accuracy gained from a large increase in the number of elements. In order to model the hydrologic response areas, two- dimensional elements are required. The elements selected for this research were the three—node triangle and the four-node quadrilateral, isoparametric elements. A computer program was written that performs the coordinate transformation and computes the partial derivatives and 99 1 other. components of the two-dimensional flow equation [109]. The nodal coordinates were obtained from the GIS overlay of the finite element grid over the watershed. The watershed is represented in the state plane coordinate system, defined by the plane table survey of 1942 of the watershed available from the USDA-ARS-Water Data Laboratory, Beltsville, Maryland. These coordinates are in Table 3 listed by node number. Each element is represented by the nodes associated with it. Thus, for Element (8) in Figure 14, the node numbers 1, 9, and 10 represent the three-node triangular element. 100 Table Hydrologic Response Area Nodal Coordinates and Slopes NODE X- Y- X- Y- . NUMBERSCOORD.COORD.ELEVATION SLOPE ASPECT SLOPE SLOPE 1 1618 1302 1946.00 0.000 0.00 0.000 0.000 2 1610 1345 1950.00 0.081 -177.30 0.004 0.081 3 1707 1467 1960.00 0.049 -119.50 0.043 0.024 4 1897 1529 1970.00 0.063 -110.90 0.059 0.022 5 2079 1538 1978.00 0.032 -76.30 0.031 0.008 6 2111 1309 1976.30 0.022 -127.00 0.018 0.013 7 1977 1309 1970.00 0.064 -113.60 0.059 0.026 8 1841 1313 1960.00 0.159 -159.90 0.055 0.149 9 1710 1320 1950.00 0.086 -75.70 0.083 0.021 10 1625 1274 1950.00 0.118 -29.50 0.058 0.103 11 1739 1188 1960.00 0.089 -75.90 0.086 0.022 12 1931 1161 1970.00 0.033 -91.80 0.033 0.001 13 2077 1160 1975.00 0.016 -91.90 0.016 0.001 1. Coordinates are state plane (feet). 2. Elevation is mean sea level (feet). 3. Aspect is in degrees measured from north: negative is counter-clockwise from 0 to 180 and positive is clockwise from 0 to 180. 4. Slopes in the x- and y-direction are the absolute values of slopes in the polar coordinate system with zero degrees due east. 101 GRID DUERLHY ON 1 9; SLOPE OR ID OUERLF—‘IY Figure 14 Finite Element Grid Representing Hydrologic Response Areas based on slope. 102 C. Hydrograph Response. Several cases were investigated to test the accuracy of the proposed methodology. The first case was a representation of the watershed using a finite element grid of arbitrary spatial form. The spatial form of this grid is shown in Figure 13. The second case was a representation of the watershed using a hydrologic response area, finite element grid. Difficulties arose due to the anisotropic nature of the hydrologic response area grid which led to consideration of an isotropic finite element grid. 1. Arbitrary Finite Element Grid Model Calibration. Several modeling runs were performed in order to produce an outflow hydrograph that matched the actual outflow hydrograph. The initial Green and Ampt parameters were estimated after Brakenseik (1983). After calculating the runoff intensities for the three soils it was determined that the soils could be treated as a single soil possessing spatial homogeneity at the watershed scale. The total volume of the resulting runoff was 1.31 cm. The actual runoff reported was 0.38 cm. This indicated that a combination of increased hydraulic conductivity and depressional storage would be needed. The range of hydraulic conductivities due to crusting was estimated as 0.15 cm to 0.06 cm. Values of hydraulic conductivity and depressional storage were selected such that the peak and volume were predicted as accurately as possible. The 103 \ depressional storage was estimated as 0.55 cm for a recently tilled silty clay loam (personal correspondence from Dr. R.J. Rawls, ARS, Beltsville). The other Green and Ampt parameters were Wetting Front suction = 31 cm, Bulk Density = 1.0 g/ccm. The resulting hydrograph shown in Figure 15 has a peak of 0.114 m3/s at 14:31 hours on May 4, 1959. This is compared to 0.13 m’/s at 14:29 hours. The calculated hydrograph had a volume of 0.21 cm compared to an actual volume of 0.55 cm. The tabulated results, outflow at node- 4 and rainfall/runoff intensities are in Table 4. The hydrograph shown in Figure 15 has a larger volume than the actual hydrograph. The volume and peak could not be matched precisely. If sufficient depressional storage to match the outflow volume was removed, then the intense part of the storm was removed resulting in very low peak rates. The final calibration was found by inspection of the rainfall intensities that resulted from a varying hydraulic conductivities in the range reported above. The difficulty in finding a set of depressional storages and hydraulic conductivities that would produce a hydrograph suggests that the Green and Ampt model does not define a unique set of parameters. The formulation of the arbitrary finite element grid as linear elements with variable width also affects the resulting calibrated parameters. The form of the finite element grid affects the across 104 the width of the element. This assumption results since it makes the assumption of uniform flow governs the form of the hydrograph as the rainfall excess is routed downstream. Table 4 Arbitrary Finite Element Grid Outflow. RUNOFF CALCULATED ACTUAL TIME RAINFALL INTENSITY OUTFLOV OUTFLOW (HR) CM/HR FT/S M3/S M3/S 14.30 0.00 0.00E+00 0.00E+00 0.00E+00 14.32 0.00E+00 0.00E+00 0.00E+00 14.33 0.00E+00 0.00E+00 0.00E+00 14.35 2.67 0.00E+00 0.00E+00 0.00E+00 14.37 0.00E+00 0.00E+00 0.00E+00 14.38 0.00E+00 0.00E+00 0.00E+00 14.40 0.00E+00 0.00E+00 0.00E+00 14.42 12.70 3.608-05 7.188-04 0.00E+00 14.43 1.30E-04 1.27E-02 3.20E-03 14.45 1.008-04 4.688-02 2.528-02 14.47 3.6OE-05 8.08E-02 8.923-02 14.48 -1.BOE-07 1.028-01 1.288-01 14.50 7.92 -1.BOE-07 1.13E-01 1.158-01 14.52 -1.808-07 1.14E-01 1.04E-01 14.53 -1.80E-07 1.088-01 8.828-02 14.55 -1.8OE-07 9.568-02 7.288-02 14.57 -1.8OE-07 8.198-02 6.04E-02 14.58 1.83 -1.80E-07 6.88E-02 4.79E-02 14.60 -1.808-07 5.73E-02 3.92E-02 14.62 -1.8OE-07 4.77E-02 3.05E-02 14.63 -1.80E-07 3.97E-02 2.548-02 14.65 -1.8OE-07 3.328-02 2.03E-02 14.67 -1.80E-07 2.80E-02 1.74E-02 14.68 0.76 -1.8OE-07 2.37E-02 1.44E-02 14.70 -1.80E-07 2.023-02 1.14E-02 14.72 -1.8OE-07 1.73E-02 9.82E-03 14.73 -1.8OE-07 1.508-02 8.21E-03 105 .pmuu ucoeodm ouwcmm xuouumnu< Eoum :Qouoochz ma whoomm ESP—.00 040.904 1. 300.— 0091—0340 AOQDOIV flak Shin 000.¢u 000.¢.— 03.1w FLL- bl'PPthllPrb Pp pr 08.0w 8.0 l «0.0 I 00.0 I 3.0 fil 00.0 r.on .l 8.0 I. 00.0 I 00.0 .I 010 I. «.10 T.N10 nu.0 08d 3 h; goquDO Gum—944504.440 Q24 iDEUd. (“ROGER/8 an!!! 011M) mo 106 COMPARISON OF OUTFLOW - mmmmnmrtmn 0.12 0.11 "‘ 0.10 ‘ III!“ 007- 0.04“ 003" 0.02“ 0.01 ‘ CALCULATID OUT-1.0" (CUBIC KIT-IC/IICOND am I I T r I I I I T r I T 0.000 0.020 0.040 0M W 0.100 0.120 107011. 001710! (CUBIC lama/mum) Figure 16 Comparison of Outflows for the Arbitrary Finite Element Grid. 107 2. Hydrologic Response Area Finite Element Grid. The finite element grid shown in Figure 14 was formed using the spatial distribution of slope as defined by the ten foot elevation contours. The parameters used for this case were the same as those used in the calibrated, arbitrary grid model. This was chosen so that differences would not arise from two different calibrated sets of parameters. The effect of better defined finite elements formed from hydrologic response areas was expected to result in a more accurate prediction of the outflow hydrograph. Difficulty in obtaining a solution for the hydrologic response area finite element grid was encountered. During the solution, a negative flow depth occurred at node number eight. This effectively precluded solution since no flow could pass this element because of the connectivity associated with the finite element grid. The cause of this difficulty was investigated through extensive analysis of the matrices associated with the Galerkin finite element formulation for the surface water flow equations. These matrices are the [C], [bx] and [by] matrices as defined in [109]. The first difficulty examined was the first derivative term represented by [bx] and [by]. These matrices are asymmetric and therefore are difficult to integrate properly when anisotropy exists in the slope term. Anisotropy in the slope results from the two dimensional representation of overland flow. The 108 I control volume in Figures 2 and 3 assume that the flow is orthogonal to the x-y coordinate system. The Manning equation relates the total flow Q to the flow depth h as a function of slope. This slope is the slope in the principal direction and therefore the flow rate Q is a vector quantity. The resolution of the flow into the respective x-y direction components is accomplished after the computation of the flow in the principal direction. The flow rates Q in [109] are the respective flow rates per unit width in the x-y directions. Due to the derivation from the control volume and the form of the finite element formulation, each element must be rotated in a local coordinate system such that the finite element nodal coordinates are orthogonal to the principal direction of the slope. This rotation is performed for each coordinate pair using the rotation matrix x' cose sine x = [110] y' -sin0 cose y where 0 = angle of the principal slope in the global coordinate system. This rotation must be done before the integration of the element matrices. The integration of the element matrices is performed 109 I numerically for the isoparametric element using the Gauss- Legendre Quadrature. The Jacobian matrix is calculated in this integration procedure using the global coordinates. If anisotropy exists, then the global coordinates must be rotated as in [111] and the local used in place of the [global coordinates. When rotation occurs for a randomly oriented four node quadrilateral finite element, it becomes unclear as to which nodal coordinate pair must be the first pair associated with node i or {-1,-1) in the E-n coordinate system. This importance can not be over emphasized since the existence of the solution depends on it. The difficulty arises when integrating the [bx] and [by] matrices in [110]. A four node quadrilateral finite element that is rectilinear and oriented with the g-axis parallel to the x-axis and the n-axis parallel to the y- axis would' be integrated such that the lower left node closest to the origin would be associated with the (-l,-l) node in the g-n coordinate system. If this principle is not adhered to then the integrated result represented by the [b] matrices is not accurate and a solution is not achieved. A four node quadrilateral finite element was investigated to demonstrate the difficulty imposed by rotation and node ordering when integrating the [b] matrices. Under a steady rainstorm intensity the equilibrium value of the outflow equals the product of 110 1 intensity and the surface area of the finite element. The solution yields unit width nodal flow rates which must be integrated over the sides of the element in order to compare the outflow with the inflow. The finite element modeled is shown in Figure 17. The outflow is 0.4 m3/s for a 0.00001 m/s intensity over the 200x200 m finite element. Continuity is achieved since the inflow is the product of intensity and the surface area or 0.4 m’/s. Note that the nodal slopes are orthogonal to the global coordinate system and that no rotation to a local coordinate system was performed. The effect of rotation to a local coordinate system was investigated by observing the effect on the outflow that the direction of slope has. By assigning slopes of 45 degrees at the nodes and assigning boundary conditions at the nodes 1,2, and 4, outflow at node 3 results. The outflow was 0.424 m3/s. The inflow was 0.40 m’ls resulting in a 6% error. This is not a large error in this case. 111 we 0;: (01 I I >‘l : I - 4 E3 3__, .002 m3 51m ' - 3 -1 -1 "3.13351? ..... 3_, .002 m s m 1 Figure 17 Four Node Quadrilateral Finite Element Outflow. 112 A mOre serious error occurs when several finite elements are combined in a system with slopes that are not orthogonal to the coordinate system. In Figure 18 two elements are shown with nodal slopes at 45 degrees. The outflow totaled 1.155 m3/s whereas the inflow totaled 0.6 m3/s. These errors arise because of the difficulty in integrating the [b] matrices. These matrices are shown below for the four node quadrilateral when successive nodes are interpreted as node i in the integration. This amounts to extreme cases of rotation such that different nodes become the closest node to the global origin. The slope angle with respect to north and the first node are shown for each case. Slope Aspect = 90 First Node = l -2 2 l -1 [bx] = 200 -l 1 2 -2 I2 -1 l 2 -2 -2 2 l -l L . Slope Aspect = 90 First Node = 2 p- q -1 l 2 -2 [bx] = 200 -2 2 1 -l ‘12 -2 2 l -l L-1 l 2 -2J 113 Slope Aspect = 90 First Node = 3 '-2 2 1 ~11 [bx] = 200 -l 1 2 -2 I2 -1 l 2 -2 -2 2 l -l Slope Aspect = 90 First Node = 4 [-1 1 2 -2 [bx] = 200 -2 2 l -1 I2 -2 2 1 -l -l l 2 '21 Because of the incompatibility of the four node quadrilateral finite element with the computation of the [b] matrices when rotation occurs, it should not be used for the solution of the surface water equations when anisotropy exists in the slope. 114 :00 E N v' '9 . . "°+\’ 1 I ' I 3 -1 -l .006 01 s "I so s—J -—-'-.6;'4Irn|39:1 -—-. .006 mas1n'I1 IF ”I -*'.11 rnasJ F Figure 18 Two Element System Outflow. 115 The three node triangle was investigated in a similar manner. Figure 19 shows the three node triangle used to demonstrate the effect of rotation on the outflow and the integrated value of the [b] matrix. Successive nodes were used as the node associated with node i or the first shape function for the triangle finite element. The effect on the [b] matrices was observed for each of the following cases. Slope Aspect = 90 First Node 8 1 -l O 1 [bx] = 200 -l 0 l 6 -1 O l Slope Aspect 8 90 First Node = 2 -1 0 1 [bx] = 200 -l 0 1 6 -l 0 1 Slope Aspect = 90 First Node = 3 -1 0 l [bx] = 200 -1 0 1 6 -l 0 l 116 0) >4 3" >. 3 —> slope x-axis g - r 1 2 Figure 19 Three Node Triangle Isotropic Outflow. 117 The result is that the node ordering and hence the rotation has no effect on the magnitude or sign of the [bx] matrix. This indicates that the triangle is less sensitive to node ordering and rotation than the four node quadrilateral finite element in the effect on the [b] matrices. This fact recommends the use of the triangle for use in anisotropic surface water flow problems. Rotation of the three node triangle finite element is done so that the x'-y' axes are orthogonal to the principal direction of slope. This is done in the same manner as in [111]. The boundary conditions are applied to the triangle as with the quadrilateral. Another difficulty arises, however, since all three nodes must be specified as boundary values of zero when the principal slope direction is in a direction away from two sides as shown in Figure 20. The solution for the triangle shown cannot be obtained. This is due to the difficulty in specifying the boundary condition. If only one or two nodes are specified, then the boundary conditions are under specified for anisotropic slopes. If three nodes are specified then the element is a null solution since all nodes are zero for all time. The use of the three node triangle finite element for anisotropic flow requires that the element be oriented in the global coordinate system such that one side is parallel to the principal direction of slope. If this is done, then the triangle may be used with any node ordering and 118 I orientation in the global coordinate system since rotation does not result in an invalid integration of the [b] matrix. 119 Figure 20 Three Node Anisotropic Triangle. 120 I 3.. Isotropic Finite Element Grid Inorder to investigate the hydrologic response of this watershed using a two-dimensional finite element configuration it was necessary to create a rectilinear grid that is oriented such that each element is orthogonal to the principal direction of slope. Due to difficulties arising from node ordering and rotation, the slopes derived from the ARC/INFO TIN slope map in Figure 11 were used with a realigned aspect of 90 degrees or due west. A new grid that possesses only isotropic slopes was used. The node ordering preserves the ordering necessary to achieve correct integraton of the [b] matrices. This ordering is preserved since no rotation is necessary due to the isotropic slopes. Figure 21 shows the isotropic finite element grid derived from the hydrologic response area grid and the TIN slope map. 121 '- 3 G \ a In 1- - N at s- p- ("I v- N F Figure 21 Isotropic Finite Element Grid. 122 .pmuo acoanm.o_wcmm UMQoLuomH scum camumOLU>m Nm ousuwm 3005900 4-4003 h8.v« 08.v« 08.¢« P D - b P F D .- F P b b _ PL b + ESLBDO 09—24003 SE85 :3... 2e: «8.: 68.3 84.3 DbbbpbPIbDFhLbhbbbP 08.v« 08.¢« r TI ' r... IIrI I r I IIirI I 304 @530 000« .v F; AquEOdw Q24 QmfijDOAdd U 0 «0.0 «0.0 00.0 «0.0 00.0 00.0 8.0 00.0 00.0 «.0 ««.0 N«.0 0«.0 *«.0 0«.0 0«.0 b«.0 0«.0 (GKOOIE an can!!! sumo) 1101.11.00 123 The hydrologic response from the isotropic finite element grid is shown in Figure 22. The isotropic finite element grid yields a hydrograph that exceeds the actual outflow hydrograph. The peak outflow of 0.18 m3/s compares to the actual of 0.13 m’ls. The time to peak is at 14:30 which compares to the actual time to peak of 14:29 hours. The rainfall excesses used were the same as those used for the arbitrary finite element grid. The tabulated results, outflow at node 1 and 2 integrated across the edges of elements (1) and (8) are presented in Table 5. Table 5 Isotropic Finite Element Grid Outflow. INTEGRATED OUTFLOW CALCULATED ACTUAL TIME NODES 1,2 OUTFLOW OUTFLOW HRS M3/S M3/S 14.400 0.00 0.00 0.00 14.417 0.03 0.00 0.00 14.433 0.60 0.02 0.00 14.450 2.44 0.07 0.03 14.467 4.78 0.13 0.09 14.483 6.31 0.17 0.13 14.500 6.62 0.18 0.12 14.517 6.16 0.17 0.10 14.533 5.31 0.14 0.09 14.550 4.35 0.12 0.07 14.567 3.45 0.09 0.06 14.583 2.71 0.07 0.05 14.600 2.14 0.06 0.04 14.617 1.70 0.05 0.03 14.633 1.37 0.04 0.03 14.650 1.12 0.03 0.02 14.667 0.93 0.03 0.02 14.683 0.78 0.02 0.01 14.700 0.66 0.02 0.01 14.717 0.57 0.02 0.01 14.733 0.50 0.01 0.01 124 ««.0 owaouuomu ozu you msoauuso no .pmuu acoeoum oumcmm commuoaeou Anzounenflflhfll 2009 ESP—.00 a<0h0< «.0 8.0 p h - Ir v0.0 8.0 b r p - «0.0 "0.0 00.0 v0.0 00.0 r. 00.0 I 50.0 I 00.0 I 00.0 I «.0 I ««.0 l N«.0 ”«.0 .- v«.0 I 0«.0 .- 0«.0 r. b«.0 IrIT T I 9504.13.30 08« .4 #4: mo ZOmZZEEOO 0«.0 OIIDO) “0'31an GILV'IDO‘TV‘D GNOOII/ CHILI)! mm ouaumm 125 I 4. Statistical Comparison of Arbitrary Grid Outflows The statistical comparison of the arbitrary finite element grid outflow hydrograph with the actual outflow is facilitated by plotting the calculated versus the actual outflows. The close agreement between the calculated and actual values would result in a line plotted along a 45° axis. Figure 16 shows close agreement in the rising limb. Regression analysis presented in the Appendix represents values for both the ascending and recession limbs and a correlation coefficient of .92. The regression analysis was performed with Lotus 1-2-3 release 2.01. The regression analysis is a measure gf goodness of fit. This technique does not address the sensitivity of the output (flow depth) to the input (slope, Manning n, and rainfall excess) parameters. Nor does it define the uniqueness of the calibrated parameters in the solution. The hydrologic response area finite element grid did not result in an improvment over the arbitrary finite element grid. The more correct representation of the slope parameter by GIS analysis may improve the predictive capability of the finite element model but is not identyfiable by the approach taken. This approach assumed that the Green and Ampt parameters calibrated for the arbitrary finite element grid would result in accurate results when used with the hydrologic response area finite element grid. The difficulty with this approach is that the hydraulic conductivity and depressional storage are not 126 4 unique nor known apriori from measurable characteristics of the soil-landuse complex. Due to the difficulty in identifying a unique set of infiltration parameters as well as others such as Manning n, further calibration or statistical comparison was not pursued. D. Discussion Three cases were investigated, the arbitrary finite element grid, a hydrologic response area finite element grid and the isotropic finite element grid for use in modeling watershed outflow under an unsteady rainfall event. The arbitrary grid consisted of one-dimensional, linear elements of variable width. These linear elements were of the same spatial form as used by Peters, Blandford and Meadows (1983). This grid was termed arbitrary because they do not necessarily conform to spatially distributed parameters. The hydrologic response areas were defined using a Geographic Information System. The finite element grid used to represent the hydrologic response areas conforms to the spatially distributed parameters. The isotropic finite element grid was investigated because of significant difficulties in modeling randomly oriented finite elements that possess anisotropic slope. 1. Input Parameters The essence of the difficulty lies in the representation of the watershed topography using finite 127 elements. Two dimensional analysis of the surface water equations was complicated by the use of a four node quadrilateral for computation of first order derivatives when anisotropy required rotation to an orthogonal set of local axes. The rotation caused errors in the integration of the [b] matrices. The other finite element investigated was the three node triangle. This finite element could be rotated to the orthogonal local axes without error in the [b] matrices integration. However, the proper assignment of boundary values requires that the triangle be oriented such that one node lies on a stream line, i.e. an isotropic finite element grid. ‘ The arbitrary grid representation of slope, i.e., a uniformly sloped plane incorrectly represents the flow paths on a curvilinear surface. The flow paths are always perpendicular to the elevation contours if inertia is insignificant as in the kinematic wave equation. This amounts to a solution of the Laplacian equation for potential and streamlines. The equipotential lines are in this case the equi-elevation lines and the streamlines or orthogonal trajectories are the flow paths. Of course, depending on the scale of the modeled flow, micro-scale topography changes may obscure a smooth trend in flow path across the watershed. Significant difficulties arise when the element grid is not orthogonal to the principal direction of the slope as it varies over the watershed surface. 128 In the present watershed, the slope was defined by the contour intervals from the topographic survey. The variation of the elevation contours sets the scale of the spatial variation of the slope parameter. Between the contour lines, nothing is known of the micro-scale topography and therefore is assumed to be purely stochastic. The average slope within this interval was used to define the slope for the elements. The ARC-INFO, TIN program was used to determine the slope and aspect at each node. The TIN program simply determines the slope of the ground surface as defined by the digitized elevation contours. As more and more finite elements are used to represent the curvilinear, spatially nonhomogeneous watershed, the nodal slopes tend toward elemental representation of slope, and in the limit, they are the same. In the case of the arbitrary grid, the slope is averaged over the plane by taking several measurements within the plane. This does not allow accurate representation of the slope by the finite element grid since the spatial nonhomogeneity is lost by averaging. The linear element is not well suited to representing the curvature of the watershed surface. However, as more and more linear, one-dimensional elements of decreasing size are used to represent the watershed, the representation of topography improves as with any other element. The difficulty arises from the use of one-dimensional elements 129 of variable width to represent two-dimensional,_ spatially variable, parameters. Difficulties were found when using a two-dimensional scheme in the modeling of two-dimensional, spatially variable parameters that are anisotropic. In such a case the finite element grid must be orthogonal to the principal direction of slope at each element. A two-dimensional finite element grid is capable of more accurately representing the spatial nonhomogeneity of the spatially variable parameters. The nodal slopes were determined by the value of the slope and aspect represented by the TIN Slope map in Figure 11. These slopes vary over the element in the case of the four-node quadrilateral because this element is capable of representing a warped surface. The average slope could have been determined by integrating the slope over the element and dividing this integrand by the area, and the consistent stress method used to obtain the nodal slope values to be used in the modeling. This procedure, however, ignores micro-scale slope variation interior to the finite element. It is though, responsive to meso-scale 'slope variation represented by the spatial form of the lO-foot contour lines. Two-foot contour intervals could have been used but this would have resulted in considerably more finite elements than the first case of the arbitrary grid. This increased fineness in the grid representation would have increased the accuracy due to the increased number of 130 elements not due, solely, to the better finite element representation of a spatially variable parameter. Soil properties were found to produce very little difference in the Green and Ampt parameters and consequently in the rainfall excess intensities. The watershed was therefore treated as one soil. The effect of the lumping of the soils is obscured by the difficulty in defining a unique set of infiltration parameters. Based on evidence from other researchers (Brakenseik and Onstad, 1977) the relative error in peak runoff to hydraulic conductivity is 2.68 percent. If a one percent error exists in the hydraulic conductivity then a 2.68 percent error in peak runoff rate results. In this case the hydraulic conductivities predicted according to the method of Rawls, Lane and Nicks (1987) method were as follows Hastings Silt Loam, K = 0.0198 cm/hr, Hastings Silty Clay Loam, K = 0.0087 cm/hr, Colby Silt Loam, K = 0.0167 cm/hr, Average, K = 0.0151 cm/hr This average was not actually used in the modeling run since the effect of crusting obscured the range of different soil hydraulic conductivities. The hydraulic conductivities are estimated to range from 0.15 to 0.06 cm/hr due to the crusting. This range is estimated by using the variation of the CBC/Clay ratio from 0.2 to 0.65. The hydraulic conductivity used in the modeling of the May 4, 1959 rainstorm event, was K = 0.056 cm/hr. This 131 corresponds to the CEC/CLAY ratio of 0.65, which was the high end of the range of data used in the WEPP Model development (cf. Rawls, Lane and Nicks). The actual CBC/CLAY ranged from 0.67-0.79. The Green and Ampt parameters used in both the arbitrary finite element grid and the hydrologic response area finite element grid are not a unique set. The rainfall intensitites used in calibration were found by inspection from a range of intensities produced by a range of hydraulic conductivities. Extreme difficulty was encountered in matching both volume and peak rate of outflow. When these same rainfall intensities were used with the hydrologic response area finite element grid a more accurate solution did not result. In fact a poorer agreement was found. This suggests that not only are the Green and Ampt parameters nonunique but the calibration is dependent on the finite element grid representation of the watershed domain. Depending on the spatial form of the finite element grid, different sets of input parameters will result from the calibration procedure. This difficulty obscures the advantages of better representation of spatially nonhomogeneous parameters. 2. Numerical Errors Numerical errors also result from the mathematical formulation of the solution. These errors result from several sources within the finite element method itself. Other errors in represetation of the physical phenomenon 132 I have been previously discussed. Numerical oscillations and stability errors arise from the size of the eigenvalues in the equation [c1{i} = {9*} . [1121 The eigenvalues are calculated by solving for the eigenvector [E]. {A} [31(2} - [113] or, {A} [3112} [114] where, {z} a solution vector of flow depths. Writing [114] as [c][a]{z} = {F*} [115] . . ‘1 and multiplying by [E] , [ul'ltc1[31{é} = [a] {F*} ' [1161‘ -1 The matrices [E] [C][E] form the eigenvector [A], For any row i, 1121 and i. or 2i 1331 f"\ ..0 Z 1 0.0.0000 z 2 000.000 {?3 . 0 0 . .o '1 0 z n . n k) R*i R*/Ai I R*/li dt = Writing [119] as “I’NN'O-‘I-j Y I A weeewww I- K. J At(R*/Ai] the recursive formulation results in aAt + z 0 aAt + 21 = 2a naAt + Z0 . At + z 0 [117] [118] [119] [120] [121] [122] 134 S Equation [122] does not amplify errors that may occur in the right-hand side term naAt since it does not contain any 2 terms. Because of this the solution of the kinematic wave equaitons in general may proceed without numerical oscillations or stability errors. This is a consequence of the partial differential equation form and the finite difference solution in the time domain. Nonlinearity in the partial differential equation however can pose difficulties in the solution algorithm. It is necessary to start with sufficiently close values of flow depth at each time step such that the iterative solution converges. If large step increases in rainfall excess intensities. occurs or large values of rainfall excesses occur, the iterative solution will not converge but rather flow depths becomes negative. If initial estimates of flow depth at each time step are used that are not the previous solution values, the solution again does not converge. Convergence was achieved by using the dry bed initial values at the first time step and the previous- time flow depths as the beginning values for iterative solution at succeeding time-steps. As evidenced in the literature, the difficulties encountered in applying a deterministic, distributed parameter model to watershed hydrology fell into three categories. These categories largely center around the types of errors encountered in the literature for the type of modeling performed for this research. These categories 135 1 are numeric errors, model equation errors, and parameter estimation errors. The approaches taken to meet the project objectives were performed in order to reach the overall goal of this research to more accurately predict the actual outflow hydrograph from a watershed of nonuniform, spatially distributed parameters. Conclusions derived from this research follow together with recommendations and future applications of the method developed. Numeric errors are those that arise from the method itself such as inaccuracies and instability in the time solution of the finite element method. These errors arise due to the form of the time dependent equation that is solved by standard finite difference techniques such as in equation [92]. Numeric errors can be subdivided into physical reality, numerical oscillations, accuracy, and stability. Physical reality is observed whenever rainfall excess added to a node causes surrounding nodes to also increase. In the solution of heat flow and groundwater flow, it is possible to obtain solutions that for initial time steps the temperature or pressure head decreases when it should be increasing (Segerlind, 1984). A physical reality error results whenever this occurs. Numerical oscillations occur when at each time step the solution is first higher then lower than the true solution thus causing oscillations about the solution. 136 \ Stability errors occur when an incremental error at each time step grows until the solution deteriorates. As a result of the equation [116], the solution of the kinematic wave equation as represented by [106] does not suffer from stability and oscillation errors prevalent in other governing equations such as the field equation for heat flow. The only other possible source of numeric errors are those arising from the accuracy of the ordinary differential equations, e.g. the time dependent equation [105] in representing the partial differential equations, i.e., the conservation of mass equation [55]. This error is accounted for by observing the Courant condition. The Courant condition is a consequence of the partial differential equation theory, the Method of Characteristics. It was observed through the choice of the time step for the maximum rainfall excess and the longest plane. This choice of time step resulted in a solution that was accurate and free of oscillations or instability. The computation of the Courant condition from known physical parameters allows selection of the time step prior to modeling using the finite element method. As seen in Figures 14 and 15, a solution free of instability and oscillations resulted. Model equation errors arise from the simplification of the full dynamic equation by the kinematic or diffusion analogy and kinematic shocks. The kinematic shocks that plagued previous researchers were avoided in the method 137 4 developed by this research. The method used here calculates the nodal flow rates using nodal slopes and Manning n values. This improvement over previous methods allowed an accurate solution of the kinematic wave equations using the finite element method. Both linear one-dimensional and two-dimensional elements were used successfully in this solution in the modeling of a two- dimensional domain of spatially nonhomogeneous slopes. This two-dimensional domain possessed complex topographical curvature as represented by the TIN slope map in Figure 11. This slope map provided the basis for the selection of nodal slope values that resulted in the solution of the kinematic wave equation free of kinematic shock. Parameter estimation errors arise from uncertainty and from spatial variation of infiltration, rainfall, roughness and other parameters over the watershed domain. The spatial variation of the watershed parameters were represented by the finite element grid by assigning to each element and 'node the appropriate parameters. Rainfall excess was assigned as an element constant as well as assumed constant over the watershed. Errors in rainfall excess arise from many sources including the assumption that infiltration, rainfall, roughness and other parameters are uniform over the entire watershed. Where the accuracy and the number of sampling points are sufficient, spatial statistics may be used to advantage to interpolate with known variance at a specified 138 location the value of a parameter. The location of the parameter is dictated by the finite element grid which in turn is dictated by the spatial variability of the parameters. Lack of knowledge apriori of the infiltration parameters and the inability to select an unique set by calibration limit the viability of calibrating. Proper representation arises from the proper selection or the use of elements that maintain faithful representation of the hydrologically homogeneous character of subareas within the watershed. Improper representation may result from an incomplete knowledge of the parameters and the variation over the watershed or with time during the modeling process. The Geographic Information System allowed proper representation of the spatially nonhomogeneous parameter, slope. By utilizing the hydrologic response area, finite-element grid proper representation occurred. Accurate modeling of the hydrologic response areas within the watershed during an unsteady rainstorm is possible if the infiltration parameters used to define rainfall excess is sufficiently well known. Better definition of the rainfall-runoff relation is needed. The finite element method holds promise as a mathematical model capable of accurately and efficiently solving the distributed, deterministic surface water equations in a watershed. The Geographic Information System holds promise as an efficient means of handling the large volume of input 139 \ data required by distributed, deterministic models. The goal of this research has been achieved by accomplishing the following three objectives using the following approaches. The objectives and approaches were: Objective 1. Define hydrologic response areas that exhibit similar soil infiltration parameters, surface roughness, and slope for a given watershed. Approach: A Geographic Information System was used to search, smooth and aggregate areas of similar soil infiltration parameters, surface roughness, and slope thus producing the specific hydrologic response areas. ' The approach to Objective 1 utilized the ARC-INFO GIS to digitize the maps of soils, landuse, and topography. Based on the infiltration parameters associated with the soils, the Manning n associated with the landuse, and the slope calculated from the 2-foot contour elevation map; hydrologic response areas were produced. Due to the homogeneity of the soils and landuse, lepe was the only spatially nonhomogeneous parameter modeled. The ARC-INFO Triangular Irregular Network (TIN) program was used to define the slope and aspect from the 2-foot contour elevation map. From the TIN Slope map, nodal slope values were determined for each node of each element. The Slope map shown in Figure 11 was produced for a 1% interval. The hydrologic response area finite element grid was based on the 10-foot contour elevations. The GIS proved to be an 140 effective means of searching, smoothing and aggregating values of the slope parameter and thus defining hydrologic response areas for the anisotropic finite element grid. The GIS was also used to determine the slopes at the nodes for the isotropic finite element grid, however, the aspects were aligned such that the slope was parallel to the sides of the elements representing streamlines. Objective 2. Apply the finite element method to the specific hydrologic response areas to compute and route the overland flow to the outlet. Approach: The rate and volume of infiltration was modeled by the Green and Ampt infiltration equation. The equation parameters were calibrated for the watershed using an arbitrary finite element grid. The rainfall excess thus defined becomes the lateral inflow for use in solution of the overland flow equations. The approach to Objective 2 utilized the Green and Ampt infiltration equation. The Green and Ampt parameters along with the rainfall excess were computed using the WEPP project programs with the aid of Dr. Walter Rawls, ARS Water Data Laboratory, Beltsville Maryland. For a description of the procedures used see Rawls, Lane and Nicks (1987). The initial estimates of the Green and Ampt parameters used measured soil properties. The rainfall excesses computed for the three soils were nearly identical during the rainstorm event. Due to the similarity, the watershed was modeled as having one soil. 141 1 Calibration of the Green and Ampt parameters was performed by varying the hydraulic conductivity and the initial abstraction until the calculated volume and peak outflow rate matched as closely as possible the actual outflow rate. The calibration was performed for the arbitrary finite element grid. The final values of initial abstraction and hydraulic conductivity were not determined directly but rather the rainfall excess intensities were determined by inspection of several sets of hydraulic conductivities and initial abstraction values. These Green and Ampt parameters appear to vary during the storm. The simply infiltration model was not capable of predicting the correct Green and Ampt parameters and resulting rainfall intensities. After calibrating the Green and Ampt parameters for the arbitrary finite element grid, the hydrologic response area grid was modeled using the calibrated rainfall excess. The purpose for this was to investigate the effect of better description of the spatial parameters such as slope by the two dimensional, hydrologic response area grid. The Galerkin finite element formulation was applied to the kinematic wave equations producing the flow rate and flow depth solution. The solution was performed in the time domain during the unsteady rainstorm event. The time solution was performed using the central difference, finite difference form. This formulation resulted in a solution free from the physical reality, numerical oscillations, 142 1 instability, kinematic shock and accuracy errors reported previous literature. It was necessary to form an isotropic finite element grid such that the solution may be obtained. The resulting finite element grid was a better representation of the spatially varying slope and the solution to a spatially varying flow depth in two-dimensions. The relative accuracy of the hydrologic response area finite element grid to the arbitrary finite element grid is unknown due to uncertainty in the infiltration parameters and resulting rainfall excesses. Objective 3. Compare the accuracy of the outflow hydrographs to the actual outflow hydrograph for a given rainstorm event for the following two cases: a. A finite element grid that is of an arbitrary spatial form. b. A finite element grid formed from hydrologic response areas defined by the Geographic Information System. Approach: The rate and volume of runoff was modeled by the finite difference/finite element method of solving the kinematic wave equation for overland flow. The excess rainfall was defined by the infiltration equation. The outflow hydrograph was calculated for the outlet of the watershed for the two cases described above. A finite 143 element model, utilizing linear, one-dimensional and two- dimensional finite elements was used to compute the overland flow equations. Lotus 1—2-3 release 2.01 was used for the solution of the linear, one-dimensional finite element grid. A FORTRAN computer program identical in algorithm to the Lotus 1-2-3 spread sheet was used for the linear, two-dimensional finite element grid. The spread sheet formulas and the computer program are contained in the Appendix. The validity of the method was checked by comparing the computed and analytical solution outflow hydrographs. The validity was checked for a single plane with the results contained in the Appendix. The finite element solution of the kinematic wave equations when solved with nodal rather than elemental flow rates produces a close approximation to the analytic solution and resulted in a close approximation of the actual outflow hydrograph for this rainstorm event. The flexibility of this method allows modeling of watershed runoff using spatially nonhomogeneous parameters. The spatial variability must be lumped interior to the finite element. Finite elements corresponding to the spatially nonhomogeneous parameter allow a much more accurate representation of the parameters. It was not determined whether the hydrologic response area, finite element grid produced a better agreement between the calculated and actual hydrograph due to the uncertainty in the infiltration parameters. 144 The method developed thfough this research calculates the two-dimensional solution to the kinematic wave equation hydrograph for a watershed of nonuniform, spatially nonhomogeneous parameters such as slope. Further, this method is free from kinematic shocks which had prevented earlier researchers from solving the problem as a continuum. V. CONCLUSIONS AND RECOMMENDATIONS The approaches taken to meet the project objectives were performed in order to reach the overall goal of this research--to more accurately predict the actual outflow hydrograph from a watershed of nonuniform, spatially distributed parameters. Conclusions derived from this research follow together with recommendations and future applications of the method developed. A. Conclusions The method developed through this research calculates the two-dimensional solution to the kinematic wave equation hydrograph for a watershed of nonuniform, spatially nonhomogeneous parameters such as slope. From this research it may be concluded that: l. The GIS was effective in searching, smoothing and aggregating values of the spatially variable slope parameter for use in the finite element model. 2. The Green and Ampt infiltration model used did not adequately predict the rainfall excess intensities in order to accurately simulate runoff from an actual storm using the finite 145 146 element model. 3. The finite element model developed through this research results in a solution free from kinematic shock for the two-dimensional kinematic wave equation in a watershed continuum. B. Recommendations The ability of the Geographic Information System and the finite element method to model and display modeled results of the watershed surface runoff recommends its use for further research. It is recommended that: 1. Further research is needed to better measure and describe the spatially variable parameters affecting - surface flow, particularly infiltration. 2. The use of geostatistics may provide valuable information as to the detail required to accurately model spatially variability in the input parameters. 3. This research should be extended to consider both overland and channel flow for larger watershed systems. This extension should include diffusion and full dynamic equation modeling. 4. Better calibration techniques should be investigated that are capable of handling the 148 flow depth (mm) Figure 24 Spatial Distribution of Peak Runoff Flow Depth. 149 the nodes were used in producing the contours. The peak ‘values occur at different times during the storm runoff period. The value of predicting the spatial distribution of surface flow is realized if the distribution of concentration of nutrients or pesticides is capable of being predicted for a watershed. This ability provides insight into the location and source of contaminants in overland flow. The flow depth shown in Figure 24 may be interpreted as the highest concentrations of contaminants. Since the flow rate over any particular location during a storm event is related to the flow depth by the Manning or Chezy equation', the concentration expressed as a mass fraction of the flow rate will correspond to the flow depth. Further research is needed to better describe the spatially variable parameters affecting surface flow, particularly infiltration. The method developed through this research provides a better description of the hydrologic processes in a two-dimensional domain. It also ' provides insight into the transport phenomena of agricultural pollution by pesticides and nutrients in surface and subsurface water as affected by overland flow and infiltration for an agricultural watershed. VI. REFERENCES Abdel-Razaq, A.Y., W. Viessman, and J.W. Hernandez. ”A Solution to the Surface Runoff Problem," Journal of the Hydraulics Division, ASCE, Vol. 93, No. HEE’ (November 1 - 0 Andrews, R.G. "The Use of Relative Infiltration Indices in Computing Runoff." Paper presented at the American Society of Agricultural Engineers, 1954. In R.E. Rallison, q.v. Band, L.E. "Topographic Partition of Watersheds with Digital Elevation Models,” Water Resources Research, Vol. 22, No. 1 (January 1986). Puinshed by the American Geophysical Union. Bartholic, J. and K. Kittleson. "Temperature and Reflectance Monitoring from Satelites as an Indication of Shift and Impact of Vegetation Change." Paper presented at the Nineteenth International Symposium on Remote Sensing of Environment, Ann Arbor, Michigan, October 1985. Biswas, A.K. "Beginning of Quantitative Hydrology," Journal of the Hydraulics Division, ASCE, Vol. 94, No. HYS ISeptemEer 1968), 1299-1316. Brakenseik, D.L. and C.A. Onstad. "Parameter Estimation of the Green and Ampt Infiltration Equation," Water Resources Research, Vol. 13, No. 6 (December 1977). Published by the American Geophysical Union. and W.J. Rawls. "Agricultural Effects on Soil Water Processes, yPart II: Green and Ampt Parameters for Crusting Soils," Transactions of the ASAE, 26:(l983), 1753-1757. Brutsaert, W. "The Initial Phase of the Rising Hydrograph of Turbulent Free Surface Flow with Unsteady Lateral Inflow," Water Resources Research, Vol. 4, No. 6 (December 1968), 1189-1192. 4_. "De Saint Venant Equations Experimentally Verified," Journal of the Hydraulics Division, ASCE, Vol. 97, No. HY9 (September 1971), 1387-1401. "Central Great Plains Experimental Watershed, A Summary Report of 30 Years of Hydrologic Research," Agricultural 150 151 I Research Service, USDA, n.d. Chen, C.L. and Ven Te Chow. ”Formulation of the JMathematical Watershed-Flow Model," Journal of the twechanics Division, ASCE, Vol. EM3 (June I971), 809-828. Chu, S.T. "Infiltration During an Unsteady Rain," Water Resources Research, Vol. 14, No. 3 (June 1978). Published by the American Geophysical Union. Crandall, S.H. Engineering Analysis, A Survey of Numerical Procedures. McGraw-Hill Book7Company, Inc., New Grace, R.A. and P.S. Eagleson. "The Modeling of Overland Flow," Water Resources Research, Vol. 2, No. 3 (1966), 393-403. Grayman, W.M. ”Land-Based Modeling System for Water Quality Management Studies," Journal of the H draulics Division, ASCE, Vol. 101, No. HYS (May I975), 567-580. Gupta, S.K. and S.I. Solomon. "Distributed Numerical Model for Estimating Runoff and Sediment Discharge of Ungaged Rivers, 1. The Information System," Water Resources Research, Vol. 13, No. 3 (June 1977). PuBIished* By the American Geophysical Union. Guymon, G.L. "Application of the Finite Element Method for Simulation of Surface Water Transport Problems,“ Institute of Water Resources Report No. IWR-21, University of Alaska, Anchorage, 1972. Handbook of Chemistry and Physics, 40th ed., Chemical Rubber Publishing Co., Cleveland, Ohio, 1972-1973. Henderson, F.M. Open Channel Flow, Macmillan Publishing Co., Inc., New York, 1966, pp. 288-324. Huggins, L.F. and J.R. Burney. "Surface Runoff, Storage and 'Routing, Chapter 5." In Hydrologic Modeling of Small Watersheds, edited by C.T. Haan, ASAE Monograph No. 5 (I962), I69-225. Huggins, L.F. and E.J. Monke. ”The Mathematical Simulation of the Hydrology of Small Watersheds," Water Resource Research Center Technical Report No. 1, Purdue University, 1966. Jayawardena, A.W. On The Application of the Finite Element Method to Catchment ModeIing. Ph.D. Thesis, University of London, King's college, 1976. 152 and J.K. White. ”A Finite Element Distributed Catchment Model, II. Application to Real Catchments," Journal of Hydrolo , Vol. 42, No. 3/4 (July 1979). PfiBlished by the E sevier Scientific Publishing Company, Amsterdam, Oxford, New York. Judah, O.M. Simulation ofRunoff Hydrographs from Natural watersheds by Finite EIement Method, Ph.D. Thesis, Virginia Polytechnic Institute and State University, Blacksburg, Virginia, 1973. , V.O. Shanholtz, and D.N. Contractor. ”Finite Element Simulation of Flood Hydrographs," Transactions of the ASAE, 1975, pp. 518-522. Kibler, D.F. and D.A. Woolhiser. The Kinematic Cascade, Colorado State University, Fort CoIlins, Publication No. 39, March 1970. Liggett, J.A. Unsteady Flow in Open Channels, Volume I, edited by Mahmood and Yevjevich, Water Resources Publications, Fort Collins, Colorado, 1975. ‘ Mein, R.G. and C.L. Larson, Modeling Infiltration Component of the Rainfall-Runoff Process. Water Resources Research Center, University of Minnesota, Graduate School, 1971. Mockus, Victor. "Estimation of Total (and Peak Rates of) Surface Runoff for Individual Storms," Exhibit A of Appendix B, Interim Survey Report Grand Neosho River Watershed, USDA, December 1, 1949. In Rallison, q.v. Morris, E.M., K. Blyth, and R.T. Clarke. "Watershed and River Characteristics, and Their Use in a Mathematical Model to Predict Flood Hydrographs,” in Remote Sensing Application in Agriculture and Hydrology: ed. by Georges Fraysse, A.A. Balkema, Rotterdam, 1980. Morris, E.M. and D.A. Woolhiser. "Unsteady One-Dimensional Flow over a Plane: Partial Equilibrium and Recession Hydrographs," Water Resources Research, Vol. 16, No. 2 (April 1980), 355-360. Peters N., G. Blandford, and M.E. Meadows. "Finite Element Simulation of Overland Flow," Proceedings of the Specialty Conference on Advances in Irrigation and Drainage: SurvivingExternaIPressures, ASCE, July 1983, pp. 466- 475. "Phase I Oconee Basin Pilot Study Trail Creek Test." Hydrologic Engineering Center, U.S. Army Corps of Engineers, Davis California, September 1975. I! {11"} 5’... . . 153 Potter, M.C. and J.F. Foss, Fluid Mechanics. Great Lakes Press, Inc., Okemos, Michigan, 1982. Rallison, R.E. "Origin and Evolution of the SCS Runoff Equation," Proceedings of the Symposium on Watershed Managgment '80, ASCE, 1980. Rawls, W. J., D.L. Brakensiek, and B. Soni. "Agricultural Management Effects on Soil Water Processes, Part I: Soil Water Retention and Green and Ampt Infiltration Parameters." Transactions of the ASAE 26 (1983). , L.J. Lane, and A.D. Nicks. "Hydrologic Components." In Compilation of Water Erosion Prediction Project (WEPP) Papers, presented at the 1987 International Winter Meeting of ASAE, Chicago, Illinois, 1987. Segerlind, L.J. Applied Fipite Element Analysis, lst ed., John Wiley and Sons, New York, 1976. . Applied Finite Element Analysis, 2nd ed., John Wiley and Sons, New York, 1984. - Smith, R.E. and D.A. Woolhiser. "Mathematical Simulation of Infiltrating Watersheds," gydrologprapers, No. 47, Colorado State University, Fort Collins, January 1971. Smith, V.E. Philosophical Physics, Harper & Brothers, New York, 1950, p. 8. Taylor, C. "A Computer Simulation of Direct Run-Off," in Finite Elements in Water Resources, edited by W.G. Gray, G.F.Pinder, and C.A. Brebbia, Pentech Press, London, 1976. , G. Al-Mashidani, and J.M. Davis. "A Finite Element Approach to Watershed Runoff," Journal of Hydrology, Vol. 21, No. 3 (March 1974), 231-246. and P.S. Huyakorn. ”A Comparison of Finite Element Based Solution Schemes for Depicting Overland Flow," Applied Mathematical Modelling, edited by C.A. Brebbia, Vol. 2, No. 3 (September 1978), 185-190. Vitruvius, Ten Books on Architecture. Translated by Morris H. Morgan, Dover Publications, Inc., New York, 1960, pp. 225-232. Woolhiser, D.A. Unsteady Flow in Open Channels, Volume II, edited by Mahmood and Yevjevich, Water Resources Publications, Fort Collins, Colorado, 1975, pp. 485-508. and J.A. Liggett. "Unsteady, One-dimensional Flow over a P1ane--the Rising Hydrograph," Water Resources 154 Research, Vol. 3, No. 3 (Third Quarter, 1967), 753-771. Yen, B.C. "Open Channel Flow Equations Revisited," Journal of the Engineering Mechanics Division, ASCE, Vol. 99, No. EMS (October 1973), 979-1009. Zobrist, A.L. "Elements of an Image Based Information System," Manual of Remote Sensing, 2nd ed., Volume 1, American Society of Photogrammetry, p. 946. APPENDIX Modeling Data 156 The following is a listing of the Lotus 1-2-3 work sheet for the arbitrary finite element grid solution. The cell formulas are printed unformatted for each section of the worksheet separated by headings describing the function of the section. ELEMENT DATA READING A1: [w12] 'ELEMENT Bl: 'DATA H1: [W12] 66 A3: [W12] Aelement B3: “nodel C3: [W12] “node2 D3: [W12] “node3 E3: [W12] Anode4 F3: [W12] Aslopel GB: [W12] “slope2 H3: [W12] “slope3 I3: [W12] “slope4 J3: [W14] “n1 K3: “n2 L3: [W12] “n3 M3: An4 N3: Awidth O3: “ROLD P3: [W12] “RNEW A4: [W12] “1 B4: 0 C4: [W12] 171 F4: [W12] 0.0406 G4: [W12] 0.0406 J4: [W14] 0.035 K4: 0.035 N4: 353 O4: 0 P4: [W12] 0 A5: [W12] A2 BS: 171 CS: [w12] 342 F5: [W12] 0.0406 GS: [W12] 0.0669 J5: [W14] 0.035 K5: 0.035 NS: 333 A6: [W12] “3 B6: 342 C6: [W12] 507 F6: [W12] 0.0669 G6: [W12] 0.0562 J6: [W14] 0.035 K6: 0.035 N6: 206 157 t ELEMENTAL EQUATIONS or THE FORM [C][A} = [F] A9: A10: C10: D10: E10: F10: 610: H10: A11: B11: D11: E11: F11: Gll: H11: Ill: A12: 812: D12: E12: F12: 612: H12: 112: A13: A14: C14: D14: E14: F14: 614: H14: A15: B15: D15: E15: F15: 615: H15: I15: A16: B16: D16: 316: F16: 616: H16: [w12] 'ELEMENT NO 1 [W12] “[C] [W12] A[Mnew [w12] '[F}= [w12] A[c] [W12] “*[AOLD) [W12] “DT/2[K]* [w12] 'QOLD+QNEW [w12] 2*($C$4-$B$4)*$N$4/6 ($C$4-$B$4)*$N$4/6 [w12] +$E$ll+$F$ll+$G$11+$H$ll+$I$11 [w12] 2*($C$4-$B$4)*$N$4*$B$27/6 [w12] ($C$4-$B$4)*$N$4*$B$28/6 [W12] -($B$37/2)*(-1)*((1-$B$38)*$D$27+$B$38*$D$33) [W12] -($B$37/2)*(1)*((1-$B$38)*$D$28+$B$38*$D$34) [W12] +$B$37*($C$4-$B$4)/2*((1-$B$38)*$N$4*$O$4+ $B$38*$N$4*$P$4) [w12] ($C$4-$B$4)*$N$4/6 2*(scs4-SBS4)*$N$4/6 [W12] +$E$12+$F$12+$6$12+$H$12+$I$12 [w12] ($C$4-$B$4)*$N$4*$B$27/6 [w12] 2*($C$4-$B$4)*$N$4*$B$28/6 [W12] -($B$37/2)*(-1)*((1-58338)*$D$27+$B$38*$D$33) [W12] -($B$37/2)*(1)*((1-$B$38)*$D$28+$B$38*$D$34) [W12] +$B$37*($C$4-$B$4)/2*((1-SBS38)*$N$4*$O$4+ $B$38*$N$4*$P$4) [W12] 'ELEMENT N02 [w12] “[c] [W12] A[Mnew [w12] '[F}= [w12] *[c] [w12] “*{AOLD) [w12] “DT/2[K]* [w12] 'QOLD+QNEW [W12] 2*($C$5-$B$5)*$N$5/6 ($C$5-$B$5)*$N$5/6 [w12] +$E$15+$F$15+$G$15+$H$15+$I$15 [w12] 2*($C$5-$B$5)*$N$5*$B$28/6 [w12] ($C$5-$B$5)*$N$5*$B$29/6 [W12] ’($B$37/2)*(-1)*((l-$B$38)*$D$28+$B$38*$D$34) [W12] -($B$37/2)*(l)*((1-$B$38)*$D$29+$B$38*$D$35) [w12] +$B$37*($C$5-$B$5)/2*((1-$B$38)*$N$5*$O$4+ $B$38*$N$5*$P$4) [w12] ($C$5-$B$5)*$N$5/6 2*($C$5-$B$5)*$N$5/6 [W12] +$E$16+$F$l6+$G$16+$H$16+$I$16 [W12] ($C$5-$B$5)*$N$5*$B$28/6 [W12] 2*($C$5-SB$5)*$N$5*$B$29/6 [W12] -($B$37/2)*(-1)*((1-$B$38)*$D$28+$B$38*$D$34) [W12] -($B$37/2)*(1)*((1-SBS38)*$D$29+$B$38*$D$35) I16: A17: A18: C18: D18: E18: F18: 618: H18: A19: B19: D19: E19: F19: 619: H19: 119: A20: B20: D20: E20: F20: 620: H20: I20: 158 [W12] +$B$37*($C$5-$B$5)/2*((1-$B$38)#$N$5*$O$4+ $B$38*$N$5*$P$4) [w12] 'ELEMENT N03 [w12] “[C] [W12] A{h}new [w12] 'irl- [w12] “[c] [w12] “*{AOLD) [w12] “DT/2[K]* [w12] 'QOLD+QNEW [w12] 2*($C$6-$B$6)*$N$6/6 ($C$6-$B$6)*$N$6/6 [W12] +$E$19+$F$l9+$6$19+$H$19+$I$19 [w12] 2*($C$6-$B$6)*$N$6*$B$29/6 [w12] ($C$6-$B$6)*$N$6*$B$30/6 [w12] -($B$37/2)*(-1)*((1-$B$38)*$D$29+$B$38*$D$35) [w12] -($B$37/2)*(1)*((l-$B$38)*$D$30+$B$38*$D$36) [W12] +$B$37*($C$6-$B$6)/2*((1-$B$38)*$N$6*$O$4+ $B$38*$N$6*$P$4) [w12] ($C$6-$B$6)*$N$6/6 2*($C$6-$B$6)*$N$6/6 [W12] +$E$20+$F$20+$6$20+$H$20+$I$20 [w12] ($C$6-$B$6)*$N$6*$B$29/6 [w12] 2*($C$6-$B$6)*$N$6*$B$30/6 [W12] -($B$37/2)*(-1)*((1-$B$38)*$D$29+$B$38*$D$35) [W12] -($B$37/2)*(1)*((1-SB$38)*$D$30+$B$38*$D$36) [w12] +$B$37*($C$6-$B$6)/2*((1-$B$38)*$N$6*$O$4+ $B$38*$N$6*$P$4) FLOW DEPTHS AND FLOW RATES AT EACH NODE, BOTH NEW AND OLD VALUES DEFINED BY THE MANNING EQUATION. INVERSE AND SOLUTION MACROS ARE PRECEEDED BY '/. A27: 827: C27: 027: F27: 627: A28: 828: C28: D28: F28: 628: A29: 829: C29: D29: A30: B30: C30: D30: [W12] Ahlold U o [w12] “Qlold [w12] l.486/$J$4*$F$4“0.5*$B$27“(5/3)*$N$4 [w12] 'INVERSE U [w12] 'SOLUTION\D [w12] “h201d U 0.0017120387 [W12] “Q201d [w12] (1.486/$K$4)*$F$4“0.5*$B$28“(5/3)*($N$4+$N$5)/2 [w12] '/DMI~~~ [w12] '/DMM~~~ [w12] “h301d U 0.0037515851 [W12] “Q301d [W12] 1.486/$K$5*$G$5“0.5*$B$29“(5/3)*($N$5+$N$6)/2 [w12] “h4old U 0.0072532181 [w12] “Q4old [w12] (1.486/$K$6)*$G$6“0.5*$B$30“(5/3)*$N$6 F32: 632: A33: B33: C33: D33: E33: F33: A34: B34: C34: D34: A35: B35: C35: D35: A36: 836: C36: D36: TIME STEP A37: 837: A38: 838: 159 U [w12] 'COPY\C U [w12] 'NEXT TIME STEP [w12] “thEW U 0 [w12] “QlNEW [w12] l.486/$J$4*$F$4“0.5*$B$33“(5/3)*$N$4 U [w12] 6 U [w12] ’/C$B$34..$B$36~$8$28..$B$30~/C$P$4~$O$4~ [w12] “hZNEW U 0.0015953261 [w12] “QZNEW [W12] (1.486/$K$4)*$F$4“0.5*$8$34“(5/3)*($N$4+$N$5)/2 [w12] “h3NEW U 0.0034417683 [w12] “QBNEW [W12] 1.486/$K$5*$G$5“0.5*$B$35“(5/3)*($N$5+$N$6)/2 [w12] “h4NEW U 0.00667484 [w12] ‘Q4NEW [w12] (1.486/$K$6)*$G$6“0.5*$B$36‘(5/3)*$N$6 IN SECONDS AND TIME WEIGHTING COEFFICIENT [w12] 'DELTA 60 [w12] 'THETA 0.5 DIRECT STIFFNESS METHOD SUMMING THE ELEMENTAL MATRIX ENTRIES TOGETHER WITH BOUNDARY CONDITIONS INTO THE GLOBAL MATRIX G38: 638: C39: D39: E39: F39: 639: C40: D40: E40: 640: C41: D41: E41: 641: C42: D42: E42: 642: U [w12] “'{F} [w12] “'[F] [w12] \- [w12] \- [w12] \- [w12] \- [w12] \- [w12] +812+A15 [w12] +315 [w12] 0 [w12] +012+015 [w12] +A16 [w12] +816+Al9 [W12] +819 [w12] +016+Dl9 [w12] 0 [W12] +A20 [W12] +820 [w12] +020 CCCCCCCCCCCCCCCCCC 160 t INVERSE OF THE GLOBAL [C] MATRIX MULTIPLIED TIMES THE {F*} VECTOR FOR THE SOLUTION OF THE NODAL FLOW DEPTHS h2,h3,h4 C45: U [w12] 0.0000279141 D45: U [w12] -0.0000096409 E45: U [w12] 0.0000048205 C46: U [w12] -o.ooooo96409 D46: U [w12] 0.0000397219 E46: U [w12] ‘0.0000198609 C47: U [w12] 0.0000048205 D47: U [w12] -0.0000198609 E47: U [w12] 0.0000981917 The worksheet solution and methodology was tested for a single plane consisting of the same input data for element No. 1. The solution was compared to the analytic solution using the Method of Characteristics. The analytic solution states that prior to equilibrium the flow depth is h = i*t where, i = rainfall intensity, t = time since start of the rainfall excess. h z flow depth. The flow rate is q = mha where, m = (1.486/n)"'Sl/2 a = 5/3 for overland flow. The following Figure Al shows the comparison of the finite element analysis and the analytic solution for a single plane. The parameters used for this plane is the upper plane in the arbitrary grid finite element solution. 161 .mcmam mamcwm m co comma scam LOW cowusaom owuaamc< 0cm LamEmHm Ouwcwm 0:» HO commLmaeou nan-a 33.. 65.3.4 + 6.8 as: 80.0 80.0 80.0 80.0 P P - p b lb) P b - 89d 80:— 8a..— b P n n n n 80.4 . H4 muamwm IRA—ID 30.: a.- D 80.0 T 8" r. 80d 83.124 232-...- umnzfl mam mdozmm (W) mm m .ocon odocwm o scum mum: Soda new :Omuadom umuaaoc< can acoaodu Oumcmm any no commuoaeou ~< «gamma .94: Ian. 9.54424 + .94: I0...- I-- D 6.3 3:... Acaduosonhv 804 03.4 8&4 084 80.0 80.0 03.0 8H0 08.0 P p b p P p p . p b — . p n r n 1 goo I 000.0 I 0u0.0 162. I 80.0 (MK) IL“ A011 Sufi-<24 Ella-AH 953.5 ijm mdwzmm 163 Conputer solution of the hydrologic response- area finite element grid was performed utilizing the following listing of FORTRAN code. 164 WWW OOOOOOOOOOOOOOOOOOOOOOOOODOOOOOOOOOOO 00 00 O on PROGRAM MAINSTRM This program finds a linear solution to a finite element overland flow problem for a global system of two-dimensional elements. Units of measurement are taken from the English System To avoid internal storage of large amounts of data. this program reads element and node information from external files as needed. In the variable listing below, the X-direction is identical to the direction of flow at a given node. DELTAT 8 Length (seconds) of time step. GCHOLD( ) = GCM( . ) * GHOLD(2) GCM( , ) = Global Capacitance Matrix f T j'bN'z’: ‘hN‘b GF( ) = Global Force vector GHNEWU) = Flow depth H (feet) at node I at end of current time step. GHOLD(I) = Flow depth H (feet) at node I at end of previous time step. GHTRYU) = Flow depth H (feet) at node I tried for solution iteration. When solution GHNEW is close to GHTRY, solution is accepted. GQRU) = Resultant of flow rate per unit width (ftflZ/sec) in X and Y directions at node I. GOX(|) = Weighted flow rate per unit width (ftHZ/sec) in X-direction at node I for given time step. GOXNEWH) 3 Flow rate per unit width (ftHZ/sec) in X-direction at node l at end of current time step. GOXOLDll) 8 Flow rate per unit width (ftflZ/sec) in X-direction at node I at tend of previous time step. GRF( ) = GF 1' 'hweighted rainfall‘h GSMX( ) = Global Stiffness Matrix for X-direction f T j'hN'lz 'lsz/dx'b GSMQX( ) = GSMX * GOX HLOLIM = Negative tolerance for flow depth H (feet). HUPLIM = Largest expected flow depth H (feet). INDEXM ) = Index for array of node numbers. NODE( ) = Array of node numbers. NP = Total number of nodal pOints being evaluated. NSOLV = Solution iteration number. MANNGN = Mannings N. MAXNOD = Maximum number of nodes assocuted With an element in the system. MAXSOL = Maximum number of solution iterations expected. Pl = 3.14159... ONEW = Flow rate per unit width at a node calculated from GHTRYU). 000000 0000000 10 20 70 60 72 165 RAINN Rate of rainfall at end of present time step (ft/sec). RAINO Rate of rainfall at end of previous time step (ft/sec). SLOPE Slope (ft/ft) at a node. SLOPEX == Slope (ft/ft) in X direction at a node. THETA Weighting coefficient. TIMEN Clock time (seconds) at end of present time step. TIMEO Clock time (seconds) at end of previous time step. Tracing definitions: ILEVEL 8 Tracing level for downstream write statements. ITRACK = Unit (TRACEOUT) to which results of intermediate calculations are written. LEVELT = Tracing level relative to the calling program. LTRACE = Tracing level for entire program. COMMON/TRACE/LTRACE.ITRACK COMMON/MATRIXIGSMX(20. 20).GSMY( 20. 20).GCM(20. 20).GF(20) COMMON/NODES/NP.NBW,NODE(20).INDEXN(20) DIMENSION GHNEW(20).GHOLD(20).GRF(20).GCHOLD(20).GSUM(20). + GQXNEW(20),GOXOLD(20).GQX(20).G$MQX(20).GHTRY(20). + GQYNEW(20).GQYOLD(20).GQY(20).GSMOY(20).GORES(20) INTEGER HOUR REAL MANNGN CHARACTERNIO TITLE DOUBLE PRECISION ASPECT,ASPRAD.PI DATA THETALSI. MANNGN/.035/. PI/3.I41592653589794/ DATA MAXNOD/4/ LEVELT/O/ HUPLIM/IOJ HLOLlMl-O.1E-06/ DATA GHNEW.GOXNEW.GOYNEWI60*O.I ITRACK = 8 OPEN(1,F|LE ='SY$TEM.IN') OPEN(3.F|LE='STORM.IN') OPEN(7,FILE='STREAM.OUT') OPEN(ITRACK.FILE='TRACE.OUT') READ(‘|,10) TITLE FORMAT(4OX.4OA) WRITE(7,20) TITLE WRITE(ITRACK.20) TITLE FORMAT(/.5X.40A./) ILEVEL = LEVELT+1 WRITE(6.70) FORMATI1X.'Enter the program tracing (debug) level: ‘) REACH”) LTRACE WRITEIITRACK.60) LTRACE FORMAT(/.10X.'The program tracing level is 212) WRITE(6.72) FORMAT(5X.'Executing MAINSTRM: Running main program.') Determine the total number of elements and nodes in the global system: CALL COUNT(ILEVEL.NELTOT.NODTOT) 166 C Determine which elements and nodes and their coordinates C are associated with the subsystem selected for analysis: CALL SUBSYS(ILEVEL.NELTOT.NODTOT) C Construct the global stiffness and capacitance matrices and the C force vector: CALL BUILDGULEVELNELTOT) C Modify the global matrices: Call MODGMflLEVEL) READ(3.10) TITLE WRITE(ITRACK,20) TITLE WRITE(7.20) TITLE WRITE(7. 722) (NODE(I),I= 1,NP) 722 FORMAT(6X. IOOI 10) ILEV1 =ILEVEL+ 1 lFllLEV1.LE.LTRACE) WRITEUTRACK, I O 2) 102 FORMAT(5X.'Executing time 'oop...') MAXSOL = NP510 Wiimw BEGIN TIME DEPENDENT CALCULATIONS 44444»qu C Read the rainfall at the beginning of the time period: READ(3. 2 3 7) HOUR.MINUTE,RAINN 237 FORMATGZ. 1X.l2.F 10.5) TSTART = (HOUR*SO.+MINUTE)*60 TIMEN ‘4 TSTART DELTAT=O WRITE(6,235) HOUR,M|NUTE.RAINN,DELTAT IF(|LEV1 + 2.LE.LTRACE) THEN WRITEllTRACK.'(1X)') WRITEIITRACK,235) HOUR,MINUTE,RAINN.DELTAT 235 FORMAT(5X.'At ',12,':'.|2.' rain intensity is ‘. 1 EIO.3.' ftlsec DELTAT s ‘, FS.1,' sec.') ENDIF WRITE(7.720) HOUR,MINUTE.(GORE$(I).I = LNP) WW START TIME LOOP annemeauuenmeuuunn 1000 DO 2000 ITlME=1,1000 TIMEO TIMEN RAINO RAINN READ(3.237,ERR=3000) HOUR,MINUTE,RAINN TIMEN = (HOUR'60+MINUTE)*60 DELTAT = TIMEN - TIMEO IF(DELTAT.LT.O) THEN WRITE(6.236) 236 FORMAT(1X,'TIME STEPS NOT CONSECUTIVE. PROGRAM STOPPED.') STOP 1030 497 1 6 7 ENDIF WRITE(6.235) HOUR,MINUTE.RAINN.DELTAT IFIILEVI + 2.LE.LTRACE) THEN WRITEilTRACK,‘(1X)‘) WRITEIITRACK.235) HOUR.MINUTE,RAINN.DELTAT ENDIF Assign the last computed values to the Cold‘ locations: DO 10301 = 1,NP GHOLD(I) 3 GHNEWII) GQXOLDU) = GQXNEWU) GQYOLDU) = GOYNEWll) CONTINUE Multiply the global capacitance matrix times the old depth values: DO 497 l=1, NP GCHOLDII) = 0.0 DO 497 K=1,NP GCHOLD(I) = GCHOLD(I) + GCMII.K)*GHOLD(K) CONTINUE NSOLV= 0 a H *4! H *I N START SOLUTION LOOP in! H a .9, H u 300 301 600 610 700 704 WRITE(6.'( 1 X)’) ILEV4=ILEV1 +3 NSOLV=N$OLV+ 1 IF(ILEV4.LE.LTRACE) THEN WRITE(|TRACK.'( 1X)’) WRITE(|TRACK.3O 1) NSOLV ENDIF WRITE(6,301) NSOLV FORMAT('+‘,4X.'Executing solution iteration number',l4) Compute the flow in each direction: IFILE= 1 DO 700 l=1,NP GHTRYII) 8 GHNEWII) NSKIP=NELTOT+NODE(I)+ 1 CALL SKIPIIFILENSKIP) READ(IFILE.610) SLOPE, WIDTH FORMATIZSX,2F 10.3) GOXNEWU) = 1.486/MANNGN 4* SLOPEflOS * GHTRY(I)*I(5./3.) lF(WlDTH.GT.0) GOXNEWU) = GQXNEWO) It WIDTH GQX(I) = (1.-THETA)*GQXOLD(I) + THETA .. GOXNEWU) GQYII) = (1.-THETA)*GQYOLD(I) + THETA 1* GOYNEWII) GRHI) = (1.-THETA)*GF(I)*RAINO + THETA if GF(l) ll! RAINN CONTINUE Write intermediate results: IF(ILEV4+ 1.LE.LTRACE) THEN WRITEIITRACK,704) FORMATUSX,’ NODE GHTRY GOXNEW GOYNEW GOX‘. 706 C 498 C 499 800 900 904 906 1400 1440 168 I ' GOY GRF') WRITE(|TRACK.706) (NODE(I).GHNEW(I).GOXNEW(I).GOYNEW(|). 1 GQX(|).GOY(I).GRF(I).I = ‘l.NP) FORMATISX.I5.6E 1 2.4) ENDIF Multiply the global X-stiffness matrix times the X-flow vector: D0 498 l=1, NP GSMOXII) 8 0.0 DO 498 K=1,NP — GSMOXO) 8 GSMQX(I) + GSMX(I.K) 1' GQXIK) CONTINUE Multiply the global Y-stiffness matrix times the Y-flow vector: DO 499 Isl. NP GSMOYII) = 0.0 00 499 K81,NP GSMQYU) : GSMOYU) + GSMYII,K) * GOY(K) CONTINUE Compute the global SUM vector: DO 900 l=1.NP GSUM(I) = GCHOLD(I) + DELTAT 4* ( GRF(I)- 1 (GSMQXII) + GSMOY(I)) ) CONTINUE Write intermediate results: lF(lLEV4.LE.LTRACE) THEN WRITEIITRACK.904) FORMATf/SX.’ NODE CHOLD GSMQX GSMOY GSUM‘) WRITE(ITRACK.906) (NODE(I),GCHOLD(I).GSMOXIDCSMOYII). 1 GSUM(|).|= 1,NP) FORMATISX.15.4E1 2.4) ENDIF Solve the system of equations 'bGCM'hit‘laGHNEW‘lFVzGSUM'h N= 20 CALL SLVSYSIGCM.GHNEW.G$UM.NP,N) Check estimated depth against result FLAG = 0. D0 1400 l=1.NP IF(ABS(GHNEW(I)-GHTRY(I)).LT.1.E-08) GO TO 1400 FLAG = 1. CONTINUE If all differences are within tolerance. go to next time step. IF(FLAG.E0.0) GO TO 1997 IF(NSOLV.GE.MAXSOL) THEN WRITE(6.1440) MAXSOL FORMAT('+',5X.'Solution conSidered convergent after ',l4, + ' iterations.) GO TO 1930 1 69 ENDIF 1 C Compare solution at each node with upper and lower limits: DO 1500 l=1.NP |F(GHNEW(I).LT.0) THEN IFIGHNEWIIIITHLOLIM) THEN WRITEIITRACKJ I 2) 7 1 2 FORMATISX.'NODE DEPTH') WRITE(ITRACK.7 1 5)(NODE(J).GHNEW(J).J= 1.NP) 7 1 5 FORMAT(5X,|4.512.4) WRITEUTRACK. 1520) WRITE(6. 1520) 1520 FORMAT(5X.'SOLUTION FOR FLOW DEPTH lS NEGATIVE. ', + 'PROGRAM EXECUTION STOPPED.') STOP ELSE GHNEW(I)=0. WRITE(6.1521) NODE(I) 1521 FORMAT(5X.'SMALL NEGATIVE FLOW DEPTH AT NODE ',I4, 1 '. 0 DEPTH ASSUMED.'/) ENDIF ELSEIHGHNEWO).GT.HUPLIM) THEN WRITEIITRACKJ 1 2) WRITEIITRACK, 7 1 5)(NODE(J).GHNEW(J).J= 1,NP) WRITE(6. 1530) HUPLIM.NODE(I) WRITE(ITRACK. 1530) HUPLIM.NODE(I) 1530 FORMATIIOX'SOLUTION FOR FLOW DEPTH 1‘: ‘.F5.1. 1 'AT NODE '.l4.'. PROGRAM EXECUTION STOPPED.') STOP ENDIF 1500 CONTINUE C lterate solution again: GO TO 300 Hanan“ ENDSOLUTIONLOOP “anaemia“ 1997 CONTINUE 1999 IF(ILEV4.LE.LTRACE) WRITE(6.1920) NSOLV 1920 FORMAT(‘+'.5X,'Solution converges after ',I3.' iterations.‘) C Compute resultant flow for each direction: 1930 D0 1940 l=1.NP 1940 GORESII) 8 SORT(ABS(GOXNEW(I))**2.+ABS(GOYNEW(I))**2.) |F(|LEVEL+ 1.LE.LTRACE) THEN WRITE(|TRACK.702) 702 FORMATI/SX.'NODE DEPTH FLOWX FLOWY TOTAL FLOW') WRITE(|TRACK,705) (NODE(I).GHNEW(I). 1 GOXNEW(|).GOYNEW(I).GORES(|),l= 1,NP) 705 FORMAT(5X,I4,4E12.4) ENDIF WRITE(7,720) HOUR. MINUTE. (GORESII).I=I,NP) 720 . FORMATIBX,12.':',12.2X, 100E 10.3) 170 2 000 CONTINUE ‘ WW END TIME LOOP W 3 000 CONTINUE 3001 WRITE(6.3001) IF(|LEVEL.LE.LTRACE) WRITE(ITRACK.300 1) FORMAT(l/.1X,'END OF PROGRAM REACHED') STOP END ”WWW 0000000000000 0 0000000 SUBROUTINE BUILDG(LEVELT,NELTOT) This subroutine reads the nodal coordinates associated with each element. rotates the element so that the local X axis lies in the direction of the element aspect angle. and sends the transformed coordinates to MATRIC for computation of element matrices. ASSMBL is called to place the coefficients of the element matrices into their respective positions in the global matrices. Subroutine arguments LEVELT and NELTOT are sent TO this subroutine. ASPECT = The aspect angle (degrees) for an element ASPRAD 8 The aspect angle (radians) for an element. GLOBALI . ) = Global coordinates associated with an element GTRANI . ) 3 Transformed global coordinates. LABELEI ) = Array of element numberals associated with subsystem. MAXNOD = Maximum number of nodes which is associated With any element in the subsystem. NELEMS = Number of elements associated with the subsystem. NELTOT = Total number of elements in the global system. NNDDELI ) 8 Array of nodes associated with an element. SLOPEI ) = Slope at a node. TRANSFI . ) = Transformation matrix. WIDTH(I) = Width of linear element at node I. COMMON/ELEMS/NELEMS.LABELE(20) COMMON/TRACE/LTRACE,ITRACK DIMENSION GLOBAL(4.2).NNODEL(4).W|DTH(4).$LOPE(4) DIMENSION TRANSF(2.2),GTRAN(4.2) DOUBLE PRECISION PI,ASPECT,ASPRAD DATA Pl/3. 141592653589 79400! DATA MAXNOD/4/ ILEVEL =LEVELT+ 1 WRITEIG. 10) IF(|LEVEL.LE.LTRACE) WRITE(|TRACK. 10) 10 11 100 120 130 140 160 520 540 560 714 716 718 720 171 FORMAT(SX.'Executing BUILDG: Building global matrices!) WRITEI6.1 I) . FORMATIIX) DO 600 INDEXE= 1.NELEMS IFILE=1 NSKIP=LABELEIINDEXE) CALL SKIPIIFILENSKIP) READ(IFILE. 1 20) NELEM.(NNODEL(J).J= 1.MAXNOD).ASPECT FORMAT(5I5.F 1 0.2) WRITEI6.130) NELEM IF(|LEVEL+2.LE.LTRACE) WRITEIITRACKJSO) NELEM FORMAT(‘+',4X,'Computing matrices for element number ',13) Determine how many nodes are associated with the element: NODSEL = MAXNOD D0 160 J=1.MAXNOD _ IF(NNODELIJIEQO) NODSEL=NODSEL- 1 Determine the element matrices: DO 560 I=1,NODSEL NSKIP = NELTOT+NNODEL(|)+1 CALL SKIPIIFILENSKIP) READ(1.540) NODNUM,(GLOBAL(|.J).J= 1.2).SLOPE(|).WIDTH(|) FORMAT(|5.4F 10.2) CONTINUE IF(NODSELEQZ) THEN CALL LINMAT(ILEVEL,NODSEL.NNODEL.GLOBALWIDTH) GO TO 599 ENDIF ASPRAD =Pl/ 180.DO*ASPECT Compute the coefficents of the transformation matrix: TRANSFI 1. 1)= DSINIASPRAD) TRANSFI 1.2)= DCOSlASPRAD) IF(ABSIASPECT).EO. 90.) TRANSH 1. 2)= 0. TRANSF(2. 1)- -TRANSF(1.2) TRAN$F(2.2)= TRANSFIIJ) IF(|LEVEL+2.LE.LTRACE) WRITE(ITRACK.714) ASPECT,ASPRAD FORMAT(SX.'Element aspect angle is ',F6.2. ' degrees or ',F7.4,‘ radians.‘) IF(|LEVEL+4.LE.LTRACE) THEN WRITE(ITRACK,716) FORMAT(ISX.'Coordinate transformation matrix:') WRITEIITRACK, 7 18)((TRAN$F(I.J).J= 1,2),I= 1,2) FORMAT(SX.ZE 10.3) ENDIF Transform the global coordinates: DO 730 l=1.NODSEL DO 726 J=1.2 726 730 510 542 550 599 600 660 172 GT RANII.J)=0. DO 726 K=1,2 I GTRANII.J)=GTRAN(I.J)+TRANSFIJ.K)OGLOBAL(I.K) CONTINUE ' CONTINUE lFIILEVEL+3.LE.LTRACE) THEN WRITEIITRACK,5 I 0) FORMAT(5X.'The following coordinates are associated with ', 1 'this element'./.5X.' NODE XGLOBAL YGLOBAL'. 2 ' Slope Xlocal Ylocal') DO 542 l=1.NODSEL WRITEIITRACK.550) NNODEL(I).(GLOBALII.J).J= 1,2).SLOPEII). 1 (GTRAN(|.J).J= 1. 2) FORMAT(SX.|5,5F10.2) WRITE(ITRACK,'( 1 X)') ENDIF Determine the element matrices: CALL MATRIC(ILEVEL,NODSEL.NNODEL.GTRAN) Add the element matrices to the global matrices: CALL ASSMBLIILEVELNODSEL,INDEXE.NELEM$,NELEM.NNODEL) CONTINUE End element loop. WRITEI6.660) IF(|LEVEL+ 1.LE.LTRACE) WRITE(ITRACK.660) FORMAT(‘+',4X.'Returning from BUILDG: ' 1 'Global matrices complete.',/) RETURN END mama" 00 00 000 SUBROUTINE COUNT(LEVELT,NELEMSNODSYS) This subroutine determines the total number of nodes and elements in the major global system data file. This permits subsequent reading from the correct lines of input data. The subroutine arguments NELEMS and NODSYS are passed FROM the subroutine. NELEMS = The total number of elements in the major global system. NODSYS = The total number of nodes in the major global system. COMMON/TRACE/LTRACE,ITRACK DATA IFILE/ 1/ ILEVEL =LEVELT+ 1 IF(|LEVEL.LE.LTRACE) THEN 1 7 3 WRITE(|TRACK. 10) WRITEI6.10) ‘ 10 FORMAT(SX.'Executing COUNT: Counting nodes and elements in '. 1 'global system...) ENDIF C Skip the record containing input column headings: NSK|P=1 CALL SKIP(IFILE. 1) 100 D0 200 IELEM=1,100 READ(1.1 20.ERR=220) NELEM 1 20 FORMAT(IS) NELEMS = IELEM 200 CONTINUE 220 NSKIP=NELEMS+2 CALL SKIP(|FILE,NSKIP) DO 400 lNODE=1,100 READ( 1. 1 20.ERR=420) NODNUM NODSYS = INODE 400 CONTINUE 420 IF(|LEVEL+2.LE.LTRACE) THEN WRITE(ITRACK.80) NELEMS.NODSYS 80 FORMAT(ISX.'TotaI number of elements in global system= ‘.l3. 1 /5X,‘Total number of nodes in global system= '.l3./) ENDIF IF(|LEVEL+1.LE.LTRACE) THEN WRITE(6.510) WRITE(|TRACK.510) 510 FORMAT(SX,‘Returning from COUNT: Counting complete.',/) ENDIF RETURN END WWW ”WWW PROGRAM MAINSTRM This program finds a linear solution to a finite element overland flow problem for a global system of two-dimensional elements. Units of measurement are taken from the English System. To avoid internal storage of large amounts of data. this program reads element and node information from external files as needed. In the variable listing below, the X-direction is identical to the direction of flow at a given node. 00 00 0 00 0 nonnnn nnonnnnnnnnnnnonnnnnOOOOOOOOOOOOODODODDDDDOO 174 DELTAT = Length (seconds) of time step. GCHOLD( ) = GCM( . ) It GHOLD(2) l GCM( . ) 8 Global Capacitance Matrix f T j‘lzN'b 'laN‘lz GF( ) = Global Force vector GHNEW(I) = Flow depth H (feet) at node I at end of current time step. GHOLD(I) 8 Flow depth H (feet) at node I at end of previous time step. GHTRY(I) 3 Flow depth H (feet) at node I tried for solution iteration. When solution GHNEW is close to GHTRY, solution is accepted. GQR(I) = Resultant of flow rate per unit width (ftflZ/sec) in X and Y directions at node I. GOX“) = Weighted flow rate per unit width (ftH2/sec) in X-direction at node I for given time step. GQXNEWO) = Flow rate per unit width (ftflZ/sec) in X-direction at node I at end of current time step. GOXOLDO) = Flow rate per unit Width (ftHZ/sec) in X-direction at node I at tend of previous time step. GRF( ) = GF it 'bweighted rainfall'lz GSMX( ) = Global Stiffness Matrix for X-direction f T j'bN'b VadN/dx'b GSMOX( ) = GSMX * GOX HLOLIM = Negative tolerance for flow depth H (feet). HUPLIM Largest expected flow depth H (feet). INDEXN( ) 3 Index for array of node numbers. NODE( ) 8 Array of node numbers. NP = Total number of nodal points being evaluated. NSOLV = Solution iteration number. MANNGN = Mannings N. MAXNOD = Maximum number of nodes associated with an element in the system. MAXSOL = Maximum number of solution iterations expected. Pl = 3.14159... ONEW = Flow rate per unit width at a node calculated from GHTRY(I). RAINN a Rate Of rainfall at end of present time step (ftlsec). RAINO = Rate of rainfall at end of previous time step (ftlsec). SLOPE = Slope (ft/ft) at a node. SLOPEX = Slope (ft/ft) in X direction at a node. THETA 8 Weighting coefficient. TIMEN = Clock time (seconds) at end of present time step. TIMEO = Clock time (seconds) at end of previous time step. Tracing definitions: ILEVEL = Tracing level for downstream write statements. ITRACK = File to which results of intermediate calculations are written. LEVELT = Tracing level relative to the calling program. LTRACE = Tracing level for entire program. WNW”§Q§”§§H§ 1 7 5 SUBROUTINE MATRIC(LEVELT.NODSEL.NSEND.GLOBAL) t This subroutine calculates the stiffness and capacitance matrices and the force vector for a triangular or 'quadrilateral element havmg a node at each corner. Subroutine arguments LEVELT.NODSEL,NSEND,XGLOBL.YGLOBL are sent TO this subroutine. Coefficients of KX,KY,C,F are sent FROM this subroutine via COMMON/ELMATS. B( . ) 8 Matrix containing the derivatives of shape function N with respect to the global coordinate system an. ) a ‘lsz/dx'h 8(2. ) = 'th/dyl‘: C( . ) =- Capacitance matrix for the element = f T j‘lzN‘l: 'bN‘b DNDNAT(I.J) = Array of derivatives of N with respect to KSI (1:1) and ETA (l=2) in Cnatural' coordinate system. F( ) = Force vector for the element GLOBAL( .J) = Array of X (J=1) and Y (J=2) coordinates for an element. ILEVEL = Tracing level used for debugging this subroutine. IJACOB(I,J) = Inverse of JACOBian matrix. JACOB(I.J) JACOBian matrix of partial derivatives of global coordinates with respect to natural coordinates ETA and KSI. Kx( . ) =- Stiffness matrix with respect to rotated local direction X. = f T j'hNVz Vsz/dx'lz KY( . ) a u 00 N I U N I N Y LEVELT = Tracing level sent from upstream program NODSEL = Number of nodes associated with the element. NSEND( ) = Node numerals assocmted with the element VN(l) = Shape function N at node I in Cnatural' coordinate system for given integration point. Dhatt, G. and G. Touzot The finite element method displayed John Wiley and Sons. NY, 1984. 509pp. (See pp. 45. 102.) no on nonnnnnnnnnnnnnnnnnnnnnn 0 non 000 Segerlind. L. Applied finite element analysis (2nd ed.). John Wiley and Sons. NY, 1984. 427pp. (See pp. 73. 365. 372-374.) COMMON/ELMATS/KX(4,4).KY(4,4).C(4,4).F(4) COMMON/TRACE/LTRACE.ITRACK DIMENSION VN(4).DNDNAT(2, 4).NSEND(4).GLOBAL(4. 2).B( 2. 4) REAL KX,KY, JACOB(2,2).IJACOB(2,2) ILEVEL=LEVELT+ 1 IF(|LEVEL.LE.LTRACE) WRITE(ITRACK. 1 2) 12 FORMAT(SX.'Executing MATRIC: Calculating elemental matrices.) 16 499 498 315 316 176 IF(|LEVEL+2.LE.LTRACE) WRITE(ITRACI<.16) NODSEL FORMAT(SX.'There are ',I2,‘ nodes associated with this element) IF(|LEVEL+3.LE.LTRACE) THEN WRITE(ITRACK. 1 7) FORMAT(I 1 2X.'NODE'.8X,'X'. 14X,'Y') WRITE(|TRACK. 18) (NSEND(I).(CSLOBAL(I.J).J= 1.2).l8 1.NODSEL) FORMAT(4(10X.15.2E15.6./)) ENDIF Initialize the element matrices: D0 200 I81.NODSEL DO 190 J81,NODSEL C(I.J) 8 0. KX(I.J) 8 0. KY(I.J) 8 0. CONTINUE F0) 8 0. CONTINUE IF(NODSELEQ 4) INTPTS IF(NODSEL.EO. 3) INTPTS II II ‘ ILEV58ILEVEL+5 DO 2000 INTPT 8 1,|NTPTS IF(NODSELEQ4) CALL QSHAPE(ILEV5,INTPT,VN,DNDNAT) IF(NODSEL.EO.3) CALL TSHAPE(VN.DNDNAT) Calculate the Jacobian matrix: JACOB=DNDNAT*GLOBAL DO 499 l=1.2 DO 499 J=1.2 JACOB(I.J) 8 0.0 DO 499 K=I.NODSEL JACOB(|.J) 8 JACOB(I.J) + DNDNATII,K) 8 GLOBAL(K,J) CONTINUE Determine the inverse of the Jacobian matrix: CALL INV2X2(JACOB,IJACOB.DETJAC) Multiply : lJACOB if DNDNAT 8 8 DO 498 181. 2 DO 498 J81. NODSEL B(l.J) = 0.0 DO 498 K81.2 B(l.J) = B(l.J) 4’ IJACOB(I.K)*DNDNAT(K.J) CONTINUE IF(|LEVEL + 5.LE.LTRACE) THEN WRITE(ITRACK,3 1 5) FORMAT(/5X."b ----- JACOBIAN ----- v: )‘2- - - -JACOBIANinv- - - m 00 3 16 i= 1,2 WRITE(ITRACK.318)(JACOB(I.J).J8 1.2).(IJACOB(I,J).J8 1.2) 177 3 18 FORMAT(SX.2E10.3.5X,2E10.3) WRITE(ITRACK.3 10) INTPT 310 FORMAT(ISX.'8 matrix at integration point‘.l2,' :) WRITE(ITRACK.31 1) (NSEND(J).J8 1.NODSEL) 31 1 FORMAT(SX.I5,5X.I5.5X,I5.5X,I5) ' DO 312 I=1,2 312 WRITE(ITRACK.314) (B(l.J).J81,NODSEL) 3 1 4 FORMAT(SX.8E10.3) ENDIF IF(NODSEL.EQ.3 .AND. ILEVEL+4.LE.LTRACE) THEN AREA8ABS(DETJAC/2.) WRITE(ITRACK.510) AREA 510 FORMAT(SX.’Area of triangle 8 '.EIO.4) ENDIF C Compute the element matrices: IF(NODSEL.EQ.4) WC 8 1.0 IF(NODSEL.EO.3) WC 8 0 5 SCALAR = wc » ABSIDETJAC) 700 D0 800 I=1,NODSEL 710 D0 790 J81.NODSEL C The capacitance matrix C is the sum over the integration C points of N'lztranspose * N multipled by a scalar value. C(I.J) 8 C(I,J) + VN(I) If VN(J) 1' SCALAR C The stiffness matrix K is determined for each direction as C the sum over the integration points of N'lztranspose If B C multiplied by a scalar value KX(I.J) 8 KX(I.J) + VN(I) It B(l.J) 9 SCALAR KY(I.J) 8 KY(I.J) + VN(I) * B(2.J) 4' SCALAR 790 CONTINUE C The force vector F is the sum over the integration points C of N multiplied by a scalar value. F0) 8 F0) + VN(I) * SCALAR 800 CONTINUE 2000 CONTINUE C Write the element matrices: 2010 lFIlLEVEL+4.LE.LTRACE) CALL WRITEM(4,NODSEL,NSEND,KX,KY,C,F) IF(|LEVEL+ 1.LE.LTRACE) WRITE(ITRACK.2 1 40) 178 2140 FORMAT(5X.'Returning from MATRIC: Elemental matrix complete.'./) RETURN ‘ END "MWWWHWWWWWM SUBROUTINE OSHAPE(LEVELT.INTPT,VN.DNDNAT) For a quadrilateral having a node at each corner. this subroutine determines the shape functions and the derivatives of the shape function with respect to the natural (KSI-ETA) coordinate system. Subroutine argument INTPT is passed TO this subroutine. while subroutine arguments VN and DNDNAT are passed FROM this subroutine. INTPT 8 Number associated with the integration point in the KSI-ETA coordinate system. VKSIN(I) 8 Value of KSI at node I in Cnatural' coordinate system. ( 8 Abscissa of node I " “ " " ) VETAN(I) 8 Value of ETA at node I “ “ " " ( 8 Ordinate of node I " " “ " ) XINTEG(I) 8 Abscissa of integration point in Cnatural‘ coord. system. YINTEG(I) 8 Ordinate of integration point in Cnatural' coord. system. 00 0000000 00 000 COMMON/TRACE/LTRACE,ITRACK DIMENSION VKSIN(4).VETAN(4).XINTEG(4).YINTEG(4).VN(4).DNDNAT(2,4) DATA VKSIN/-1.0. 1.0. 1.0, -1.0/. VETAN/-1.0. -1.0. 1.0. 1.0/ Calculate the location of the integration points in the Cnatural' KSI-ETA coordinate system for the linear quadrilateral. 000 ILEVEL 8LEVELT+ 1 IF(|LEVEL.LE.LTRACE) WRITE(ITRACK. 10) 10 FORMAT(SX.'Executing OSHAPE: Determining shape functions for ‘. 1 ‘an integration point.) WHERE 8 .577350 00 100 l=1.4 XINTEG(I) 8 WHERE it VKSIN(I) YINTEG(I) 8 WHERE * VETAN(I) 100 CONTINUE C At each integration point : DO 450 181,4 C Compute the shape function N at each node I: VN(I) 8 0.25 * (1.+VKSIN(I) If XINTEG(INTPT)) + * (1.+VETAN(I) It YINTEG(INTPTl) C Compute the derivatives of N wrt natural coordinates: 179 DNDNAT(1.I) = 0.25 1: VKSIN(I) a (1+VETAN(I) a XINTEG(INTPT)) DNDNAT(2.!) = 0.25 i VETAN(I) a Ii+Vi60 - a Mooemrsl SANXTARY racerrres cowsrnucvxou HAYERIAL i-bx: sLxCMt l-lsx: FAIR-LOU STRENGTH SEPTIC tau: s-xsx: MODERATE-SLOPE ns-zss: FAIR-SLOPE.LOU steamer. .esoarwxoa isos: sevens-Stove aoaorxLL zsos: soon-store FlELOS 1-2s: MODERATE-SEEPAGE IMPROBAILE-EXCESS FINES ssv.cs 2-7s: NOOERAYE-SEEPAGE.SUOFE LAGOON 1oz: sevgnc-SLovC SAND .ness l-bx: SL1GHT IMPROSABLE-EXCESS FINES ssuxtsnv e-xss: MODERATE-SLOPE LANDFILL isos: SEVERE-SLOPE GRAVEL (teeucn) ' l L-ox: sLxth I L-ax: 0000 I SANITARY s-nsx: MODERATE-SLOPE a-nsx: FAIR-SLOPE I LAMO‘ILL x5os: SEVERE-SLOPE TOPSOIL Lsox: POOR-SLOPE I Lane.) I 1 I l-dz: GOOD DAILY s-Lss: FAIR-SLOPE HATER HANAGEHENT coves roe isos: FOOR-SLOPE i-Jx: MODERATE-SEEPAGE I LAMostLL no.0 J-as: HOOERAte-SEEPAGE.SLO9E | nesenvoxn sot: SEVERE-SLOPE I AREA I EQILOING sxrc DEVELOPNENT l l-es: SLIGHT SEVERE-PIPING I SNALLOH a-iss: HOOERATE-SLOPE ENOANKHENTS l Excav.txc~s 15oz: SEVERE-SLOPE oxnes AMo l l LEVEES I i-as: SLICMt SEVERE-NO HATER l cutLLIMcs s-zss: MODERATE-SLOPE ExCavstEO I I vxtMout 15os: SEVERE-SLOPE acuos l aASCMC~ts AOUIFER can I i-sx: SLIGHI 0559 10 wsiER I oucLLxucs s-xss: HOOERATE-SLOPE I I HITM isos: stveac~5Loo£ I oaaxuncc l IASEHEutS : l-tl: SLIGHT n-Js SIL.L: enooes EASILV I s~ALL a-es: HOOEIAVE°SLOP£ so: st.L: SLOPE.EROO£S EASILV l I CONMfRCiAt .os: SEVERE-SLOPE II lflRIGAilON I-Js v:$L: canoes CASIIV.SOIt .Lowxuc I IIiIItnififi5 I II 3.: vrSL: SLOPE.tHOU(S tASltY.501t BlOHING I 195 'ropography, soils, and landuse were taken from "THE CENTRAL GREAT PLAINS EXPERIMENTAL WATERSHED, A Summary Report of 30 'Years of Hydrologic Research". The following map and data are excerpts. from this report and from Plane Table 'ropographic Surveys performed in 1942 by the Soil conservation Service-USDA under Hugh Hammond Bennet. The report is available from the Water Data Laboratory, Agriculutral Research Service-USDA, Beltsville Maryland. The folio of maps is available for loan from the same Laboratory. 196 WATERSHED l-H 3.62 acres , WATERSHED 3-H 3.95 acres ‘ £543.53“). 836R . rue (0:: 7 Area atter Ill/59. 3.77 ac. Dlstam and direction to neareet rain gage. Gaqinq station. Original watershed boundary. Watershed boundary after Ill/$9. Contoure. Sail boundariee. Numeratar it coil type; denamlnatar in range at top eat! depth.in tacbee. WATERSHED Z-H . 3.40 acres WATERSHED 4-H 3.84 acres . Area after tit/59. 3.64 ac. SCALE using-— 30 0 w m FEET Watersheds l-H,2-H, B-H, B 4-H. 197 Tabla 3 .—Crop and trantn-mt plan for 4-oara watersheds for the period of 1958 through 1967 Crop and treatment by your: y fl I watershed 1958 1959 1960 1961 1962 1963 1904 1965 1966 1967 1.n Mn Mn Mn Mn Mn ' Mn 3: Sm nn un 2-H Mn Mn Mn Mn Mn Mn ' P P P P 22.3 . . Mr M: M: Mr Mr Mr 215-}! M:- M: Mr Mr M: Mr 254! Mn Mn ' Mn Mn Mn 18-8 P P P P P P P P ' P P 341’ Wm Sm m Mn Sm Fm Hm - Sm I-‘m Hm B-E Wm Sm Rn H m Sm Fm Hm Sm fin Um- 4-3’ Sm Fm Wm Sm I-‘m Wm Sm Fm Hm Sm 7-3’ Sm Fm an Sm n: Hm Sm' m m; Sm 54: Pm Wm 82: Pm wn Sm Pm m: Sm an 6—8 m Hm Sm I'll Um Sm Eh 9m Sm F: 3.] Symbols used in columns are: [In . native meadow} Hr - meadow of reseeded native grasses; St . forage sorghum; S - sorghum in rows for grain; H . wheat; F - fallow; P . native pasture; m - stubble mulch farmed. SIIZCTID .0307? EVENTS Aflroceosur C0001T10fl Date 0 (1a.! ‘-10 .02 ‘-17 .J‘ ‘-1, 051 5-2 .25 S-J .1‘ watershed conditions: fellow good residue cover. 0 (1a.) .00 .00 .01 .00 .00 198 A25 IAIN'ALL Onto 71.. Rat. ten/ht! Event of HA1 5, 1959 RC D-JS-R S-‘ 1‘10 .00 1‘20 .90 1‘22 2.10 1‘2‘ 1.60 1‘26 S.‘0 1‘20 3.60 1‘10 2.‘0 1‘12 1.20 1507 .19 1527 .0) 1510 .00 1750 .06 Ace. (Ln.t .00 ‘ 0°, .22 .‘0 .52 .60 .6‘ .75 .70 .26 .I‘ HATEISHED O-U (40.001 21n- 1425 1‘27 1‘29 13‘! 1‘15 1‘39 1‘82 1“5 1‘52 120’ 1589 11" IUNO'? Iot- (In/hr! .3100 .2‘1 1.2) .997 0“: .295 .110 .0517 .0215 .0005 .0051 .0000 Ace. (tn.1 .0000 Do) .07 .11 .1, 01‘ .15 .15 .10 .15 .16 199‘ Regression analysis of Table 4 outflows. This analysis was performed with Lotus 1-2-3 for the actual and calculated outflows. The pairs of values are at constant time for which the analysis was performed. Regression Output: Constant 2.9OE-01 Std Err of Y Est 3.968—01 R Squared 9.268—01 No. of Observations 2.7OE+01 Degrees of Freedom 2.SOE+01 X Coefficient(s) 3.43E+Ol Std Err of Coef. 1.94E+00 "llllllll'llflfllllfllITS