WATER T A B L E D E P T H SIMULATION FOR FLAT AGRICULTURAL L A N D UNDER SUBSURFACE DRAINAGE A N D SUBIRRIGATION PRACTICES by E N A C.Y. C H A O B.ASc , The University of British Columbia, 1983 A THESIS SUBMITTED IN PARTIAL F U L F I L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF MASTER OF APPLIED SCIENCE in THE F A C U L T Y O F G R A D U A T E STUDIES DEPARTMENT OF BIO-RESOURCE ENGINEERING We accept this thesis as conforming to the required standard THE UNTVERSITY O F BRITISH COLUMBIA APRIL, 1987 e ENA C.Y. CHAO, 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of BIO-RESOURCE ENGINEERING The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date APRIL, 1987 ABSTRACT Drainable porosity as a function of water table was investigated to replace the common practice of treating it as a constant A continuous function in the form of a negative exponential equation relating drainable porosity to water table depth was developed by three methods: (1) laboratory core-sample analysis; (2) rainfall rate and water table depth analysis; (3) drainage rate and water table depth analysis. Furthermore, this function was derived for four different water table regimes: (1) subsurface drainage; (2) low subirrigation and subsurface drainage; (3) high subirrigation and subsurface drainage; (4) no drainage and no subirrigation. The drainable porosity function was incorporated into a water balance model which simulated the soil moisture profile and the water table depth on a daily basis. Major modification of the previous model was the elimination of separate falling and rising water table equations since discrete porosity values were no longer assigned to particular soil depth intervals. A subroutine program which computed the total maximum transient storage and the transient storages to each of the four successive soil zones was also incorporated. The 'maximum drainable porosity' and the 'rate constant' parameters in the negative exponential equation were found to be different among the three methods of analysis and among the four water table regimes. Good agreement between simulated and actual water table depths of each regime for 1984 and 1985 was found. The modified water balance model could be used to generate different water table depths by changing the input parameter of design drainage rate. From these outputs, a appropriate drainage rate which gives the desired water table depth could be selected for the purpose of horizontal subsurface drainage system design. ii TABLE OF CONTENTS ABSTRACT ii LIST OF TABLES v LIST OF FIGURES vi LIST OF SYMBOLS A N D ABBREVIATIONS vii A C K N O W U 3 X 3 E M E N T x 1. INTRODUCTION 1 1.1. STATEMENT A N D N A T U R E O F P R O B L E M 1 1.2. OBJECTIVES 2 1.3. SCOPE 3 2. LITERATURE REVIEW 4 2.1. SUBSURFACE D R A I N A G E 4 2.1.1. Physical Properties of the Soil 6 2.1.2. Trafficability and Timeliness 8 2.1.3. Plant Growth and Crop Production 10 2.2. SUBIRRIGATION 12 2.3. CAPILLARY FLOW 13 2.4. DRAINABLE POROSITY 18 2.5. WATER B A L A N C E M O D E L 19 3. MATERIALS A N D METHODS 23 3.1. EXPERIMENTAL FIELD A N D MEASUREMENTS 23 32. PHYSICAL CHARACTERIZATION O F SOU 27 3.3. F O R M U L A T I O N OF D R A I N A B L E POROSITY EQUATION 29 3.3.1. In Situ Drainable Porosity 29 3.3.2. Curve Fitting 30 3.3.3. Corrected f-z Function 30 3.4. WATER B A L A N C E M O D E L 31 3.5. STATISTICAL ANALYSES 32 4. RESULTS A N D DISCUSSION 33 4.1. SOIL PHYSICAL ANALYSES 33 4.2. SOIL H Y D R A U L I C ANALYSES 34 4.3. LABORATORY DETERMINED I N SITU D R A I N A B L E POROSITY 37 4.4. IN SITU D R A I N A B L E POROSITY 38 4.4.1. Rainfall Rate Method 38 4.4.2. Drainage Rate Method 41 4.5. THE DRAINABLE POROSITY EQUATION 41 4.5.1. Curve Fitting 41 4.5.2. Corrected f-z function 46 4.6. T H E MODIFIED WATER B A L A N C E M O D E L 47 4.6.1. Assumptions to the Model 48 4.6.1.1. The Water Balance Model 48 i i i 4.6.1.2. The Water Table Model 50 4.6.1.3. The Computer Model 51 4.7. WATER TABLE SIMULATION 52 5. SUMMARY AND CONCLUSIONS 64 6. RECOMMENDATIONS 67 REFERENCES 68 APPENDIX 76 iv LIST OF TABLES TABLE I: Opportunity Days Gained From Drainage Improvement (Chieng et al., 1987) 10 TABLE II: Characteristic grain size distributions of sandy subsoils of various geo-genetical groups, in percentages of the mineral part of the soil (Bloemen, 1980) 17 TABLE III: Nonlinear Regression Equations (Bhattacharya and Broughton, 1979) 20 TABLE IV: Particle Size Distribution (USDA) 34 TABLE V: Textural Classification (USDA) 34 TABLE VI: Dry Bulk Density 35 TABLE VII: Satiated Hydraulic Conductivity 35 TABLE VIII: Regression Statistics to the Negative Exponential Curve Fitting 43 TABLE IX: Regression Statistics to the Polynomial Curve Fitting 43 TABLE X : Correlation Index Statistics to the Negative Exponential Curve Fitting 43 TABLE XI: Parameters to the Negative Exponenetial Drainable Porosity Equation 45 TABLE XII: Parameters to the Corrected Drainable Porosity Equation 47 TABLE XIII: Maximum Design Drainage Rate R M A X 62 v LIST OF FIGURES FIGURE 1: Subsurface Drainage System Layouts (Chieng, 1986) 5 FIGURE 2: Relationships between capillary rise and water table depth for sandy subsoils of various geo-genetical groups (Bloemen, 1980) ..15 F IGURE 3: Height of capillary rise in relation to soil moisture suction in marine clays at increasing clay contents (Bloemen, 1980) 16 F IGURE 4: Experimental Field Layout 24 F I G U R E 5: The Water Control System 25 FIGURE 6: Soil Moisture Retention Curves 36 F IGURE 7: Drainable Porosity Curves From Laboratory Analysis 36 F IGURE 8: Drainable Porosity Based on Rainfall Rate Method For Regime A 39 F IGURE 9: Drainable Porosity Based on Rainfall Rate Method For Regime B 39 F IGURE 10: Drainable Porosity Based on Rainfall Rate Method For Regime C 40 FIGURE 11: Drainable Porosity Based on Rainfall Rate Method For Regime D 40 FIGURE 12: Drainable Porosity Based on Drainage Rate Method For Regime A 42 F I G U R E 13: Drainable Porosity Based on Drainage Rate Method For Regime B .— 42 F I G U R E 14: Generalized Flowchart- of the Water Balance Model 53 F IGURE 15: Actual Water Table Curves For Regimes A, B, C, D of 1984 54 F IGURE 16: Actual Water Table Curves For Regimes A, B, C, D of 1985 54 F IGURE 17: Laboratory-Based Water Table Simulation For Regime A of 1984 55 F I G U R E 18: Laboratory-Based Water Table Simulation For Regime B of 1984 55 F IGURE 19: Laboratory-Based Water Table Simulation For Regime C of 1984 56 F IGURE 20: Laboratory-Based Water Table Simulation For Regime D of 1984 56 F I G U R E 21: Laboratory-Based Water Table Simulation For Regime A of 1985 57 F I G U R E 22: Laboratory-Based Water Table Simulation For Regime B of 1985 57 F IGURE 23: Laboratory-Based Water Table Simulation For Regime C of 1985 58 F I G U R E 24: Laboratory-Based Water Table Simulation For Regime D of 1985 58 F IGURE 25: Drainage-Rate-Based Water Table Simulation For Regime A of 1984 ..59 F IGURE 26: Drainage-Rate-Based Water Table Simulation For Regime B of 1984 ..59 F IGURE 27: Drainage-Rate-Based Water Table Simulation For Regime A of 1985 ..60 F IGURE 28: Drainage-Rate-Based Water Table Simulation For Regime B of 1985 ..60 F I G U R E 29: Simulated Water Table Heights For Regimes A , B, C of 1984 (Based on Laboratory Analysis) 63 F IGURE 30: Simulated Water Table Heights For Regimes A , B, C of 1985 (Based on Laboratory Analysis) 63 vi LIST OF SYMBOLS AND ABBREVIATIONS A area A E actual evapotranspiration b0,bi,bj,...,b6 constants to curve fitting equations C governing fact, 4K./L 2 C.V. coeffient of variance CI A E / P E ratio for first soil zone C2 A E / P E ratio for second soil zone C3 A E / P E ratio for third soil zone C4 A E / P E ratio for fourth soil zone D vertical distance from drain center to the impermeable soil layer D l or D x depth from the soil surface of the first zone D Hooghoudt's equivalent depth of soil, below the cirain center, through which flow to chains occurs df degree of freedom DP chainable porosity Drain. drainage rate DW depth of zone from soil surface that R M A X occurs E X excess water F or f drainable porosity f* f, . or f. . drain lab ^diain determined o V drainage rate method vii f^k f determined by laboratory method f • f determined by rainfall rate method ram f. drainable porosity of the ith zone fj drainable porosity of the first zone f2 drainable porosity of the second zone H or H M A X total drain depth P index of correlation h height of water table at midspacing, above subdrain center (I) counter for the current day K. hydraulic conductivity (K) counter for the previous day m number of independent variables n number of samples L horizontal spacing distance between subdrain centers LAB. laboratory PE potential evapotranspiration PPT or PRE precipitation R drainage rate (drainage coefficient), volume of outflow per unit area of land drained per unit time R* coefficienct of determination R E G regime R M A X maximum drainage rate or design drainage rate RO surface runoff S.D. standard deviation viii SEWJO sum of total daily water table height recorded above the depth of 30cm at midpoint between drains SMCF total available soil moisture content SMC1 available soil moisture content for first soil zone SMC2 available soil moisture content for second soil zone SMC3 available soil moisture content for third soil zone SMC4 available soil moisture content for fourth soil zone SSreg sum of squares of regression SSres sum of squares of residuals SSy total sum of squares Syx standard error TR transient storage TR1 transient storage for^ first soil zone TR2 transient storage for second soil zone TR3 transient storage for third soil zone TR4 transient storage for fourth soil zone T R M A X maximum total transient storage Y R year Z or z water table depth from soil surface A H change in water table height At change in time A(TR) or ATR change in transient storage A Z change in water table depth ix ACKNOWLEDGEMENT The author wishes to extend her sincere appreciation to Dr. S.T. Chieng for his guidance and advice during the course of this research, and for his valuable criticisms and comments during the preparation of this thesis. The author is also deeply grateful to Drs. S.O. Russell and J.L. Smith for reviewing this work. Special thanks go to Dr. J. DeVries and Mr. M . Bonsu for their assistance with the laboratory analyses, to Mr. N . Jackson and Mr. J. Pelke for their assistance with photography, and to all the graduate colleagues who gave the much needed encouragement. Above all, the author wishes to extend her love and thanks to her family, especially her sister and mother for their patience, understanding and support. This study was financially supported by a postgraduate scholarship from the Natural Science and Engineering Research Council of Canada awarded to the author. x 1. INTRODUCTION 1.1. S T A T E M E N T AND N A T U R E OF P R O B L E M The control of moisture level in a soil that is otherwise water-saturated during excess precipitation periods or water depleted during drought periods is agriculturally beneficial. The implementation of a dual water management system, which serves to drain or irrigate, can bring previously uncultivable waterlogged land into production or enhance the production of existing arable lands. In the Lower Fraser Valley of British Columbia, Canada, where land is relatively flat and the climate is such that alternate seasons of water surplus and deficit occur, a combined subsurface drainage and subirrigation system is well suited to controlling the water table (Goldstein, 1978). The predominant hydraulic condition exhibited by the local soils is moderately poor to poor drainage. It is estimated that an area of 40,000 to 50,000 ha requires some form of water management scheme (Driehuyzen, 1983). In arriving at a practical water management system, the design drainage rate or drainage coefficient (R) must be determined first The general design criterion is to have minimal water control system costs with maximal crop yield and positive soil conditioning. The drainage coefficient is a function of the soil's drainable porosity as well as the soil's water holding capacity, regional climate, water table depth requirements of particular crops and field machine operations. In order to attain a realistic drainage coefficient the drainable porosity itself must be properly described. The common definition of drainable porosity as the volume of water that can be drained from a unit of soil by a certain water table drawdown implies it is a constant independent of pressure head. Such constants are commonly assigned to 1 INTRODUCTION / 2 particular soil types or soil depths. This, however, is not correct Childs (1960) and Taylor (1960) have reported the dependence of drainable porosity upon water table depth and the rate of change of water table depth. The pressure gradient created by a falling water table and the time interval over which this event occurs affect the volume of water released from the soil profile. Thus a continuous relationship between drainable porosity and water table position is needed. This replaces the inadequate approach of using a constant value to represent drainable porosity. A function that expresses drainable porosity in terms of water table position can then be applied in a water balance model to derive an appropriate drainage coefficient 1.2. OBJECTIVES The main objectives of this research project are: 1. to derive a functional relationship between drainable porosity (f) and water table depth (z) for the local region by: a. laboratory analysis on disturbed and undisturbed soil samples; b. analysis of in situ rainfall rate and water table depth data; c. analysis of in situ drainage rate and water table depth data; 2. to incorporate the derived f-z relationship into an existing water balance model that simulates soil moisture profile and water table movements for flat agricultural lands; 3. to verify the modified water balance model by using recorded field data; 4. to use the verified model to generate water table depths for a given set of conditions. INTRODUCTION / 3 1.3. SCOPE The scope of this thesis is restricted to flat agricultural land whose water regime can be managed by a horizontal subsurface drainage system. This requires relatively uniform soil properties over an area of a few hectares and seasonal drainage characteristics of alternating high and low water table levels during the year. The refined water balance model will take the local conditions into account when the local weather data are used. 2. LITERATURE REVIEW 2.1. SUBSURFACE DRAINAGE Proper drainage of agricultural land is of major importance to increasing the efficiency and reliability of crop production. Subsurface drainage is a means of implementing this needed water management It is a system of artificial channels existing beneath the soil surface at a specifc depth which collects and conveys excess water away from the root zone. The system consists of buried perforated plastic tubings or clay tiles commonly arranged in three possible layouts: random, herringbone, or gridiron as shown in Figure 1 (Donnan and Schwab, 1974). The general topographic and soil features govern the appropriate system to be installed (Chieng, 1986): 1. Parallel pattern: A series of drains placed parallel to each other and usually spaced at regular intervals. 2. Herringbone pattern: This arrangement consists of one main line and a series of laterals which enter the main line at an angle. This system is most suited to drainage of relatively narrow and sloping depressions or valleys where the main line is placed in the lowest areas. 3. Random pattern: This system is generally used where the topography is undulating or soils vary and the field contains isolated wet areas. Drains are laid out in a manner to connect relatively small, isolated, poorly drained depressions in the most economic and effective way. Subsurface drainage serves primarily two functions: (a) to improve the trafficability of the land for field operations such as seed bed preparation, planting, and harvesting, (b) to provide a growth-maintaining environment within the root zone during high rainfall events. Subsurface drainage therefore modifies certain properties of 4 FIGURE 1: Subsurface Drainage.System Layouts (Chieng, 1986) 3 LITERATURE REVIEW / 6 the soil matrix and processes within the soil-water-atmosphere regime. Three areas of particular significance are: physical properties of the soil; trafficability and timeliness; and plant crop production. 2.1.1. Physical Properties of the Soil Hundel et al. (1976) reported a number of measurable soil properties that are affected by drainage: drainable porosity and pore size distribution, hydraulic conductivity, crust density and penetration, and compressible strength. The combined analyses of these properties indicate the soil's behaviour or vulnerability to external forces — structural and aggregate stability. The porosity is an index of the relative pore volume in a soil profile. Both total porosity and pore size distribution characterize the process of aggregation. Soils with a high volume of air-filled pores are indicative of good soil structure. Hundel et al. (1976) found subsurface drainage to improve soil structure by increasing the density of large and medium sized pores. The effectiveness of both air and water pathways is enhanced due to the increased cohesiveness of soil aggregates resulting from a drying soil structure (Hillel, 1980). Closely related to the porosity is the soil's hydraulic conductivity. By definition, the hydraulic conductivity is the ratio of the specific discharge to the hydraulic gradient This soil property is characteristically constant for saturated soil of stable structure or for rigid porous medium; the hydraulic conductivity is 10"2 to 10"3cm/sec for sandy soils and 10"4 to 10"'cm/sec for clayey soils (Hillel, 1980). In addition to texture, soil structure and in particular, the number and geometry of the pores affect the hydraulic conductivity (Doering, 1965). It then follows that since subsurface drainage alters the pore distribution favourably, the hydraulic conductivity is correspondingly LITERATURE REVIEW / 7 enhanced. Leyton and Yadar (1960), and Hundel et al. (1960) have demonstrated this observation between subsurface drained and undrained fields. Soil properties which are important to the root development of seedlings are penetration resistance and compressive strength. Because individual root branches encounter mechanical obstructions within the soil matrix, it is beneficial to minimize these two factors. Fausey and Schwab (1969) found them to be less in subsurface drainage treatments than those with surface or no drainage treatments. The problem of crusting also arises in the latter two cases. Upon drying, these fields exhibit high crust density, wide crack spacings, and large, thick crust units (Hundel et al., 1976). The crust creates poor tillage for seed bed preparation. In addition to causing obstructive and distorted pathways to root penetration, crusting also damages seedling roots by tearing them apart Subsurface drainage reduces surface crust strength which is a function of the drying rate (Hillel, 1980). Two thermal properties that are affected by soil moisture variations are specific heat capacity and thermal conductivity. The specific heat capacity (C) is defined as the change in heat content of a unit mass per unit change in temperature. It is a function of the soil's composition; according to De Vries (1975), C is given by the arithmetic sum of the heat capacities of its various constituents. Therefore an increase in the soil moisture capacity causes an increase in the soil's heat capacity. Also with wet soils, more radiant energy is absorbed than radiated back to the atmosphere, and a longer warming time is required when compared with dry soils (Unesco/FAO, 1973). Much of the incident energy is utilized by the increased rate of evaporation. Similarily, soil with a high water content exhibits high thermal conductivity. It is the amount of heat transferred through a unit area in a unit time under a unit temperature gradient Thermal conductivity is dependent upon the soil's solid, water and LITERATURE REVIEW / 8 air fractions. Hillel (1980) states that the thermal conductivity of water is one order of magnitude greater than that of air. When soil pores become constricted or saturated with water, air is displaced thus increasing the soil's thermal conductivity. 2.1.2. Trafficability and Timeliness Farm operations are time dependent and cannot usually be postponed for any length of time without incurring further risks of reduced yield and increased cost of the production system. Timeliness of farm operations involves the consideration of the soil's condition and response to any external forces subjected to i t The soil structure is affected directly by the shrinking and swelling due to varying water table level and indirectly by the weather and biotic elements. Introducing a subsurface drainage system further influences the soil structure to the farmer's advantage. By reducing the soil moisture content and correspondingly increasing the matric suction, the soil strength is improved. By their field experiments, Paul and De Vries (1979) showed a linear relationship between soil strength and water table depth. Trafficability is the soil's capacity to provide traction and to allow unobstructed overland passage of vehicles without receiving serious structural damage (Hillel, 1980). Amount of soil moisture removal and, hence, degree of trafficability are a function of the antecedent moisture content, type of soil, and farm activities. Bornstein and Hedstrom (1981) reported moisture level in the top. 15cm of silty clay soils to have the greatest effect on trafficability. Subsurface drainage shortens the soil drying time which is dependent upon evapotranspiration. Land preparation for seeding may therefore start earlier in the growing season. Fausey and Schwab (1969) reported that subsurface drainage allowed spring cultivation to commence one month ahead than was possible with surface drainage. Paul and De Vries (1983) reported earlier trafficability of two to LITERATURE REVIEW / 9 four weeks in silty clay loam fields and by up to fiv; weeks in a muck field when comparing subsurface drainage with no drainage. Chieng et al. (1987) further demonstrated this measure of time available for field work or opportunity days. The benfits of a well designed subsurface drainage system was found to be an increase in the number of opportunity days by 18 to 78 days from March 1st to May 31st during an average year. Table 1 shows the actual gain in opportunity days during the spring from 1983 to 1985 in the Boundary Bay area. Prasher et al. (1985) made theoretical applications of this dependence of soil trafficability on the water table depth. They investigated a method of deriving the probability of obtaining a satisfactory workable period. The probability was assumed to follow a Markov chain process and was constrained to the criterion of at least 12 consecutive days in the month of March when the water table is 60cm or more below the soil surface. The 12 consecutive days represented the minimum number of necessary opportunity days. It was proposed that this probabilistic method be used as an alternative tool for assessing the performance of a proposed drainage system. Costly computer runs of drainage simulation models could thus be avoided. When soil structure is poor due to high moisture content, long term losses in crop production may result Excessive traffic on plastic-behaving soil causes puddling and eventual compaction upon drying. This is particularly important to soils high in clay content (Reeve and Fausey, 1974). Subsurface drainage reduces this risk and also lengthens considerably the time available for tilling. Steinhardt and Trafford (1974) observed effective reduction in wheel submergence and lateral compaction when a subsurface drainage was installed. The absence of a compacted layer allows for easier root penetration and water infiltration. Also, traction pull is reduced during tillage operation. Thus, the interplay between timely drainage of the field and subsequent LITERATURE REVIEW / 10 Table I: Opportunity Days Gained From Drainage Improvement (Chieng et al., 1987) YEARS R E G I M E 1983 1984 1985 A 83 90 89 B 75 59 80 C 80 54 74 D 19 4 29 Ave. opportunity days gained (A-D) 64 86 60 production operations is vital. 2.1.3. Plant Growth and Crop Production Proper management of the water table improves the growth condition of the root zone. High water table results in insufficient soil aeration. The infiltrating water not only displaces the air, but also hinders gaseous diffusivity by decreasing air phase and increasing discontinuities between air-filled pores. (Harris and van Barel, 1967). Clark and Kemper (1967) observed oxygen to diffuse through air-filled pores 10000 times faster than through water-filled pores. One immediate growth-limiting effect of deficient aeration is reduced root respiration (Wesseling, 1974). This is crucial since energy produced from respiration is required in all the metabolic processes. A properly functioning cell is able to control solute absorption. As a cell dies, its membrane permeability increases, allowing the leakage of ions and other small molecules (Salisbury and Ross, 1978). Consequently, a poorly aerated root zone would lower crop production. A number of field experiments has been conducted to evaluate the effectiveness LITERATURE REVIEW / 11 of subsurface drainage system in improving the soil for crop production. Carter and coworkers (1973, 1983) have reported an increase in annual sugar cane yield of 24 to 32% and in the number of single crop grown per year from three to five. Sugar cane is very susceptible to prolonged high water table conditions. During winter, cane seed stubbles decay. In the spring, early root development is inhibited by low soil temperature and anaerobic soil conditions (Carter and Camp, 1981). Increased crop yields from lowered water table were attributed to: 1) increased depth of root development, 2) improved soil aeration, 3) warmer soil temperature earlier in the growth season, and 4) improved trafficability and timeliness of farm operations. This is supported by the studies of Yang et al. (1977) which demonstrated that higher root density and deeper root growth enabled the root system to explore a larger volume for nutrients and, during drought, for water in the capillary fringe. In terms of efficient water, land, and energy usages, subsurface drainage may be an option to other systems such as conventional water furrows. Campbell et al. (1978) found that potatoe plants grown under this type of water management produced 10% more with 12% less water. The overall yield increase averaged 40% in a year. In this case, up to 12.5% of productive land area was made available for cultivation. Where land is scarce, areas occupied by open drains need to be minimized. This also lessens the difficulty in maneovering farm equipments. Other crops benefitting from subsurface drainage are corn, soybeans and oats (Erickson, 1965; Schwab et al., 1966, 1974; Carter and Camp, 1978). LITERATURE REVIEW / 12 2.2. SUBIRRIGATION Given the proper physical conditions, an existing subsurface drainage system can be operated as a "controlled and reversible" drainage system. Termed as subirrigation, this is the process of upward water movement from the water table and capillary fringe as a result of an artificially generated hydraulic gradient Incorporating subirrigation replenishes the soil profile with needed water during periods of drought due to lack of irrigation, rainfall or overdrainage. The feasibility of subirrigation is determined by several existing natural conditions. Early experiments by Fox et al. (1956) indicated the importance of a high natural water table or a shallow impermeable layer to prevent excessive seepage losses. The dual drainage/subirrigation system also requires a gently sloping topography of 0.5% or less (Massey et al., 1983). It is crucial that the water table rise to the proper height at an adequate rate. Thus the soil should provide high hydraulic conductivity similar to those exhibited by fine sandy soils and coarse silty soils. In the case of heavier or finer soil structure, in which water rises slowly, two alternatives have been studied. Working with silty clay and silty loam soils, Carter et al. (1983) used concrete and plastic borders to confine the field plot and restrict lateral water movement Doty and Parsons (1979) experimented with a water mound developed above the drain tiles. The latter method enhanced the precision of subirrigation by reducing the large water table fluctuation found within the sandy clay loam profile. The use of subirrigation is well suited to crops whose yield response is maximum at a certain minimum water table depth. Doty et al. concluded from regression analysis that the number of days the water table was less than 106.7cm from the soil surface influenced production yield. Their second-order regression function shows that for each additional day between 25 and 55 days that the water table was LITERATURE REVIEW / 13 106.7cm or less from the surface, silage yield could be increased 110 to 220 kg/ha. Other advantages to this type of water management system include: 1) reduced initial cost as compared to separate drainage and irrigation set ups, 2) low labour and maintenance requirements, 3) simultaneous operation of subirrigation with other farm activities, and 4) reduced nutrient leaching due to reduced downward water flow. Massey et al. (1983) made a comparative study between subirrigation and sprinkler irrigation while focussing on the energy and water requirements of the systems. They showed that a gravity-operated subirrigation system was more energy efficient than a sprinkler nozzle operation. 2.3. CAPILLARY FLOW Capillary flow is an important factor affecting subsurface drainage and subirrigation. It occurs in the capillary fringe where the moisture is continuous. The capillary fringe is commonly defined as the soil zone in which the pores are saturated, but the water pressure is less than that of atmospheric. It is also the region of uniform moisture above the water table where the hydraulic conductivity is essentially the same as the saturated hydraulic conductivity (Childs, 1957). Thus the capillary fringe provides an additional water pathway. The phenomenon of soil water capillarity has been extensively modelled by the behaviour of water rise and retention in a capillary tube. The migration of capillary water occurs in response to equilibrating water pressure differences. When expressed as a stress term of capillary potential, capillary flow is similar to electric or heat current in that the direction of flow is always from high to low potential levels. Physical factors determining the extent of capillary flow include: soil pore size and conformation, degree of saturation, surface tension and pressures in the water relative LITERATURE REVIEW / 14 to atmospheric pressure (Spangler, 1960). • The distribution of capillary potential fluctuates simultaneously with changes in soil moisture content Capillary potential itself is dependent on the radii of menisci between soil grains, which is determined by the moisture level. A reduction in overall soil moisture content causes water to recede farther into the interstices between grains, and therefore decreases the radii and the capillary potential. This results in an increased attraction for water in order to re-establish a state of equilibrium. The nonequilibrium that causes capillary flow can occur under several conditions. Water is drawn upward in response to water loss from the capillary zone by transpiration of crops and by evaporation at the soil surface. Stuff and Dale (1978) have found capillary water to supply an average of 27% of the evapotranspiration during periods of little or no precipitation. Thomas et al. (1977) demonstrated the differing capillary potential distribution below and above subirrigation laterals for barley and corn crops grown in two differently textured soils. They accounted the distribution patterns to varying degrees of root activity and the nonlinearity of the hydraulic conductivity. The actual height of capillary rise from the water table can vary widely depending on the textural and humus content differences in soil profile and on the depth of the water table (Bloemen, 1980). This is illustrated in Figures 2, 3 and Table II. The direction of flow is reversed when rain or irrigation water falls on the ground surface, inducing downward seepage. This reversed capillary potential gradient is enhanced by the presence of subsurface drainage tiles which act as a sink. LITERATURE REVIEW / 15 80 - i 0 @3 i i i i i i i i | 1 1—i—i—i 'i i i | r——i—i—i—i i i i | 1 1—i—i—i i i i | 0.V 1 10 100 1000 SUCTION h (cm) FIGURE 2: Relationship between c a p i l l a r y r i s e and water table depth for sandy subsoils of various geo-genetical group (Bloemen, 1980) UTERATURE REVIEW / 16 2 0 0 n 0 n 1 1 1 1—i— i—i— i—i 1 1 1 1—i—i—i—i—i 100 1000 10000 SUCTION h (cm) • FIGURE 3: Height of c a p i l l a r y r i s e i n r e l a t i o n to s o i l moisture suction i n marine clays at increasing clay contents (Bloemen, 1980) Table II: C h a r a c t e r i s t i c grain s i z e d i s t r i b u t i o n of sandy: s u b s o i l s of various geo-genetical groups, i n percentages of the mineral p a r t of the s o i l (Bloemen, 1980) GROUP 0 - 2 2 - 1 6 1 6 - 5 0 5 0 - 7 5 7 5 - 105- 150 2 1 0 - 3 0 0 - 4 2 0 - 6 0 0 - 850+ 105 150 2 1 0 3 0 0 4 2 0 6 0 0 8 5 0 (/jm) (yum) (yum) (jam) (jum) (jum) 0-jrn) (;um) (yum) (jam) (jjm) Oum) Young marine 2 . 1 1 .0 2 . 3 - 9 . 4 2 9 . 4 2 7 . 6 2 4 . 0 3 . 8 0 . 3 0 . 1 O l d marine 6 . 2 2 . 6 1 3 . 8 3 2 7 . 5 3 3 . 8 1 4 . 7 1 .6 - - -F l u v i a l 1 1 0 . 5 2 . 5 2 . 0 5 . 0 1 8 . 0 2 4 . 0 1 6 . 0 1 7 . 5 9 . 5 3 . 0 P l e i s t o c e n e a e o l i a n - 2 . 5 4 . 0 6 . 0 1 4 . 0 3 6 . 0 2 1 . 5 1 0 . 5 4 . 0 1 .0 0 . 5 LITERATURE REVIEW / 18 2.4. DRAINABLE POROSITY Drainable porosity (0, or specific yield, is one of the important parameters used in subsurface drainage design. Drainable porosity is the volume of water that can be drained from a unit of soil by a unit depth of water table drawdown in the absence of a natural source or sink. Conversely, f is the volume of water per unit area taken up by the unsaturated soil above the water table for a unit rise in the water table (Prasher, 1978). Drainable porosity has not always been included in the analysis of agricultural drainage, especially when steady state flow into drainage tiles is assumed (Bouwer and van Schilfgaarde, 1963; Toksoz and Kirkham, 1971). In the area of analytical and water balance modelling of non-steady subsurface drainage, f is commonly assigned a constant value particular to the soil under study. The analytical model of Young and Lignon (1972) uses a logarithmic relationship between drainable porosity (f) and hydraulic conductivity (K) developed originally by the U.S. Bureau of Reclamation. Skagg (1976,1981) has further replaced this basic parameter by a ratio expression of K / f in the water management model, DRAINMOD. The weakness of these analyses arises from the erroneous assumption of f being a constant for the entire soil profile considered. Early investigation by Childs (1960) and Taylor (1960) have shown f to be affected by both the water table depth and the rate of depth change. A functional relationship between f and water table depth was proposed; Taylor (1960) suggested an exponential relationship. Dos Santos and Young (1969) derived the functional relationship as "the sum of the air content at the soil surface and the ratio of the rate of change of the volume of water held in the moisture profile above the water table to the rate of rise or fall of the water table". Inherent in their calculation is the effects of rainfall and evaporation upon the soil profile. Another approach taken by Duke (1972) considers the hydraulic properties of capillary LITERATURE REVIEW / 19 pressure, bubbling pressure, and effective saturation as defined by Brooks and Corey (1964). To account for this characteristic of variability, researchers at MacDonald College of McGil l University (Montreal, Canada) developed a water balance model which assigned different values of f in relation to varying depth zones in the soil profile (Foroud, 1974; Chieng, 1975; Bhattacharya et al., 1977; Chieng et al., 1978). In their case, the soil profile was divided into four successive zones of interest Each zone was assumed to have its available soil moisture storage depleted by evapotranspiration after its transient storage was depleted. Bhattacharya and Broughton (1979) subsequently derived a curvilinear relationship between drainable porosity (f) and water table depth (h). A nonlinear least-square curve fitting computer program was applied to their extensive laboratory data. The resulting nonlinear regression equations shown in Table 3 were generated for the sand and the clay soils. These f-h relations are applicable within the constraints of a maximum suction of 2m of water and a maximum drain depth of 1.5m, which are conditions generally found in the arable regions of the Ottawa-St Lawrence Lowland in Canada. 2.5. WATER BALANCE MODEL In the present study, a water balance approach is used to analyze water movement within a defined soil volume. Simply stated, the concept of water balance defines the water content of a given soil volume at any time as the balance between water inflow and outflow and the change of the soil water storage. Mathematically, a water balance model can be expressed as: Inflow = Outflow ± Change in Storage (2.1) Such water balance models have been developed to simulate soil moisture changes due LITERATURE REVIEW / 20 Table III: Nonlinear Regression Equations (Bhattacharya and Broughton, 1979) SOIL TYPE TYPE OF EQUATION VALUES OF FITTED PARAMETERS Upland sand f - ( l/{ 1 + 2 ~ B l ( h " °' 5 )} - l/{ 1 + 2 ° - 5 B l } )6 2 0, = -11.6466 Bz = 3.2993 f = h ( e _ e ' h - e" 6 , , h ) + B s ( 1 - e" B e h ) Ba - 0.4154 B » = 0.3484 Bs - 0.3305 3 6 - 0.6480 to the combined effects of climate, crop growth, and multicomponent water management systems, which facilitate subsurface drainage, surface drainage, subirrigation, and sprinkler irrigation (Baier and Robertson, 1966; Jensen et al.. 1971; Hillel et al., 1975; Anderson et al., 1978; Skaggs, 1978;. Broughton and Foroud, 1978; Foroud and Broughton, 1978). The water balance model applied in this study was originally developed by Foroud (1974) and further modified by Chieng (1975) to optimize the subsurface drainage design. The model computed the daily available soil moisture content at the midspacing between adjacent parallel drains along a profile of four successive soil zones extending from the soil surface to the drain depth. In addition, the model computed the water table depth to be expected in a certain soil for certain drainage system capacities with varying drain depths. An appropriate drainage coefficient could then be selected to determine the spacing between laterals from the Hooghoudt's relationship (Equation (4.7)). This water balance model has been described and applied by Chieng et al. (1978). The simulation model was designed to determine the daily change in water Ste. Rosalie clay LITERATURE REVIEW / 21 table depth: Az = A(TR)/f. = R A t / f (2.2) where: Az = change in water table depth A(TR) = change in transient storage f. = drainable porosity of the ith zone R = drainage rate At = time interval The result was then used to obtain the revised water table depth : \ = \_i ± Az (2.3) Since each successive soil zone was assigned a fixed drainable porosity value, separate falling and rising water table equations were needed to describe the movement across the boundaries of any two zones. These are shown below: 1) the falling water table equation: z t = Dj_ - [ (D : - z ^ ) ^ + A T R ] / f 2 (2.4) 2) the rising water table equation: z t = D x + [ (z t _ 1 - D : ) f 2 + ATR]/ f j (2.5) where: = depth of the 1st soil zone f-^ = drainable porosity of the 1st soil zone ?2 = drainable porosity of the 2nd soil zone In this model, the drainable porosities for the third and fourth zones were assumed to be the same as that of the second zone. The modifications made in the present model re-expresses the drainable porosity as an explicit function of water table depth. This replaces the unrealistic representation of an average porosity value for a particular depth interval. The derived function is LITERATURE REVIEW / 22 based on Taylor's (1960) observation as possibly being an exponential one with parameters specific to the local conditions. More recent work by Bhattacharya and Broughton (1979) has supported this assumption. For this study, effort was made to keep the function simple and parameters few. 3. MATERIALS AND METHODS 3.1. EXPERIMENTAL FIELD AND MEASUREMENTS The experimental site of 3.4 ha was located at the northwest corner of the Boundary Bay Airport in the Lower Fraser Valley of British Columbia, Canada. Common to this part of the province is the low-lying, flat alluvial landform. The soil was classified as Humic Luvic Gleysol developed on moderately fine to fine textured deposits with moderately poor to poor drainage (Driehuyzen, 1983). The soil was also found to belong to the Ladner Series. Four different types of crops were planted on individual treatment plots. These included corns, grass forage, strawberry, and potatoes. The crop layout in the experimental field was randomly designed with the corns planted on the west side to reduce shading effects on adjacent crops. Each water-controlled treatment plot size was 42m wide by 100m long, giving an area of 4200m2. The drainage set up for each plot consisted of three parallel drainage pipes installed at an average depth of 1.1m from the soil surface. These 100mm-diameter laterals spaced at 14m apart were perforated and corrugated polyethylene tubings. They ran a length of 100m in the east-west direction before emptying into a nonlined earth ditch at the east end as shown in Figure 4. An open ditch was excavated and used as a collector. For the purpose of the study, the ditch was subdivided into three sections and interconnected by plastic pipes at the ditch bottom. A pump was installed to empty the collector ditch into the main ditch. The highway ditch along 72nd Street then carried all the water from the main ditch and away from the drained area. For the purpose of subirrigation, overflow stand pipes of proper height controlled the different water levels in the experimental plots. Water 23 Open Ditch r A (See Figure 5) Laterals 100mm DIA REGIME A REGIME B REGIME C V 20m 14m 100m • • A" j — j ^ - R a i n Gauge ^Collector Ditch L, REGIME D LEGEND • A Automatic water table recorder Sampling s i t e • Drain t i l e :Open d i t c h FIGURE 4: Experimental F i e l d Layout MATERIALS A N D METHODS / 25 SECTION A - A (NOT TO SCALE) FIGURE 5: The Water Control System from the domestic supply line was introduced into the collector ditch from which it flowed above the drain lines and into the soil The water control system is diagrammatically illustrated in Figure 5. In operating the water control system, certain criteria were considered. One criterion for subsurface drainage was to keep the water table low enough such that the SEW 3 0 value did not exceed 200cm during the 1st of of November to the 31st of March. SEW J 0 is the sum of daily water table recorded above the 30cm depth at midspacing between laterals. Defined originally by Seiben (1965), the SEW value is an index to crop yield in relation to the water regime in the upper soil layer. Above a certain minimum value, crop yield is negatively affected by the presence of excessive water. A second criterion was to obtain trafficability of the land during early spring from the 31st of March to the 1st of May. Different subirrigation operations were to be maintained throughout the summer and early fall months when periods of water deficit were prevalent Four water table regimes were tested at the experimental set up: 1. REGIME A: free subsurface drainage at all times; water table controlled at or below subdrains MATERIALS A N D METHODS / 26 2. REGIME B: subsurface drainage during high precipitation periods; subirrigation by maintaining water table depth at 60cm from soil surface during periods of water deficit 3. REGIME C: subsurface drainage during high precipitation periods; subirrigation by maintaining water table depth at 30cm from soil surface during periods of water deficit 4. REGIME D: no drainage or subirrigation at any time The subirrigation levels of 60cm and 30cm water table depth were based on published sources. In Holland, the water table is to be kept at a minimum of 50cm below soil surface for arable land and 40cm for grassland (Luthin, 1978). In the B.C. Agricultural Drainage Manual, the required water table depth is at 50cm or more below soil surface 48 hours after a storm rainfall event ceases so that crop and soil structure damage is minimized. The height of capillary rise in the local field have been observed to be approximately 20cm (Chieng, 1987). It was concluded that the two proposed subirrigation levels represented the lower and upper limits from the reported optimum of approximately 50cm depth. Subsurface drainage began during the winter of 1981 when the subdrains and sub-ditches were installed at the experimental site. The pump was installed in early 1983 and put into operation for subirrigation in the summer of 1983. The water table height at the midpoint between laterals were recorded using automatic water table chart recorders. Three were installed in the drained area (Regimes A, B and C) and one in the adjacent unchained area (Regime D). In situ rainfall readings were not collected until 1984 when a rain gauge was placed next to the pump station. For this reason, rainfall measurements recorded at the Delta Ladner South weather station, which was about 5km away from the site, were used for analyses prior to this date. The MATERIALS A N D METHODS / 27 locations of the recording instruments are indicated in Figure 4. The daily values of potential evapotranspiration were computed from Equation 1 of the technique proposed in 1965 by Baier and Robertson (Russelo et al., 1974). This method took into account the local maximum tempertature, temperature range, and incident solar energy on a daily basis. 3.2. PHYSICAL CHARACTERIZATION OF SOIL Laboratory analyses of the soil in each treatment plot were conducted to determine the physical and hydraulic properities of: satiated hydraulic conductivity, drainable porosity curve, bulk density, particle size distribution, and moisture retention curve. Three undisturbed core samples were randomly collected from each of the four treatment plots. Their approximate locations are indicated in Figure 4. The core dimensions were 7.3cm in diameter and 7.5 cm in height These were removed from a depth halfway between maintained water table level and soil surface with a cylindrical metal sampler. The sampling depth for Regimes A, B, C and D were 50, 30, 15, and 50 cm, respectivley. At the same time, loose soil samples were obtained at successive depth intervals from each core sampling site for the analyses of soil moiture content particle size distribution, and moisture retention curve. The satiated hydraulic conductivity was determined for each core sample. Since the cores were extracted above the water table level where the soil pores were essentially saturated but with the water held at less than atmospheric pressure, the satiated and not the saturated hydraulic conductivity was representative of the field condition. The falling head method described by Klute (1965) was applied with a few modifications. The core bottom was wrapped with cheesecloth and placed on top of a MATERIALS A N D METHODS / 28 steel grid, thus exposing it to the atmosphere. An accompanying steel cylinder taped to its top end acted to contain the draining water. In place of the standard stand pipe, a screw dial gauge was set up to measure the hydraulic head drop. At least four runs were performed on each core sample. The statistics of mean, standard deviations, and coefficient of variance were computed for analysis. The drainable porosity curve consisted of discrete porosity values for the tension range of 0 to 80cm at 10cm intervals. The unconsolidated porous plate was constructed of saturated, seived sand high in ferrous and aluminum oxide content The plate was large enough to hold all twelve core samples. The sand was initially passed through a 0.246mm seive. To ensure complete saturation, the plate was vibrated temporarily using a pneumatic vibrator set at 70kPa. This step caused close-packing configuration to take place, thus removing any trapped air. Preparation of the core samples involved initial drainage at 60cm tension for 24 hours, followed by saturation from below for 24 hours. Tension measurements were taken from the midpoint of the core height Each sample was weighed at saturation and again at each successive tension level when drainage ceased to continue. The plate was moistened prior to replacement of the core to prevent discontinuity of water at the interface. The apparatus was observed to fail due to air entry at 90cm tension. Bulk density determination was the final analysis performed on these core samples. The standard core method (Blake, 1965) was applied. For speed and convinience, the hydrometer method for particle size analysis was chosen. Pretreatment of the loose soil samples consisted of screening through a 2mm seive followed by removal of organic matter using a 6% hydrogen peroxide solution. Soil dispersion was aided by the addition of a Calgon solution. The entire experiment was carried out in a temperature-controlled room. The step-by-step procedure is MATERIALS A N D METHODS / 29 detailed by Day (1965). The soil moisture retention curve for each sample was obtained with a standard pressure-plate and pressure-membrane apparatus (Richards, 1965). Tensions of 3, 9, 30, and 150m of water were applied on all samples. Due to the low holding capacity of the pressure chambers, samples were tested only in duplicates. 3.3. F O R M U L A T I O N OF DRAINABLE POROSITY EQUATION 3.3.1. In Situ Drainable Porosity The in situ drainable porosity values were derived in two ways using the equation given by Broughton and Foroud (1978): f = R(At/Az) (3.1) One method was to define R as the daily rainfall intensity measured in the field. From the water table charts, the total change in water table level from start to end of the related rainfall event was measured to give Az. This dealt with the rising limb of the water table curve and hence represented the subirrigation aspect of the water table control. The time lag between start of a rainfall event and response of the water table was observed to be approximately 12 hours. Only distinct rainfall events and water table peaks were chosen to calculate drainable porosity values. The other method was to define R as the in situ drainage rate or specific discharge. From measurements of total discharge in the field and the known area of discharge, the drainage rate was calculated. These field data were made available by Kerr-Wilson (1985). Since the instantaneous drainage rate was recorded, the change in both time and water table position were re-expressed as the slope on the water table charts at the time of measurement It should be noted that the drainable porosity MATERIALS A N D METHODS / 30 derived this way concerned with the falling limb of the water table curve and hence the drainage behaviour of the water table movement Again, only drainage rate which indicated a clear drop in water table level were used. 3.3.2. Curve Fitting Two possible types of regression functions were considered in the curve fitting process: (1) the polynomial equation f = b 0 + bjz1 + b 2z 2 + b 3z 3 + ... + b n z n (3.2) (2) the negative exponential equation f = b„( 1 - e " b l Z ) (3.3) Both functions consisted of the independent variable, tension (z) and the dependent variable, drainable porosity (0- The initial regression analyses were performed on the laboratory-derived f data to determine the more suitable function. These data were used since they provided an orderly range of drainable porosity and tension values. Also the drainable porosity values determined in the laboratory were considered to be correct over those determined from rainfall intensities and drainage rates. In this study, it was assumed the conditions in the laboratory setting were more controlled than those found in the field. The appropriate function was then fitted to the in situ data for each of the four water regimes. 3.3.3. Corrected f-z Function The drainable porosity curve derived from laboratory procedures on core samples was considered to be an accurate representation of what happens in the field, especially when given an adequate sample size. However, due to the time and labour MATERIALS A N D METHODS / 31 consumptiveness of this method, it is oftentimes impractical. In this research, one of the objectives was to use in situ data of rainfall intensity, drainage rate, and water table depth, to generate the drainable porosity curve. Such data could be conviniently collected using meteorological instruments that are easy to install and operate on site. Also, these data inherently take into account the actual effects of evapotranspiration, rainfall, and hysteresis on water table movement Although use of meteorological data is emphasized, the measurement of a drainage event was found to be difficult compared to that of a rainfall event Thus, the conclusion was to focus on the rainfall data to generate the drainable porosity curve. However, since a rainfall event represented only the rising limb on a water table curve, the derived f-z function needed to be corrected in order to account for water table drawdown. 3.4. WATER BALANCE MODEL The main and subroutine programmes of the water balance model was written in the F O R T R A N IV G language for an IBM 360/370 computer. For practicality, the incorporation of the f-z function was best facilitated by a subroutine programme. Its usage was then easily invoked by a C A L L command within the main programme. The parameters passed were the water table depth and the constants of the f-z equation which were specific to each water regime. Another subroutine programme added dealt with the transient storage of the soil moisture. The original water balance model described soil moisture storage by two components: (1) soil moisture storage which was the difference between the field capacity and the wilting point; and (2) transient storage which occurred in the drainable pore space. The soil moisture storage was assumed to drain as quickly as MATERIALS AND METHODS / 32 the specified drainage coefficient would allow. When soil moisture of the whole soil profile reached its field capacity, excess moisture was assumed to enter transient storage causing the water table to rise. Since the drainable porosity across the soil profile was defined by a continuous function of depth, the maximum total and zonal transient storages must be calculated instead of assigning constant values as previously. 3.5. STATISTICAL ANALYSES A nonparametric sign test was used to test the degree of accuracy by which the water balance model and hence the incorporated f-z function could simulate the behaviour of the water table over a period of one year. The Wilcoxon's matched-pairs signed-rank test, or simply, the Wilcoxon paired-sample test, was chosen to test the null hypothesis that there is no difference between actual and simulated water table levels. This test utilized both direction and magnitude of the difference between paired observations (Siegel, 1956). The assumptions made were (1) a random sample of pairs of observations have been taken, and (2) the absolute difference in the paired observations could be ranked. As with any nonparametric method, no assumptions were made about the form of the population probability distribution. Since the direction of the difference was not predicted in advance, a two-tailed region of rejection was required. 4. RESULTS A N D DISCUSSION 4.1. SOIL PHYSICAL ANALYSES The results of the physical analyses of the four experimental plots (Regimes A, B, C and D) are tabulated in Table IV to Table VI. Each analysis was performed on a sample size of three (n=3). A l l plots were found to have a textural classification of silt loam under the standards of the United States Department of Agriculture (USDA). In the analysis of particle size distribution, Regime B exhibited the highest silt content while Regime C exhibited the lowest The reverse trend was observed in their sand contents. The clay contents were similar among the plots with their sand contents varying the greatest The dry bulk density values of 1.46, 1.32, 1.29 g/cm3 for Regimes A, B, C, respectively, must be explained with respect to their sampling depths — 50, 30, 15cm from the soil surface. Thus, although the drainage practice in Regime A with no subirrigation should lead to a less dense soil matrix, its associated core analysis at a much greater depth revealed a naturally increasing dry bulk density with soil depth. The overall range of 1.29g/cm3 to 1.46g/cm3 was realistically agreeable with the value of 1.40g/cm3 reported by Driehuyzen (1983) given the small sample size and the inherent variability in evaluating this soil property. In addition, a concurrent study by Kodsi (1987) using 12 composite soil samples per plot resulted in a characteristic dry bulk density of 1.47g/cm3. 33 RESULTS AND DISCUSSION / 34 Table IV: Soil Particle Size Distribution (USDA) REGIME %CLAY %SILT %SAND n DEPTH (cm) A 21 69 10 3 50 B 23 73 4 3 30 C 26 58 16 3 15 D 22 64 14 3 50 Table V: Textural Classification (USDA) REGIME CLASSIFICATION n DEPTH (cm) A silt loam 3 50 B silt loam 3 30 C silt loam 3 15 D silt loam 3 50 4.2. SOIL HYDRAULIC ANALYSES Hundel et al. (1976) reported a number of soil properties affected by drainage. Particularly important to the hydraulic behaviour of the soil profile are the hydraulic conductivity, pore size distribution and drainable porosity. The freely drained plot of Regime A displayed the highest satiated hydraulic conductivity among the treatment plots as shown in Table VII. Conversely, the satiated hydraulic conductivity decreased progressively with increased subirrigation practices. This was illustrated by the near 50% reduction between Regimes B and A and by the reduction of one order of magnitude between Regimes C and A. From Figure 6, Regime A also correspondingly exhibited RESULTS A N D DISCUSSION / 35 Table VI: Dry Bulk Density G I M E M E A N S.D. C.V. R A N G E n DEPTH (cm) (g/cm3) (g/cm3) (g/cm3) A 1.46 0.04 0.03 1.43 - 1.51 3 50 B 1.33 0.14 0.11 1.03 - 1.50 3 30 C 1.29 0.06 0.05 1.23 - 1.35 3 15 D 1.30 0.31 0.24 0.94 - 1.49 3 50 Table VII: Satiated Hydraulic Conductivity C.V. R E G M E A N (m/s) S.D. (m/s) R A N G E (m/s) n DEPTH (cm) A 2.6(10-6) 0.3(10-6) 0.13 1.6(10"6) - 3.8(10"6) 12 50 B 1.5(10"6) 0.5(10"6) 0.35 2.4Q0"7) - 4.2(10"6) 15 30 C 3.9(10-7) 0.4(10-7) 0.19 9.8(10-8) - 9.5(10"7) 12 15 D 1.7(10"6) 0.2(10"6) 0.13 4.0(10"7) - 3.2(10"') 15 50 the lowest moisture retention behaviour at all water tension levels. Similar works by McWhorten and Duke (1976). Skagg (1981) and Shih (1983) support this observation. In the presence of a shallower water table resulting from poor drainage or subirrigation, more moisture is brought into the surface zone by the interactive forces of evapotranspiration and capillary rise. Between Regimes B and C, however, Regime B consistently showed greater moisture content over the entire tension range of 0m to 158m. Due to the inadequate sample size of 6 per regime per tension levels of 0, 3, RESULTS AND DISCUSSION / 36 -O o 50 40 30 H O O 20-| UJ (Z Z) X 10 20 —I— 40 60 ~80~ 100 TENSION (m) Legend A SCHEME A X SCHEME B • SCHEME C B SCHEME D 120 i 140 160 FIGURE 6: S o i l Moisture Retention Curves 10 n 1 4 Q TENSION (cmh zo) FIGURE 7 : Drainable Porosity Curves From Laboratory Analysis RESULTS A N D DISCUSSION / 37 9, 29 and 158m, no statistical analysis was performed to determine if this difference is significant From simple observation, DeVries (1984) concluded that the effects of the two subirrigation practices would not be apparent after only one summer of actual operation. 4.3. LABORATORY DETERMINED IN SITU DRAINABLE POROSITY The drainable porosity curves for Regimes A, B, C and D in Figure 7 show clearly the effect of subirrigation on drainable porosity. The approximate drainable porosity at 120cm tension were 9, 7, 5 and 9% by volume for Regimes A, B, C and D, respectively. These differences were due to the change in soil structure resulting from reversed drainage flow. It was speculated that possibly the transport and deposition of soil particles within soil pores reduced the soil's capacity to release or absorb water in reponse water table movement However, some researchers have argued that this convective process would be minimal for a stable soil where the extent and frequency of water table fluctuations were not excessive. Further field work would be needed to confirm this. Referring again to Figure 7, the drainable porosity of Regime A was consistently greater than that of Regime D between the tension range of 0cm to 100cm. According to Hundel et al. (1976), subsurface drainage increases the density of large and medium sized pores. The improved hydraulic behaviour of Regime A over Regime D was also made evident in the hydraulic conductivity and moisture retention analyses. RESULTS A N D DISCUSSION / 38 4.4. IN SITU DRAINABLE POROSITY 4.4.1. Rainfall Rate Method Drainable porosity computed from water table charts and rainfall readings for the years 1983 and 1984 are shown graphically in Figures 8 to 11. A typical hydrograph followed the pattern of sharp rise and then slow recession with the receding limb extending much with respect to time before reaching a plateau. Symmetry was preserved only when a subsequent rainfall event occurred prior to completion of recession of the present hydrograph event. Another point of importance was the time interval the meteorological data represented. Since the rainfall data was collected on a daily basis, the accumulated value was considered to occur over the entire 24-hour period. This limited the accuracy with which the porosity could be computed although the scale on the water table charts was more refined at 8 hours per axis division. A general comparison between these porosities and those determined in the laboratory indicated a trend of lower values by this approach. This was especially true for porosity values falling within the laboratory tension range. The best agreement was found with the higher subirrigation of Regime C. Regime D showed the poorest agreement although there was a lack of valid data. With the undrained treatment, it was observed that during initially high water table conditions, a rainfall event did not induce a large rise in water table level. From Equation (2.2), this condition resulted in large porosity values. Physically, this implied large pore sizes and well distributed pore spaces within the upper 10cm depth in the undrained plot R E S U L T S A N D D I S C U S S I O N / 39 o > IS >-C O o o 0_ Ld _ l CD < < cr 16-, 14-12-10 6-2-20 I 40 60 i 80 100 TENSION (cmh 2 o) 120 140 FIGURE 8: Drainable Porosity Based on R a i n f a l l Rate Method For Regime A 16-j > 14--O 12-TY 10-o 8-o CL. 6-Ld j C D 4-< < 2-Q 0-0 20 40 60 80 100 120 140 TENSION (cmh 2 o) FIGURE 9: Drainable Porosity Based on R a i n f a l l Rate Method For Regime B RESULTS AND DISCUSSION / 40 16 —i o > 14-s 12-ITY 10-(/) o 8-cz o Q_ 111 6-BLI 4-< < 2-cz Q 0-20 40 60 80 100 TENSION (cmh 2 o) 120 14 0 F I G U R E 10: D r a i n a b l e P o r o s i t y B a s e d o n R a i n f a l l R a t e M e t h o d F o r R e g i m e C 16-j o > 14--O 12-i— 10-00 O 8-cz o Q_ i. i 6-BLE 4-< <c Z — cz o 0-0 I 20 40 60 80 100. TENSION (cmh ? o) —i 1 120 140 F I G U R E 11: D r a i n a b l e P o r o s i t y B a s e d o n R a i n f a l l R a t e M e t h o d F o r R e g i m e D RESULTS A N D DISCUSSION / 41 4.4.2. Drainage Rate Method The drainage rates computed from the division of collected total discharge by area drained were representative mainly of the summer months in 1984 (Kerr-Wilson, 1985). Since the total discharge was collected from the middle lateral of a treatment plot, the drainage area was 1400m2. Data collection was performed only on the treatment plots of Regimes A and B. These data were limited due to the absence of drain flow resulting from either lack of rainfall events or random sampling schedules. Porosities derived from this method are shown in Figures 12 and 13. By general inspection, these porosities tended to be higher than those derived by either the laboratory or rainfall rate method. 4.5. T H E DRAINABLE POROSITY EQUATION 4.5.1. Curve Fitting A MIDAS (Fox and Guire, 1976) program and a SAS (Ray, 1982) program were used to determine the suitable curve fitting equation between the polynomial and the negative exponential regressions. The sample coefficients of determination, R 2 , for all regimes under each regression are shown in Tables VIII and IX. The R 2 value was found to be quite high in either cases when considered separately. To compare the two regressions, the correlation index, F , was necessary. This statistic was computed for the negative exponential regression, and, from Table X , again a high correlation value was found with each regime. The negative exponential equation was the logical choice for the following reasons: 1. simplicity in expression, fewer number of parameters; 2. inherent asymptotic characteristic. RESULTS AND DISCUSSION / 42 16-i O > CO O o Q_ CO < z; < o 14 12 10-8-6-4-2 H 0- l 20 r~ 40 60 80 100 TENSION (cmruo) 120 14 0 FIGURE 12: Drainable Porosity Based on Drainage Rate Method For Regime A ^ 16 n > 14 1* 1 2 >- 10 CO g 8 o U J _ l m 4-< z: < 2-ct: Q 0-20 40 60 80 100 TENSION (cmh-o) 120 140 FIGURE 13: Drainable Porosity Based on Drainage Rate Method For Regime B RESULTS A N D DISCUSSION / 43 Table VIII: Regression Statistics to the Negative Exponential Curve Fitting R E G n m df SSreg SSres SSy R 2 Syx A 39 1 37 2107.594 28.130 2135.724 0.987 0.872 B 39 1 37 1329.257 15.242 1344.499 0.989 0.642 C 39 1 37 710.533 11.127 721.660 0.985 0.548 D 39 1 37 1941.974 10.423 1952.397 0.994 0.531 Table IX: Regression Statistics to the Polynomial Curve Fitting R E G n m df SSreg SSres SSy R 2 Syx A 39 6 32 233.660 25.351 259.011 0.902 0.890 B 39 6 32 146.160 14.842 161.002 0.901 0.681 C 39 6 32 80.899 10.828 91.727 0.882 0.582 D 39 6 32 252.120 7.631 259.750 0.971 0.488 Table X : Correlation Index Statistics to the Negative Exponential Curve Fitting R E G n m df SSreg SSres SSy I 2 Syx A 39 1 37 2107.571 28.153 2135.724 0.993 0.872 B 39 1 37 1329.117 15.382 1344.499 0.994 0.645 C 39 1 37 712.560 9.100 721.660 0.994 0.496 D 39 1 37 1941.818 10.578 1952.396 0.997 0.535 RESULTS A N D DISCUSSION / 44 With the polynomial regression, a higher R 2 value was always possible given enough parameters. However, the expression would become unnecessarily cumbersome. In this case, a polynomial of the sixth power was needed. Also the resulting polynomial fit would have to be constrained within the Ocm to 120cm tension range since it was not asymptotic. In Table XI, the parameters b 0 and bj to the chosen equation f = b„( 1 - e~ b l Z ) (4.1) are given with the R 2 statistic. Focussing on the drainable porosity derived by laboratory analysis, the maximum porosity, b 0, confirmed the apparent effect of subirrigation by its successively reduced value when comparing the regimes in the order of D, A, B and C. The rate at which the maximum porosity was reached with respect to increased tension was slowest for Regime B. In the drainage rate analysis of f, the available data was limited to the summer period of April 1st to September 30th. By this method of analysis, the resulting b 0 values were highest, 13.8% for Regime A and 11.7% for Regime B. This observation was caused by the nature of the soil to form cracks and channels, and also by the channeling effects of plant root and earthworm activities. In the rainfall rate analysis, the noticeable difference was the consistently lower f values in the tension range up to 80cm when compared to the corresponding f values derived from the previous two methods. The low sample coefficient of determination for Regimes B, C and D were due to the lack of data for tensions greater than 80cm, 50cm and the entire tension range, respectively. The drainable porosity derived by all three methods represented a cumulative value with respect to the direction in which the water table travels. With the rainfall rate analysis, the water table was rising; the wetting front was moving through a soil profile of small pores to larger pores. The incremental change in drainable porosity RESULTS A N D DISCUSSION / 45 Table XI : Parameters to the Negative Exponential Drainaible Porosity Equation POROSITY REGIME n b.(%) biU/cm) R 2 flab A 39 8.6023 0.04978 0.987 B 39 6.7775 0.05248 0.989 C 39 5.0225 0.04808 0.985 D 39 8.7246 0.03706 0.995 f . ram A 30 12.8716 0.01460 0.915 B 24 10.9392 0.01245 0.887 C 20 7.2742 0.01755 0.872 D 8 7.5398 0.03495 0.859 f drain A 13 13.7657 0.01925 0.910 B 6 11.7094 0.01956 0.971 increased with height Drainable porosity itself increased sharply upon nearing the soil surface. The water table was receding in the drainage rate analysis. Conversely, the direction of flow was from larger pores to smaller pores. The incremental change in drainable porosity near the soil surface was always greatest and decreased with each successive depth interval. Thus the incremetal increase in drainable porosity decreased as the water table dropped farther. Test for variability was performed on the two parameters b 0 and bi for the laboratory analysis. There was insufficient evidence to indicate a significant difference in bi among the four regimes. By the Duncan's Multiple Range Test, a significant difference was found to exist in b0 between Regimes C (subirrigation at 30cm) and D RESULTS AND DISCUSSION / 46 (undrained). The effects between drained and subirrigated treatments, and between the two levels of subirrigation were found not to be statistically significant in the research in view of the regression parameters. However, the trend in the differing treatment effects were present 4.5.2. Corrected f-z function An underlying objective of the present research was to utilize readily obtainable meteorological data collected on site and at weather stations, in particular water table charts and rainfall data. The application of these data to derive the drainable porosity curve was considered more convinient than that by direct laboratory measurement method and drainage rate method. The latter two methods entailed a greater amount of physical expenditure. Before the drainable porosity function based on rainfall analysis alone could be incorporated into the water balance model, it was corrected by the two regression equations derived from laboratory and drainage rate analyses. This took into account the reverse direction with which drainable porosity was determined, a rising water table event corrected by a falling water table event in order to reduce the hysteresis effect Firstly, ratio values of f . /f, , and f . /f, . were computed at various tension rain lab ram drain levels. For all regimes, a constant value, denoted by b4, existed for tensions less than 0.1cm. Above this level, the curves all followed the general negative exponential form. It was above this point on each curve that a regression equation was fitted and the constant b4 appended afterwards. The correction term was of the general mathematical form: f/f* = b2( 1 - e" b 3 Z ) + b4 (4.2) Thus the final corrected f-z equation was of the form: RESULTS A N D DISCUSSION / 47 Table XII: Parameters to the Corrected Drainable Porosity Equation POROSITY RATIO f . ram flab R E G I M E A B C D 1.08616 1.26521 0.93948 0.04984 0.009976 0.009192 0.001214 0.001942 0.438 0.383 0.529 0.815 R2 0.998 0.995 0.995 0.998 f . ram f drain A B 0.23289 0.35408 0.008327 0.007193 0.709 0.595 0.998 0.999 e " b j Z ) + b« ] (4.3) The parameters b 2, b 3 and b 4 were obtained for all regimes under the laboratory-based correction. However with the drainage-rate-based correction, only those for Regimes A and B were obtainable. Those for Regimes C and D were not obtainable due to lack of available data. These parameters are tabulated in Table XII. With these parameters, Equation (4.3) was used in the water balance model to compute the drainable porosity for a given water table depth. This was facilitated by a subroutine program. corrected = b„( 1 - e " b l Z ) / [ b2( 1 -4.6. T H E MODIFIED WATER B A L A N C E M O D E L RESULTS A N D DISCUSSION / 48 4.6.1. Assumptions to the Model The following assumptions were made for: 4.6.1.1. The Water Balance Model 1. The basic principle of water balance in the form of: Inflow = Outflow ± Change in Storage (4.4) is applied with inflow of precipitation and outflows of surface runoff, subsurface runoff and actual evapotranspiration. No inflow of irrigation is considered since artificial water replenishment was through subirrigation. 2. Two distinct storages are considered: a. soil moisture storage which is the diffience between field capacity and wilting point; b. transient storage which is the drainable pore space. 3. The soil profile is divided into four zones. Depths of first, second, third and fourth zones are 30cm, 30cm, 40cm and the remaining distance to drain depth, respectively. 4. The water balance equation to each zone is given by: SMC(I) = SMC(K) + PPT(I) -AE(I) -EX(I) (4.5) The variable counters I and K represent the current day and previous day, respectively. 5. Excess Water: EX(I) a. When the soil profile is at field capacity before a rainfall, EX(I) is equivalent to the net precipitation for that day. It is the difference between precipitation, PPT(I) and actual evapotranspiraton, AE(I). EX(I) goes into transient storage causing a rise in water table. RESULTS A N D DISCUSSION / 49 b. When the soil moisture content is less then field capacity before a rainfall, the net precipitation replenishes the soil moisture of the first, second, third and fourth zone, respectively, up to field capacity. The remaining amount of water becomes EX(I) and goes into transient storage. Actual Evapotranspiration: AE(I) When an insufficient rainfall occurs such that soil moisture is not brought up to field capacity, EX(I) is zero and actual evapotranspiration, AE(I), takes place. The total evapotranspiration from surface retention, bare soil and plant transpiration is assumed as AE(I). Water extracted as AE(I) occurs from the uppermost zone firstly and then progresses down to the next lower zone. Therefore, AE(I) is the product of the potential evapotranspiration and the A E / P E value of the zone from which water is withdrawn. a. AE(I) occurs directly from transient storage at the maximum rate of PE(I) x (AE/PE) if the water table is at the surface or in the top 25cm of the profile. b. On days with precipitation, AE(I) takes place from that rainfall event at a maximum rate of PE(I) x (AE/PE). c. On days when precipitation exceeds AE(I), the remainder replenishes the moisture storage of the first, second, third and fourth zones, respectively, followed by transient storage, subsurface runoff and surface runoff. Al l subsurface runoff occurs through the subsurface drainage system. Surface runoff occurs only when the soil profile is saturated to the surface. This is not true when the rainfall rate exceeds the infiltration rate so that surface runoff occurs before the water table reaches the surface. However, this is a conservative assumption since the event of an earlier surface runoff results in RESULTS A N D DISCUSSION / 50 less transient storage and a faster lowering of the water table then predicted by the model. 9. Water held in the drainable pore space between saturation and field capacity percolates to the drains at the design drainage rate with the drain tube center located at the bottom of the fourth zone. 10. Deep seepage is considered as part of the subsurface drain outflow. 4.6.1.2. The Water Table Model 1. Assumption of a nearly flat water table allows for the expression of daily change in water table at midspacing between drains as: Az = A(TR)/fj = R A t / f (4.6) This implies that the subsurface drainage flow is restricted by drain tube capacity. 2. When the water table is in the top 40cm of the soil profile, drainage rate is restricted by the drain tube capacity and is at the maximum design rate, R M A X . Below this depth, drainage rate, R, follows the Hooghoudt's relationship: R = 4K( 2d g H + H 2 ) / S2 (4.7) where R = drainage rate K = hydraulic conductivity d g = equivalent depth H = water table height above drain S = lateral spacing RESULTS A N D DISCUSSION / 51 4.6.13. The Computer Model 1. The water table is at the soil surface at the start of each water balance year (January 1st to December 31st). Thus the initial water table depth is assigned zero and the inital soil moisture for each zone is assigned their respective maximum values. Although in reality, this is not always the case, by middle of January the water table is usually at or very near the soil surface at least once in this region and hence, the accuracy of the model is restored there afterward. 2. In determining the drainable porosity from the f-z function, the associated water table depth is necessary. This poses a problem since the current water table depth, Z(I), to be predicted is itself needed to compute f(I). To overcome this problem, a procedure of a first crude prediction followed by refinement is used. This involves the computation of an initial crude Z'(I) based on previous day's Z(K) and f(K). Then a second drainable porosity f(I) is obtained from this Z'(I). The final predicted water table depth is derived by the use of f(I) which is an average of f(K) and f(I). Due to the shape of the f-z equation given by the negative exponential form, the iterative procedure is convergent 3. Given the condition that the water table is at the drain level and change in transient storage is positive due to positive EX(I), drainable porosity based on previous day's water table depth, Z(K), is revised to allow for the depletion of transient storage before the water table reaches the drain level. It is assumed no current rainfall event will result in a change in transient storage greater than the maximum transient storage for the fourth zone. 4. Given the conditions that the water table is low and a sudden large rainfall event occurs, only f(I), the second drainable porosity resulting from Z'(I), is used. This ensures the selection of a larger and more realistic drainable porosity RESULTS A N D DISCUSSION / 52 value and the result is a more subdued rise in water table. The basic daily algorithm to the water balance model is illustrated in Figure 14 as a generalized flowchart. The F O R T R A N program itself is listed in the appendix. 4.7. WATER T A B L E SIMULATION Weather conditions for the two consecutive years of 1984 and 1985 were quite different The low precipitation and thus dry conditions of 1985 were apparent in the lower water table levels and fewer sharp peaks compared to those of 1984. The actual water tables for each regime of these two years are shown in Figures 15 and 16. In the simulation, the maximum drainage rate, R M A X , was the varying parameter used to arrive at a visually acceptable agreement between the model and actual water table depths. Figures 17 to 24 give the simulation results to Regimes A, B, C and D for the two years. These were based on the laboratory-corrected drainable porosities. The results based on drainage-rate-corrected drainable porosities are given in Figures 25 to 28. The model set the initial position of the water table at the soil surface. There was no minimum depth to which the daily water table could be computed. The water table height was determined by the Hooghoudt's relationship, Equation (4.7). This equation was valid only for the prediction of water table height at or above the drain level. In the simulation, whenever the water table was at or below the drain level, the water table height was set to zero. This assumption would not affect the continuity principle since the mass balance was performed for the soil profile above the drain lines only. From all the simulatons, the water table tended to reach and stay at the drain level during the period of June 11th to November 10th in 1984 and May 10th to November 3rd in 1985. These results corresponded well with the Q START^ S INPUTS: 1. B t a r t l n g date 2. ending date INPUTS: 1. maximum drainage rate, EHAX 2. maximum drain depth 3. depth of impermeable layer 4. depth to where RHAX occurs 5. total available s o i l moisture 6. available s o i l moisture of each zone Inputs: 1. parameters to f-t equation 2, AE/PE ratio of each tone COMPUTE: 1. total maximum transient storage 2. transient storage of each zone 3. governing factor C CD-CD-/ INPUTS: 1. vesther station number 2. year, month, day 3. precipitation PRE 4. potential evapotranspiration PE t YES ^ PRINT OUT: 1. day, month, year 2. PE, AE, PRE, DP 3. RO, R, TR 4. SMC1, SMC2, SMC3, SMC4 5. ATR, az, Z, H 1 COMPUTE: 1. actual evapotranspiration AE 2. SMC1, SMC2, SMC3, SMC4 COMPUTE: 1. actual evapotranspiration AE 2. SMC1, SMC2, SMC3, SMC4 3. excess EX COMPUTE: 1. actual evapotranspiration AE 2. assign SMC as previous day's values © COMPUTE: 1. surface runoff 2. subsurface drainage rate 1 COMPUTE: 1. transient storage 2. change in transient storage 3. f(K) based on previous day's Z(K) 4. f i r s t water table depth Z'(I) 5. f'(I) based on Z'(I) 6. averaged drainable porosity f(I) 7. change In water table depth 8. final water table depth for the current day Z(I) © 00 G 5 00 O c 00 oo O z FIGURE 14: Generalized Flowchart of the Water Balance Model RESULTS AND DISCUSSION / 54 1.4 i 0 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (doy) FIGURE 15: Actual Water Table Curves For Regimes A, B, C, D of 1984 '1 1.2 J i 1 1 1 1 1 1 1 1 1 1 1 1 1 0 30 60 90 120 150 180 210 240 270 3 0 0 3 3 0 3 6 0 3 9 0 TIME (day) FIGURE 16: Actual Water Table Curves For Regimes A, B, C, D of 1985 RESULTS A N D DISCUSSION / 55 x UJ I UJ i m cn UJ -0.6 Legend A ACTUAL X MODEL l 1 1 1 1 1 1 1 1 1 1 1 1 0 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 17: Laboratory-Based Water Table Simulation For Regime A of 1984 FIGURE 18: Laboratory-Based Water Table Simulation For Regime B of 1984 RESULTS AND DISCUSSION / 56 -0.6 Legend A ACTUAL X MODEL - 1 1 1 1 1 1 1 1 1 1 1 1 1 0 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 19: Laboratory-Based Water Table Simulation For Regime C of 1984 x o CD i< IX U J Legend A ACTUAL X MODEL —i 1 i 1 1 1 r 1 1 1 1 1 1 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 20: Laboratory-Based Water Table Simulation For Regime D of 1984 RESULTS AND DISCUSSION / 57 x UJ x UJ _i CO oc UJ - | 1 1 1 1 1 1 1 1 1 r~. 1 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 21: Laboratory-Based Water Table Simulation For Regime A of 1985 1.4-j 1.2-Legend A ACTUAL X MODEL -1.2--1.4-H 1 1 1 1 1 1 1 1 1 1 1 1 1 ,0 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 22: Laboratory-Based Water Table Simulation For Regime B of 1985 RESULTS AND DISCUSSION / 58 1-3-1 Legend A ACTUAL X MODEL § -0.9--1.1--1.3--1-5 H 1 1 1 1 1 1 1 1 1 1 j- 1 1 0 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 23: Laboratory-Based Water Table Simulation For Regime C of 1985 FIGURE 24: Laboratory-Based Water Table Simulation For Regime D of 1985 RESULTS AND DISCUSSION / 59 x o m UJ Legend A ACTUAL X MODEL —1 1 1 1 1 1 1 1 1 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 25: Drainage-Rate-Based Water Table Simulation For Regime A of 1984 FIGURE 26: Drainage-Rate-Based Water Table Simulation For Regime B of 1984 RESULTS AND DISCUSSION / 60 x UJ n UJ i m ,< or -T 30 60 T 120 150 180 210 TIME (day) r 240 ~i r 270 300 330 360 FIGURE 27: Drainage-Rate-Based Water Table Simulation For Regime A of 1985 1.4-i 1.2-Legend A ACTUAL X MODEL -1.2--1.4--1.6 H 1 1 1 1 1 1 1 1 1 1 1 1 1 0 30 60 90 120 150 180 210 240 270 300 330 360 390 TIME (day) FIGURE 28: Drainage-Rate-Based Water Table Simulation For Regime B of 1985 RESULTS A N D DISCUSSION / 61 observed water table levels as shown in Figures 15 and 16. In general, the model simulated the trend of the actual water table in that peaks and valleys did coincide. This was especially true with Regimes B and C for both years, and with Regime D for the earlier part of both years. In the latter part of the year with Regime D, the trend was present but at a consistently lower depth. With Regime A, the model tended to overshoot with peaks and undershoot with valleys before June but simulated well for the last three months of the year. The water table was observed to fluctuate most pronouncely with Regime A. Referring to Table 13, the model confirmed the apparent effects of the different water control regimes. The R M A X values decreased in the order of Regimes A, B, C and D. This followed correctly the respective water regimes of "drainage", "low subirrigation and drainage", "high subirrigation and drainage", "no drainage and no subirrigation". The different drainable porosity characterizing each regime was also made evident in Figures 29 and 30. In both years, the water table heights of Regime A were consistently lower than those of Regimes B and C. Due to the high subirrigation operation, Regime C exhibited the highest levels. These trends were agreeable with those of Figures 15 and 16 which showed the actual field water table heights. Another observation made was the similarity of R M A X between laboratory-corrected simulations and drainage-rate-corrected simulations. Thus, it is possible to work with field data of rainfall and drainage rates in place of laboratory core samples to derive the drainable porosity equation. Application of the Wilcoxon's paired-sample test was attempted. However due to the nature of the test, what was visually a good simulation needed not to be a statistically significant agreement Since both direction and magnitude of the difference between the paired observations were considered, the test was not able to optimize the RESULTS A N D DISCUSSION / 62 Table XIII: Maximum Design Drainage Rate R M A X R M A X (mm/d) M E T H O D R E G I M E 1984 1985 A 4.5 3.0 B 2.5 2.0 Lab. Based C 2.5 2.0 D 0.6 1.5 Drain. A 4.0 3.0 Rate Based B 3.0 2.5 absolute minimum differences summed over the analysis period. Thus a large overshoot could be compensated by a large undershoot occurring anywhere within the one year period. At the end. the model simulations were tested visually. In being able to reasonably simulate the water table, the usefulness of the model is in its relationship to the SEW 3 0 criterion. For this research field and the Lower Fraser Valley, the sum of daily water table height recorded above the depth of 30cm at the midpoint between drains is not to exceed 200cm (SEW J 0<200cm) during November 1st to March 31st A simple accompanying computer program could conveniently compute the SEW 3 0 value after each simulation. Furthermore, the graphical plots of the simulation could provide the information of approximately when the water table is at or above the 30cm depth during the year. RESULTS AND DISCUSSION / 63 Legend & R E G I M E A X R E G I M E B • R E C I M E C 1 1 1 1 1 1 1 1 1 1 1 1 1 0 30 60 90 120 150 180 210 240 270 300 330 3 6 0 3 9 0 TIME (day) FIGURE 29: Water Table Simulation For Regimes A, B and C of 1984 (Based on Laboratory Method) 1 . 4 - i 1.2-0 30 60 90 120 150 180 210 240 2 7 0 3 0 0 3 3 0 3 6 0 390 TIME (day) : FIGURE 30: Water Table Simulation For Regimes A, B and C of 1985 (Based on Laboratory Method) 5. SUMMARY AND CONCLUSIONS 1. The experimental field at Boundary Bay, Delta, was placed under four types of water control regimes: drainage, low subirrigation with drainage, high subirrigation with drainage, and no subirrigation with no drainage (Regimes A, B, C, D, respectively). Physical and hydraulic soil properties of each treatment plots were analyzed using undisturbed core samples and loose soil samples. All the regimes were classified as silt loam with high percentage content of silt The bulk density ranged from 1.29 g/cm3 to 1.46 g/cm3 which agreed well with the results of other independent studies on the same field. The satiated hydraulic conductivity of Regime A was the highest followed by those of Regimes B, C and D. Regime A also exhibited the lowest water retention curve among the four treatments. 2. Drainable porosities derived from laboratory analysis showed clearly the negative effect of subirrigation on this property. The approximate drainable porosity at 120cm tension for Regimes A, B, C and D were 9, 7, 5 and 9% by volume, respectively. The decrease in drainable porosity was probably attributed to the transport and deposition of fine soil particles brought about by subirrigation. The reduced pore spaces led to the reduction in drainable porosity. 3. Drainable porosities determined from the recorded water table charts and measured drainage rates were generally higher than those determined in the laboratory. This conclusion, however, was not substantiated by an adequate sample size. 4. Drainable porosities determined from water table charts and rainfall data were generally lower than those determined by either laboratory or drainage rate method. This was especially true for the tension range of 0cm to 120cm. Since the drainable porosity is a cumulative value affected by the direction of flow, these observations are reasonable. In the rainfall rate analysis, the water table was rising, 64 S U M M A R Y A N D CONCLUSIONS / 65 moving through a profile of smaller to larger pores. The incremental change in drainable porosity increased with the height of the water table and was greatest near the soil surface. Conversely with the drainage rate analysis, the water table was receding from the surface, moving through larger to smaller pores. The incremental change in drainable porosity decreased progressively with water table depth. However, the initial drainable porosities near the surface depth was always greater than those derived by a rising water table reaching the soil surface. 5. Drainable porosities derived in the laboratory were used for the curve fitting procedure. The negative exponential form was chosen over the polynomial form since it was simple in expression and was inherently asymptotic. The polynomial fit was equally acceptable but only when fitted to the sixth power. 6. Test for variability was performed on the two parameters, b 0 and bi, to the negative exponential drainable porosity equation: f = b„( 1 - e " b l Z ) (5.1) b 0 represented the maximum drainable porosity while bi represented the rate constant No significant difference was found for bj among treatments. For b 0, a significant difference existed between Regimes C and D. At this point in the research, only the apparent trends indicating the effects of different water treatments were reportable. 7. Since an underlying objective of the study was to utilize readily available meteorological data in defining drainable porosity, the drainable porosity function based on rainfall analysis was corrected by those based on laboratory analysis and drainage rate analysis. The corrected f-z expression was of the general mathematical form of: Corrected = b0( 1 - e " b l Z ) / [ b,( 1 - e " b 3 Z ) + b 4 ] (5.2) 8. Equation (5.2) was incorporated into the computer water balance model by a subroutine program, DPORE(....). In the simulation process, Z(I) itself is needed in S U M M A R Y A N D CONCLUSIONS / 66 advance for use in the f-z equation. This problem was solved by using Z(K), the previous day's water table depth, to obtain f(K) which in turn gave a first approximated Z'(I) value. A second f(I) was obtained from Z'(I). The final f(I) was derived from the average of f(K) and f(I). This f(I) value was used to compute the present day's water table depth, Z(I). 9. The results generated by the modefied water balance model were visually tested for good agreement with actual water table depths. Although weather conditions for 1984 and 1985 were quite different, the model did follow the trend of actual water table patterns. Better simulation was found for Regimes B, C and D over Regime A. Water table fluctuation was most pronounced for Regime A such that the model tended to correspondingly overshoot or undershoot 10. The Wilcoxon's paired-sample test was found to be inappropriate for this study in comparing simulated and actual water table levels. It did not optimize the absolute minimum differences summed over the analysis period. 11. One of the additional usefulness of this water balance model lies with the SEW 3 0 information that it is able to provide. Furthermore, the graphical representation of its simulation results easily shows when the water table is likely to be at or above this critical depth of 30cm given a wet or a dry year. 6. RECOMMENDATIONS As a result of this study, the following recommendations are suggested for further investigation: 1. To determine the variability of drainable porosity within a regime, drainable porosity curves obtained from laboratory analysis on core samples extracted from varying depths are required. 2. Further observation of the effects the four types of water regimes have on drainable porosity is needed along with obtaining information of possible statistical significance on the negative exponential parameters. 3. Determination of the applicability of the negative exponential function to other soil types, and, if so, the determination of the parameters for these soil types. 4. Continued drainage rate data collection is needed to supplement this part of the drainable porosity analysis. 5. A more suitable statistical analysis is necessary to test the accuracy with which the water balance model is able to simulate water table depth. 6. The implementation of field trials of the drainage coefficient determined from model simulation should be performed to confirm the validity of RMAX and thus the practicality of the model itself. 67 REFERENCES Anderson, C.E., H.P. Johnson and W.L. Powers. 1978. A Water Balance Model for Deep Loess Soils. Transactions of ASAE. Vol.21(2):314-320. Baier, W. and G.W. Robertson. 1966. A New Versatile Soil Moisture Budget Can. J. Plant Science, Vol.46:299-315. Bhattacharya, A .K. and R.S. Broughton. 1979. Variable Drainable Porosity in Drainage Design. J. Irrigation and Drainage Division, ASCE. Vol.l05(IRl):71-86. Bhattacharya, A.K. , N . Foroud, S.T. Chieng and R.S. Broughton. 1977. Subsurface Drainage Cost and Hydrologic Model. J. Irrigation and Drainage Division, ASCE. Vol.l03(IR3):299-308. Black, C.A. (Ed.) 1965. Methods of Soil Analysis. Part I. Physical and mineralogical properties, including statistics of measurement and sampling. Agronomy(9). Am. Soc. Agron. Inc., Madison,Wisconson. Blake, G.R. 1965. Bulk Density. In: C.A. Black (Ed.) Methods of Soil Analysis, Part I. Agronomy (9): 374-390. Am. Soc. of Agronomy. Inc. Madison, Wisconson. Bloemen, G.W. 1980. Calculation of Steady State Capillary Rise from the Groundwater Table in Multi-Layered Soil Profiles. Z. Pflanzenernaehr. Bodenkd. Vol.l43:701-718. Bornstein, J. and W.E. Hedstrom. 1981. Trafficability Factor in Efficient Forage Production. ASAE Paper No.81-2068. Bouwer, H. and J. van Schilfgaarde. 1963. Simplified Method of Predicting Fall of Water Table in Drained Land. Transactions of ASAE 6(4): 288-296. Brooks, R.H. and A.T. Corey. 1964. Hydraualic Properties of Porous Media. Colorado State University. Hydrology Paper No.3. Fort Collins, Colorado State University. Colorado. 27pp. Broughton, R.S. and N . Foroud. 1978. A Model to Predict Water Table Depths for Flatlands. Canadaian Agri. Eng. Vol.20(20:81-86. 68 / 69 Campbell, K . L , J.S. Rogers and R R . Hensel. 1978. Water Table Control for Potatoes in FLorida. Transactions of ASAE. Vol.21)40:701-705. Carter, C.E. and C R . Camp. 1978. Subsurface Drainage of an Alluvial Clay Soil For Soybeans. ASAE Paper No.78-2040. Carter, C.E. and C R . Camp. 1981. Subsurface Drainage of an Alluvial Soil Increased Sugar Cane Yield. ASAE Paper No.81-2536. Carter, C.E. and C R . Camp. 1983. Subsurface Drainage of an Alluvial Soil Increased Sugarcane Yield. Transactions of ASAE. Vol.l6(2):426-429. Carter, C.E., J.E. Irvine and C R . Camp. 1983. A System for Managing Water in Fine Textured Soils. Transactions of ASAE. Vol.l6(2):767-771. Chieng, S.T. 1975. The Effects of Subsurface Drain Depths and Drainage Rates on Water Table Levels. Thesis presented to McGil l University, at Montreal, Canada in 1975, in partial fulfillment of the requirements for the degree of Master of Science. Chieng, S.T. 1986. A D A L D : An Interim Report to CIDA. Chieng, S.T. 1987. Personal communications. Bio-Resource Engineering Dept University of British Columbia. British Columbia, Canada. Chieng, S.T., R.S. Broughton and N . Foroud. 1978. Drainage Rates and Water Table Depths. J. Irrigation and Drainage Division, ASCE. Vol.l04(IR4): 413-433. Chieng, S.T., J. Keng and M . G . Driehuyzen. 1987. Effects of Subsurface Drainage and Subirrigation on the Yields of Four Crops. Canadian Agri. Eng. Vol.29(l):21-26. Childs, E.C. 1957. The Physics of Land Drainage. Agronomy Monograph No.7:l-78. J.N. Luthin (Ed). American Soc. Agronomy. Childs, E.C. 1960. The Nonsteady State of the Water Table in Drained Land. J.Geophysical Research. Vol.65(2):780-782. Clark, F.E. and W.D. Kemper. 1967. Microbial Activity in Relation to Soil, Water and / 70 Soil Aeration. In Irrigation of Agricultural Lands. R .M. Hagen, H.R. Haise, T.W. Edminister (Eds.). Agronomy 11. American Soc. Agronamy. Day, P.R. 1965. Particle Fractionation and Particle Size Analysis. In: C.A. Black (Ed). Methods of Soil Analysis. I. Physical and minerological properties, including statistics of measurement and sampling. Agronomy (9): 545-566. DeVries, D.A. 1975. Heat Transfer in Soils. In: Heat and Mass Transfer in the Biosphere. D.A. DeVries, N .H. Afgen (Eds.). Washington, D.C.:Scripta Book Co. DeVries, J. 1985. Personal communications. Soil Science Department, University of British Columbia, B.C., Canada. Doering, E.J. 1965. Change of Permeability With Time. In: Drainage For Efficient Crop Production. Conference Proceedings. ASAE Publication. Donnan, W.W. and G.O. Schwab. 1974. Current Drainage Methods in the U.S.A. In: Drainage For Agriculture. Jan Van Schilfgaarde. (Ed.). Agronomy Monograph 17. American Soc. Agronomy, pp.93-114. Doty, C W . and G.D. Christenbury. 1979. Controlled and Reversible Drainage Past, Present, and Future. ASAE Paper No.79-2545. Doty, C.W., S.T. Currin and R.E. McLin. 1975. Controlled Subsurface Drainage For Southern Coastal Plains Soil. J. Soil and Water Conservation. Vol.30(2):82-87. Doty, C W . and J.E. Parsons. 1979. Water Requirements and Water Table Variations for a Controlled and Reversible Drainage System. Transactions of ASAE. Vol.22(3):532-536. Driehuyzen, M . G . 1983. Boundary Bay Water Control Project December, 1983. B.C. Ministry of Agriculture and Food. Cloverdale, B.C. Duke, H.R. 1972. Capillary Properties of Soils - Influence Upon Specific Yield. Transactions of ASAE. Vol.l5(4):688-691. Erickson, A.E. 1965. Short Term Oxygen Deficiencies and Plant Responses. In: Drainage / 71 for Efficient Crop Production Conference Proceedings. ASAE Publications. Fausey, N.R. and R.D. Brehn. 1976. Shallow Subsurface Drainage Field Perfomance. Transactions of ASAE. Vol.l9(6): 1082-1084, 1088. Fausey, N.R. and G.O. Schwab. 1969. Soil Moisture Content, Tilth, and Soybean (Glycine Max) Response with Surface and Subsurface Drainage. J. Agronomy. Vol.61:554-557. Foroud, N . 1974. A Hydrologic Study for Determining Subsurface Drainage Coefficients for the S L Lawrence Lowland Region. Thesis presented to McGil l University, at Montreal, Canada, in 1974, in partial fulfillment of the requirements for the degree of Master of Science. Foroud, N . and R.S. Broughton. 1978. Effects of Different Drainage Rates on the Duration of High Water Tables in South Western Quebec. Canadian Agric. Eng. Vol.20(2):71-75. Fox, D.J. and K . E Guire. 1976. Documentation For MIDAS. Statistical Research Laboratory. The University of Michigan. Fox, R.L., J.T. Phelan and W.D. Criddle. 1956. Design of Subirrigation Systems. Agricultural Engineering. Vol.37(2):103-108. Gardner, W.H. 1965. Water Content In: C.A. Black (Ed.). Methods of Soil Analysis. I. Physical and minerological properties, including statistics of measurement and sampling. Agronomy (9):83-125. American Soc. Agronomy. Goldstein, M.J. 1978. Evaluation and Prediction of Subirrigation Using a Subsurface Drainage System. Thesis submitted to the University of British Columbia, Vancouver, B.C., Canada, in partial fulfillment of the requirements for the degree of Bachelor of Science in Agricultural Sciences in the Dept of Soil Science. Harris, D.G. and C.H. Van Barel. 1967. Root Respiration of Tobacco, Corn and Cotton Plants. J. Agronomy. Vol.49:182-184. Hillel, D. 1980(a). Fundamentals of Soil Physics. New York: Academic Press Inc. Ltd. / 72 Hillel, D. 1980(b). Applications of Soil Physics. New York: Academic Press Inc. Ltd. Hillel, D.I., C.H.M. Van Barel and H . Talpaz. 1975. Dynamic Simulation of Water Sorage in Fallow Soil asAffected by Mulch of Hydrophobic Aggregates. Soil Sc. Soc. America Proceedings. Vol.39:826-833. Hundel, S.S., G.O. Schwab and G.S. Taylor. 1976. Drainage System Effects on Physical Properties of a Lakebed Clay Soil. American J. Soil Sc. Vol.40:300-305. Jensen, M.E., J.L. Wright and B.J. Pratt. 1971. Estimating Soil Moisture Depletion from Climate, Crop and Soil Data. Transactions of ASAE. Vol.l4(5): 954-259. Kerr-Wilson, G . 1985. A Study of the Characteristics of Subsurface Drain Flow. Thesis submitted in partial fulfillment of the requirements for the degree of Bachelor of Applied Science in Bio-Resource Engineering, University of British Columbia, Vancouver, B.C., Canada. Klute, A. 1965. Laboratory Measurement of Hydraulic Conductivity of Saturated Soil. In: C.A. Black (Ed.). Methods of Soil Analysis. I. Physical and minerological properties, including statistics of measurement and sampling. Agronomy (9):210-220. American Soc. Agronomy. Kodsi, E 1987. Influence of Sewage Sludge Appliction on Hydraulic and Physical Properties of a Ladner South Silty Clay Loam Soil. A thesis submitted in partial fulfillment of the degree of Master of Science in Bio-Resource Engineering. University of British Columbia, Vancouver, B.C., Canada. Lagace, R. and A. St-Yves-Desilets. 1979. A Water Table and Water Balance Model for the Quebec's Conditions. ASAE Paper No.79-2072. Leyton, L. and J.S.P. Yadar. 1960. Effect of Drainage on Certain Physical Properties of a Heavy Clay Soil. J. Soil Sc. Vol.ll(2):305-313. Luthin, J.N. 1978. Drainage Engineering. Huntington, N.Y. : Robert E. Krieger Publishing Company. Massey, F.C., R.W. Skaggs and R . E Sneed. 1983. Energy and Water Requirements for / 73 Subirrigation vs Sprinkler Irrigation. Transactions of ASAE. Vol.26(l):126-133. McWhorten, D.B. and H.R. Duke. 1976. Transient Drainage with Nonlinearity and Capillarity. J. Irrigation and Drainage Division, A S C E Vol.l02(IR2): 193-204. Paul, C L . and J. DeVries. 1979. Effect of Soil Water Status and Strength on Trafficability. Canadian J. Soil Sc. Vol.59:313-324. Paul, C L . and J. DeVries. 1983. Soil Trafficability in Spring. 2. Prediction and the Effect of Subsurface Drainage. Canadian J. Soil Sc. Vol.63:27-35. Prasher, S.O. 1982. Examination of the Design Procedure for Drainage/Subirrigation Systems in the Lower Fraser Valley, British Columbia. Thesis presented to the University of British Columbia at British Columbia, Canada, in 1974, in partial fulfillment of the requirements for the degree of Doctor of Philosophy. Prasher, S.O., S.O. Russell, S.T. Chieng and T.H. Podmore. 1985. A Probabilistic Approach to Agricultural Drainage Design. Transactions of the A S A E Vol.28(5): 1494-1498. Reeve, R.C. and N.R. Fausey. 1974. Drainage and Timeliness of Farming Operations. In: Jan Van Schilfgaarde (Ed.). Drainage for Agriculture. Agronomy Monograph 17. American Soc. Agronomy. Ray, A.A. (Ed) 1982. SAS User's Guide: Statistics. SAS Institute Inc. Cary, N . C : SAS Institute Inc. Richards, L.A. 1965. Physical Conditions of Water in Soil. In: C A . Black (Ed.). Methods of Soil Analysis. I. Physical and minerological properities, including statistics of measurement and sampling. Agronomy (9): 128-157. American Soc. Agronomy. Russelo, D., S. Edey and J. Godfrey. (Eds.) 1974. Selected Tables and Conversions Used in Agrometeorology and Related Fields. Canada Department of Agriculture. Publication 1522, 1974. Salisbury, F.B. and C.W. Ross. 1979. Plant Physiology. 2nd Edition. Belmont, / 74 California: Wadsworth Publication Company Inc. Schwab, G.O., N.R. Fausey and D.W. Michener. 1974. Comparision of Drainage Methods in a Heavy Textured Soil. Transactions of A S A E Vol.l7(3): 424-425. Schwab, G.O., G.S. Taylor, J.L. Fouss and E. Stibbe. 1966. Crop Response from Tile and Subsurface Drainage. American Soc. Soil Sc. Proceedings. Vol.30:634-637. Shih, S.F. 1983. Soil Surface Evaporation and Water Table Depth. J. Irrigation and Drainage Division, ASCE. Vol.l09(IR4):366-376. Sieben, W.H. 1965. From: G.A.W. Van de Goor. 1972. Plant Growth in Relation to Drainage. In: Drainage Principles and Applications I. Introductory Subjects. Publication 16, Vol 1. International Institute For Land Reclamation and Improvement Wageningen, Netherlands. Siegel, S. 1965. Nonparametric Statistics For the Behavioral Sciences. N.Y.: McGraw-Hill Books Company. Skaggs, R.W. 1976. Determination of Hydraulic Conductivity-Drainable Porosity Ratio from Water Table Measurements. Transactions of ASAE. Vol.l9(l):71-86. Skaggs, R.W. 1978. A Water Management Model for Shallow Water Table Soils. Technical Report No. 134. Water Resources Research Institute of the Univ. of N . Carolina. N.C. State University, Raleigh, N . Carolina. Skaggs, R.W. 1981. Water Movement Factors Important to the Design and Operation of Subirrigation Systems. Transactions of ASAE. Vol.24(6): 1553-1561. Spangler, M.G . 1960. Soil Engineering. 2nd Edition. Scranton, Pennsylvania: International Textbook Company. Steinhardt R. and B.D. Trafford. 1974. Some Effects of Subsurface Drainage and Ploughing on the Structure and Compatibility of a Clay Soil. J. Soil Sc. Vol.25(2):138-152. Stuff, R.G. and R.F. Dale. 1978. A Soil Moisture Budget Model Accounting for / 75 Shallow Water Table Influences. J. American Soc. Soil Sc. Vol.42(4): 637-644. Taylor, G.S. 1960. Drainable Porosity Evaluation From Outflow Measurements and Its Use in Drawdown Equations. Soil Sc. Vol.90:338-343. Thomas, A.W., H.R. Duke, and E G . Kruse. 1977. Capillary Potential Distribution in Root Zones Using Subsurface Irrigation. Transactions of ASAE. Vol.20(l):62-67,75. Toksoz, S. and D. Kirkham. 1971. Steady Drainage of Layered Soils I. Theory. J. Irrigation and Drainage Division of ASCE. Vol.97(IRl):l-18. Proceeding Paper No.7985. Unesco/FAO. 1973. Irrigation, Drainage and Salinity. An International Source Book. London: Hutchensen & Co. Ltd. Van der Gulik, T.S. Chieng and F. Smith. 1986. B.C. Agricultural Drainage Manual. B.C. Ministry of Agriculture and Food. Agricultural Engineering Branch. Vomocil, J.A. 1965. Porosity. In: C.A. Black (Ed.). Methods of Soil Analysis. I. Physical and minerological properties, including statistics of measurement and sampling. Agronomy (9): 299-314. American Soc. Agronomy. Wesseling, J. 1974. Crop Growth and Wet Soils. In: Drainage for Agriculture. Jan Van Schilfgaarde (Ed.). Monograph 17. American Soc. Agronomy. Wesseling, J. and W.R. Van Wijle. 1957. Soil Physical Conditions in Relation to Drain Depth. In. Drainage of Agricultural Lands. J.N. Luthin (Ed.). American Soc. Agronomy Publisher. Yang, S.J., P.L. Wang, Y.T. Chang and T.C. Chang. 1977. The Performance of Subsurface Drainage and Its Effect on Sugarcane Grown in Fine Textured Soil in Taiwan. Taiwan Sugar. Vol24(4):378-384. Young, T.C. and J.T. Lignon. 1972. Water Table and Soil Moisture Probabilities with Tile Drainage. Transactions of ASAE. Vol.l5(3): 488-451. APPENDIX The water table model is accessible on the MTS System at the University of British Columbia from the Department of Bio-Resource Engineering. 1 (;*•»•»«***»»***•***»»*****»«*»*•**»•*««***»*«*»*»»*•«**»•»»*•»»*»»««»»*»»»***»»* 2 C MAIN PROGRAM : MOD.DRAIN 3 A H 5 c ACORR O MAXIMUM POROSITY OF CORRECTION EQUATION. 6 c AE n ACTUAL EVAPOTRANSPIRATION. 7 c ALAB s MAXIMUM POROSITY OF LABORATORY f EQUATION. 8 c AFIELD = MAXIMUM POROSITY OF IN SITU f EQUATION. 9 c C GOVERNMENT FACTOR = 4K/L**2 10 c CI a VALUE OF AE/PE FOR THE 1ST ZONE. 11 c C2 s VALUE OF AE/PE FOR THE 2ND ZONE. 12 c C3 m VALUE OF AE/PE FOR THE 3RD ZONE. 13 c C4 = VALUE OF AE/PE FOR THE 4TH ZONE. 14 c D - DEPTH OF IMPERMEABLE LAYER FROM THE DRAIN CENTER IS c D1 B DEPTH FROM SOIL SURFACE OF FIRST SOIL ZONE. 16 c DP.DPK.DPI » DRAINABLE POROSITY. 17 c DPORE S SUBROUTINE TO COMPUTE DRAINABLE POROSITY. 18 c DP4 B SUBROUTINE TO ASSIGN f (z IS AT DRAIN AND TR IS POSITIVE). 19 c DELZ B CHANGE IN WATER TABLE DEPTH. 20 c DELTR S CHANGE IN TRANSIENT STORAGE. 21 c DW a DEPTH OF ZONE FROM SOIL SURFACE THAT RMAX OCCURS. 22 c FACTOR a SUBROUTINE TO GENERATE PARAMETER INPUT TABLE. 23 c H B WATER TABLE HEIGHT ABOVE DRAIN CENTER AT MIDSPACING. 24 c HMAX = MAXIMUM DRAIN DEPTH ( TOTAL DRAIN DEPTH ). 25 c I = COUNTER FOR CURRENT DAY. 26 c IDAY s DAY. 27 c I END = ENDING YEAR. 28 c IMONTH = MONTH. 29 c IRD = INPUT DEVICE (9). 30 c ISTART s STARTING YEAR. 31 c OYEAR s YEAR. 32 c K = COUNTER FOR PREVIOUS DAY. 33 c KCORR 3 RATE CONSTANT OF CORRECTION EQUATION. 34 c KLAB B RATE CONSTANT OF LABORATORY f EQUATION. 35 c KFIELD B RATE CONSTANT OF IN SITU f EQUATION. 36 c MEND S ENDING MONTH. 37 c MSTART S STARTING MONTH. 38 c OPE B PRECIPITATION READ FROM THE TAPE. 39 c OPR S POTENTIAL EVAPOTRANSPIRATION READ FROM THE TAPE. 40 c PE=PET = POTENTIAL EVAPOTRANSPIRATION. 41 c PRE=PPT = PRECIPITATION. 42 c REM B REMAINDER OR EXCESS. 43 c RMAX B ALLOWABLE MAXIMUM DRAINAGE-RATE (DESIGN FLOW-RATE). 44 c RO a SURFACE RUNOFF. 45 c SD B SUBSURFACE DRAINAGE-RATE. 46 c SMCF B AVAILABLE SOIL MOISTURE AT FIELD CAPACITY OF WHOLE PROFILE. 47 c SZMCF1 S AVAILABLE SOIL MOISTURE AT FIELD CAPACITY OF 1ST ZONE. 48 c SZMCF2 8 AVAILABLE SOIL MOISTURE AT FIELD CAPACITY OF 2ND ZONE. 49 c SZMCF3 m AVAILABLE SOIL MOISTURE AT FIELD CAPACITY OF 3RD ZONE. 50 c SZMCF4 B AVAILABLE SOIL MOISTURE AT FIELD CAPACITY OF 4TH ZONE. 51 c SMC1 B AVAILABLE SOIL MOISTURE OF THE 1ST ZONE. 52 c SMC2 S AVAILABLE SOIL MOISTURE OF THE 2ND ZONE. 53 c SMC 3 B AVAILABLE SOIL MOISTURE OF THE 3RD ZONE. 54 c SMC4 S AVAILABLE SOIL MOISTURE OF THE 4TH ZONE. 55 c STA S WEATHER STATION NAME. 56 c STAT • S WEATHER STATION READ FROM TAPE. 57 c TITLE 1 s SUBROUTINE TO GENERATE HEADING TO EACH PAGE OF*OUTPUT. 58 c TRANS B SUBROUTINE TO COMPUTE TOTAL & ZONAL TRANSIENT STORAGE. 76 / 77 59 C TRMAX «= MAXIMUM TRANSIENT STORAGE FOR WHOLE PROFILE. 60 C TR1.TRISUM » MAXIMUM TRANSIENT STORAGE FOR 1ST ZONE. 61 C TR2,TR2SUM « MAXIMUM TRANSIENT STORAGE FOR 2ND ZONE. 62 C TR3,TR3SUM = MAXIMUM TRANSIENT STORAGE FOR 3RD ZONE. 63 C TR4,TR4SUM » MAXIMUM TRANSIENT STORAGE FOR 4TH ZONE. 64 C WTH = MAX. HEIGHT FROM DRAIN CENTER THAT HOOGHOUT'S EQUATION APPLIES. 65 C Z « WATER TABLE DEPTH FROM SOIL SURFACE. 66 C 67 C TABLES : 68 C 1. DISTRIBUTION OF DAILY WATER TABLE DEPTH. 69 C 2. MONTHLY SUMMARY WATER TABLE DEPTH DISTRIBUTION FOR EACH 70 C WATER BALANCE YEAR. 71 C 3. MONTHLY SUMMARY WATER TABLE DEPTH DISTRIBUTION FOR N YEARS. 72 C 4. MONTHLY FREQUENCY DISTRIBUTION OF WATER TABLE DEPTH FOR N YEARS. 73 C 5. SEASONAL FREQUENCY DISTRIBUTION TABLE. 74 C 6. SUMMARY OF FREQUENCY DISTRIBUTION FOR N YEARS. 75 C 76 C NOTE : N IS THE TOTAL NUMBER OF YEARS FOR THAT RUN. 77 C 78 C 79 C MISSING = WHEN " MISSING " APPEARS, THIS MEANS THE DAILY DATA FOR 80 C PRECIPITATION OR/AND PE WAS/WERE NOT AVAILABLE FOR THESE DAYS. 81 C CALCULATIONS WILL NOT BE DONE FOR THESE DAYS BUT CONTINUE RIGHT 82 C AFTER MISSING DAYS WHEN DAILY DATA IS AVAILABLE BY USING THE 83 C INFORMATION FROM THE DAY BEFORE THE MISSING DAYS. 84 C 85 C ICOND = COUNTER FOR DAILY WATER TABLE DISTRIBUTION. 86 C ICONT = COUNTER FOR MONTHLY SUMMARY OF WATER TABLE DEPTH. 87 C IFDOWT = COUNTER FOR FREQUENCY DISTRIBUTION OF WATER TABLE DEPTH. 88 C IC0NF1 = COUNTER FOR TOTAL FREQUENCY OF WATER TABLE DEPTH. 89 C IC0NF2 = COUNTER FOR SEASONAL FREQUENCY DISTRIBUTION OF WATER TABLE 90 C DEPTH ( APRIL-MAY ). 91 C IC0NF3 = COUNTER FOR SEASONAL FREQUENCY DISTRIBUTION OF WATER TABLE 92 C DEPTH (JUNE-SEPTEMBER ). 93 C IC0NF4 = COUNTER FOR SEASONAL FREQUENCY DISTRIBUTION OF WATER TABLE 94 C DEPTH (OCTOBER-NOVEMBER). 95 C IC0NF5 = COUNTER FOR MONTHLY FREQUENCY DISTRIBUTION OF WATER TABLE 96 C DEPTH FOR N YEARS. 97 C 98 C NOTE : WHEN THE COUNTER IS ASSIGNED TO 0, NO OUTPUT OF THE 99 C APPROPRIATE TABLES ARE WANTED. IF THE COUNTER IS ASSIGNED 100 C TO 1, THEN THE APPROPRIATE TABLES WILL BE PRINTED. 101 C 102 DIMENSION MONTH(380).PRE(380),PE(380),AE(380),REM(380),DP(380), 103 * R0(380),TR(380),SMC 1(380),SMC2(380),SMC3(380).DELTR(380),Z(380), 104 * DELZ(380),SD(380),IDAY(380),KM0NT(12,5),DEP(10),IYEAR(380), 105 * MAT(150.10),MJAF(62,10),MANF(124.10),M4F(31.10),M5F(31,10), 106 * MATT(12.14) , ITOTAL(12. 14) ,ITOTF(380,14),H(380).MMJF(62.10). 107 * M6F(31.10),M7F(31,10),M8F(31,10),M9F(31,10),M10F(31.10), 108 * M11F(31.10),PERI0D(5.18) ,STA(10).SOIL(10),MNDF(62,10),MS0F(62,10) 109 *,M1F(31.10),M2F(31,10),M3F(31,10),M12F(31,10) 110 REAL KCORR,KLAB.KFIELD 111 INTEGER STAT 112 READ(5.10) IRD,ISTART,IEND,MSTART.MEND 113 10 F0RMAT(5I2) 1 14 READ(5, 11) RMAX,HMAX,D,DW,D1,SMCF.SZMCF1,SZMCF2.SZMCF3,SZMCF4 115 11 FORMAT(F5.1.2F5.0.7F5.1) 1 16 READ(5, 12) ALAB,KLAB,AFI ELD,KFI ELD.ACORR.KCORR,CONST,C1.C2.C3.C4 / 78 117 12 FORMAT(3(F8.4,F8.5),F6.3,4F5.2) 1 18 READ(5,13)( (KMONT(I .J) ,J-1,5) , I - 1.12) 1 19 13 F0RMAT(8(5A2)) 120 READ(5,14)((PERI0D(I,d),J=1.18),I=1,5) 121 14 F0RMAT(3(18A1,2X)) 122 READ(5,15) (DEP(I).1=1,10) 123 15 FORMAT(10F5.0) 124 READ(5.20)ICOND.ICONT,IFDOWT,ICONF1,IC0NF2,IC0NF3,IC0NF4.IC0NF5 125 *,IC0NF6.IC0NF7 126 20 FORMAT(1012) 127 WTH=HMAX-DW 128 C=RMAX/(2*D*WTH+WTH**2) 129 CALL TRANST(D1,KFIELD,KCORR.AFIELD,ACORR,CONST,HMAX, 130 *TR1 SUM,TR2SUM,TR3SUM.TR4SUM,TRMAX) 131 TR1=TR1SUM 132 TR2=TR2SUM 133 TR3=TR3SUM 134 TR4=TR4SUM 135 CALL FACTOR(RMAX,DW,D,HMAX,C,SMCF,SZMCF1.SZMCF2,SZMCF3,SZMCF4, 136 •TRMAX,TR1,TR2.TR3.TR4.C1,C2 , C3,AFIELD,KFIELD,ALAB,KLAB,ACORR, 137 •KCORR,CONST) 138 TEMP 1=0.99 139 TEMP2=99.99 140 DO 25 N=1,366 141 DO 25 M=1,14 142 25 ITOTF(N,M)«0 143 DO 30 N=1,12 144 DO 30 M=1,14 145 IT0TAL(N,M)=O 146 30 MATT(N,M)=0 147 DO 50 1=1.31 148 DO 50 J=1.10 149 M1F(I,J)=0 150 M2F(I,d)=0 151 M3F(I ,<J)=0 152 M4F(I,d)=0 153 M12F(I,J)=0 154 IF(I .EQ.31) GO TO 35 155 35 M5F(I,J)=0 156 IF ( I .EQ.a i ) GO TO 40 157 M6F(I,d)=0 158 40 M7F(I,d)=0 159 M8F(I,J)=0 160 IF(I .EO.31) GO TO 45 161 M9F(I,J)=0 162 45 M10F(I,J)=0 163 IF( I .E0.31) GO TO 50 164 M1 1F(I,J)=0 165 50 CONTINUE 166 DO 55 KJM=1,124 167 DO 55 MLK*1.10 168 55 MANF(KJM.MLK)=0 169 DO 60 IJN=1,62 170 DO 60 NJI=1,10 171 60 MJAF(IdN.NJI)=0 172 DO 65 1=1,62 173 DO 65 J=1,10 174 MNDF(I,J)=0 / 79 175 MSOF(I,J)=0 176 65 MMJF(I.d)=0 177 MAX 1=0 178 MAX2=0 179 MAX3=0 180 MAX4=0 181 MAX5=0 182 MAX6=0 183 IJK=0 184 C 185 1=2 186 READ(5,120) ( STA( IJ) , Id =1,10) 187 90 READ(IRD,100,END=130) STAT,JYEAR,IMONT,JDAY,PPT,PET 188 100 FORMAT(17,312,16X.F4 . 1 , 11X.F4.1 ) 189 110 IF(STAT.EO.9999999.OR.STAT.EO.1111) GOTO 1470 190 120 F0RMATO0A2) 191 IF(JYEAR.LT.ISTART) GO TO 90 192 IF(UYEAR.GT.IEND) GO TO 1470 193 IF(IMONT.LT.MSTART) GO TO 90 194 IF(IMONT.GT.MEND) GO TO 130 195 IF(IMONT.EO.12)G0 TO 131 196 GO TO 132 197 131 IdK=IJK+1 198 MONTH(I)=IMONT 199 IDAY(I)=JDAY 200 PRE(I)=PPT 201 PE(I)=PET 202 1=1+1 203 IF(IJK.EQ.31)G0 TO 130 204 GO TO 90 205 132 CONTINUE 206 MONTH(I)=IMONT .207 IDAY(I)=JDAY 208 PRE(I)=PPT 209 PE(I)=PET 210 1=1+1 211 GO TO 90 212 130 CONTINUE 213 11=1-1 214 IIK=II 215 C 216 140 IF(ICOND.EO.O) GO TO 150 217 CALL TITLE1(STA,JYEAR,RMAX,SMCF,HMAX) 218 150 K=.1 219 ZKEEP=0 220 Z(K)=0.0 221 CALL DPORE(Z(K).HMAX,KFIELD,AFIELD,KCORR.ACORR,CONST,DP(K)) 222 DELTR(K)=0.0 223 DELZ(K)=0.0 224 SMC1(K)=SZMCF1 225 SMC2(K)=SZMCF2 226 SMC3(K)=SZMCF3 227 TR(K)=TRMAX 228 H(K)=HMAX-Z(K) 229 SD(K)=RMAX 230 ICHECK=0 231 IDAYC=1 232 ICOM=0 / 80 233 IF(ICONO.EQ.O) GO TO 170 234 WRITE(6.160) SD(K),TR(K) ,SMC 1(K).SMC2(K),SMC3(K).DELTR(K). 235 1 DELZ(K),Z(K),H(K) 236 160 F0RMAT(58X,5F6.2,4F8.2) 237 170 DO 1020 I-2.11 238 IF(PE(I).EQ.TEMP1.OR.PRE(I).EO.TEMP2) GO TO 980 239 ICOM-0 240 IF(ICHECK.EQ.O) GO TO 180 241 K=dSTAY+1 242 ICHECK-0 243 GO TO 190 244 180 K » I - 1 245 190 IF(PRE(I)-C1*PE(I)) 650.640,200 246 200 REM(I)=SMC1(K)+SMC2(K )+SMC3(K)-SMCF+PRE(I)-C1*PE(I) 247 NA-0 248 IF(REM(I)) 640,210.210 249 210 IF(TR(K)+REM(I)-TR2-TR3-TR4) 560,550.220 250 220 CALL DPORE(Z(K),HMAX.KFIELD,AFIELD,KCORR,ACORR,CONST,DP(K)) 251 DPK-DP(K) 252 SD(I)=RMAX 253 IF(TR(K)+REM(I)-RMAX-TRMAX) 540.530,230 254 230 R0(I)=TR(K)+REM(I)-RMAX-TRMAX 255 240 TR(I)=TRMAX 256 250 IF(NA.EO.1) GO TO 690 257 260 SMC1(I)«SZMCF1 258 270 SMC2(I)«SZMCF2 259 280 SMC3(I )=SZMCF3 260 290 AE(I)=C1*PE(I) 261 300 DELTR(I)=TR(I)-TR(K) 262 C FOR THE CONDITION WHERE THE WATER TABLE IS AT DRAIN LEVEL AND 263 C DELTR(I) IS POSITIVE DUE TO POSITIVE EXCESS, DP(K) IS ASSIGNED 264 C A VALUE WITH THE ASSUMPTION THAT NO CURRENT PRECIPITATION EVENT 265 C WILL RESULT IN DELTR(I) GREATER THAN TR4. 266 IF (DELTR(I).GT.0.0.AND.Z(K).GE.HMAX) GOTO 301 267 GOTO 302 268 301 CALL DP4(DELTR(I),D1,HMAX,TR4,Z(K).KFIELD,KCORR,AFIELD,ACORR,CONST,DP(K) ) 269 DPK=DP(K) 270 302 IF (DP(K).GT.O.O) GOTO 304 271 DP(K)=0.0 272 DELZ(I)=0.0 273 GOTO 306 274 304 DELZ( I ) «DELTR( I ) /DP(K) 275 306 Z(I)=Z(K)-DELZ(I) 276 IF (Z(I).GE.HMAX) Z(I)=HMAX 277 310 REM(I)=SMC1(K)+SMC2(K)+SMC3(K)-SMCF+PRE(I)-C1*PE( I) 278 320 CALL DPORE(Z(I),HMAX,KFIELD,AFIELD,KCORR,ACORR,CONST,DP(I)) 279 DPI=DP(I) 280 330 DP(I)=(DPK+DPI)/2.0 281 IF (DP(I) .GT.O.O) GOTO 326 282 DP(I)=0.0 283 DELZ(I)=0.0 284 GOTO 328 285 326 DELZ(I)=DELTR(I)/DP(I) 286 328 Z(I)=Z(K)-DELZ(I) 287 C GIVEN THE CONDITION THAT THE WATER TABLE IS LOW AND A 288 C SUDDEN LARGE RAINFALL EVENT OCCURS, DRAINABLE POROSITY 289 C IS CORRECTED (INCREASED) BY USING THE MOST CURRENTLY 290 C PREDICTED WATER TABLE DEPTH OF TOOAY. / 81 291 IF (H(K).LT.200.0.AND.PRE(I) .GT.15.0) GOTO 329 292 GOTO 435 293 329 CALL DPORE(Z(I),HMAX.KFI ELD.AFI ELD.KCORR,ACORR.CONST,DP(I)) 294 IF (DP(I).GT.O.O) GOTO 432 295 DP(I)=0.0 296 DELZ(I)=0.0 297 GOTO 433 298 432 DELZ(I)=DELTR(I)/DP(I) 299 433 Z(I)=Z(K)-DELZ(I) 300 435 H(I)=HMAX-Z(I) 301 430 IF(Z(I).GE.(HMAX-0.01)) Z(I)=HMAX 302 IF(H( I ) .LE.0.01) H(I)=0.0 303 IF(Z( I ) .LT.O) Z(I)=0.0 304 IF(H(I).GT.HMAX)H(I)=HMAX 305 IF(ICOND.EQ.O) GO TO 470 306 IF(M0NTH(I).GE.2.AND.M0NTH(I).LE.12) GO TO 440 307 GO TO 450 308 440 IF(IDAY(I).NE.1) GO TO 450 309 C 310 CALL TITLE 1(STA.JYEAR,RMAX.SMCF,HMAX) 311 C 312 450 WRITE(6,460) IDAY(I),MONTH(I),<JYEAR,PE(I),AE(I),PRE(I),RO(I),DP(I) 313 1,SD(I) .TR(I) ,SMC1(I) .SMC2(I) ,SMC3(I) ,DELTR(I) .DELZ(I) .Z(I) ,H(I) 314 460 FORMAT(' ' ,19X,3( I2 ,1X) ,F5 .2 ,9F6.2 .4F8.2) 315 WRITE(7.455) IDAY(I),MONTH(I),JYEAR,H(I) 316 455 F0RMAT(3(I2,1X),F8.2) 317 470 IF(ICONT.EO.O) GO TO 1020 318 490 <Jd=MONTH( I ) 319 IF(Z(I) .GE.999.99) GO TO 1020 320 DO 500 IM=1 ,10 321 IMM=10-IM+1 322 IK=100*IM 323 IF(Z( I ) .LT. IK.AND.Z( I ) .GE.( IK-100)) GO TO 510 324 500 CONTINUE 325 GO TO 1020 326 C 327 510 KL=IMM 328 DO 520 IM-1.KL 329 520 MATT(dd, IM)«=MATT(JJ. IM) + 1 330 GO TO 1020 331 530 RO(I)=0. 332 GO TO 240 333 540 RO(I)=0. 334 GO TO 570 335 550 CALL DPORE(Z(K).HMAX.KFIELD,AFIELD,KCORR.ACORR.CONST,DP(K)) 336 DPK=DP(K) 337 RO(I)=0.0 338 GO TO 570 339 560 CALL DPORE(Z(K).HMAX.KFIELD.AFIELD.KCORR,ACORR.CONST,DP(K)) 340 DPK=DP(K) 341 RO(I)=0.0 342 IF(TR(K)+REM(I)-TR3-TR4) 610,570,570 343 570 IF(H(K).GE.(HMAX-DW)) GO TO 590 344 580 IF (TR(K).GT.O.O.AND.H(K).LE.O.O) GOTO 582 345 SD(I)=C*(2*D*H(K)+(H(K))**2) 346 GOTO 600 347 582 SD(I)=C*(2*D*0.005*HMAX+(0.005*HMAX)**2) 348 GOTO 600 / 82 349 590 SD(I)=RMAX 350 600 T R ( I ) » T R ( K ) + R E M ( I ) - S D ( I ) 351 R0(I)=0.0 352 GO TO 250 353 610 IF(TR(K)+REM(I )-TR4) 620,580,580 354 620 IF (TR(K).GT.O.O.AND.H(K).LE.O.O) GOTO 625 355 SD(I)=C*(2*D*H(K)+(H(K))**2) 356 GOTO 626 357 625 SD(I)=C*(2*D*0.005*HMAX+(0.005*HMAX)**2) 358 626 IF(SD(I).GE.(TR(K)+REM(I))) GO TO 630 359 GO TO 600 360 630 TR(I)=0.0 361 RO(I)=0.0 362 GO TO 250 363 640 REM(I)=0. 364 NA=1 365 GO TO 210 366 650 IF(TR(K)-TR4-TR3-TR2-RMAX) 651,660,660 367 651 IF (TR(K)-TR2-TR3-TR4) 652,660,660 368 652 IF (TR(K)-TR3-TR4) 653,660.660 369 653 IF (TR(K)-TR4) 654,660,660 370 654 IF (TR(K)+PRE(I)-C1*PE(I)) 656,655.655 371 655 REM(I)=PRE(I)-C1*PE(I) 372 NA=0 373 GOTO 620 374 656 TR(I)=0.0 375 REM(I)=TR(K)+PRE(I)-C1*PE(I) 376 NA=1 377 IF (H(K).GT.O.O) GOTO 657 378 SD(I)=0.0 379 GOTO 658 380 657 SD(I)=C*(2*D*H(K)+(H(K))**2) 381 658 RO(I)=0.0 382 CALL DPORE(Z(K),HMAX,KFIELD,AFIELD,KCORR.ACORR.CONST.DP(K)) 383 DPK=DP(K) 384 GOTO 250 385 660 IF(H(K).GE.(HMAX-DW)) GO TO 670 386 IF (TR(K).GT.O.O.AND.H(K).LE.O.O) GOTO 672 387 SD(I)=C*(2*D*H(K)+(H(K))**2) 388 GOTO 680 389 672 SD(I)=C*(2*0*0.005*HMAX+(0.005*HMAX)* *2) 390 GO TO 680 391 670 SD(I)=RMAX 392 680 TR(I)=TR(K)-SD(I)+PRE(I)-C1*PE(I) 393 R0(I)=0.0 394 CALL DPORE(Z(K),HMAX,KFIELD,AFIELD,KCORR.ACORR,CONST,DP(K)) 395 DPK=DP(K) 396 GO TO 840 397 690 IF(PRE(I)-C1*PE(I)) 870,840,700 398 700 IF(SMC2(K)+SMC3(K)-SZMCF2-SZMCF3) 720,710,710 399 710 SMC1(I)=SMC1(K)+PRE(I)-C1*PE(I) 400 GO TO 270 401 720 IF(SMC3(K)-SZMCF3) 790,730,730 402 730 IF(SMC1(K)+PRE(I )-C1*PE(I)-SZMCF1)740,770,780 403 740 SMC 1(1)=SMC1(K)+PRE(I)-C1*PE(I) 404 750 SMC2(I)=SMC2(K) 405 760 SMC3(I)=SMC3(K) 406 GO TO 290 / 83 407 770 SNICK I )-SZMCF1 408 GO TO 850 409 780 SMC1(I)-SZMCF1 410 SMC2(I)-SMC2(K)+SMC1(K)+PRE(I)-C1*PE(I)-SZMCF1 411 GO TO 760 412 790 IF(SMC1(K)+SMC2(K)+PRE(I)-C1*PE(I)-SZMCF1-SZMCF2) 730,820.800 413 800 SMC3(I)-SMC3(K)+SMC1(K ) + SMC2(K)+PRE(I)-C1*PE(I)-SZMCF1-SZMCF2 414 810 SMC 1(I)"SZMCF1 415 SMC2(I)-SZMCF2 416 GO TO 290 417 820 SMC3(I)-SMC3(K) 418 GO TO 810 419 840 IF (TR(I) .GT.O.O) GOTO 842 420 SMC1( I )«SMC1(K)+TR(I ) 421 T R ( I ) » 0 . 0 422 GOTO 850 423 842 SMC1(I)«SMC1(K) 424 850 SMC2(I)»SMC2(K) 425 860 SMC3(I)-SMC3(K) 426 GO TO 290 427 870 IF(SMC1(K)+SMC2(K)+SMC3(K)-SMCF) 890,880,880 428 880 SMC1(I)=SMC1(K)+PRE(I)-C1*PE(I)+TR(K) 429 GO TO 850 430 890 IF(SMC2(K)+SMC3(K)-SZMCF2-SZMCF3) 920,900.900 431 900 SMC1(I)=SMC1(K)+PRE(I)-C1*PE(I) 432 IF(SMC1(I)) 910.850,850 433 910 SMC 1(1)-0.0 434 SMC2(I)*SMC2(K)+C2*(SMC1(K)+PRE(I)-C1*PE(I)) 435 AE(I)=SMC1(K)-C2*(SMC1(K)+PRE(I)-C1*PE(I)) 436 SMC3(I)=SMC3(K) 437 GO TO 300 438 920 SMC2(I)=SMC2(K)+PRE(I)-C1*PE(I) 439 SMC1(I)=0.0 440 IF(SMC2(I)) 830,860,860 441 930 S M C 1 ( I ) « 0 . 0 442 SMC2( I )«SMC2(K)+C2*(SMC1(K)+PRE( I ) -C1*PE( I ) ) 443 IF(SMC2(I)) 940,970,970 444 940 S M C 2 ( I ) « 0 . 0 445 IF(SMC2(K).E0.O.O)G0 TO 950 446 SMC3( I )«SMC3(K)+C3*(SMC2(K)+C2*(SMC1(K)+PRE( I ) -C1*PE( I ) ) ) 447 IF(SMC3(I).LT.0.0)G0 TO 960 448 AE(I)=SMC1(K)+SMC2(K)-C3*(SMC2(K)+C2*(SMC1(K)+PRE(I)-C1*PE(I))) 449 GO TO 300 450 950 SMC3(I)-SMC3(K)+C3*(SMC1(K)+PRE(I)-C1*PE(I)) 451 IF(SMC3(I).LT.O.O) GO TO 960 452 AE(I)=SMC1(K)-C3*(SMC1(K)+PRE(I)-C1*PE(I)) 453 GO TO 300 454 960 SMC3(I)=0.0 455 AE(I)=SMC1(K)+SMC2(K)+SMC3(K) 456 GO TO 300 457 970 AE(I)=SMC1(K)-C2*(SMC1(K)+PRE(I)-C1*PE(I)) 458 SMC3(I)=SMC3(K) 459 GO TO 300 460 C OUTPUT DATA CHECK FOR TABLE 1 461 980 IC0M=IC0M+1 462 JSTAY=K-(IC0M-1) 463 ICHECK=1 464 IF(ICOND.EQ.O) GO TO 1020 / 84 465 IF(IDAY(I) .E0 .1 .AND.M0NTH(I) .NE .1) CALL TITLE 1(STA.JYEAR, 466 1 RMAX,SMCF,HMAX) 467 IF(PE(I).EQ.TEMPI.AND.PRE(I).EQ.TEMP2) WRITE(6,990) IDAY(I),M0NTH( 468 11),JYEAR 469 990 F0RMAT(' ' ,19X,3( I2,1X), 'MISSING',5X. 'MISSING') 470 IF(PE(I).E0.TEMP1.AND.PRE(I).NE.TEMP2) WRITE(6,1000)IDAY(I).M0NTH( 471 11) ,JYEAR,PRE(I) 472 1000 F0RMAT(' ' ,19X,3( I2 .1X) , 'MISSING' .4X,F6.2) 473 IF(PRE(I).E0.TEMP2.AND.PE(I).NE.TEMP 1) WRITE(6,1010)IDAY(I),M0NTH( 474 11) ,JYEAR,PE( I ) 475 1010 F0RMAT(' ' ,19X,3( I2 ,1X) ,F5.2,7X, 'MISSING') 476 1020 CONTINUE 477 C END LOOP FOR YEAR 478 C OUTPUT TABLE 2 479 IF(ICONT.EO.O) GO TO 1090 480 WRITE(6,1030)(STA(I),I=1,10).JYEAR 481 1030 FORMAT( '1 ' ,5( / ) ,27X, 'TABLE : MONTHLY SUMMARY OF THE NUMBER OF 482 1DAYS THAT THE PREDICTED WATER', / ,42X, 483 2'TABLE DEPTH IS LESS THAN THE INDICATED WATER TABLE DEPTH' , / / , 45X , 484 3'STATION : ' ,10A2.9X, 'YEAR : 19',12) 485 WRITE(6,1040) (DEP(I),I=1.10) 486 1040 F 0 R M A T ( / , 2 1 X , 8 7 ( ' * ' ) , / , 2 1 X . ' * WATER TABLE * ' . 7 0 X . ' * ' , / , 487 121X,'* DEPTH * ' . 1X.F5.0.9(F7.0) ,1X. ' * ' . / , 2 1 X , 488 2 ' * (MM) * ' , 7 0 X , ' * ' , / . 2 1 X , 8 7 ( ' * ' ) ) 489 DO 1050 J1=1,12 490 1050 WRITE(6,1060) (KMONT(J1,J2),J2=1,5), (MATT(J1.J2),J2=1,10) 491 1060 F0RMAT(21X,'* ' , 5 A 2 , ' * ' . 14,3X,8(14,3X), 14,2X, 492 1 ' * ' , / , 2 1 X , ' * ' , 1 4 X . ' * ' , 7 0 X , ' * ' ) 493 WRITE(6,1070) RMAX,SMCF,HMAX 494 1070 FORMAT(21X,87( ' * ' ) , / / ,22X, 495 1'RMAX = ' . F 4 . 0 , ' MM/DAY',/ ,22X,'AVAILABLE SOIL MOISTURE', 496 2' FOR TOP THREE ZONES = ' . F 5 . 0 , ' MM',/,22X,'SUBSURFACE DRAIN', 497 3' DEPTH = ' . F 6 . 0 , ' MM') 498 1080 DO 1085 IJ=1,12 499 DO 1085 JI = 1 , 14 500 ITOTAL(IJ,JI)=ITOTAL(IJ,JI)+MATT(IJ,JI) 501 1085 MATT(IJ,JI ) =0 502 1090 IF(IFDOWT.EO.O) GO TO 1450 503 C 504 C . 505 C 506 C CHECKING FOR LEAP YEAR 507 C 508 C 509 MZ=1900+JYEAR 510 MM=MZ/4 511 MM=4*MM 512 IF(MM.EQ.MZ)IIK=367 513 IIK=366 514 DO 1440 11=1.10 515 Z1=DEP(II) 516 IR=II 517 IF(Z1-1OO0.0) 1100,1100,3000 518 1100 IC=1 519 JC=1 520 INMJ=1 521 INJA=1 522 INAN=1 523 INGS-1 524 INNS-1 525 1110 DO 1430 1*2,IIK 526 IF(PE(I).EO.TEMP 1.OR.PRE(I).EO.TEMP2) GO TO 1430 527 IF (Z( I ) .LE.Z1) GO TO 1120 528 IC=1 529 <JC=1 530 INMd=1 531 INJA=1 532 INAN-1 533 INGS-1 534 INNS=1 535 GO TO 1430 536 1120 IF(IC0NF1.EQ.O) GO TO 1140 537 DO 1130 IN-1.IC 538 INF=IN 539 1130 I T 0 T F ( I N , I R ) « I T 0 T F ( I N , I R ) + 1 540 IF(MAX1.LE.INF) MAX1-INF 541 IC=IC+1 542 IF(IC.EQ.367) GO TO 4000 543 1140 IF(IC0NF2.EQ.O) GO TO 1160 544 C ACCOUNTING FOR LEAP YEAR 545 C 546 MZ=1900+JYEAR 547 MM=MZ/4 548 MM=4*MM 549 IF(MM.E0.MZ)G0 TO 1 550 MN=60 551 KKL=MN 552 GO TO 2 553 1 MN=61 554 KKL=MN 555 2 CONTINUE 556 IF( I .LT.2.0R. I .GT.MN) GO TO 1170 557 DO 1150 IMJ=1,INMJ 558 IMF2=IMd 559 1150 MMJF(IMJ,IR)=MMJF(IMd,IR)+1 560 IF(MAX2.LE.IMF2) MAX2=IMF2 561 INMd«INMJ+1 562 IF(INMd.EQ.63) GO TO 5000 563 1160 IF(IC0NF3.EQ.0) GO TO 1190 564 1170 CONTINUE 565 N1=KKL+61 566 N2=KKL+1 567 IF(I .LT.N2.OR.I .GT.N1)G0 TO 1200 568 DO 1180 IMA=1,INJA 569 IMF3=IMA 570 1180 MJAF(IMA,IR)=MJAF(IMA,IR)+1 571 IF(MAX3.LE.IMF3) MAX3=IMF3 572 INUA=INJA+1 573 IFONJA.EQ.63) GO TO 6000 574 1190 IF(IC0NF4.E0.O) GO TO 1220 575 1200 CONTINUE 576 N3=N1+123 577 N4=N1+1 578 IF(I .LT.N4.OR.I .GT.N3)G0 TO 1201 579 DO 1210 IAN=1,INAN 580 IMF4=IAN 581 1210 MANF(IAN, IR)»MANF(IAN, IR)+1 582 IF(MAX4.LE.IMF4) MAX4-IMF4 583 INAN=INAN+1 584 IF(INAN.EO.125) GO TO 7000 585 1220 IF(IC0NF5.E0.O) GO TO 1204 586 1201 CONTINUE 587 N5«N3+61 588 N6-N3+1 589 IF( I .LT.N6.OR.I .GT.N5)G0 TO 1202 590 DO 1203 ING-LINGS 591 IMF5-ING 592 1203 MS0F(ING,IR)-MS0F(ING.IR)+1 593 IF(MAX5.LE.IMF5)MAX5»IMF5 594 INGS-INGS+1 595 IF(INGS.EQ.63)G0 TO 7001 596 1204 IF(IC0NF6.E0.0)G0 TO 1206 597 1202 N7-N5+61 598 N8-N5+1 599 IF( I .LT.N8.0R.I .GT.N7)G0 TO 1230 600 DO 1205 INN»1,INNS 601 IMF6-INN 602 1205 MNDF(INN,IR)-MNDF(INN,IR)+1 603 IF(MAX6.LE.IMF6)MAX6*IMF6 604 INNS-INNS-H 605 IF(INNS.E0.63)G0 TO 7003 606 1206 IF(IC0NF4.EQ.O)G0 TO 1430 607 1230 IF(IDAY(I).GT.1) GO TO 1240 608 dC=1 609 1240 DO 1250NM-1.12 610 IF(MONTH(I).GT.NM) GO TO 1250 611 NN-NM-3 612 GO TO (3,4,5,1260,1280,1300,1320.1340.1360,1380,1400,6),NM 613 1250 CONTINUE 614 3 DO 44 M1=1,dC 615 44 M1F(M1,IR)*M1F(M1,IR)+1 616 GO TO 1420 617 4 DO 56 M2-1,dC 618 56 M2F(M2. IR)«M2F(M2, IR)+1 619 GO TO 1420 620 5 DO 57 M3-1.JC 621 57 M3F(M3, IR)«M3F(M3, IR)+1 622 GO TO 1420 623 1260 DO 1270M4-1.dC 624 1270 M4F(M4. IR)«M4F(M4, IR)+1 625 GO TO 1420 626 1280 DO 1290M5-1,dC 627 1290 M5F(M5,IR)-M5F(M5,IR)+1 628 GO TO 1420 629 1300 DO 1310 M6-1,dC 630 1310 M6F(M6,IR)-M6F(M6,IR)+1 631 GO TO 1420 632 1320 DO 1330M7-1.dC 633 1330 M7F(M7,IR)=M7F(M7,IR)+1 634 GO TO 1420 635 1340 DO 1350 M8=1,dC 636 1350 M8F(M8.IR)=M8F(M8,IR)+1 637 GO TO 1420 638 1360 DO 1370M9=1,dC / 87 639 1370 M9F(M9,IR)=M9F(M9,IR)+1 640 GO TO 1420 641 1380 DO 1390 M10=1,JC 642 1390 M10F(M10,IR)=M10F(M10,IR)+1 643 GO TO 1420 644 1400 DO 1410 M11=1,JC 645 1410 M11F(M11,IR)=M11F(M11,IR)+1 646 GOTO 1420 647 6 DO 58 M12=1,JC 648 58 M12F(M12,IR)=M12F(M12.IR)+1 649 1420 JC=JC+1 650 1430 CONTINUE 651 1440 CONTINUE 652 C 653 1450 READ(IRD,1460.END=1470)STAT,JYEAR,IMONT,JDAY,PPT.PET 654 1460 FORMAT(17,312, 1GX,F4.1,1IX.F4.1) 655 C 656 1=2 657 IJK=0 658 GO TO 110 659 C .• 660 1470 IF(ICONT.EQ.O) GO TO 1500 661 WRITE(6.1480) (STA(I) ,1 = 1,10),ISTART,IEND 662 1480 FORMAT( '1 ' ,5 ( / ) ,27X, 'TABLE : MONTHLY SUMMARY OF THE NUMBER OF 663 1DAYS THAT THE PREDICTED WATER',/ .42X, 664 2'TABLE DEPTH IS LESS THAN THE INDICATED WATER TABLE DEPTH ' , / / , 45X , 665 3'STATION : ' ,10A2,4X, 'YEAR : 19 ' ,12 , ' -19 ' ,12 ) 666 WRITE(6,1040) (DEP(IJ),IJ=1,10) 667 00 1490 ITIT=1 , 12 668 1490 WRITE(6,1060)(KMONT(ITIT,JJJ),JJJ=1,5),(ITOTAL(ITIT,JK),JK=1,10) 669 WRITE(6,1070) RMAX,SMCF,HMAX 670 1500 IF(IFDOWT.EO.O) GO TO 2000 671 IF(IC0NF4.E0.O) GO TO 1760 672 WRITE(6,1600) (KM0NT(1,K1),K1 = 1,5), ISTART,I END,(STA(I),1 = 1,10) 673 WRITE(6,1650) (DEP(I),I=1,10) 674 DO 1691 1=1,31 675 1691 WRITE(6,1670) I , (M1F(I,J),J=1,10) 676 WRITE(6,1680) RMAX,SMCF,HMAX 677 WRITE(6,1600) (KM0NT(2,K1),K1 = 1,5) , I START.I END,(STA(I) ,1=1.10) 678 WRITE(6,1650) (DEP(I),1=1,10) 679 DO 1692 1=1,31 680 1692 WRITE(6,1670) I .(M2F(I,d).J=1,10) 681 WRITE(6,1680) RMAX,SMCF,HMAX 682 WRITE(6.1600) (KM0NT(3.K1),K1 = 1,5), I START.I END,(STA(I) .1=1.10) 683 WRITE(6,1650) (DEP(I),I=1.10) 684 DO 1693 1=1,31 685 1693 WRITE(6,1670) I ,(M3F(I,J),J=1,10) 686 WRITE(6,1680) RMAX,SMCF,HMAX 687 WRITE(6,1600) (KM0NT(4,K1),K1 = 1.5), I START,IEND,(STA(I ) ,1=1.10) 688 1600 FORMAT( '1 ' ,5 ( / ) ,32X, 'TABLE : FREQUENCY DISTRIBUTION OF PREDICT 689 1 ED WATER TABLE DEPTHS', / ,46X, 'FOR THE MONTH O F ' , 5X, 5A2, 3X, 690 2'YEAR : 1 9 ' , 1 2 , ' - 19 ' ,12 , / ,50X, 'STATION : ' , 10A2) 691 WRITE(6,1650) (DEP(I),I=1.10) 692 1650 F0RMAT(31X,72( ' * ' ) , / , 31X, ' * W.T.D. * ' . 6 1 X , ' * ' , / , 3 1 X , ' * ' ,8X, ' * ' , 693 1 10F6.0.1X , ' * ' , / , 3 1 X , ' * DAY * ' , 6 1 X , ' * ' , / , 3 1 X , 7 2 ( ' * ' ) ) 694 DO 1660 1=1,30 695 1660 WRITE(6,1670) I ,(M4F(I,J).J=1,10) 696 1670 F0RMAT(31X, ' * ' , 15, ' * ' , 1 X , 10(14,2X) , ' * ' ) / 88 697 WRITE(6,1680) RMAX,SMCF,HMAX 698 1680 F0RMAT(31X,72( ' * ' ) , / / , 34X, 'W.T .D . IS WATER TABLE DEPTH IN M M ' , / , 699 134X,'DAY REPRESENTS THE NUMBER OF SUCCESSIVE DAYS IN A MOIST PE 700 2RI0D.' , / ,34X,'NUMBERS IN THE TABLE ARE THE NUMBER OF OCCURRENCES F 701 30R WHICH THE', / ,41X, 'WATER TABLE DEPTH IS EQUAL TO OR LESS THAN TH 702 4E PRESCRIBED',/ ,41X,'DEPTH FOR THE NUMBER OF SUCCESSIVE DAYS INDIC 703 5ATED.' , / ,34X, 'RMAX= ' , F 4 . 0 , ' MM/DAY',/ ,34X,'AVAILABLE SOIL ' , 704 6' MOISTURE FOR THE TOP THREE ZONES = ' , F 5 . 0 . / , 3 4 X , ' S U B S U R F A C E ' , 705 7' DRAIN DEPTH = ' . F 5 . 0 , ' MM') 706 WRITE(6.1600) (KMONT(5,K1),K1 = 1,5). ISTART,I END,(STA(I), I = 1,10) 707 WRITE(6,1650) (DEP(I),I = 1,10) 708 DO 1690 1=1,31 709 1690 WRITE(6,1670) I,(M5F(I,d).d=1.10) 710 WRITE(6,1680) RMAX,SMCF,HMAX 711 WRITE(6.1600) (KMONT(6,K1),K1=1,5), ISTART,I END,(STA(I), I = 1,10) 712 WRITE(6,1650) (DEP(I),I=1,10) 713 DO 1700 1=1,30 714 1700 WRITE(6.1670) I,(M6F(I.d),d=1.10) 715 WRITE(6,1680) RMAX,SMCF,HMAX 716 WRITE(6,1600) (KMONT(7,K1),K1 = 1,5), I START,I END,(STA(I), I = 1,10) 717 WRITE(6,1650) (DEP(I),I=1,10) 718 DO 1710 1=1,31 719 1710 WRITE(6,1670) I, (M7F(I,d) ,d =1,10) 720 WRITE(6,1680) RMAX,SMCF,HMAX 721 WRITE(6,1600) (KMONT(8,K1),K1 = 1,5), I START.I END,(STA(I) ,1 = 1,10) 722 WRITE(6,1650) (DEP(I),I=1,10) 723 DO 1720 1=1,31 724 1720 WRITE(6,1670) I,(M8F(I,d).d=1,10) 725 WRITE(6,1680) RMAX,SMCF,HMAX 726 WRITE(6,1600) (KMONT(9,K1),K1=1,5), ISTART,IEND,(STA(I),I=1,10) 727 WRITE(6,1650) (DEP(I),I=1,10) 728 DO 1730 1=1,30 729 1730 WRITE(6,1670) I,(M9F(I,d).d=1,10) 730 WRITE(6,1680) RMAX.SMCF,HMAX 731 WRITE(6,1600) (KMONT(10,K1),K1 = 1,5), I START,IEND,(STA(I),I = 1,10) 732 WRITE(6,1650) (DEP(I),I=1,10) 733 DO 1740 1=1,31 734 1740 WRITE(6,1670) I,(M10F(I,d),d=1,10) 735 WRITE(6,1680) RMAX,SMCF,HMAX 736 WRITE(6,1600) (KMONT(11,K1),K1 = 1,5), I START,I END,(STA(I), I = 1,10) 737 WRITE(6,1650) (DEP(I),I=1,10) 738 DO 1750 1=1,30 739 1750 WRITE(6,1670) I,(M11F(I,d),d=1,10) 740 WRITE(6.1680) RMAX,SMCF,HMAX 741 WRITE(6,1600) (KMONT(12,K1),K1 = 1,5), ISTART,I END,(STA(I), I = 1.10) 742 WRITE(6,1650) (DEP(I),I = 1,10) 743 DO 1694 1=1,31 744 1694 WRITE(6.1670) I,(M12F(I,d),d=1,10) 745 WRITE(6,1680) RMAX,SMCF,HMAX 746 1760 IF(ICONF2.EQ.O) GO TO 1800 747 WRITE(6,1780)(PERIOD(1,I),I=1,18),ISTART,IEND,(STA(I),I=1,10) 748 1780 FORMAT( '1 ' ,5( / ) ,32X, 'TABLE : FREQUENCY DISTRIBUTION OF PREDICT 749 1 ED WATER TABLE DEPTHS', / ,45X, 'FOR THE PERIOD: ' ,18A1, ' YEAR : 19 750 2 ' , 1 2 , ' - 1 9 ' , I 2 , / , 5 1 X , 'STATION : '.10A2) 751 WRITE(6,1650) (DEP(I),I=1,10) 752 IMA2=MAX2+1 753 DO 1790 d=1,IMA2 754 1790 WRITE(6,1670)d,(MMdF(d,d1).d1=1,10) / 89 755 WRITE(6.1680) RMAX,SMCF,HMAX 756 1800 IFUC0NF3.EQ.0) GO TO 1820 757 WRITE(6,178O)(PERI0D(2,I),I=1.18),ISTART.I END,(STA(I) , I = 1,10) 758 WRITE(6,1650) (DEP(I ) ,1 = 1,10) 759 IMA3=MAX3+1 760 DO 1810 L=1,IMA3 761 1810 WRITE(6,1670)L,(MJAF(L,L1),L1=1, 10) 762 WRITE(6.1680) RMAX.SMCF.HMAX 763 1820 IF(IC0NF5.E0.O) GO TO 1840 764 WRITE(6,178O)(PERI0D(3,I),I=1,18),ISTART.I END,(STA(I) , I * 1,10) 765 WRITE(6,1650) (DEP(I),1=1,10) 766 IMA4=MAX4+1 767 DO 1831 K=1,IMA4 768 1831 WRITE(6,1670) K,(MANF(K,K1),K1=1,10) 769 WRITE(6,1680) RMAX,SMCF,HMAX 770 WRITE(6,178O)(PERI0D(4,1),1 = 1,18),ISTART,I END.(STA(I) , I = 1,10) 771 WRITE(6,1650) (DEP(I).I = 1,10) 772 IMA5=MAX5+1 773 DO 1832 K=1,IMA5 774 1832 WRITE(6,1670) K,(MSOF(K,K1),K1 = 1 , 10) 775 WRITE(6,1680) RMAX,SMCF,HMAX 776 WRITE(6,178O)(PERI0D(5,I),I = 1,18).ISTART , I END,(STA(I) , I = 1,10) 777 WRITE(6,1650) (DEP(I).1 = 1,10) 778 IMA6=MAX6+1 779 DO 1830 K-1.IMA6 780 1830 WRITE(6,1670) K,(MNDF(K,K1),K1• 1, 10) 781 WRITE(6,1680) RMAX,SMCF,HMAX 782 1840 IF(IC0NF1.EO.O) GO TO 9000 783 WRITE(6,1850) (STA(I ) ,I = 1,10),ISTART. IEND 784 1850 FORMAT( '1 ' ,5( / ) ,32X, 'TABLE : FREQUENCY DISTRIBUTION OF PREDICT 785 1ED WATER TABLE DEPTHS', / ,46X, 786 2'STATION : ' ,10A2,4X, 'YEAR : 19 ' , 12 , ' - 19 ' . 12 ) 787 WRITE(6,1650) (DEP(I),I=1.10) 788 IMA1=MAX1+1 789 DO 1870 N=1,IMA1 790 1870 WRITE(6,1670)N,(IT0TF(N,M),M=1,10) 791 WRITE(6,1680) RMAX,SMCF,HMAX 792 1880 WRITE(6,1890) 793 1890 FORMAT('1' ,10X, ' PROGRAM HAS FUNCTIONED VERY WELL BYE BYE ') 794 GO TO 9000 795 2O00 WRITE(6,2100) 796 2100 FORMAT('0',10X,'FREQUENCY DISTRIBUTION IS NOT WANTED') 797 GO TO 9000 798 3000 WRITE(6,3100) 799 3100 FORMAT(' ' . 2X,'PROGRAM WAS TERMINATED SINCE Z IS GREATER THAN 100 800 10. - - - CHECK STATEMENTS BETWEEN 1012 AND 502 . ' ) 801 GO TO 9000 802 4OO0 WRITE(6,4100) 803 4100 FORMAT('0',5X,'YOU HAVE EXCEEDED THE SIZE OF THE ITOTF MATRIX') 804 GO TO 9000 805 5000 WRITE(6,5100) 806 5100 FORMAT('0',5X,'YOU HAVE EXCEEDED THE SIZE OF THE MMJF MATRIX') 807 GO TO 9000 808 6000 WRITE(6,6100) 809 6100 FORMAT('0'.5X,'YOU HAVE EXCEEDED THE SIZE OF THE MJAF MATRIX') 810 GO TO 9000 811 7000 WRITE(6,7100) 812 7100 FORMAT('0',10X,'YOU HAVE EXCEEDED THE SIZE OF MANF MATRIX') / 90 813 GO TO 9O0O 814 7001 WRITE(6,7002) 815 7002 FORMAT('0',10X,'YOU HAVE EXCEEDED THE SIZE OF MSOF MATRIX') 816 GO TO 9000 817 7003 WRITE(6,7004) 818 7004 FORMAT!'0'.10X,'YOU HAVE EXCEEDED THE SIZE OF MNDF MATRIX') 819 9000 STOP 820 END 821 ^***»********»*******»*************+*********+*************************** 822 SUBROUTINE TRANST(D1,KFI ELD.KCORR,AFIELD,ACORR,CONST, 823 •HMAX,TR1 SUM,TR2SUM,TR3SUM,TR4SUM,TRMAX) 824 REAL KFIELD,KCORR 825 C MAXIMUM TRANSIENT STORAGE CALCULATED BY DECREMENTS OF DELTZ FROM SOIL 826 C SURFACE TO DRAIN DEPTH (HMAX). 827 DELTZ=10.0 828 C CALCULATE MAXIMUM DRAINABLE POROSITY AT SOIL SURFACE (DPMAX) 829 Z=0.0 830 CALL DPORE(Z,HMAX,KFI ELD,AFI ELD,KCORR,ACORR,CONST,DPMAX) 831 c MAXIMUM TRANSIENT STORAGE FOR ZONE 1 832 D11=0.0 833 TR1SUM=0.0 834 CALL DPORE(011,HMAX,KFIELD,AFIELD,KCORR.ACORR,CONST.DP) 835 DP1=DP 836 D11=D11+DELTZ 837 10 IF (D11.GT.D1) GOTO 20 838 CALL DP0RE(D11,HMAX,KFIELD,AFIELD,KCORR,ACORR,CONST,DP) 839 DP2=DP 840 WATER 1=DELTZM(DP2+DP1 )/2.0) 841 TR1SUM=TR1SUM+WATER 1 842 DP1=DP2 843 D11=D11+DELTZ 844 GOTO 10 845 c MAXIMUM TRANSIENT STORAGE FOR ZONE 2 846 20 D2=D11-DELTZ 847 TR2SUM=0.0 848 30 IF (D2.GT.(D1+300.)) GOTO 40 849 CALL DPORE(D2,HMAX,KFI ELD,AFI ELD,KCORR,ACORR,CONST,DP) 850 DP2=DP 851 WATER2=DELTZ^((DP2+DP1)/2.0) 852 TR2SUM=TR2SUM+WATER2 853 DP1=DP2 854 D2=D2+DELTZ 855 GOTO 30 856 C MAXIMUM TRANSIENT STORAGE FOR ZONE 3 857 40 D3=D2-DELTZ 858 TR3SUM=0.0 859 50 IF (D3.GT.(D1+300.+400.)) GOTO 60 860 CALL DPORE(D3,HMAX,KFI ELD,AFI ELD,KCORR,ACORR,CONST,DP) 861 DP2=DP 862 WATER3=DELTZ*((DP2+DP1)/2.0) 863 TR3SUM=TR3SUM+WATER3 864 DP1=DP2 865 D3=D3+DELTZ 866 GOTO 50 867 C MAXIMUM TRANSIENT STORAGE FOR ZONE 4 868 60 D4=D3-DELTZ 869 TR4SUM=0.0 870 70 IF (D4.GT.HMAX) GOTO 80 / 91 871 CALL DPORE(D4.HMAX,KFI ELD,AFI ELD,KCORR,ACORR,CONST,DP) 872 DP2-DP 873 WATER4=DELTZ*((DP2+DP1)/2.0) 874 TR4SUM=TR4SUM+WATER4 875 DP1=DP2 876 D4-D4+DELTZ 877 GOTO 70 878 C TOTAL MAXIMUM TRANSIENT STORAGE 879 80 TRMAX-TR1SUM+TR2SUM+TR3SUM+TR4SUM 880 C MAX. REMAINING TRANSIENT WATER STORAGE FOR ZONE 1 881 TR1*TRMAX-TR1 SUM 882 C MAX. REMAINING TRANSIENT WATER STORAGE FOR ZONE 2 883 TR2-TR1-TR2SUM 884 C MAX. REMAINING TRANSIENT WATER STORAGE FOR ZONE 3 885 TR3-TR2-TR3SUM 886 C MAX. REMAINING TRANSIENT WATER STORAGE FOR ZONE 4 887 TR4=TR3-TR4SUM 888 IF (TR4.LE.0.0) TR4=0.0 889 RETURN 890 END 891 SUBROUTINE DP4(DELTR,D1,HMAX,TR4,ZK,KFIELD,KCORR,AFI ELD,ACORR, 892 * CONST,DPK) 893 REAL KFIELD,KCORR 894 C ASSIGNED DRAINABLE POROSITY GIVEN THE CONDITION THAT THE WATER 895 C TABLE IS AT THE DRAIN LEVEL AND POSITIVE TRANSIENT STORAGE 896 C AT THE SAME TIME. 897 D4 =HMAX-(D1+300.0+400.0) 898 IF (DELTR.LT.TR4*2.0/3.0) GOTO 10 899 IF (DELTR.LT.TR4/3.0) GOTO 20 900 Z=ZK-D4*5.0/6.0 901 GOTO 30 902 10 Z=ZK-D4/6.0 903 GOTO 30 904 20 Z=ZK-D4/2.0 905 30 CALL DPORE(Z,HMAX.KFI ELD,AFI ELD,KCORR,ACORR,C0N6T.DPK) 906 RETURN 907 END 908 SUBROUTINE DPORE(Z,HMAX,KFIELD,AFI ELD,KCORR,ACORR,CONST,DP) 909 REAL KFIELD,KCORR 910 C DRAINABLE POROSITY CALCULATED BY AN EXPONENTIAL RELATIONSHIP 911 C BETWEEN POROSITY AND TENSION. 912 C TRANSFORM WATER TABLE DEPTH TO WATER TABLE HEIGHT 913 HT=HMAX-Z 914 C CONVERT WATER TABLE HEIGHT HT FROM MM TO CM 915 ZCM=HT/10.0 916 C DETERMINE POROSITY BASED ON RAINFALL AND W.T. DEPTH DATA 917 PRAIN=AFIELD*(1-EXP(-KFIELD*ZCM)) 918 C CORRECT ABOVE POROSITY ON LABORATORY-POROSITY BASIS 919 PCORR=PRAIN/(ACORR*(1-EXP(-KCORR*ZCM))+CONST) 920 C CONVERT ADJUSTED POROSITY FROM PERCENTAGE TO FRACTION 921 DP=PCORR/100.0 922 RETURN 923 END 924 SUBROUTINE FACTOR(RMAX,DW,D,HMAX,C,SMCF,SZMCF1,SZMCF2,SZMCF3, 925 1SZMCF4,TRMAX,TR1.TR2,TR3,TR4,C1,C2,C3,AFI ELD.KFI ELD,ALAB,KLAB,ACORR 926 *KCORR,CONST) 927 REAL KFIELD,KLAB,KCORR 928 WRITE(6,10) / 92 929 10 FORMAT( '1 ' , 15 ( / ) , 30X ,68 ( ' * ' ) , / , 30X , 930 1'* INPUT DRAINAGE PARAMETERS FOR PROGRAMMING : ' , 2 1 X , ' * ' ) 931 WRITE(6,20) RMAX.DW,D,HMAX,C 932 20 F 0 R M A T ( 3 0 X , ' * ' , 6 6 X , ' * ' . / , 3 0 X , 933 1'* MAXIMUM DRAINAGE RATE (RMAX)' ,23( ' - ' ) ,F5.1. 'MM/DAY * ' f / , 3 0 X , 934 2 ' * DEPTH OF ZONE FROM SOIL SURFACE THAT RMAX OCCURS (DW)- - - ' , 935 3F4.0,'MM * ' , / , 3 0 X , 936 4 ' * DEPTH OF IMPERMEABLE LAYER FROM THE DRAIN ( D ) ' , 1 0 ( ' - ' ) , F 5 . 0 , 937 5'MM * ' , / , 3 0 X , 938 6 ' * MAXIMUM DRAIN DEPTH (HMAX) ' ,29 ( ' - ' ) , F5.0,'MM * ' , / , 3 0 X , 939 7 ' * DRAINAGE COEFFICIENT ( C»4K/L**2 ) ' , 1 8 ( ' - ' ) , F 1 0 . 8 . ' * ' ) 940 WRITE(6,25) KFIELD,KLAB,AFI ELD,ALAB 941 25 F0RMAT(30X,'* RATE CONSTANT OF FIELD MODEL ( K F I E L D ) ' , 1 6 ( ' - ' ) , 942 1F7.5 , ' /CM * ' . / , 3 0 X , ' * RATE CONSTANT OF LAB MODEL (KLAB) ' , 943 2 1 9 ( ' - ' ) , F 7 . 5 . ' / C M * ' , / , 3 0 X , ' * MAXIMUM POROSITY OF ' , 944 3'FIELD MODEL (AFIELD)' , 1 3 ( ' - ' ) , F 8 . 4 , ' % * ' . / , 3 0 X . ' * MAXIMUM', 945 4' POROSITY OF LAB MODEL ( A L A B ) ' , 1 6 ( ' - ' ) , F 8 . 4 , ' % * ' ) 946 WRITE(6.35) KCORR,ACORR,CONST 947 35 F0RMAT(30X,'* RATE CONSTANT OF CORRECTION EQUATION (KCORR)', 948 1 8 ( ' - ' ) , F 7 . 5 , ' / C M * ' , / , 3 0 X , ' * MAXIMUM RATIO OF CORRECTION ' , 949 2'EQUATION (ACORR ) ' , 8 ( ' - ' ) , F 6 . 4 , ' %/% * ' , / , 3 0 X , ' • CORRECTION', 950 3' CONSTANT ( C O N S T ) ' . 3 0 ( ' - ' ) , F 5 . 3 , ' * ') 951 WRITE(6,30) SMCF,SZMCF1,SZMCF2 952 30 FORMAT(30X,'* AVAILABLE SOIL MOISTURE FOR TOP THREE ZONES ' , 953 1 ' ( S M C F ) ' , 4 ( ' - ' ) , F 5 . 1 . ' MM * ' , / , 3 0 X , 954 2 ' * AVAILABLE SOIL MOISTURE FOR 1ST ZONE ( S Z M C F 1 ) ' , 9 ( ' - ' ) , 955 3F5. 1 , ' MM * ' , / , 3 0 X , 956 4 ' * AVAILABLE SOIL MOISTURE FOR 2ND ZONE (SZMCF2) ' , 9 ( ' - ' ), 957 5F5 . 1 , ' MM * ' ) 958 WRITE(6,40) SZMCF3,SZMCF4 959 40 FORMAT(30X, 960 1'* AVAILABLE SOIL MOISTURE FOR 3RD ZONE ( S Z M C F 3 ) ' , 9 ( ' - ' ) . 961 2 F 5 . 1 , ' MM * ' , / , 3 0 X , 962 3 ' * AVAILABLE SOIL MOISTURE FOR 4TH ZONE ( S Z M C F 4 ) ' , 9 ( ' - ' ) , 963 4 F 5 . 1 , ' MM * ' ) 964 WRITE(6.50) TRMAX,TR1,TR2,TR3,TR4 965 50 F0RMAT(30X,'* MAXIMUM TRANSIENT STORAGE FOR WHOLE PROFILE (TRMAX) 966 1 ' , F 5 . 1 , ' MM * ' , / , 3 0 X , 967 2 ' * MAXIMUM TRANSIENT STORAGE FOR 1ST ZONE ( T R 1 ) ' , 1 0 ( ' - ' ) ,F5 .1 , 968 3' MM * ' , / . 3 0 X , 969 4 ' * MAXIMUM TRANSIENT STORAGE FOR 2ND ZONE ( T R 2 ) ' , 1 0 ( ' - ' ) , F 5 . 1 , 970 5' MM * ' . / , 3 0 X , 971 6 ' * MAXIMUM TRANSIENT STORAGE FOR 3RD ZONE ( T R 3 ) ' , 1 0 ( ' - ' ) . F 5 . 1 , 972 7' MM * ' , / , 3 0 X , 973 8 ' * MAXIMUM TRANSIENT STORAGE FOR 4TH ZONE ( T R 4 ) ' , 1 0 ( ' - ' ) , F 5 . 1 , 974 9' MM *') 975 WRITE(6,60) C1.C2.C3 976 60 F0RMAT(30X,'* COEFFICIENT OF AE/PE FOR THE 1ST ZONE (C1 ) ' , 977 1 1 3 ( ' - ' ) , F 4 . 2 , ' MM * ' , / , 3 0 X , 978 2 ' * COEFFICIENT OF AE/PE FOR THE 2ND ZONE (C2 ) ' , 979 3 1 3 ( ' - ' ) . F 4 . 2 , ' MM * ' , / , 3 0 X , 980 4 ' * COEFFICIENT OF AE/PE FOR THE 3RD ZONE (C3 ) ' , 981 513( ' - ' ) , F 4 . 2 , ' MM * ' , / , S O X , 982 6 ' * It SINCE NO AVAILABLE SOIL MOISTURE IN 4TH ZONE C4=0 . ' , 11X , ' * ' 983 7 , / , 3 0 X , 6 8 ( ' * ' ) ) 984 RETURN 985 END 986 SUBROUTINE TITLE 1 (STA , JYEAR.RMAX,SMCF,HMAX) / 93 987 DIMENSION STA(10) 988 WRITE(6,10) STA,JYEAR,RMAX,SMCF,HMAX 989 10 FORMAT('1' ,30X,'TABLE PREDICTED DAILY WATER TABLE DEPTHS A 990 1ND WATER BALANCE PARAMETERS'./ / ,48X, 991 2'STATION : ' ,10A2,12X,'YEAR : 19' .12, / ,48X, 'RMAX * ' . F 4 . 1 , 992 3' MM/DAY',4X,'R=4K(2DH+H**2)/L**2 FOR 400<Z<HMAX',/,48X, 993 5'AVAILABLE WATER IN TOP THREE ZONES = ' . F 4 . 0 . ' MM' , / ,48X, 994 6'SUBSURFACE DRAIN DEPTH = ' . F 5 . 0 , ' MM' , / , 20X ,100 ( ' - ' ) ) 995 WRITE(6,20) 996 20 FORMAT('+' , / ,21X, 997 1'D M Y ' , / , 2 1 X , ' A 0 E 998 2 ' . / , 2 1 X , ' Y N A PE AE PRE RO DP SD 999 3 TR SMC 1 SMC2 SMC3 DELTR DELZ Z H ' , / . 2 4 X , ' T 1000 4 R MM MM MM MM MM MM MM MM MM 1001 5MM MM MM MM ' , / , 2 4 X , ' H ' , / . 2 0 X , 1 0 0 ( ' - ' ) ) 1002 RETURN 1003 END
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Water table depth simulation for flat agricultural...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Water table depth simulation for flat agricultural land under subsurface drainage and subirrigation practices Chao, Ena C. Y. 1987
pdf
Page Metadata
Item Metadata
Title | Water table depth simulation for flat agricultural land under subsurface drainage and subirrigation practices |
Creator |
Chao, Ena C. Y. |
Publisher | University of British Columbia |
Date Issued | 1987 |
Description | Drainable porosity as a function of water table was investigated to replace the common practice of treating it as a constant A continuous function in the form of a negative exponential equation relating drainable porosity to water table depth was developed by three methods: (1) laboratory core-sample analysis; (2) rainfall rate and water table depth analysis; (3) drainage rate and water table depth analysis. Furthermore, this function was derived for four different water table regimes: (1) subsurface drainage; (2) low subirrigation and subsurface drainage; (3) high subirrigation and subsurface drainage; (4) no drainage and no subirrigation. The drainable porosity function was incorporated into a water balance model which simulated the soil moisture profile and the water table depth on a daily basis. Major modification of the previous model was the elimination of separate falling and rising water table equations since discrete porosity values were no longer assigned to particular soil depth intervals. A subroutine program which computed the total maximum transient storage and the transient storages to each of the four successive soil zones was also incorporated. The 'maximum drainable porosity' and the 'rate constant' parameters in the negative exponential equation were found to be different among the three methods of analysis and among the four water table regimes. Good agreement between simulated and actual water table depths of each regime for 1984 and 1985 was found. The modified water balance model could be used to generate different water table depths by changing the input parameter of design drainage rate. From these outputs, a appropriate drainage rate which gives the desired water table depth could be selected for the purpose of horizontal subsurface drainage system design. |
Subject |
Water table |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0097202 |
URI | http://hdl.handle.net/2429/26685 |
Degree |
Master of Applied Science - MASc |
Program |
Biomedical Engineering |
Affiliation |
Applied Science, Faculty of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1987_A7 C42_8.pdf [ 5.52MB ]
- Metadata
- JSON: 831-1.0097202.json
- JSON-LD: 831-1.0097202-ld.json
- RDF/XML (Pretty): 831-1.0097202-rdf.xml
- RDF/JSON: 831-1.0097202-rdf.json
- Turtle: 831-1.0097202-turtle.txt
- N-Triples: 831-1.0097202-rdf-ntriples.txt
- Original Record: 831-1.0097202-source.json
- Full Text
- 831-1.0097202-fulltext.txt
- Citation
- 831-1.0097202.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0097202/manifest