Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Examination of the design procedures for drainage/subirrigation systems in the lower Fraser Valley, British… Prasher, Shiv Om 1982

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1982_A1 P73.pdf [ 8.85MB ]
Metadata
JSON: 831-1.0076973.json
JSON-LD: 831-1.0076973-ld.json
RDF/XML (Pretty): 831-1.0076973-rdf.xml
RDF/JSON: 831-1.0076973-rdf.json
Turtle: 831-1.0076973-turtle.txt
N-Triples: 831-1.0076973-rdf-ntriples.txt
Original Record: 831-1.0076973-source.json
Full Text
831-1.0076973-fulltext.txt
Citation
831-1.0076973.ris

Full Text

EXAMINATION OF THE DESIGN PROCEDURES FOR DRAINAGE/SUBIRRIGATION SYSTEMS IN THE LOWER FRASER VALLEY, BRITISH COLUMBIA by SHIV OM^PRASHER M.Tech., Punjab A g r i c u l t u r a l University, 1978 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES In t e r d i s c i p l i n a r y Studies We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1982 © Shiv Om Prasher, 1982 In presenting this thesis in p a r t i a l fulfilment of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t 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 f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of I n t e r d i s c i p l i n a r y Studies The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: October 9,1981 Abstract Techniques for designing drainage/subirrigation systems in the Lower Fraser Valley of B r i t i s h Columbia are examined in the present study. An attempt was made to formally define the ov e r a l l goal in the design of such systems. Various d i f f i c u l t i e s were encountered in proceeding formally in the design. Therefore, a "branch and bound" approach was used in which a series of studies are conducted to focus our attention on the key issues of the problem and to regroup or eliminate other issues of secondary importance. A simple study showed that drainage designs seemed to be the l i m i t i n g factors in drainage/subirrigation system designs for the Lower Fraser Valley. Therefore, subirrigation design was not considered in further analyses. The drainage requirements for d i f f e r e n t seasons were discussed. It was suggested that a drainage system designed to meet workability requirements in early spring should be more than s u f f i c i e n t to meet other seasonal requirements of interest from a drainage point of view. It was suggested that these requirements w i l l be met by designing a drainage system that ensures at least one workable period of twelve days in March. A Markov chain model was proposed that can simulate the t r a n s i t i o n s in the water table elevations in response to weather. Design curves were presented for some l o c a l s o i l s that can aid designers to perform drainage designs that s a t i s f y requirements of the individual farmers. A study was undertaken to investigate the importance of uncertainty in s o i l parameters on the drainage system design. F i r s t and second order methods of analyzing uncertainty were applied to Hooghoudt's equation of designing drainage systems. The a p p l i c a b i l i t y of the uncertainty approach was extended to the numerical model of designing drainage systems based on the Boussinesq equation. An example problem was solved to i l l u s t r a t e how a drainage design c r i t e r i o n can be formulated when uncertainties due to both the climate and the s o i l parameters are present at the same time. Conclusions were drawn from the present study and recommendations were made for future work. The environmental impacts of a g r i c u l t u r a l drainage were discussed. They are included in the appendix because the main thrust of the thesis was on the design of drainage/subirrigation systems. Also, a design methodology was proposed in the appendix that can be used in designing drainage systems on a steady state basis in the absence of any knowledge about the saturated hydraulic conductivity of the s o i l or the location of the impermeable layer. iv Table of Contents Abstract i i L i s t of Figures v i L i s t of Tables . v i i Acknowledgement v i i i L i s t Of Symbols x INTRODUCTION 1 LITERATURE REVIEW 9 2.1 Drainage/Subirrigation Theory 9 2.1.1 Steady saturated flows 11 2.1;2 Unsteady saturated flow 13 2.1.3 Saturated/unsaturated flow 17 2.2 Design C r i t e r i a 20 2.2.1 Providing workable days 21 2.2.2 Avoiding excessively wet conditions 24 2.2.3 Avoiding excessively dry conditions 27 2.3 Stochastic Approaches to Drainage Design ....28 CLIMATE AND DRAINAGE/SUBIRRIGATION DESIGN 39 3.1 Drainage/Subirrigation System Design 40 3.2 Drainage Design C r i t e r i a for the Lower Fraser Valley 48 3.2.1 Early spring requirements 50 3.2.2 Growing period requirements 53 3.2.3 Early f a l l requirements 53 3.2.4 Significance of early spring design c r i t e r i o n ...54 3.3 The Mathematical Drainage Model 56 3.4 Layout and Design of the Experiment 60 V 3.5 Measurement of Model Input Parameters 62 3.5.1 Pr e c i p i t a t i o n 63 3.5.2 Drainable porosity 63 3.5.3 Hydraulic conductivity 66 3.5.4 Depth to the impermeable layer 71 3.6 Model Validation 73 3.7 Stochastic Model 77 3.7.1 Markov processes i. 79 3.7.2 Stochastic model for drainage design 85 3.7.2.1 Transition matrix 86 PARAMETER UNCERTAINTY IN DRAINAGE DESIGN 107 4.1 Methods of Analysing Uncertainty 108 4.1.1 F i r s t and second order analysis 109 4.1.1.1 Steady state theory 114 4.1.1.2 Nonsteady state drainage design 126 LINKING PARAMETER AND CLIMATIC UNCERTAINTIES IN DRAINAGE SYSTEM DESIGN 136 5.1 The Stochastic Drainage Model 136 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 150 6.1 Conclusions 155 6.2 Recommendations 158 BIBLIOGRAPHY 161 APPENDIX A 171 APPENDIX B : 173 APPENDIX C 184 v i L i s t of Figures 1. Annual Water Balance of the Lower Fraser Valley from Hare and Thomas ( 1 974) 51 2. Schematic of a Drainage/Subirrigation System 57 3. Experimental Plot 61 4. Comparison of the Observed and Calculated Midspan Water Table Height Below the S o i l Surface in Spring 1982 76 5. Cumulative Probability D i s t r i b u t i o n Functions of Net R a i n f a l l 88 6. Probability of Obtaining at least One Workable Period in March as a Function of Drain Spacing for some Lowland So i l s of the Fraser Valley 104 7. Relationship between the Number of Consecutive Workable Days and their P r o b a b i l i t y of Occurrence in March for a Fine Textured S o i l 105 8. Water Table Drawdown at the Midspacing Computed from the Boussinesq Equation 131 9. Growth of Uncertainty in the Output as a Function of Time 132 10. Standard Deviation of Water Table Height above the Drain Centre Line (two days after the beginning of Drawdown ) 133 11. Probability of Obtaining at least One Workable Period in March as a function of Drain Spacing 148 v i i L i s t of Tables 1. Comparison of Drainage and Subirrigation System Designs for the Lower Fraser Valley 47 2. Results of the Auger Hole Tests and the Indirect Method to Measure Hydraulic Conductivity ....70 3. Input Data Used in the Model Validation 75 4. Output from the Drainage Model based on the Boussinesq Equation 87 5. Winter Transition Matrix 90 6. Spring Transition Matrix 91 7. Steady State P r o b a b i l i t i e s for the Winter Transition Matrix 94 8. Input Data Used in the Comparison of the Chieng Model and the Proposed Model 102 9. Results from the Parameter Uncertainty Analysis of the Hooghoudt Equation 125 10. Net R a i n f a l l Matrices 139 11. Transition Matrices 141 12. Performance of Various Drainage System Designs in the Lower Fraser Valley from a Workability Point of View 146 v i i i Acknowledgement I wish to express my appreciation and feeling of indebtedness to Dr. S. 0 . Russell, Professor of C i v i l Engineering, for his guidance, keen interest, constructive c r i t i c i s m and ever available help throughout the course of t h i s investigation. Thanks are due to Dr. T. H. Podmore, Associate Professor of A g r i c u l t u r a l Engineering in Colorado State University, for introducing me to t h i s problem and his guidance and encouragement during the i n i t i a l phase of the present invet igat ion. I am thankful to Drs. J. DeVries, T. A. Black, R. A. Freeze and S. T. Chieng for very thought-provoking discussions at various stages of the study. I would l i k e to single out Drs. J. DeVries and S. T. Chieng who were always ready to discuss the various aspects of the problem and my discussions with them were always very rewarding. Thanks are due to Ms. Evelyn Tischer for her thought-provoking discussions on the new methodology of designing drainage systems reported in Appendix B. I must thank the National Hydrology Research Institute of Environment Canada, Ottawa, for giving me the opportunity to look at some of the possible environmental impacts of a g r i c u l t u r a l drainage on watershed hydrology. I am especially thankful to Mr. A. Vandenburgh, Drs. F. Morton, V. Klemes and J . G i l l i l a n d for their suggestions and support. Sincere ix thanks are due to the staff of NHRI who made my short stays at Ottawa very pleasant. The f i n a n c i a l support provided by the B r i t i s h Columbia Ministry of Agriculture and Food and the University of B r i t i s h Columbia Graduate Fellowship are gr a t e f u l l y acknowledged. Thanks are also due to a special friend of the family, Dr. P. S. Sinha, for his constant support and words of encouragement during the course of the present study. I w i l l not be doing j u s t i c e to myself i f I do not acknowledge the hard times and s a c r i f i c e s which my wife, Sunita, and my daughter, Preeti, had to make during the course of the present study. This thesis would not have seen this day without their whole-hearted support and their love for me. I w i l l also remain indebted to my parents, s i s t e r s and brother for their good wishes and continued moral support in completing t h i s work. In the end, I am feeling very sad that my grandparents could not l i v e t h i s long to see me f u l f i l l i n g t heir dream. It was due to their blessings and endless a f f e c t i o n for me that I have completed t h i s work. I dedicate the present work in their unperishable memory. X LIST OF SYMBOLS A c r o s s - s e c t i o n a l a r e a p e r p e n d i c u l a r t o the f l o w d i r e c t i o n ( L 2 ) C(p) s p e c i f i c m o i s t u r e c a p a c i t y ( L 3 / L 3 . L ) CS c r o p s u s c e p t i b i l i t y f a c t o r d d e pth of the impermeable l a y e r below the d r a i n c e n t r e (L) D d r a i n depth below the s o i l s u r f a c e (L) E[ ] e x p e c t e d v a l u e f d r a i n a b l e p o r o s i t y ( L 3 / L 3 ) h h e i g h t of the water t a b l e above the d r a i n c e n t r e (L) I h y d r a u l i c g r a d i e n t (L/L) K(p) u n s a t u r a t e d h y d r a u l i c c o n d u c t i v i t y f u n c t i o n of the s o i l (L/T) K s a t u r a t e d h y d r a u l i c c o n d u c t i v i t y of the s o i l (L/T) MFLP m a t r i c f l u x p o t e n t i a l ( L 2 / T ) p.. t r a n s i t i o n p r o b a b i l i t i e s of the water t a b l e (%) P p r o b a b i l i t y v e c t o r (%) p p r e s s u r e head (L) q s p e c i f i c d i s c h a r g e (L/T) Q d i s c h a r g e r a t e (L'/T) R net s o u r c e or s i n k (L/T) S d r a i n s p a c i n g (L) SDI s t r e s s day i n d e x x i SD s t r e s s day f a c t o r S.D. S t a n d a r d d e v i a t i o n V a r ( ) v a r i a n c e V v o l u m e t r i c f l u x per u n i t a r e a (L./T) Z g r a v i t a t i o n a l head (L) A S change i n water s t o r a g e per u n i t s u r f a c e a r e a ( L 3 / L 2 ) n 's s t e a d y s t a t e p r o b a b i l i t i e s of the Markov c h a i n s (%) A t time s t e p (T) *e i n i t i a l midspan h e i g h t of the water t a b l e above t h d r a i n s (L) Z change i n water t a b l e d e pth (L) 1 CHAPTER I INTRODUCTION Ag r i c u l t u r a l drainage involves orderly removal and disposal of excess water from both above and below the land surface by means of surface ditches and/or subsurface drains. It i s mainly used as a means of providing good s o i l and water management in humid regions and to reclaim saline s o i l s in semi-arid regions. Drainage and i r r i g a t i o n were regarded in the Mar del Plata Action Plan (UN Water Conference, 1977) as a major component of the Action Programme on Water for Agriculture. The conference recommended that i r r i g a t i o n and drainage schemes be implemented with great p r i o r i t y in a co-ordinated manner at both national and international l e v e l s . It i s estimated (Bulavko, 1971) that there are over 3.5 m i l l i o n square kilometres of marshland and other waterlogged lands in the world, that i s 2.9% of the t o t a l land area. To s a t i s f y the future world food needs, more and more of this presently uncultivable land w i l l have to be drained in addition to a more intensive and e f f i c i e n t use of land presently under a g r i c u l t u r a l production. Under certain circumstances, a drainage system can also be used for i r r i g a t i o n (subirrigation) by a r t i f i c i a l l y supplying water below the s o i l surface through a system of p a r a l l e l drains or ditches. This method of i r r i g a t i o n i s recommended in areas where the land i s r e l a t i v e l y f l a t and where alternate seasons of water d e f i c i t s and water surplus occur. This method of i r r i g a t i o n may have some li m i t a t i o n s in certain fine textured 2 s o i l s due to very low water conducting capacity. Such an approach where the same system i s used for both drainage and subi r r i g a t i o n could be i d e a l l y suited for the Lower Fraser Valley of B r i t i s h Columbia (Goldstein, 1978). C l i m a t i c a l l y , there i s an excess of p r e c i p i t a t i o n in winter and often substantial water d e f i c i t s in the summer. Therefore, under natural conditions during winter and early spring the water table i s near the s o i l surface and in summer i t recedes into the s o i l . Both of these phenomena are undesirable from the point of view of optimal crop production, but both could be avoided by a single drainage/subirrigation system. At present, the Engineering Branch, B r i t i s h Columbia Ministry of Agriculture and Food, Abbotsford makes recommendations on drainage design based on steady state theory (Hooghoudt's method) which assumes that over longer periods of time the position of the water table remains unchanged and that t o t a l inflow into the system equals t o t a l outflow from the system. The steady state water table i s a special condition of the real s i t u a t i o n in the f i e l d where the water table is continually f l u c t u a t i n g . A more r e a l i s t i c design procedure should take th i s into account and aim at the prevention of undesirable patterns of water table fluctuations. There are some models available such as those of Skaggs(1975), Skaggs and Tang(l976), Skaggs(1978), Van SchiIfgaarde(1974) and Chieng(l975) which do consider water table fluctuations with time. However, because of the s i m p l i c i t y of the steady state theory and the fact that i t seems to give reasonable results as 3 compared with the t r a n s i e n t case which r e q u i r e s r e l a t i v e l y complex a n a l y s i s and has grea t e r need f o r data, the En g i n e e r i n g Branch f a v o r s the steady s t a t e theory. When i t comes to the combined d r a i n a g e / s u b i r r i g a t i o n system design, the E n g i n e e r i n g Branch does not have any design c r i t e r i a except that i t i s g e n e r a l l y f e l t by v a r i o u s experts i n the Fra s e r V a l l e y (Baehr, 1980) that the d r a i n spacing should be made smaller i f the system i s to work f o r both drainage and s u b i r r i g a t i o n r a t h e r than j u s t d rainage. In d e s i g n i n g a d r a i n a g e / s u b i r r i g a t i o n system, the key parameters are d r a i n spacing and d r a i n depth. In formal d e c i s i o n theory terms, the o v e r a l l goal i n the design of a system where u n c e r t a i n t y i s i n v o l v e d i s to maximize the d e c i s i o n maker's expected u t i l i t y . A u t i l i t y f u n c t i o n d e s c r i b e s the r e l a t i v e d e s i r a b i l i t y t h a t an i n d i v i d u a l a s s i g n s to the p o s s i b l e outcomes i n a given s i t u a t i o n . Expected u t i l i t y i s obtained by summing the product of u t i l i t i e s and t h e i r l i k e l i h o o d s over a l l p o s s i b l e s t a t e s of the proc e s s . In the case of a l a r g e w e l l - f i n a n c e d e n t e r p r i s e and where there i s no r i s k to l i f e and limb, u t i l i t y i s adequately represented by monetary b e n e f i t s - and over time the most s u c c e s s f u l business w i l l be the one which always seeks to maximize i t s expected monetary b e n e f i t s . However, i n a small underfinanced business where a l a r g e l o s s i n one year c o u l d endanger i t s e x i s t e n c e , the primary o b j e c t i v e may be to s u r v i v e , that i s , to minimize the p r o b a b i l i t y of unacceptable l o s s e s . A farmer may wish to i n s t a l l a d r a i n a g e / s u b i r r i g a t i o n 4 s y s t e m i n o r d e r t o m a x i m i z e h i s e x p e c t e d b e n e f i t s o r t o m i n i m i z e t h e r i s k o f l a r g e c r o p l o s s e s i n t h e e v e n t o f a d v e r s e w e a t h e r o r i t may be a c o m b i n a t i o n o f t h e t w o . T h u s i t i s d i f f i c u l t t o e v e n d e f i n e t h e o b j e c t i v e f o r a t y p i c a l s y s t e m . A c c e p t i n g f o r t h e moment t h a t t h e o b j e c t i v e i s t o m a x i m i z e e x p e c t e d b e n e f i t s , t h e r e a r e s t i l l d i f f i c u l t i e s . The n e t b e n e f i t f r o m a f a r m i n g o p e r a t i o n i s t h e d i f f e r e n c e b e t w e e n r e c e i p t s f r o m t h e s a l e o f p r o d u c e a n d t h e c o s t s o f p r o d u c t i o n . The b e n e f i t f r o m a d r a i n a g e / s u b i r r i g a t i o n s y s t e m i s t h e e x t r a b e n e f i t o b t a i n e d b e c a u s e o f t h e s y s t e m l e s s i t s c o s t . M a x i m i z i n g t h e e x p e c t e d b e n e f i t s i s e q u i v a l e n t t o m i n i m i z i n g t h e sum o f t h e c o s t o f t h e s y s t e m p l u s t h e e x p e c t e d c r o p l o s s e s . B u t e s t i m a t i n g e x p e c t e d c r o p l o s s e s i s d i f f i c u l t b e c a u s e a n i n d e x o f s y s t e m p e r f o r m a n c e a n d a r e l a t i o n s h i p b e t w e e n t h e i n d e x a n d c r o p l o s s e s i s n e e d e d . Some i n d i c e s o f s y s t e m p e r f o r m a n c e a n d t h e i r r e l a t i o n s h i p s w i t h c r o p l o s s e s a r e a v a i l a b l e i n t h e l i t e r a t u r e ( H i l e r , 1 9 6 9 ; H i l e r a n d C l a r k , 1 9 7 1 ) b u t t h e i r l a r g e d a t a r e q u i r e m e n t s r u l e o u t t h e p o s s i b i l i t y o f u s i n g t h e m i n m o s t p r a c t i c a l s i t u a t i o n s . F o r t h e L o w e r F r a s e r V a l l e y , s u c h d a t a i s n o t y e t a v a i l a b l e a n d p e r h a p s w i l l n o t be a v a i l a b l e i n t h e n e a r f u t u r e . M o r e o v e r t h e d e s i g n o f a d r a i n a g e / s u b i r r i g a t i o n s y s t e m h a s t o i n c o r p o r a t e u n c e r t a i n t i e s i n a r a t i o n a l w a y . t h a t s t e m f r o m b o t h w e a t h e r a n d s o i l p a r a m e t e r s . G i v e n t h e d i f f i c u l t y o f p r o c e e d i n g f o r m a l l y i n t h e d e s i g n o f a d r a i n a g e / s u b i r r i g a t i o n s y s t e m a n e n g i n e e r i n g a p p r o a c h c a l l e d " b r a n c h a n d b o u n d " may be u s e d . I t w o r k s on t h e p r i n c i p l e t h a t i f a n u n d e r e s t i m a t e o f t h e c o s t o f a l t e r n a t i v e A 5 exceeds an o v e r e s t i m a t e of a l t e r n a t i v e B then a l t e r n a t i v e B has t o be cheaper than A. O f t e n a s i m p l e study can show t h a t one a l t e r n a t i v e i s always i n f e r i o r t o a n o t h e r and t h u s i t can be q u i c k l y e l i m i n a t e d from f u r t h e r c o n s i d e r a t i o n . A l s o i t may be found t h a t one r e q u i r e m e n t i s of such importance t h a t i t o v e r r i d e s a l l o t h e r s . The g e n e r a l p r o c e d u r e i s t o conduct a s e r i e s of s t u d i e s t o f o c u s our a t t e n t i o n on key i s s u e s of the problem and t o r e g r o u p or e l i m i n a t e o t h e r i s s u e s of secondary impo r t a n c e . We d e a l w i t h o n l y a p o r t i o n of the problem at one t i m e , r e v i e w a l l p o s s i b l e a l t e r n a t i v e s a t our d i s p o s a l and take a c t i o n i n f a v o r of the one t h a t i s f e a s i b l e under l o c a l c o n d i t i o n s . F o l l o w i n g a l o n g t h e s e l i n e s a sequence of i n v e s t i g a t i o n s was c o n ducted as f o l l o w s : 1) To check whether d r a i n a g e - or s u b i r r i g a t i o n - s y s t e m d e s i g n i s the l i m i t i n g f a c t o r i n a d r a i n a g e / s u b i r r i g a t i o n d e s i g n 2) To f o r m u l a t e a r e a l i s t i c d e s i g n g o a l f o r a d r a i n a g e / s u b i r r i g a t i o n system f o r the Lower F r a s e r V a l l e y 3) To check on the f e a s i b i l i t y of u s i n g the B o u s s i n e s q e q u a t i o n t o s i m u l a t e d r a i n a g e and s u b i r r i g a t i o n p r o c e s s e s (a r e a l i s t i c nonsteady s t a t e model w i t h o u t the c o m p l e x i t i e s and the d a t a l i m i t a t i o n s of a s a t u r a t e d / u n s a t u r a t e d f l o w a n a l y s i s ) 4) To d e v e l o p a methodology f o r t a k i n g i n t o 6 account u n c e r t a i n t y due t o p r e c i p i t a t i o n and/or e v a p o t r a n s p i r a t i o n when f o r m u l a t i n g d e s i g n c r i t e r i a 5 ) To a p p l y . f i r s t o r d e r a n a l y s i s t o a s s e s s the e f f e c t s of u n c e r t a i n t y i n s o i l p arameters i n d r a i n a g e / s u b i r r i g a t i o n system d e s i g n 6 ) To d e v e l o p a methodology f o r t a k i n g i n t o account u n c e r t a i n t i e s due t o both weather and s o i l parameters when f o r m u l a t i n g d e s i g n c r i t e r i a To meet the above aims, a m a t h e m a t i c a l model based on the B o u s s i n e s q e q u a t i o n i s a p p l i e d t o s i m u l a t e the d r a i n a g e / s u b i r r i g a t i o n p r o c e s s . A f t e r c h e c k i n g t h a t d r a i n a g e d e s i g n i s the g o v e r n i n g f a c t o r i n a d r a i n a g e / s u b i r r i g a t i o n d e s i g n , a d r a i n a g e system d e s i g n i s proposed t o meet w o r k a b i l i t y r e q u i r e m e n t s i n the s p r i n g f o r the Lower F r a s e r V a l l e y . I t i s argued t h a t o t h e r d r a i n a g e r e q u i r e m e n t s such as an upper l i m i t t o h i g h water t a b l e s d u r i n g the growing p e r i o d of c r o p s and t r a f f i c a b i l i t y r e q u i r e m e n t s i n the f a l l a r e i m p l i c i t l y c o n s i d e r e d i n such a d e s i g n . A Markov c h a i n model i s proposed t o d e t e r m i n e the p r o b a b i l i t y of o b t a i n i n g any number of c o n s e c u t i v e workable days i n s p r i n g f o r a g i v e n s i t e and from a g i v e n d r a i n a g e system. A p r a c t i c a l way t o i n c o r p o r a t e parameter u n c e r t a i n t y i n d r a i n a g e d e s i g n i s d i s c u s s e d . Examples are g i v e n f o r d e s i g n s based upon both s t e a d y and nonsteady s t a t e t h e o r y . F u r t h e r i t 7 is shown how a drainage/subirrigation system can be designed that can consider uncertainties due to both weather and s o i l parameters in a rational way. In Appendix A, a new si m p l i f i e d design methodology i s proposed for designing drainage systems without any pre-knowledge about the s o i l hydraulic conductivity or the location of the impermeable layer. Also during the course of the study, the writer became aware of the importance of environmental impacts of a g r i c u l t u r a l drainage on watershed hydrology. Therefore, although the main thrust of the thesis i s on the "micro-level" a c t i v i t i e s on our watersheds, an investigation into some of their impacts i s included in Appendix B for Canadian conditions. The subject matter in t h i s thesis has been arranged in the following way. Chapter 2 reviews various approaches that have been proposed to design a drainage/subirrigation system. The design c r i t e r i a for such a system are discussed in d e t a i l for humid regions. Various approaches that incorporate uncertainty due to climate and s o i l parameters are reviewed and their incorporation in such designs i s discussed. In Chapter 3, i t i s shown that for l o c a l c limatic conditions, a drainage system design performed in the conventional way i s more than s u f f i c i e n t to f a c i l i t a t e s u b i r r i g a t i o n for many s o i l s selected from the l i t e r a t u r e . The design c r i t e r i a for designing drainage systems are discussed for humid regions. It i s shown that i f a drainage system i s designed to s a t i s f y workability requirements in the spring, i t should be more than s u f i c i e n t to meet drainage 8 requirements during the periods of crop growth and harvesting. A transient mathematical model based on the Boussinesq equation is validated for simulating the drainage process in early spring. Problems of characterizing input data for such a model are discussed. F i n a l l y , a Markov chain model i s proposed for the water table fluctuations that incorporates climatic uncertainty into the drainage model in a rat i o n a l way. This model can be used to determine the p r o b a b i l i t y of occurrence of any number of consecutive workable days in early spring (March) for a given location and drainage system. Chapter 4 deals with the uncertainty in drainage designs due to uncertainties in the s o i l properties. It i s shown that a f i r s t order analysis of uncertainty can be applied with l i t t l e e f f o r t to both steady and nonsteady state drainage design methods. Chapter 5 suggests a way to incorporate uncertainties due to both climate and s o i l parameters in the design. It is shown that confidence l i m i t s can be assigned to a given design. It i s hoped that the procedure may enable designers to perform designs that w i l l meet the requirements of individual farmers depending upon the risk they are w i l l i n g to accept in their farming enterpr i se. Chapter 6 presents a summary and the conclusions from t h i s study. Recommendations are made for future research in the area. 9 CHAPTER II LITERATURE REVIEW The goal of thi s chapter i s to summarize developments in drainage/subirrigation system design. Although most of the work c i t e d in the l i t e r a t u r e has focussed on drainage design, i t can be conveniently employed in subirrigation system design with some minor changes in the i n i t i a l and boundary conditions. Therefore, most discussion of drainage theory in the f i r s t section holds true for su b i r r i g a t i o n . Whenever drainage theory is the topic of discussion, subirrigation theory w i l l be i m p l i c i t l y considered in i t . Considering the vast amount of l i t e r a t u r e available on thi s subject, only those works that are d i r e c t l y related to the present study w i l l be treated in depth. 2.1 Drainage/Subirrigation Theory The design of a drainage/subirrigation system i s generally performed by using one of the numerous spacing formulae that relate drainage system parameters to s o i l properties, hydrologic conditions and water table fluctuations. Many of these formulae are si m p l i f i e d versions of the general p a r t i a l d i f f e r e n t i a l equation(PDE) of groundwater flow and therefore, are applicable for idealized conditions. The governing PDE i s obtained by combining Darcy's law with the law of conservation of matter (continuity equation). Ideally, the Navier-Stokes equations should be used for the study of groundwater flow problems (Van Schilfgaarde, 1970). These equations are derived using kinematic, dynamic and 10 thermodynamic relations for f l u i d flow. The f i r s t of these implies a statement of conservation of mass; the second can be derived from a momentum balance and ; the t h i r d takes the form of an equation of state r e l a t i n g the f l u i d ' s pressure, temperature and density. But the solution of these equations requires i d e n t i f i c a t i o n of the appropriate i n i t i a l and boundary conditions. For groundwater flow, i t implies s p e c i f i c a t i o n of the actual geometry of void spaces. This i s not p r a c t i c a l under natural conditions. Therefore, we substitute Darcy's law in place of the momentum balance in the derivation of the governing equation of groundwater flow. In common usage, i t i s said that Darcy's law leads to the consideration of a f i c t i t i o u s flux, c a l l e d s p e c i f i c discharge, in contrast to the actual s p e c i f i c discharge obtainable from Navier-Stokes equations. Darcy's law is defined as v = Q/A =-KI (1 ) where v i s the s p e c i f i c discharge, Q i s the discharge rate, A i s the cross-sectional area normal to the dire c t i o n of flow, K is the hydraulic conductivity and I i s the hydraulic gradient. However, i t is worth noting that Hubbert (1956) derived an expression for Darcy's law from the Navier-Stokes equations which can provide an estimate of the actual flux. Therefore, Darcy's law by i t s e l f does not lead to either f i c t i t i o u s or true values of s p e c i f i c discharge. Equation (1) for Darcy's law has been used in groundwater flow theory because i t allows the s h i f t 11 from the actual s p e c i f i c discharge in each pore to the flux across a unit area of the medium. The development of drainage theory has proceeded along a variety of widely divergent paths. Some researchers have concentrated on the development of mathematically exact solutions to problems that have severe p r a c t i c a l constraints. Others have employed rather s i m p l i s t i c assumptions to get some useful answers to p r a c t i c a l situations. Again, there are some who have preferred a middle path between the two. In the following pages, the above solutions w i l l be described in more d e t a i l . 2.1.1 Steady saturated flows A steady state flow condition exists when the rate of change of mass storage in the flow zone is zero. The boundaries and flow rates of such a system are time invariant. For incompressible flows in homogeneous and isotropic s o i l s , the governing PDE of groundwater flow becomes an e l l i p t i c equation, commonly known as Laplace's equation. This equation is the basis of most drainage analyses based on the potential theory (Kirkham, 1958). Its solution not only provides an estimate of the drain spacing but also the d i s t r i b u t i o n of pressure head in the flow zone. These mathematically exact solutions, although accurate, are seldom applied in p r a c t i c a l designs due to their complex nature and because they do not lend themselves readily to d a i l y use by p r a c t i c i n g engineers (Ravelo, 1978). 12 This has led to the search for approximate solutions, most of which incorporate Dupuit-Forchheimer (D-F) assumptions. These assumptions imply that the flow i s horizontal and that the flow v e l o c i t y i s proportional to the slope of the water table but independent of depth (Bear, 1979). In addition, i t i s normally assumed that the water table i s the streamline bounding the flow region. Therefore, the application of these approximate solutions should be limited to flow systems that are shallow r e l a t i v e to their areal extent and that have small water table curvature. Convergence losses and surfaces of seepage can not be determined by using these s i m p l i f i e d solutions (because the loss of hydraulic gradient in the v e r t i c a l d i r e c t i o n is ignored). Kirkham (1967) showed that there i s no difference between the potential theory and D-F theory i f the soil-has i n f i n i t e hydraulic conductivity in the v e r t i c a l d i r e c t i o n . For other s o i l s , solutions w i l l be approximate. I n f i n i t e v e r t i c a l hydraulic conductivity may be approximated by a s o i l that has numerous worm holes, root channels and cracks (preferred pathways). Hooghoudt (1940) was probably the f i r s t to obtain an approximate solution where r a d i a l flow effects were considered. This theory does not give the shape of the water table; only the midspan water table height i s obtained. However, the Hooghoudt formula considers both the r a d i a l flow near the drains and the nearly horizontal flow at greater distances from the drains. It may be expressed as: 1 3 S 2 = 4 K h ( 2 d + h ) / R where S i s the drain spacing, K i s the saturated hydraulic conductivity of the s o i l , h i s the permissible midspan water table height above the drains, d i s the depth to the impermeable layer below the drains and R i s the drainage c o e f f i c i e n t . Hooghoudt suggested that convergence of flow l i n e s near the drains can be considered i m p l i c i t l y in the above formula i f the parameter d i s replaced by some . equivalent depth of the impermeable layer below the drains. Several decades after i t s f i r s t appearance, the Hooghoudt's formula i s s t i l l favoured by pra c t i c i n g engineers for humid areas because of i t s s i m p l i c i t y and ease of application. Wesseling (1964) carried out extensive calculations using potential theory and compared the results with Hooghoudt's formula. He showed that results obtained using Hooghoudt's formula did not vary more than about 5% from those of potential theory for midspan water table heights. It should be mentioned that the above formulae are applicable only when steady state conditions occur. Since steady state conditions rarely apply, these formulae should be used with care. 2.1.2 Unsteady saturated flow The solutions to unsteady flow problems are more d i f f i c u l t than those to steady flow problems. Yet, in nature, steady 14 state conditions for the time spans of interest are seldom encountered. The solution of nonsteady saturated flow problems requires some qua n t i f i c a t i o n as to the volume of water stored in or drained from the s o i l p r o f i l e with the ri s e or f a l l of water table. A common practice i s to consider a time invariant s p e c i f i c y i e l d or drainable porosity. Childs (1960a) has pointed out that this assumption can lead to s i g n i f i c a n t errors. But, in view of the complicated nature of the unsteady problems, i t s t i l l remains the most p r a c t i c a l assumption available for the development of theory useful in p r a c t i c a l designs (Van Schilfgaarde, 1970). Laplaces's equation derived from potential theory can be applied to nonsteady flow problems by introducing time-dependent boundary conditons. However, the c a l c u l a t i o n of the location of the water table with such solutions is quite d i f f i c u l t . Also, the boundary condition along the water table i s quadratic in terms of the derivatives of hydraulic head (Davis and De Wiest, 1966). Moreover, the equation i t s e l f i s not t h e o r e t i c a l l y correct because i t ignores any water storage in the flow zone. Therefore, the other approach that has been used extensively i s the one based on the D-F assumptions. Both of these approaches can be employed to study the rate of r i s e or f a l l of water table from a known i n i t i a l condition in the presence or absence of pr e c i p i t a t i o n (Van Schilfgaarde, 1970 and 1974). The PDE resulting from incorporating the D-F assumptions in the general groundwater flow equation results in what i s 15 commonly known as the Boussinesq equation f l i l * K _ J _ (h 3h) ^ R 3t 3x 3x " where h i s the height of the water table above the impermeable layer, X is the horizontal distance, t i s time, f i s drainable porosity, K i s saturated hydraulic conductivity and, R is the rate at which water enters or leaves the saturated region and is negative for evapotranspiration or deep seepage and positive for recharge. This i s a nonlinear PDE and has been solved a n a l y t i c a l l y in the absence of a source or a sink for the case when drains are lying on the impermeable layer (Boussinesq, 1903; Dumm, 1954). Recently closed form solutions have been obtained in the presence of recharge (Basak, 1979; Singh and Rai, 1980). When drains do not l i e on top of an impermeable layer, the above a n a l y t i c a l solutions break down because one of the boundary conditions for drainage can no longer be s a t i s f i e d (Van Schilfgaarde, 1970). To circumvent this d i f f i c u l t y , the Boussinesq equation has been l i n e a r i z e d (Polubarinova-Kochina, 1962). It has been l i n e a r i z e d in many ways by making further simplifying assumptions (Glover as c i t e d in Dumm, 1954; Werner, 1957; Krayenhoff van de Leur, 1958; Massland, 1959; T e r z i d i s , 1968; Singh and Jacob, 1977; and Prych, 1980). Van Schilfgaarde (1963) found an a n a l y t i c a l solution to equation (3) in the absence of a source or a sink and for the case when drains do not l i e on top of an impermeable layer 16 without resorting to any l i n e a r i z a t i o n . Although his solution s t i l l f a i l s to s a t i s f y one of the boundary conditions t h i s problem can be a l l e v i a t e d by considering the water table recession as a sequence of small steps (Van Schilfgaarde, 1964). In the presence of a source or sink at the s o i l surface and for the case when drains do not l i e on top of the impermeable layer, there are very few closed form solutions available (Van SchiIfgaarde, 1965). In most cases, numerical techniques are employed to arrive at approximate solutions (Moody, 1966; Skaggs, 1975). An obvious flaw in the saturated flow theory i s that the contribution of the unsaturated flow region i s not accounted for in c a l c u l a t i o n s . The source or sink i s assumed to be located at the water table instead of at the s o i l surface. A correction i s sometimes made by adjusting the height of the water table or the impermeable layer by an amount equal to the height of the c a p i l l a r y fringe or a i r entry value (Childs, 1960a and 1960b; Bouwer, 1964a). Sometimes a variable drainable porosity i s used in the calculations rather than a constant one to account for the unsaturated flow region (pseudo saturated/unsaturated flow). Skaggs and Tang (1976) show that when using this approach, the solutions from the Boussinesq equation are not d i f f e r e n t from the solutions obtained by solving the governing PDE for saturated/unsaturated flow. However, the above methods do require determination bf the water retention c h a r a c t e r i s t i c s of the s o i l . Although these c h a r a c t e r i s t i c s are comparatively easier to determine than the unsaturated hydraulic conductivity 17 c h a r a c t e r i s t i c s , nevertheless i t i s a s o i l property that is not routinely measured and therefore needs extra e f f o r t . 2.1.3 Saturated/unsaturated flow The separation of flow in the unsaturated region from saturated flow i s arbitrary (Van Schilfgaarde, 1970). Physically, the flow system i s continuous in the s o i l (continuous flow lines and head d i s t r i b u t i o n ) . The only j u s t i f i c a t i o n for making the d i s t i n c t i o n can be found in the mathematical complexities and data lim i t a t i o n s associated with the unsaturated flow. The appropriate PDE for a transient system that includes both saturated and unsaturated flow in two dimensions i s Ix IMP) S • !l ««PHH 4 i » • % H-"c(p>H where K(p) i s the unsaturated hydraulic conductivity as a function of pressure head p, C(p) i s the s p e c i f i c moisture capacity of the s o i l as a function of pressure head. The above equation can be written in another form making moisture content the independent variable instead of pressure head (Freeze and Cherry, 1979). The solution of this equation along with the appropriate i n i t i a l and boundary conditions for drainage or sub i r r i g a t i o n is quite d i f f i c u l t . For drainage/subirrigation design there are very few solutions available in the l i t e r a t u r e and they are 18 mainly of academic importance (Wind and Van Doorne, 1975; Tang and Skaggs, 1977; Vauclin et a l . , 1979). This i s because f i e l d - e f f e c t i v e values of soil-water properties such as K(p) and C(p) are d i f f i c u l t to measure. Moreover, the f i e l d v a r i a b i l i t y of the s o i l hydraulic properties w i l l negate the additional s e n s i t i v i t y obtained from the more precise methods (Skaggs, 1975). Also, the computational time required for these analyses is quite prohibitive and therefore, these solutions do not lend themselves e a s i l y for use in p r a c t i c a l designs, es p e c i a l l y i f they are to be based on long-term weather records. This has led to a d i f f e r e n t type of modeling in drainage/subirrigation design, namely the water balance approach (Chieng, 1975; Skaggs, 1978; Bhattacharya and Broughton, 1979). It i s a mathematical book-keeping operation to compute changes in s o i l moisture in response to p r e c i p i t a t i o n , evapotranspiration, drain outflow, surface runoff and deep seepage. A water balance model i s generally expressed by the equation: inflow = outflow + change in storage More s p e c i f i c a l l y , the water balance model of Chieng (1975) uses the following equation: SMC(t) = SMCU-1) + PRE (t) - AE (t) - EX(t) 19 where SMC(t) i s the s o i l moisture content at time t, SMC(t-1) i s the s o i l moisture content at time t-1, PRE(t) i s the recharge (p r e c i p i t a t i o n and/or i r r i g a t i o n ) at time t, AE(t) i s the actual evapotranspiration at time t and EX(t) i s the excess water at time t. The fluctuations in the water table height at the midspacing are calculated from the equation: At T At Where A s i s the change in water storage per unit surface area over timefct, f is the drainable porosity and Z i s the change in water table depth in time A t . The term A s / A t is the rate of change water storage that can be drained under gravity. Replacing i t by q, the drainage c o e f f i c i e n t , we get q * f A t At or Az - q J-20 The above equation was used in Chieng's computer model to simulate d a i l y changes in the midspan water table height. The water balance models can not give any indication about the shape of the water table p r o f i l e (water budgeting usually takes place at one point between the drains - usually the midspan). Also, although these models can work for long-term weather records, they are not as accurate as models based upon the p a r t i a l d i f f e r e n t i a l equation governing the saturated/unsaturated groundwater flow due to various assumptions that are made to keep the analysis both tractable and v e r s a t i l e (Chieng, 1975; Skaggs, 1978). 2.2 Design C r i t e r i a The a g r i c u l t u r a l function of a drainage/subirrigation system is to help maximize the benefit/cost r a t i o of the farm enterprise. A proper design c r i t e r i o n should enable the designer to decide on the optimum subsurface drain spacing such that a favourable benefit/cost relationship holds for farmers. Therefore, the design c r i t e r i o n i s an indispensable part of a drainage/subirrigation system design (and, s u r p r i s i n g l y , i t is the least understood component of the design). Ideally, a design c r i t e r i o n should consider various effects that a drainage/subirrigation system may have on various crops, s o i l s , farm operations and environment. There i s a d e f i n i t e lack of knowledge about these e f f e c t s due to the complexity of geo-hydrologic processes and the paradoxical nature of their interactions. Although many designs done in the past have 21 withstood the test of time without precise knowledge about these e f f e c t s , they have been done "more often than not based on guides derived from experience rather than on a n a l y t i c a l formulations" (Van Schilfgaarde, 1979). The aims of a drainage/subirrigation system in humid regions can be divided into three categories: 1) Providing workable days in early spring and during harvesting 2) Avoiding conditions that are too wet during the crop growth per iod 3) Avoiding conditions that are too dry during the crop growth period 2.2.1 Providing workable days More harm i s often done to the farm enterprise due to the interference of excess water with the timeliness of farming operations than by the direct adverse effects of excess water on plant growth (Smedema, 1979). The losses may vary from complete crop f a i l u r e , i f planting i s delayed too long, to reduced yie l d s i f t i l l a g e , weed control, harvesting are not performed on time (Reeve and Fausey, 1974). Under the tight calendar of modern farming, a delay in one operation can cause delay in subsequent operations. When the operation i s not delayed but performed 22 under non-workable conditions, the effectiveness of the operation i s quite low (machinery getting stuck in s o i l or lack of t r a c t i o n ) . When an operation i s performed under wet conditions the damage to s o i l structure may be serious and repercussions may be f e l t in subsequent years. It i s mainly the formation of a t r a f f i c pan that not only r e s t r i c t s downward movement of water but also hinders root growth. Once formed, i t s e f f e c t s could be f e l t in years to come unless i t i s broken up by sub-soiling. Workable conditions exist when the moisture content of the surface s o i l i s not so high as to impede farm t r a f f i c , seed bed operations or harvesting. A d i s t i n c t i o n should be made between workable and t r a f f i c a b l e s o i l s . Whereas t r a f f i c a b i l i t y relates to the a b i l i t y of the s o i l to support t r a f f i c without being s t r u c t u r a l l y damaged beyond l i m i t s for good crop growth (Paul and deVries, 1979), workability includes the above d e f i n i t i o n as well as the f r i a b l e state of the s o i l (when s o i l w i l l f a l l into small crumbs, not. into clods and d e f i n i t e l y not be smeared and puddled). A workable s o i l i s also t r a f f i c a b l e but the converse may not be true for some s o i l s . Most of the farming operations can be c a r r i e d out i f the tension in the top s o i l is 100 cm of water or more (Oosterbaan and Wind, 1979). Skaggs (1978) defines a workable day as a day when the a i r volume in the s o i l p r o f i l e exceeds some l i m i t i n g value; the r a i n f a l l occurring on that day i s less than a minimum value; and that a minimum number of days have elapsed since that amount of r a i n f a l l occurred. For Wagram loamy sand, the 23 suggested values are 3.7 cm3/cm2, 1.2 cm and 1 day respectively. Smedema (1979) recommends the use of moisture content at the lower p l a s t i c l i m i t of the s o i l as a basis for estimating workability. Further, he shows that i f the moisture content at the lower p l a s t i c l i m i t value of the s o i l i s lower than i t s f i e l d capacity moisture content, no s i g n i f i c a n t improvement in the s o i l workability can be achieved by narrowing the drain spacings. However, increasing drain depths can make the s o i l workable much e a r l i e r (Wind, 1976). This i s mainly due to the corresponding increase in hydraulic gradients that result in increased drainage rates and therefore, the top s o i l becomes workable quite early. Paul and deVries (1983) report on the development of two procedures for the prediction of t r a f f i c a b i l i t y of ti l e - d r a i n e d lowland s o i l s in the Lower Fraser Valley. They related the tension of water in the top s o i l to the s o i l strength and related the l a t t e r to water table depth for two di f f e r e n t s o i l s (Hallart s i l t y clay loam and Lumbum muck). C r i t i c a l water table depths for t r a f f i c a b i l i t y were found to be 53, 45, and 60 cm for Lumbum muck, Hallart s i l t y clay loam (grassland), and Hallart s i l t y clay loam (cultivated land) respectively. Most of the above studies are based on saturated/unsaturated flow theory in predicting workable days. Yet, as explained e a r l i e r , the use of thi s theory in p r a c t i c a l designs is very limited. Therefore, p r a c t i c i n g engineers w i l l have a tendency to select simple indicators for workable days such as depth to the water table. However, more research i s 24 required to establish t h i s r e l a t i o n s h i p for many s o i l s . 2.2.2 Avoiding excessively wet conditions Excessively wet conditions as such are not generally detrimental to plant growth. They should be avoided during the growth period because they cause lack of aeration (poor gaseous exchange), degradation of s o i l structure, loss of nitrogen by d e n i t r i f i c a t i o n . Too wet conditions can cause excessive overland flow that may result in s o i l erosion and in flooding in the downstream reach of the basin. Moreover, s o i l temperature increases more slowly in spring when the moisture content of the s o i l is high. Low s o i l temperature reduces a v a i l a b i l i t y of nutrients and water uptake. Germination and seedling growth are also affected by low s o i l temperature. S o i l moisture conditions in the top s o i l are related to the water table depth below the s o i l surface. Therefore, a good approach to judge the effectiveness of a drainage system would be by studying fluctuating water tables under the l o c a l weather conditions in drained lands. S t a t i s t i c a l l y , i t i s then possible to compute frequency of exceedence of certain water levels for a given drainage system. This information can then be used to determine good drainage/subirrigation design for a given s i t e (Van Schilfgaarde,1965; Chieng, 1975). An optimal drainage system could be selected i f the frequency tables of water l e v e l exceedences could be translated into appropriate loss functions. This s t i l l remains the basic problem in drainage system design (Oosterbaan and Wind, 1979) 25 and hence constitutes a weak component of drainage/subirrigation system design (Williamson and K r i z , 1970; Bouwer, 1974). It should be mentioned that the use of the water table as an indicator of excess moisture stress i s born out of convenience. Actually, the flow of water into the plant roots from the s o i l i s governed by the water potential difference between plant roots and s o i l . Our i n a b i l i t y to model these potentials has led to the search for alternative (and preferably simpler) c r i t e r i a , and water table elevation seems to be a sa t i s f a c t o r y indicator (Bouwer, 1974). Sieben(l964) studied the effects of water table fluctuations on crop y i e l d s . He recommended the use of a SEW30 concept which i s a measure of high water table above the c r i t i c a l depth of 30 cm. Water table heights that are within 30 cm of the s o i l surface are summed over a period (growing period, f u l l year). Larger SEW values indicate excesssive s o i l moisture conditions and should be prevented. A flaw commonly pointed out in t h i s concept is that i d e n t i c a l SEW values may not necessarily imply i d e n t i c a l conditions (Wesseling, 1974). H i l e r (1969) refined t h i s idea and proposed the use of a stress-day index (SDI) concept to characterize crop water requirements. The SDI i s determined by multiplying a crop s u s c e p t i b i l i t y factor and a stress day factor. Then the product is summed over a l l the crop growth stages. The crop s u s c e p t i b i l i t y factor i s a function of the species and stage of development of crops, e. g., loss to a corn crop may be d i f f e r e n t from loss to a pea crop for the same le v e l of excess 26 water stress; loss to corn during the reproduction stage w i l l be more than at any other stage for the same l e v e l of excess water in the root zone. The stress day factor i s a measure of the intensity and duration of the stress to plants. Mathematically, SDI can be expressed as ° N SDI = I (CS. * SD.) i = l 1 1 where CS^ i s the crop s u s c e p t i b i l i t y factor for the i t h growth stage, SD£ i s the stress day factor for the i t h growth stage, and N i s the number of growth stages. By minimizing the SDI value, the crop y i e l d would be maximized (Ravelo, 1978). This i s due to the inverse r e l a t i o n between the crop y i e l d and SDI. The stress day factor may stand for oxygen deficiency, an osmotic stress, or SEW30. The crop s u s c e p t i b i l i t y factor can be determined only from f i e l d studies (Hiler and Clark, 1971). The above method of characterizing crop losses due to a non-optimal environment in the crop root zone requires a substantial amount of l o c a l data which, in many cases, i s seldom av a i l a b l e . Therefore, research e f f o r t s are required either to look for another appropriate c r i t e r i o n that i s not too demanding with respect to data requirements or to f i n d crop response data for l o c a l conditions. If researchers could specify what data were required, i t might be possible for an organizations such as 27 the Research Branch of Agriculture Canada or the USDA to begin to measure the necessary parameters on a long-term basis. 2.2.3 Avoiding excessively dry conditions In humid regions, y i e l d decreases can occur due to water d e f i c i t s in the summer. Less water i s available to plants in such conditions. This results in lower a v a i l a b i l i t y of nutrients to growing plants. Stomata close, causing a decline in the transpiration rate. Photosynthesis, protein dehydration, c e l l enlargement and d i v i s i o n are affected. The end result i s that there i s a-decline, in the t o t a l dry matter production. Most of the crop loss functions in this case are s t a t i c , that i s , they f a i l to consider the dynamics of the continuous growth process of plants (Morgan et a l . , 1980). Some of the dynamic models(Tanner, 1981; Bierhuizen and Slatyer, 1965; Wallis, 1981; Feddes and Van Wijk, 1977; Federer, 1979; H i l e r , 1969) are very d i f f i c u l t to apply in p r a c t i c a l designs due to dearth of appropriate data for l o c a l conditions. An assumption common to most models i s that maximum crop y i e l d occurs when the plants transpire at the potential rate (Singh, 1981). An obvious flaw to thi s assumption i s that crop loss due to factors other than water d e f i c i t s w i l l also be assigned as crop loss due to water d e f i c i t s . This makes data, where avail a b l e , somewhat questionable for use in p r a c t i c a l designs. 28 2.3 Stochastic Approaches to Drainage Design Stochastic-conceptual models of physical systems are trying to provide answers to questions being raised in the minds of s c i e n t i s t s as to why dispersion in the predicted values occurs when conceptually their analyses seem to be correct. This leads us to believe that there are some hidden but inherent sources of error in our formulations even though an appropriate mathematical model is used. If such a deterministic model i s available, then the future state of the system i s completely predictable provided the i n i t i a l and boundary conditions, the system parameters, and the forcing functions can be sp e c i f i e d exactly. However,in practice these quantities are seldom completely known and are estimated from the available data which i s usually limited. Assuming no other sources of error such as instrumental errors, reading errors and round off errors, the estimates would s t i l l be uncertain to some extent. Therefore, an ensemble of values can be defined, each having a certain p r o b a b i l i t y d i s t r i b u t i o n . A deterministic model is equally applicable to any point in t h i s ensemble and therefore, the solution too w i l l have a corresponding ensemble. In a purely deterministic case, one of the various points in the ensemble is considered as the best to the exclusion of a l l others and analysis i s based on i t (Sagar, 1978a). This choice of a best value relates to the concept of a possible equivalence between the deterministic and stochastic models though i t i s not yet clear i f such equivalence exists 29 (Freeze, 1975; Bakr, 1976; Dagan, 1979). However, most of the deterministic models are based on t h i s concept and they have been quite successful in the past. There may be a problem in characterizing this equivalence (Freeze, 1975; Sagar, 1978a) but i t seems appropriate to assume that such an equivalence does exist and e f f o r t s must be made to quantify i t rather than pursuing t r u l y stochastic models (Smith and Freeze, 1979; Gutjahr and Gelhar, " 1981) which seem to have more academic than p r a c t i c a l values at the present time. Sagar (1978a) has discussed f i v e types of uncertainties in mathematical models according to the form in which random elements enter into the system: (1) random forcing function, (2) random i n i t i a l condition, (3) random boundary conditions, (4) random c o e f f i c i e n t s or parameters and (5) a mixture of the above. The method of solution of a stochastic PDE depends on the form in which the random elements enter into the equation. The finding of a n a l y t i c a l solutions to these equations is s t i l l at an early stage because, although the equation i s linear in a determinstic sense, i t i s nonlinear in stochastic quantities. Therefore, only comparatively simple problems have been solved a n a l y t i c a l l y (Sagar and K i s i e l , 1972; Bakr, 1976; Sagar, 1978a). In drainage/subirrigation design, the fluctuation of the water table is d i r e c t l y linked to the recharge from natural p r e c i p i t a t i o n and to discharge through evapotranspiration. Both of these quantities are highly dependent on weather and therefore, the risk and uncertainty associated with the weather is also associated with the net recharge to the water table. 30 Van Schilfgaarde (1965) developed a relationship between the water table and intermittent r a i n f a l l . He studied the risk of f a i l u r e of a drainage system due to the uncertainty in r a i n f a l l by running a drainage model with long h i s t o r i c a l records. He argued that since s o i l s may not be at f i e l d capacity when a r a i n f a l l event occurs, some portion of the t o t a l r a i n f a l l i s held in the s o i l and i t might not reach the water table. Again, some portions may contribute towards surface runoff and evapotranspiration. Therefore, only a portion of the to t a l r a i n f a l l reaches the water table. He used a water balance approach to get net recharge to the water table. This net recharge replaced the r a i n f a l l term in his model. Extensions have been made to this approach by Young and Ligon (1972), Wiser et a l . (1974), Chieng (1975) and Skaggs (1978). However, such an approach i s performed using h i s t o r i c a l data only and therefore, the randomness associated with the natural processes of p r e c i p i t a t i o n and evapotranspiration i s not considered p r o b a b i l i s t i c a l l y (Sagar and P r e l l e r , 1980). Musy and Duckstein (1976) applied Bayesian s t a t i s t i c a l techniques to solve the problem of designing subsurface drainage systems under hydrologic uncertainties. Their method requires a probability density function of extreme r a i n f a l l and an e x p l i c i t crop loss model both of which are very tedious to estimate for lo c a l conditions. Further they assumed that only one submersion event can occur per year and_that such an event always occurred at the worst possible growth stage. Later, Fogel et a l . (1979) proposed the use of a stochastic model of r a i n f a l l which no 31 longer requires the above assumptions to be incorporated. Again, their methodology presupposes a crop response model which should be s i t e - and crop- s p e c i f i c . Cornell (1972) proposed the use of a f i r s t order uncertainty analysis to model hydrologic uncertainties. This method requires the s p e c i f i c a t i o n of the f i r s t two moments of the random variable only and can y i e l d these two moments for the solution of the problem quite e a s i l y . However, i t gives accurate r e s u l t s only i f the c o e f f i c i e n t of v a r i a t i o n i s small. This method w i l l be discussed in more d e t a i l later when discussing parameter uncertainty where i t w i l l be compared with f u l l d i s t r i b u t i o n analysis methods of incorporating parameter uncertainty in p r a c t i c a l designs. Sagar (1979) showed that for the l i n e a r i z e d Boussinesq equation, the expected value of the hydraulic head i s the same as would have been obtained had the problem been posed as a deterministic one employing expected values of the i n i t i a l condition, boundaries and net recharge. This was to be expected because he used a linear equation. The same conclusion w i l l not hold for a nonlinear equation. He obtained the solution in the form of two expressions, one for the expected value of the head and the other for i t s covariance function. In another study, Sagar and P r e l l e r (1980) solved the case for drainage under uncertain weather conditions. They started with the l i n e a r i z e d form of the Boussinesq equation and assumed di f f e r e n t stochastic functional forms for r a i n f a l l and evapotranspiration. They were able to solve t h i s 32 initial-boundary value problem and obtained the f i r s t two moments for the hydraulic head(mean and the variance). They suggested that standard deviation can be considered as a measure of the uncertainty in the water table l e v e l and therefore, of. the r i s k involved. A drainage/subirrigation design should also include uncertainties due to s p a t i a l and temporal variations of s o i l parameters. As stated e a r l i e r , equivalent s o i l properties and their incorporation in p r a c t i c a l designs should be investigated. This holds a lot of promise in such designs as compared to the research e f f o r t s done to quantify small scale v a r i a t i o n s . The writer is not aware of any study in s o i l and water engineering where temporal variations in s o i l properties are considered. This concept can be very important in drainage/subirr-igation designs because from an observation point of view, i t i s known that hydraulic conductivity of the top s o i l increases after drainage because of the improvement in s o i l structure. Also, due to increased worm a c t i v i t i e s in the top one metre of the s o i l in drained f i e l d s the hydraulic conductivity of the sub-soil may also increase. This observation can not be considered o r d i n a r i l y in such designs because the system has to perform i t s function in the f i r s t few years after i n s t a l l a t i o n . This concept may be considered in designs by imposing a mole drainage system, where possible, over a t i l e drainage system. Although the l i f e of mole drainage systems is not f u l l y known at t h i s time, t h i s system should be ef f e c t i v e for at least a few years from i t s i n s t a l l a t i o n . By 33 that time, the hydraulic conductivity may have increased and the drains which were i n i t i a l l y i n s t a l l e d at a spacing wider than the normal w i l l not become an under-designed system as may have been perceived e a r l i e r . However, such a drainage system design is only possible i f some qua n t i f i c a t i o n with regard to the percentage increase in hydraulic conductivity as well as the l i f e of a mole drainage system can be made. A ra t i o n a l drainage/subirrigation design procedure which includes the effects of uncertainty due to climate and s o i l parameters has not been developed to date. Though as reported e a r l i e r , d i f f e r e n t techniques ranging from purely a n a l y t i c a l (Bakr, 1976) to f u l l y numerical (Freeze, 1975) have been employed in analyzing t h i s uncertainty in initial-boundary value problems of groundwater flow. The various approaches dealing with parameter uncertainty in groundwater flow can be divided into two main groups: f u l l d i s t r i b u t i o n analyses, and f i r s t and second moment methods. Dettinger and Wilson (1981) provide a c r i t i c a l appraisal of these techniques, underscoring their advantages and disadvantages. The discussion given below summarizes their findings. F u l l d i s t r i b u t i o n methods require a complete s p e c i f i c a t i o n of the p r o b a b i l i s t i c properties of a l l stochastic inputs and parameters of the flow system and y i e l d (or attempt to yield) a complete pro b a b i l i t y d i s t r i b u t i o n of the resulting flow. F i r s t and second moment methods assume that the f i r s t two moments of a random variable are s u f f i c i e n t to characterize i t . They 34 consider only the mean and variance-covariance of the dependent as well as independent variables. The method of derived d i s t r i b u t i o n s and Monte Carlo simulations are the two most important f u l l d i s t r i b u t i o n techniques. The derived d i s t r i b u t i o n approach i s an a n a l y t i c a l method of deriving the p r o b a b i l i t y d i s t r i b u t i o n of a random function given the d i s t r i b u t i o n s of i t s independent variables (Benjamin and Cornell, 1970). The analysis becomes very complicated unless applied to simple systems with r e l a t i v e l y simple p r o b a b i l i s t i c properties. An example of this type of analysis i s Sagar and K i s i e l ' s (1970) examination of parameter uncertainty for aquifer pump tests. The Monte Carlo method employs numerous rep l i c a t i o n s of flow system simulations, with the parameters and inputs of each simulation generated at random from their respective p r o b a b i l i t y d i s t r i b u t i o n s . This method can be readily automated and i s applicable to r e l a t i v e l y complex problems. The disadvantages of t h i s method are high computational costs and the fact that results are never in the closed form which a derived d i s t r i b u t i o n s t r i v e s for, and therefore results are not readily transferable to a new s i t u a t i o n . Monte Carlo methods have been applied to study the ef f e c t of s p a t i a l v a r i a b i l i t y of physical properties in flow through porous media including Warren and Price (1961), Freeze (1975) and Smith and Freeze (1979). It i s often quite d i f f i c u l t to obtain f u l l p r o b a b i l i t y d i s t r i b u t i o n s that are input to a f u l l d i s t r i b u t i o n method and i t i s s l i g h t l y easier to estimate the moments. Because the 35 results obtained using f u l l d i s t r i b u t i o n methods are highly dependent on the d i s t r i b u t i o n selected for the input parameters they can be quite deceptive i f the assumed d i s t r i b u t i o n happens to be wrong. It may be worth noting that the results obtained by using f u l l d i s t r i b u t i o n analyses are generally shown in terms of the f i r s t two moments of the dependent variable (Sagar and K i s i e l , 1972; Freeze, 1975; Eagleson, 1978). The f i r s t and second moment method assumes that the information about random variables (or functions) can be characterized by the mean representing the central or expected tendency of the variable (or function) and the variance-covariance representing the amount of scatter around the mean. A function completely defined by these moments i s the one which i s normally d i s t r i b u t e d . It has zero skewness and other higher moments can be calculated from the variance (Benjamin and Cornell, 1970). The disadvantages of th i s method are that i t i s impossible to test the assumption that the mean and variance-covariance f u l l y describe the function and that certain relationships of interest (e.g., Y=max[x]) do not lend themselves to this analysis (e.g., are not d i f f e r e n t i a b l e ) . In many problems, these disadvantages are more than offset by several major advantages. The f i r s t and second moment method involves the same familiar tools and procedures of algebra and calculus that are commonly used in the more deterministic analysis of. the same problem (e. g. f i n i t e difference methods, matrix algebra e t c . ) . The p r a c t i c a l f e a s i b i l i t y of analyzing more thoroughly stochastic models encourages such modeling to 36 take place where i t might not have otherwise been taken. It i s better to have an approximate model of the whole problem than an exact model of only a portion of i t . In addition to the above advantages, an approach based on means and variance-covariances may be a l l that i s j u s t i f i e d when one notes (Cornell, 1972): (1) that the data and physical arguments are often i n s u f f i c i e n t to e s t a b l i s h the f u l l d i s t r i b u t i o n ; (2) that most engineering analyses include an important component of the real (but d i f f i c u l t to measure), professional uncertainty(due, for example, to imperfect physical theories and to engineering approximations); and (3) that the f i n a l output (decision or design parameter) i s often not sensitive to moments higher than the mean and var iance-covar iance. F i r s t and second moment methods can be applied in two ways: perturbation and/or Taylor series expansions. In perturbation analyses, the PDE governing the piezometric heads i s perturbed s l i g h t l y , y i e l d i n g a new equation that governs the random component of the piezometric surface. Bakr et a l . (1978) used t h i s technique in a wave number domain to solve the problem of s p a t i a l v a r i a b i l i t y . They represented the p r o b a b i l i s t i c moments of the logarithm of hydraulic conductivity by F o u r i e r - S t i e l t j e s integrals and produced second moment properties of heads in one-, two- and three-dimensional flows. • These results s t a t i s t i c a l l y describe flow in aquifers whose properties are continuously variable, rather than aquifers with properties that vary in discrete blocks of aquifers described by lumped 37 parameters. This continuous description is c e r t a i n l y the most r e a l i s t i c . Also these wave number domain analyses are quite rigourous mathematically. Further, they require estimates of input covariances which are very d i f f i c u l t to estimate even i f extensive data are available (Gutjhar and Gelhar,1981). Sometimes the terms f i r s t order and second order analysis are used for the f i r s t and second moment methods. F i r s t order analysis i s the analysis of mean and variance-covariance of a random function based on i t s f i r s t order Taylor series expansion. Second order analysis i s the analysis of the mean based on a second order Taylor series expansion (the variance-covariance analysis is s t i l l r e s t r i c t e d to the use of f i r s t order Taylor series expansion). Therefore, the mean estimated from f i r s t and second order analysis may be d i f f e r e n t but the variance-covariance w i l l not. F i r s t and second moment methods based on Taylor series have been employed by a number of authors. Cornell (1972) presented applications of thi s approach to a wide variety of simple hydrologic and water resources problems and suggested much wider applications. Sagar (1978b) used th i s approach in analyzing one-dimensional flow in confined aquifers by the f i n i t e element method. In view of the above facts, the writer has employed f i r s t and second moment methods to study effects of parameter uncertainty in drainage/subirrigation design. This approach i s easy to apply and i t can y i e l d valuable answers, especially in those situations where other methods of incorporating 38 uncertainty may be impractical. The uncertainty due to weather has been incorporated by using a Markov chain model that estimates the probability of success of a drainage/subirrigation design. 39 CHAPTER III CLIMATE AND DRAINAGE/SUBIRRIGATION DESIGN In previous chapters the d e s i r a b i l i t y of a more r e a l i s t i c approach to drainage/subirrigation design than the steady state design was discussed. However, d i f f i c u l t i e s inherent in a more comprehensive analysis were also noted. Also d i f f i c u l t i e s in defining an objective function for a drainage system were discussed. In t h i s chapter i t is argued and supporting evidence i s presented that for a drainage/subirrigation system design for the Lower Fraser Valley, drainage requirements govern the system design. It i s suggested that drainage requirements during the c u l t i v a t i o n and planting period in early spring override a l l other seasonal requirements. It i s suggested that a suitable index for the performance of a drainage system i s the pr o b a b i l i t y of providing a s a t i s f a c t o r y period with workable s o i l conditions in early spring. A d e f i n i t i o n of satisfactory workable s o i l conditions in quantifiable p r a c t i c a l terms i s presented and a model which has been developed to estimate the p r o b a b i l i t y of achieving such conditions i s described. The model includes the use of the Boussinesq equation to simulate nonsteady state saturated flow towards the drains, a Markov chain model to handle v a r i a b i l i t y in weather conditions and a computer simulation model to provide estimates of the p r o b a b i l i t i e s of obtaining a workable period in early spring. An experimental study to check on the Boussinesq equation i s also described. 40 3.1 Drainage/Subirrigation System Design A drainage/subirrigation system design involves the design of a system which can be used for both drainage and subi r r i g a t i o n as and when required by climatic conditions. This is a r e l a t i v e l y new practice but holds promise because the high i n i t i a l cost of a drainage system can be j u s t i f i e d more eas i l y i f i t can also be used for i r r i g a t i o n . More farmers can be persuaded to adopt such systems because of the increased returns from investment. It i s d i f f i c u l t to design such systems. Although drainage system design is quite tractable, the subir r i g a t i o n system design is s t i l l in i t s infancy. Not much is known at thi s time with regard to the hydraulic e f f i c i e n c y of such systems. Moreover, unsaturated flow plays a more s i g n i f i c a n t role in subi r r i g a t i o n than in drainage. It i s the unsaturated region in the s o i l p r o f i l e through which water moves upwards to meet the demands of the plant root system. But as discussed e a r l i e r , the incorporation of the unsaturated flow in a p r a c t i c a l design procedure i s quite impractical at thi s time. The complexities increase exponentially and data l i m i t a t i o n s are, for a l l p r a c t i c a l purposes, insurmountable in such an approach. Although not much is known about the performance of a drainage/subirrigation system in subir r i g a t i o n mode, i t has been held that in the Lower Fraser Valley, subirrigation system design w i l l be a c o n t r o l l i n g factor in the design of a drainage/subirrigation system. These assumptions seem to be based more on i n t u i t i o n than actual observations in the f i e l d . 41 Drainage systems are designed in the valley on the basis of steady state theory (Hooghoudt's method) to meet a drainage c o e f f i c i e n t of 1.27 cm per day for most locations (Anon., 1972). A subi r r i g a t i o n system for the valley should be able to s a t i s f y a s u b i r r i g a t i o n c o e f f i c i e n t of 5 mm per day which i s the maximum potential evapotranspiration rate expected to occur over a period of several days (Black and Goldstein, 1977 as c i t e d by Oke, 1978). They also observed that for a clay loam s o i l , about 70% of the crop water requirements were met by the drainage system when used for subi r r i g a t i o n by routine r a i s i n g of water lev e l in the outlet. However, they did not measure the unsaturated hydraulic conductivity of the s o i l as a function of pressure head that would have enabled one to check the s u i t a b i l i t y of using a drainage system for su b i r r i g a t i o n when the former i s designed in the conventional way. They suggested that for l o c a l conditions a subirrigation system should be able to control the water table in such a way that a steady upward flux of c a p i l l a r y water from the water table can be maintained at a rate of 5 mm per day to the bottom of the active root zone of the plants (a starting value of 30 cm may be adequate for most crops in the v a l l e y ) . Drainage requirements thus appear to be more than two and a half times the subirrigation requirements. A common practice for designing subirrigation systems i s to base them on steady state theory (Skaggs, 1978). The design requirements for such a system may be expressed as the a v a i l a b i l i t y of water at the bottom of an active root zone from 42 a water table at some distance below i t . Water i s made available to the plant roots by c a p i l l a r y flow and therefore the unsaturated flow properties of the s o i l play a dominant role in subirr i g a t i o n system design. However, the unsaturated flow properties are seldom known for l o c a l s o i l s in many locations. The Lower Fraser Valley i s no exception. The water retention c h a r a c t e r i s t i c s for some lo c a l s o i l s have been measured recently (Anon., 1980). However, without knowing the unsaturated hydraulic conductivity of s o i l s at some low tension, the above data can not provide estimates of the required unsaturated hydraulic conductivity - pressure head relationships. Therefore, to test the hypothesis similar s o i l s were selected from the l i t e r a t u r e whose unsaturated flow c h a r a c t e r i s t i c s are well documented. The analyses on drainage/subirrigation system design in this section were performed assuming these s o i l s were located in the Valley. The model validation for subir r i g a t i o n was not possible though the model has been found to be quite accurate by various researchers as i t is based on Darcy's law (Gardner, 1958; Gardner and Fireman, 1958). Cap i l l a r y r i s e from the water table is governed by Darcy's law for unsaturated flow, which i s v = -K(p) dh/dz (1 ) Where v i s the volumetric flux per unit area perpendicular to the flow, K(p) i s the unsaturated hydraulic conductivity, h i s 43 the h y d r a u l i c head being the sum of pressure head (p) and the g r a v i t a t i o n a l head ( z ) . Pressure head i s negative i n the unsaturated zone. The g r a v i t a t i o n a l head i s zero at the water t a b l e and i s assumed p o s i t i v e i n upward d i r e c t i o n . R e w r i t i n g Darcy's law as v = K(p) { (dp/dz) - 1} (2) For steady s t a t e flow i n the v e r t i c a l d i r e c t i o n , v reamins c o n s t a n t . S o l v i n g f o r Z (p=o at the water t a b l e ) , we have Z = J { dp/(l+v/K(p)) } (3) Equation (3) was suggested by Gardner (1958) f o r steady e v a p o r a t i o n from a s t a t i o n a r y water t a b l e . I n t e g r a t i o n of equation (3) depends mainly upon the f u n c t i o n a l r e l a t i o n s h i p between K and p. Closed-form s o l u t i o n s to equation (3) are a v a i l a b l e f o r some f i x e d f u n c t i o n a l forms of K(p) r e l a t i o n s h i p s (Gardner, 1958; Gardner and Fireman, 1958; C i s l e r , 1969). However, numerical s o l u t i o n s to equation (3) can be obtained without much e f f o r t (Anat et a l . , 1965). In numerical s o l u t i o n s , the K(p) r e l a t i o n s h i p s can e i t h e r be .provided i n f u n c t i o n a l form or i n simple t a b u l a r form. The advantages of the numerical s o l u t i o n s are obvious: there i s no assumption with regard to the f u n c t i o n a l form of K(p) r e l a t i o n s h i p . Whatever i s observed i s used i n the a n a l y s i s . S o l u t i o n s to equation (3) can y i e l d the maximum a l l o w a b l e upward f l u x (vmax) that i s p o s s i b l e from a given s o i l f o r a give n depth. A l s o , the above a n a l y s i s can be made to y i e l d that 44 depth of the water table which, for a given s o i l , can provide a given value of flux at the bottom of the active root zone without l e t t i n g the tension in the root zone exceed some threshold value. This threshold value i s a function of the crop species and may vary with the stage of growth of the crop. Taylor and Ashcroft (1972) provide a table of such values for di f f e r e n t crops. The maximum allowable tension at the bottom of the active root zone plays a very s i g n i f i c a n t role in subi r r i g a t i o n design. A somewhat d i f f e r e n t approach to determine the upward flux can be formulated by using the concept of matric flux potential introduced by Gardner (1958) and further discussed by Raats (1971) and Shaykewich and Stroosnijder (1977). The matrix flux.potentia1 i s defined as MFLP(p) = ( K(p) dp (4) The integration is performed between some lower l i m i t of pressure head |^ at which moisture content i s © 0 and the desired pressure head p at which the water content i s Q. The MFLP combines both parameters c o n t r o l l i n g water flow (hydraulic conductivity and hydraulic head) into a single parameter. Writing Darcy's law as equation (2) v = K(p) ( (dp/dz) - 1} 45 For steady state flows v = K(p) (dp/dz) - K(p) P i v + K(p) = (1/dz) J K(p) dp = (1/dz) {MFLP(p2)-MFLP(p1 ) } (5) where dz i s small. The extreme right-hand side of equation (5) i s an exact representation of the i n t e g r a l . Therefore, the procedure to get a given vmax at the bottom of the root zone f i r s t involves the cal c u l a t i o n of MFLP at d i f f e r e n t suction heads. Then to get a given vmax for a l i m i t i n g pmax, one has to substitute for MFLP(pmax) and MFLP(O) and solve equation (5). It was shown by Shaykewich and Stroosnijder (1977) that the MFLP method is a more exact method of simulating water movement than using the Gardner approach in which . hydraulic conductivity and pressure head are treated separately. However, when extra care was taken to make depth increments very small in the case of equation (3), no noticeable difference was found between the two approaches. Therefore, in the present study, equation (3) was used to give the depth to the water table for a l i m i t i n g pressure head at the bottom of the e f f e c t i v e root zone for a given upward flux. Depth increments of 1 cm were used in a l l cases. The s o i l s used in the study were Ha l l a r t s i l t y clay loam (Paul and deVries, 1979); Pachappa fine sandy loam, Yolo l i g h t clay, Chino clay (Gardener and Fireman, 1958); a s i l t y clay loam (Bouma et 46 a l . , 1972) and a hypothetical medium textured s o i l used by Pikul et a l . (1974) as suggested by Freeze. The results of the drainage and subirrigation designs, had these s o i l s been located in the Lower Fraser Valley, are shown in Table 1. The maximum depth to the water table was calculated for the condition that the steady upward flux i s 5 mm per day and pmax at the bottom of the active root zone (30 cm) i s less than or equal to 400 cm of tension. Drain spacings required for drainage and subirrigation systems using steady state theory (Skaggs, 1979) are shown in Table 1. In most cases, spacings for drainage systems are smaller than those for subir r i g a t i o n systems. Goldstein (1978) also found the above to be true for two s o i l s selected from l i t e r a t u r e . Using drainage systems for sub i r r i g a t i o n , the height to which the water table in the drainage outlet should be raised was found to be very p r a c t i c a l in most cases. These heights are also shown in Table 1 for d i f f e r e n t s o i l s . Unless the s o i l texture is very fine, the head required to move water back into the s o i l p r o f i l e e f f i c i e n t l y was smaller than was perceived. In almost a l l cases, i t was found that for l o c a l conditions a suitable drainage system would be more than adequate for subirrigation purposes. Although the s o i l s used in this analysis were not l o c a l s o i l s but similar s o i l s selected from l i t e r a t u r e , they are believed to be reasonably representative of l o c a l s o i l s . Thus i t seems reasonable to conclude that for the Lower Fraser Valley drainage requirements w i l l govern in the design of a dual purpose 47 TABLE 1 COMPARISON OF DRAINAGE AND SUBIRRIGATION SYSTEM DESIGNS FOR THE LOWER FRASER VALLEY S o i l Type Dralnage (m) S u b ! r r i g a t ion (m) Head f o r Subi r r i gat ion (m) H a l l a r t Si 1ty c l a y loam 53-0 68.00 0.39 Pachhapa f i n e sandy loam 10.0 12.00 0.-»3 Y o l o 1i ght c l a y 3.0 2.09 0.70 Chino c l a y 3.8 3.92 0.68 S i l t y c l a y loam 32.0 32 .00 0.70 S o i l ( f r e e z e ) 32.0 32.00 0.39 NOTE: d r a i n a g e c o e f f i c i e n t - 1.2 cm/day s u b i r r i g a t i o n c o e f f i c i e n t - 0.5 cm/day 48 drainage/subirrigation system. However, some further studies using actual data on l o c a l s o i l s would be desirable to provide more confidence in the conclusion which i s so important to drainage/subirrigation system designs for the Valley. From th i s point onwards, t h i s study concentrates only on drainage design. As discussed above, s u b i r r i g a t i o n design i s i m p l i c i t in such analysis. However, i t i s reminded that topography of the f i e l d has an important bearing on the subirrigation design. The conclusion reached in this section should hold only for the r e l a t i v e l y f l a t lands in the Valley. 3.2 Drainage Design C r i t e r i a for the Lower Fraser Valley Drainage systems are designed in the valley mostly on the basis of steady state theory (Hooghoudt's ' method). The major assumption in this theory i s that inflows equal outflows in a given time (zero storage). Therefore, the water table i s assumed to be stationary on a long-term basis and i s in balance with recharge. The design c r i t e r i o n for a steady state condition can be described by a certain combination of drainage c o e f f i c i e n t and a minimum permissible depth to the water table midway between drains. The drainage c o e f f i c i e n t i s the rate of water removal by a drainage system in one day. Generally, a value of 1.30 cm/day i s recommended for the Lower Fraser Valley (Anon., 1972). The minimum permissible midspan water table depth i s a function of both workability requirements of the s o i l s and the crops to be grown in the area. This parameter i s the least understood 49 factor in the design for l o c a l conditions and generally some rules of thumb are used in characterizing i t for actual designs. Steady state conditions seldom, i f ever, p r e v a i l in the Valley. Although we do get low-intensity, long duration p r e c i p i t a t i o n storms in the Valley, the real s i t u a t i o n in the f i e l d i s a fluctuating water table mostly due to variable intensity r a i n f a l l and varying lengths of time between rainstorms. The steady state design method i s used for convenience and in the absence of any other suitable p r a c t i c a l method. In general, designs based upon the steady state theory seem to have performed s a t i s f a c t o r i l y from a drainage point of view. But i t must be borne in mind that the assumption of steady state i s not physically correct for l o c a l conditions. Therefore, the p o s s i b i l i t y of using a more r e a l i s t i c nonsteady state design method was investigated. In a nonsteady state design, fluctuations of the water table are considered e x p l i c i t l y . The design would normally be based upon the prevention of harmful sequences of water table fluctuations. What these undesirable patterns are i s s t i l l not known. There have been some studies done elsewhere which report on some possible c r i t e r i a that could be used in drainage design (Bouwer, 1974). But most of them are either at a rudimentary stage or are based on l o c a l experience. Therefore, their transposition to other than the l o c a l conditions for which they were developed is somewhat unclear. 50 For the Valley, the drainage requirements can be divided into three groups., depending upon the d i f f e r e n t seasons in a year namely, early spring requirements, growing period requirements and early f a l l requirements. The winter requirements have not been l i s t e d because there are two schools of thought as to whether a zone of unsaturated flow should be maintained at a l l times in the winter (De Vries, 1981) or not (Baehr,1980). Those in favour of maintaining such a zone argue that aerobic conditions in winter w i l l help maintain or improve the s o i l structure due to worm a c t i v i t i e s . Those who see no value in maintaining such a zone argue that occasional flooding in the winter may k i l l some harmful nematodes in the s o i l which would otherwise be harmful to crops. Although the writer i s in favour of maintaining t h i s unsaturated zone on a continuing basis, t h i s question is not addressed in the present work. 3.2.1 Early spring requirements The annual water balance of the Lower Fraser Valley i s shown in Figure 2. Winter i s usually very wet in the Valley with the most of the annual p r e c i p i t a t i o n f a l l i n g in t h i s period. Normally because of the wet winter, s o i l s are quite wet in early spring. To a large degree an a g r i c u l t u r a l enterprise is at the mercy of the weather. The more a cropping season gets delayed, the greater the risk involved. An important a g r i c u l t u r a l function of a drainage system i s to provide workable conditions on the farm to meet crop planting needs in the early spring. Deficient workability in spring can 51 Precipitation Potential evapotranspiration Change in soil storage Actual evapotranspiration Surplus Deficit o E \ E E -40 T r T 1 1 1 1 1 r J i i L A-zo J F M A M J J A S O N D r « -H — Drainage Irrigation Drainage o E Figure 1 . Annual water balance of the Lower Fraser Valley from Hare and Thomas (1974). 52 cause delays in planting and hence can result in either a poor quality crop or reduction in crop yields (Reeve and Fausey, 1974). Farming operation schedules may get disrupted due to poor workability conditions .thereby causing a jump in production costs. For some crops l i k e potatoes, the losses due to the ov e r a l l delay may be very s i g n i f i c a n t because early potatoes can fetch higher prices than late ones. If planting i s delayed, the other follow-up operations are also delayed and these repurcussions may sometimes prove f a t a l to crops. Paul and deVries (1983) have developed procedures for the prediction of t r a f f i c a b i l i t y of t i l e - d r a i n e d lowland s o i l s in the Lower Fraser Valley. They report c r i t i c a l water table depths for t r a f f i c a b i l i t y as 53, 45, and 60 cm for Lumbum muck, Hall a r t s i l t y clay loam (grassland) and Ha l l a r t s i l t y clay loam (cultivated land) respectively. Wind (1976) used an average tension of 300 cm in the top 5 cm as the l i m i t i n g value for workability. As discussed in chapter II, workability and t r a f f i c a b i l i t y do not imply the same thing. However, both of them do depend upon water tension of the top layer of s o i l which i s highly dependent upon the weather. Because of the p r a c t i c a l d i f f i c u l t i e s of applying unsaturated flow analysis in actual designs, as discussed in chapter II, i t was decided to employ the reported values of Paul and deVries as indicators of workability in the Valley. A 60 cm depth to the water table on a p a r t i c u l a r day was assumed to ensure workable s o i l conditions on that day. 53 3.2.2 Growing period requirements The function of a .drainage system does not stop at providing workable conditions in early spring. These conditions should p r e v a i l during the whole crop growing period to enable one to perform other farming operations such as weed and pest control as and when the need a r i s e s . Another important function of drainage in the growing period is the prevention of harmful fluctuations of the water table in the active root zone. Wet conditions during th i s period can prove completely f a t a l i f excess water i s allowed to stay in t h i s zone for long durations or to reduced yie l d s i f allowed to stay for shorter periods of time. The exact form of the crop loss functions for many crops of interest is seldom known in most locations. The nature of this function i s very complex and poorly understood. However, i t i s known that crop loss i s highly dependent upon the magnitude of the excess water stress and the timing of such a stress in the growing period. At the present time, these functions are not known for any crop in the Lower Fraser Valley. 3.2.3 Early f a l l requirements The a g r i c u l t u r a l function of a drainage system does not end in ensuring workable conditions and a good plant environment in spring and during the growing season. T r a f f i c a b i l i t y during harvesting is also an important c r i t e r i o n . Not much w i l l be gained i f conditions were good throughout planting and growing period but wet in the f a l l . If harvesting i s not done on time, 54 there may be a poor quality crop i f i t i s somewhat delayed or an almost t o t a l loss i f i t i s very delayed. A good drainage system design may be looked upon as a way of reducing the risk in farming enterprise. An e f f i c i e n t drainage system w i l l increase the p r o b a b i l i t y that t r a f f i c a b l e conditions w i l l p r e v a i l in f a l l . The results obtained by Paul and deVries (1983) can be used to check t r a f f i c a b l e conditions in f a l l for the l o c a l s o i l s . 3.2.4 Significance of early spring design c r i t e r i o n The climate in the Lower Fraser Valley i s such that the p r e c i p i t a t i o n i s greater in spring than in the other seasons of interest from a drainage system design point of view (Baier et a l . , 1969). This implies that the drainage requirements are more stringent in spring than in any other season. Therefore, i f a drainage system is designed to provide satisfactory workable conditions in early spring, i t should be more than s u f f i c i e n t to meet drainage requirements for summer and early f a l l . This leads to great s i m p l i f i c a t i o n in drainage design. The most complex and the least understood c r i t e r i o n , that i s , the crop loss function, does not need to be e x p l i c i t l y considered in such a design. Although there may be some problems in doing a benefit/cost analysis in drainage design with t h i s approach, such a design w i l l help i d e n t i f y the best design for a given location from a designer's point of view. A benefit/cost analysis can be done for both early spring and early f a l l periods i f benefits of providing workable and 55 t r a f f i c a b l e days can be expressed in terms of some monetary value. We are now at a point where we must define satisfactory workability requirements, that i s , what i s a workable day, when and how many workable days would be enough to perform necessary operations. For example seed bed preparations with usual tools (tractor-mounted ploughs, c u l t i v a t o r s ) require s o i l strength to be s u f f i c i e n t to support farm t r a f f i c without experiencing any appreciable damage to i t s structure (a function of the tension in the plough layer) and when i t w i l l f a l l into small crumbs, not into clods and d e f i n i t e l y not be smeared and puddled. Skaggs (1978) has defined a workable day as a day when the a i r volume exceeds a l i m i t i n g value; the r a i n f a l l occuring on that day i s less than a maximum value and; that a minimum number of days have elapsed since that amount of r a i n f a l l occurred. This d e f i n i t i o n requires, among other things, some sort of qua n t i f i c a t i o n with regard to the status of a i r in the s o i l p r o f i l e . This is possible only in an unsaturated flow analysis and once again i t must be ruled out from such a study because of p r a c t i c a l l i m i t a t i o n s as discussed in section 2.1.3. The findings of Paul and deVries (1983) with regard to depth of the water table below the s o i l surface offer a possible c r i t e r i o n for designing drainage systems in the Valley. A workable day i s the one on which water table stays below 60 cm from the s o i l surface. It should be noted that the above c r i t e r i o n was arrived after a detailed analysis that involved both unsaturated flow analysis and s o i l strength - pressure head 5 6 analysis. The next hurdle i s the consecutive number of workable days required to perform farming operations in the spring. The word consecutive i s important because not much w i l l be gained i f the required number of days are spread a l l over the month of March. This w i l l result in increased production costs and may, in fact, be equivalent to having no workable period at a l l . The consecutive number of days w i l l depend, among other factors, on the crop to be grown. Ten to twelve days appear to be s u f f i c i e n t for most crops in the Valley (Driehuyzen, 1982; deVries, 1982; Skaggs, 1978). In t h i s study, a twelve-day workable period was selected as the design c r i t e r i o n . This includes a small factor of safety since the depth to the water table i s only an indi r e c t indicator of workability. However, two or more periods of six days are also very good for seedbed preparation and planting. The methodology can be e a s i l y modified to estimate the p r o b a b i l i t y of such events. A methodology i s developed to estimate the p r o b a b i l i t y of any number of consecutive workable days. A period of twelve days i s used in the present work for i l l u s t r a t i o n purposes. Any other number could be used without a f f e c t i n g the methodology. 3.3 The Mathematical Drainage Model A schematic of the drainage/subirrigation system considered herein was shown in Figure 2. Under the D-F assumptions and no storage in the unsaturated zone, the Boussinesq equation may be written as f »h . K. 3 (h j h ) (6) dt dx dx 57 Precipitation or evapotranspiration 1 1 1 t t t t Figure 2. Schematic of a drainage/subirrigation system. 5 8 If a source/sink i s present, then equation(6) may be written as f 3h » K _JL_ (h 3h) * R (7) 3t 3x 3x " where h i s the distance of the water table above the impermeable layer, t i s time, K i s the saturated "hydraulic conductivity of the s o i l , f i s the drainable porosity, X i s the horizontal distance and R i s the rate at which water i s added or subtracted to the water table by r a i n f a l l or evapotranspiration respectively. Equations(6) or (7) can be solved numerically for drainage or s u b i r r i g a t i o n depending upon pertinent i n i t i a l and boundary conditions. For drainage, these are u A ^ a 3x 2 . 4x s 2x*\ . h « d + 8 y Q ( § - 5 7 + -gr ) . t > o ; o < x < S (7a) h = d ; t > 0 ; x=0 and x=S (7b) 3h = 0 ; t > 0 ; x=S/2 (7c) 3x 59 For s u b i r r i g a t i o n , these are ' d + 8 V s " s2" s 1 - r - ' t = o ; o < x < S/2. (7d) h = h ; t > o ; x=o and x=S o (7e) where i n i t i a l conditions (7a) and (7b) imply a fourth degree polynomial shape for the water table (Dumm, 1964). There are variety of ways in which equations (6) or (7) can be solved on the computer using numerical methods such as f u l l y i m p l i c i t method, Crank-Nicholson i m p l i c i t scheme and predictor-corrector method. A l l of these methods have merits and demerits (Remson et a l . , 1971). The alternate predictor-corrector method given by Douglas and Jones (1963) was used in the present study to solve equation (7) subject to appropriate i n i t i a l and boundary conditions for drainage and s u b i r r i g a t i o n . This method i s the commonly used predictor corrector method with an i m p l i c i t predictor (Remson et a l . , .1971). The use of th i s method results in a linear tridiagonal system of equations at each step which can be solved very e f f i c i e n t l y on the computer. Further, there i s no need to ite r a t e to get solutions at the next time step. This i s in contrast with the f u l l y i m p l i c i t scheme (Skaggs, 1975) where 60 convergence problems can be encountered for some cases. 3.4 Layout and Design of the Experiment Confidence in the prediction by mathematical models can be gained only by comparing how well their predictions match r e a l i t y . Therefore, a f i e l d experiment was planned to validate the mathematical model proposed in the last section. A 5 ha f i e l d experimental s i t e was selected in the Lower Fraser Valley in the d i s t r i c t of Delta. Details of the experimental layout are i l l u s t r a t e d in Figure 3. The s o i l in this area i s Ladner s i l t y clay loam. The f i e l d was not previously under any a g r i c u l t u r a l crop. Although the land was r e l a t i v e l y f l a t , i t was not l e v e l l e d before i n s t a l l i n g drains due to lim i t e d resources. Drain l a t e r a l s of 10 cm diameter were spaced at 14 m adjacent to an area with no subsurface drainage. The average depth of l a t e r a l s was 1.2 m. The s o i l p r o f i l e was quite homogeneous down to a depth of - about 2.5 m, mainly consisting of clay. The s o i l p r o f i l e description i s given in Appendix A. The p l a s t i c l a t e r a l s drained into a c o l l e c t o r ditch which consisted of three separate surface ditches connected with one another by p l a s t i c pipes. The above scheme was devised to have better control over the water table in various parts of the f i e l d . The water table can be raised, for example, in the southern part of the f i e l d by preventing the water flowing from the southern c o l l e c t o r sub-ditch to other sub-ditches. Similar controls can be exercised to raise the water table in the central or northern 19 62 part of the f i e l d . The northern c o l l e c t o r sub-ditch i s connected to the main ditch going along the north-western boundary of the drained area by means of a p l a s t i c pipe. The main di t c h meets the highway di t c h along 72nd street and c a r r i e s water away from the area. Sampling locations to measure hydraulic conductivity by auger hole test are shown in Figure 3. Two batteries of piezometers are also shown that were used to check the performance of automatic water table recorders. Three automatic water table recorders were i n s t a l l e d in the drained area and one was i n s t a l l e d in the adjacent undrained area. Their approximate location in the plot are shown in Figure 3. 3.5 Measurement of Model Input Parameters The input data requirement for the mathematical model used in the study is more or less the same as i s required for designs that are based upon steady state theory. However, being a transient model, a time dependent source or sink is required in the simulations. Further, since fluctuating water tables are considered in such designs, a knowledge of drainable pore space of the s o i l i s also required. The rest of the input data i s the same as required for steady state theory. It includes information on the hydraulic conducting capacity of the s o i l (saturated hydraulic conductivity), the resistance to flow towards drains (depth to the impermeable layer below drains, e f f e c t i v e drain diameter), and the hydraulic forcing functions ( p r e c i p i t a t i o n , evapotranspiration). During the v a l i d a t i o n 63 period, the water le v e l in the main ditch was lower than the drain l e v e l . Therefore, no seepage problems were encountered. 3.5.1 Pr e c i p i t a t i o n P r e c i p i t a t i o n was measured using a standard tipping bucket gauge in combination with an event recorder. The gauge was cali b r a t e d f i r s t in the laboratory and i t was found that when 15.1 cm3 of water i s funneled through a 20.0 cm diameter receptor, the tipping action takes place. This i s equivalent to saying that the bucket t i p s when i t rains 0.481 mm. The event recorder was p a r t i c u l a r l y sensitive to high humidity and low tempratures. Therefore, the recorder and battery were housed in a wooden box. The wooden box was insulated and s i l i c a gel c r y s t a l s t i e d in cheese cloth were enclosed to control temperature and humidity respectively. 3.5.2 Drainable porosity A l l nonsteady saturated flow methods need some qua n t i f i c a t i o n as to the volume of water stored in or drained from the s o i l column when the water table fluctuates. The drainable porosity (f) can be defined as the volume of water per unit area released from storage when the water table f a l l s by a unit distance in the absence of a natural source or sink. It can be defined for r i s i n g water tables also (subirrigation) as the volume of water per unit area taken up by the unsaturated s o i l above the water table for a unit r i s e in the water table in the absence of a natural source or sink. 64 These d e f i n i t i o n s l e a d t o the n o t i o n t h a t f i s a c o n s t a n t ( t h a t i s , independent of p r e s s u r e head). But i t i s known t h a t s o i l s w i t h deep water t a b l e s have a l a r g e r c a p a c i t y t o s t o r e incoming r a i n water than i f the water t a b l e s were s h a l l o w . I t means t h a t a g i v e n r a i n f a l l w i l l cause a v a r i a b l e r i s e i n the water t a b l e , depending upon the l o c a t i o n of the water t a b l e p r i o r t o the r a i n f a l l e v e n t. In d r a i n a g e / s u b i r r i g a t i o n d e s i g n s , a l t h o u g h i t i s a common p r a c t i c e t o t r e a t f as a c o n s t a n t , C h i l d s (1960a) and T a y l o r (1960) have shown t h a t f i s a f u n c t i o n of the water t a b l e h e i g h t . The assumptions of a c o n s t a n t f may not be c o r r e c t e s p e c i a l l y i n s h a l l o w water t a b l e s f o r f i n e t e x t u r e d s o i l s (Duke, 1972). Duke p r e s e n t e d a c l o s e d form e x p r e s s i o n f o r s p e c i f i c y i e l d ( d r a i n a b l e p o r o s i t y ) d e r i v e d from a s o i l m o i s t u r e r e t e n t i o n c u r v e . A composite d r a i n a b l e p o r o s i t y was o b t a i n e d by T a y l o r and s u b s e q u e n t l y used by Hoffman and Schwab (1964) i n d r a i n a g e d e s i g n . Skaggs (1976) has proposed an i n d i r e c t method of o b t a i n i n g K/f r a t i o s from water t a b l e r e c e s s i o n s i n d r a i n e d or about t o be d r a i n e d f i e l d s . B h a t t a c h a r y a and Broughton (1979) p r e s e n t e d a method t o compute water t a b l e depths t h a t t a k e s i n t o account the v a r i a b l e n a t u r e of d r a i n a b l e p o r o s i t y . They r e p o r t e d f u n c t i o n a l forms of f f o r two s o i l s (Upland sand and S t e . R o s a l i e c l a y ) . The outflow-measurement methods f o r c h a r a c t e r i z i n g f a r e q u i t e d i f f i c u l t i n p r a c t i c e . T h e r e f o r e , e f f o r t s were made t o get an e s t i m a t e of the d r a i n a b l e p o r o s i t y from the s o i l m o i s t u r e r e t e n t i o n c u r v e s . However, the s o i l m o i s t u r e r e t e n t i o n c u r v e s 65 are generally not e a s i l y obtainable and for most l o c a l s o i l s , their complete characterization was not a v a i l a b l e . In most p r a c t i c a l designs drainable porosity i s treated as a constant due to the d i f f i c u l t i e s in characterizing the exact nature of the drainable porosity/pressure head realtionship (Van SchiIfgaarde, 1970). Moreover the drainage designs performed in the past treating f as a constant have withstood the test of time and unless some more p r a c t i c a l methods of obtaining i t s functional form are made ava i l a b l e , t h e i r use in the present designs should continue. Dumm (1968) has shown a general rel a t i o n s h i p between hydraulic conductivity and drainable porosity. Sometimes an e f f e c t i v e f can be estimated from the water table recession measurements for the times when i t rains. It i s assumed that the t o t a l r i s e in the water table caused by r a i n f a l l i s due to the net effect of simultaneous drainage and r a i n f a l l , that i s f * D Z = R- q - & SMC - ET where £ SMC is the change in moisture content during the rain period, ET i s the evapotranspiration during that period, q i s the drainage rate that took place during that period, R i s the t o t a l rain that f e l l during that period and DZ i s the actual r i s e in the water table over that period. The observations of water table height were taken when ET was small and the s o i l had previously been at f i e l d capacity. Therefore, the above equation s i m p l i f i e s to 66 f * DZ = R - q (8) Using the above method, an average value of f was found to be approximately 0.03 which i s very close to the value recommended for such clay s o i l s . Therefore, a constant f value of 0.03 was used in the subsequent analysis. 3.5.3 Hydraulic conductivity Hydraulic conductivity (K) i s a measure of the a b i l i t y of the s o i l to conduct water. Its e f f e c t i v e value i s one of the most important parameters in drainage/subirrigation design. Ideally, i t should represent the saturated hydraulic conductivity of that part of the s o i l p r o f i l e where most of the flow occurs during drainage and su b i r r i g a t i o n . There are many methods available for estimating hydraulic conductivity (Bouwer and Jackson, 1974). They are mostly point measurement methods and they are available for the case when the water table i s present or when i t i s absent. In the presence of a water table, methods are the auger hole method (Ernst, 1950; Kirkham and Van Bavel, 1948; Van beers, 1970); the tube method (Frevert and Kirkham, 1948); the piezometer method (Luthin and Kirkham, 1949; Young, 1968); the two-well method (Childs, 1952; Childs et a l . , 1953); the multiple well technique (Smiles and Young, 1963) and the four well method (Kirkham, 1954; Snell and Van Schilfgaarde, 1964). In the absence of the water table, methods are the shallow well pump-in method (Winger, 1960; Boersma, 1965); the i n f i l t r a t i o n gradient technique (Bouwer, 67 1964); the a i r entry permeameter (Bouwer, 1966) and the double tube method (Bouwer, 1961 and 1962). Shady et a l . (1977) reported on the use of the auger hole method to estimate f i e l d scale hydraulic conductivity in Canada. They reported that on the basis of these estimations, wider drain spacings are used now than was the practice before the tests were performed. Most of the above methods can be readily applied and can y i e l d answers in a r e l a t i v e l y short time. But the resulting K values represent only a single point in the f i e l d . To determine an e f f e c t i v e (equivalent) K, numerous tests are required at d i f f e r e n t locations in the test area. This i s often not feasible due to the time and costs involved in such an exercise. Moreover i t i s unclear at the present time how many tests are required to quantify the f i e l d hydraulic conductivity. The only thing that we know with certainty i s that the greater the number of tests, the greater i s our confidence in the resultant mean hydraulic conductivity. S t a t i s t i c a l approaches hold promise in such cases but there have been few applications in p r a c t i c a l drainage designs. A better way to determine the equivalent K i s to estimate i t from the phenomenon that w i l l be ultimately predicted by i t s use. Therefore, an indirect way to estimate the e f f e c t i v e K i s to use water table recession hydrographs (Hooghoudt, 1938; Hoffman and Schwab, 1964; Skaggs, 1976 and 1979). The e f f e c t i v e K values are estimated by back c a l c u l a t i n g from the observed water table drawdown in a given time period. The process i s repeated for d i f f e r e n t drawdown events and an e f f e c t i v e K value 68 is thus estimated. These methods take into consideration, though i m p l i c i t l y , the e f f e c t s of s o i l p r o f i l e heterogeneities, nonuniformities and anisotropies (Skaggs, 1979). The resultant K value i s a lumped parameter that takes into account the above anomalies in a p r a c t i c a l way. Their characterization in p r a c t i c a l designs is otherwise very d i f f i c u l t , i f not impossible. Further, errors due to unknown location of the depth to the impermeable layer are also reduced in such methods. The only catch i s that t h i s method needs more time and expense than do point measurement techniques. However, such methods are p r a c t i c a l i f they are included in the advanced planning of drainage/subirrigation projects. The required measurements of the water table fluctuations can then be planned and subsequently taken after laying one or two l a t e r a l s or the main l i n e (Skaggs, 1976). The benefits drawn from such an exercise may be considerable due to the savings possible with improved designs. For the present study, both the auger hole method and the indirect method were used to determine e f f e c t i v e hydraulic conductivity. The auger hole method (Van Beers, 1970) was used to determine point K values. These tests were performed at four locations in the experimental f i e l d . Their locations are marked in Figure 3. Holes were dug into the s o i l with an auger to a depth approximately 1.2 m below the s o i l surface (equal to drain depth). When water in the hole reached equilibrium with the surrounding groundwater, part of i t was bailed out quickly with a b a i l e r . Ground water began to seep into the hole and the rate 69 of i t s r i s e in the hole was recorded. The hydraulic conductivity of the s o i l was computed from a graph describing the rel a t i o n s h i p among the rate of r i s e , the equilibrium groundwater l e v e l , and the hole geometry. This method y i e l d s an average K of a s o i l column about 30 cm in radius and extending from the water table to about 20 cm below the bottom of the hole. The computed values are shown in Table 2. The l a s t three values in the Table are more than double the f i r s t value. On r e - v i s i t i n g the f i e l d , i t was found that the last three holes were located in s l i g h t l y depressed areas and therefore, seepage from surrounding areas probably perturbed the observations to some extent. The f i e l d has an undulating top surface owing to the fact that i t has not been under a g r i c u l t u r a l production before. At the time of i n s t a l l i n g the drainage/subirrigation system in the f i e l d the land surface was not graded and . l e v e l l e d due to limited resources available at that time. Therefore, i t was concluded that not much weight can be attached to the last three entries in Table 2. Only the hole in location "A" was fortunately located in a r e l a t i v e l y f l a t area and may represent the K value for that point. It should be noted that variations of up to one order of magnitude or more are quite common with the auger hole tests. The experience described in the above discussion points towards the importance of c a r e f u l l y selecting r e l a t i v e l y f l a t areas for performing these tests and i l l u s t r a t e s the d i f f i c u l t y in .obtaining representative r e s u l t s . 70 TABLE 2 RESULTS OF AUGER HOLE TESTS AND INDIRECT METHOD TO MEASURE HYDRAULIC CONDUCTIVITY Descr i p t ion K cm/day Loca t ion #1 Loca t ion #2 L o c a t i o n #3 L o c a t i o n #A From the water t a b l e drawdown (assuming f = 0.03) 10.00 22.00 32.00 27.00 2.26 71 In view of the above factors, i t was decided to use the indi r e c t method of Skaggs (1976 and 1979) to calculate a f i e l d e f f e c t i v e value of K from water table recession measurements. A few r a i n l e s s periods were selected and an e f f e c t i v e K was calculated using nondimensional plots of water table drawdown given by Moody (1966). Similar plots given by Skaggs (1975) were not used because the water table was not f l a t i n i t i a l l y as i s required in his method. The above periods were dif f e r e n t from the ones used for model v a l i d a t i o n . However, the average e f f e c t i v e value of K can be computed only i f drainable porosity is known. For drainable porosity of 0.03 the average e f f e c t i v e K value i s 2.26 cm/day. Whereas the indirect method lumps together various unknowns in the drainage design such as the location of the impermeable layer, the r a d i a l flow near the drains, the entry resistance to flow near the drains and the saturated hydraulic conductivity of the s o i l into a single parameter, the auger hole method provides only point estimates for the saturated hydraulic conductivity of the s o i l . As discussed e a r l i a r , the e f f e c t s of s o i l p r o f i l e heterogeneities, nonuniformities and anisotropies are also i m p l i c i t in an indi r e c t method. 3.5.4 Depth to the impermeable layer The impermeable layer i s a layer that has a saturated hydraulic conductivity approximately one-tenth of the saturated hydraulic conductivities of the other layers. The effect of an impermeable layer below the drains on a drainage system is very 72 pronounced. As depth to the impermeable layer increases below the drains (d), the height of the water table above the drains decreases. This i s due to lesser resistance to the flow as d increases. Although present drainage/subirrigation theory i s s u f f i c i e n t to give p r a c t i c a l answers when d i s zero or at i n f i n i t y , i t i s inadequate when d i s not at one of these extremes. However, i t s q u a n t i f i c a t i o n is indispensable in drainage/subirrigation system design because i t forms one of the boundary conditions in the design. As stated e a r l i e r , most of the available p r a c t i c a l design tools incorporate D-F assumptions. In these formulations, i t is often suggested to use some equivalent depth to the impermeable layer instead of d. It i s argued that convergence losses near the drains are accounted for by t h i s substitution (Van Schi1fgaarde, 1963). It should be noted that the equivalent depth to the impermeable layer becomes independent of the actual depth above a certain threshold value for a given drainage system (Van SchiIfgaarde, 1974; Sakkas, 1975; Boumans, 1979; Chieng, 1980; Sakkas and Antonopoulos, 1981). Therefore, i f d is greater than the threshold value, i t s effect on the drainage system design is minimal. I n t u i t i v e l y , this seems correct also because an impermeable layer at a distance more than 5 m below the drains may not af f e c t the flow net at a l l in most cases of p r a c t i c a l interest. Boumans (1979) recommends the use of an equivalent depth of one-quarter of the drain spacing for those cases where i t i s unknown. It may be noted that the above substitution i s to be done when the position of the impermeable 73 layer i s unknown due to i t s distance below the drains. During the course of thi s investigation, a new method of designing drainage systems based upon the steady state theory has been developed. This method does not require any pri o r knowledge about s o i l hydraulic conductivity and depth to the impermeable layer. For most of the cases tested to date, i t seems to lead to the same system design as would have resulted from standard drainage models (Kirkham*s method, Hooghoudt's method) with e x p l i c i t knowledge of the above parameters. Since, the main thrust of the thesis i s with nonsteady state design, this new method i s discussed in Appendix B. To detect the location of d for the experimental plot, a hole was dug to a depth of 3 m. There was no indication of the presence of an impermeable layer. Therefore, for the system in question, an equivalent depth of 1 m was selected (Chieng, 1980). This i s the maximum equivalent depth that can af f e c t flows for the drainage system in the experimental pl o t . 3.6 Model Validation In the vali d a t i o n phase of mathematical modelling, model predictions are compared with actual f i e l d observations on the underlying physical process. A statement can then be made as to the appropriateness of the model in hand and the data employed by the model. If results are unsuitable, some of the input data can be varied and the model i s rerun. This exercise i s repeated u n t i l an acceptable reproduction of the observed behaviour i s achieved. 74 The s o i l input data used to validate the model is l i s t e d in Table 3. P r e c i p i t a t i o n data was also required in the model va l i d a t i o n . F i r s t l y , i t was thought that the p r e c i p i t a t i o n data recorded at the Vancouver International Airport (about 10 Km from the experimental plot) could be transposed to the research plo t . But soon i t became apparent that p r e c i p i t a t i o n patterns vary considerably from the sea shore inland. This variation was further confirmed by comparing the data recorded at the airport with another weather station in the Ladner south region , about 10 Km from both the airport and the reserch pl o t . The two sets of data vary considerably in terms of the i n t e n s i t i e s of p r e c i p i t a t i o n . Therefore, a tipping bucket rain gauge was set up in the plot to record p r e c i p i t a t i o n representative of the ' area of the pl o t . The recording covered the period of March 5 to March 25, 1982. The model predictions are compared with the observed data in Figure 4. It seems that the model is appropriate for simulating the response of a drainage system at least for that time of the year. This i s to be expected because the Boussinesq equation i s applicable to saturated flows and in early spring the s o i l conditions j u s t i f y the use of a saturated flow analysis. It should be noted that the methodology to be developed in later sections of t h i s chapter i s independent of the mathematical model used. What is required i s some model which can relate drainage/subirrigation system design with appropriate inputs. The intent of the model va l i d a t i o n was to f i n d out the appropriateness of the Boussinesq equation for simulating the 75 TABLE 3 INPUT DATA USED IN THE MODEL VALIDATION Parameter Val ue H y d r a u l i c c o n d u c t i v i t y D ra inab le p o r o s i t y E q u i v a l e n t depth to the impermeable l a y e r below the d r a i n cent re Dra in spac ing Dra in depth Dra in d iameter 2.26 cm/day 0 .03 100.00 cm 1400.00 cm 120.00 cm 10.00 cm Measured Calculated On 10H 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 March 1982 Figure 4. Comparison of the observed and calculated midspan water table depth below the soil surface in spring 1982. 77 drainage process in spring. This i t does s a t i s f a c t o r i l y . 3.7 Stochastic Model Workability i s d i r e c t l y related to moisture tension in the s o i l . Moisture tension, in turn, i s dependent, among others, on depth of the water table from the s o i l surface. Further, water table depth i s greatly affected by a r t i f i c i a l drainage. Therefore, workability i s related to the degree of drainage. However, workability in spring i s especially related to the moisture tension in the plough layer of the s o i l . Moisture tension in t h i s layer i s d i r e c t l y related to weather. The p r o b a b i l i t y of obtaining workable days becomes very high i f the weather i s dry in the late winter and early spring (that i s , i f i t does not rain too often and there i s some evapotranspiration during that period). Dependence upon weather can be reduced considerably i f the f i e l d i s a r t i f i c i a l l y drained by means of drainage systems. Timely removal of excess water from the root zone by drainage increases the probability of obtaining the required number of workable days that are essential for performing f i e l d work. Therefore, concepts of r i s k and uncertainty associated with weather should become an integral part of a ra t i o n a l drainage system design. Van Schilfgaarde (1965) was perhaps the f i r s t to consider v a r i a b i l i t y of the weather in drainage design in a p r o b a b i l i s t i c way. Using a water budget model, he was able to calculate lengths of periods during which the water table exceeded a given height for a given return period. Other drainage models that 78 have been used to ob t a i n s i m i l a r i n f o r m a t i o n are those of Chieng (1975) and Skaggs(1978). The above models are j u s t a few of the many conceptual models of s i m u l a t i n g the drainage p r o c e s s . Most of these models i n c o r p o r a t e d i f f e r e n t c o n c e p t u a l i z a t i o n s with r e s p e c t t o the u n d e r l y i n g drainage process but the treatment of c l i m a t i c u n c e r t a i n t y i n almost a l l of them i s the same. Each model i s made to run f o r the p e r i o d of h i s t o r i c r e c o r d on recorded weather data and the performance of a given drainage system i s analysed to f i n d how that system would have performed over the p e r i o d of r e c o r d . Normally, the s i m u l a t i o n s are done at d a i l y i n t e r v a l s s i n c e one day i s the most common i n t e r v a l f o r which most of the r e q u i r e d weather data such as r a i n f a l l i s a v a i l a b l e . That i s , a model i s made to run f o r t h i r t y or more years on a d a i l y b a s i s and the performance of a given drainage system i s then s t u d i e d . T h i s process i s repeated f o r v a r i o u s combinations of d r a i n spacing and depths ( d i f f e r e n t system d e s i g n s ) . A drainage system that f u l f i l s the design c r i t e r i o n s a t i s f a c t o r i l y in an average year i s then s e l e c t e d f o r that l o c a t i o n . The above procedure i s repeated i f e i t h e r the l o c a t i o n or s o i l type changes. Flaws i n t h i s type of approach are q u i t e obvious. The s t o c h a s t i c nature of the weather i s not co n s i d e r e d i n such s i m u l a t i o n s . Models are run f o r recorded weather p a t t e r n s o n l y . The u n d e r l y i n g assumption i s that weather p a t t e r n s i n f u t u r e w i l l be the same as they were i n the p a s t . However, i f a design i s done using the l o c a l r e c o r d and that r e c o r d happens to be of 79 s u f f i c i e n t length, the design w i l l probably be quite safe. Another problem with the conventional way of designing drainage systems on long-term basis i s that a new simulation i s required every time a new drainage system design needs testing, or the location changes with a corresponding change in s o i l type. It can cause a considerable computational burden in terms of the costs involved and the exercise of finding an optimal design may run into f i n a n c i a l constraints. The magnitude of these costs w i l l vary with the model used for the design but in a l l cases, they w i l l remain s i g n i f i c a n t . In view of the above observations, a f i r s t order Markov chain model for water table fluctuations was formulated to consider weather uncertainty in drainage design in a rati o n a l way. This approach eliminates the need for repeated use of long-term simulations as required by most models and therefore a s i g n i f i c a n t reduction in computational costs i s achieved. Since Markov chains w i l l be repeatedly referred to in the following discussion, the concept i s discussed below in d e t a i l . 3.7.1 Markov processes A random process i s described as Markovian when dependence exi s t s between sequential events and the dependence i s such that the next step depends only on the present state and not on other past states of the system. This step-by-step dependence , c a l l e d a f i r s t order process, has lag-one autocorrelation. However, the theory can be extended to multiple-lag autocorrelations among sequential events (higher order Markov 80 processes). We w i l l be discussing only the f i r s t order Markov process in t h i s section. The following discussion follows c l o s e l y the description of Markov processes given by H i l l i e r and Lieberman (1981). The Markovian property implies that the conditional p r o b a b i l i t y of any future event, given any past event and present state of the system, is independent of the past event and depends only upon the present state of the process. The conditional p r o b a b i l i t i e s P { x(t+1)= j | X(t)= i } are c a l l e d " t r a n s i t i o n " p r o b a b i l i t i e s . Transition p r o b a b i l i t i e s (one-step) are stationary i f P { X(t+1)= j | X(t)= i } = P { X(1)= j | X(0)=i } for a l l t=0,1,... Stationary p r o b a b i l i t i e s are generally denoted by p^•. S t a t i o n a r i t y refers to their independence with respect to time, that i s , for each i , j , and n (n=0,1,2,.... ) P { X(t+n)= j | X(t)=i }= P { X(n)= j | X(0)=i } for a l l t=0,1,... 81 00 These conditional p r o b a b i l i t i e s are denoted by p.^ and are c a l l e d n-step t r a n s i t i o n matrices. Therefore, p ^ i s the conditional p r o b a b i l i t y that a random process, s t a r t i n g in state i , w i l l be in state j after exactly n t r a n s i t i o n s . If we deal with discrete processes, the Markov process i s c a l l e d a Markov chain. A f i n i t e state Markov chain includes the following: (1) A f i n i t e number of states (2) The Markovian property (3) Stationary t r a n s i t i o n p r o b a b i l i t i e s (4) I n i t i a l p r o b a b i l i t i e s , P { X(0)= i } The behaviour of such Markov chains is described by square t r a n s i t i o n matrices which have the following general form ( n ) ( n ) (n) POO ( n ) ( n ) ( n ) P10 ( n ) ( n ) ( n ) PmO Where (n) P.. > 0, for a l l 1 and j , ^ 3 ™" and n = 0 , 1, 2. 82 and (n) m j=0 1, for a l l 1 and j , and n = 0 , 1» 2. The i t h row of the stochastic matrix (p. # , P., . P( i s a p r o b a b i l i t y vector of a l l possible outcomes of the next state, given i t i s i at present. The Chapman-Kolmogorov equations are used to compute n-step t r a n s i t i o n p r o b a b i l i t i e s , that i s , the pro b a b i l i t y that the process w i l l be in state j after n periods, given i t i s in i at present. These equations are (n) m (v) ( n - v ) Plk B z Pik Pkj • f o r a 1 1 1 » J » n k = 0 and 0 < v < n (v) ( n - v ) The above equations imply that Pik P^j i s the conditional p r o b a b i l i t y that, s t a r t i n g from state i , the process reaches state k after v steps and then to state j in (n-v) steps. The importance of these equations l i e s in the fact that n-step t r a n s i t i o n p r o b a b i l i t i e s can be obtained from the one-step t r a n s i t i o n p r o b a b i l i t i e s r e c u r s i v e l y . For example, i f n=2, these equations become 83 (2) m P i j * klQ P 1k p k j , for a l l 1 , j . It may be noted that p i j are elements of matrix P These elements are obtained by multiplying the matrix of one-step t r a n s i t i o n p r o b a b i l i t i e s by i t s e l f , that i s P = P.P = P: Therefore, i t follows that the matrix of n-step t r a n s i t i o n p r o b a b i l i t i e s can be obtained as (") n n - 1 n - 1 P • P. P. P P - P « PP « p p that i s , the n-step t r a n s i t i o n matrix i s obtained by computing the nth power of one-step t r a n s i t i o n matrix. Another property of Markov chains that has a great p r a c t i c a l significance i s the steady state p r o b a b i l i t y of Markov chains. This property i s exhibited by a l l Markov chains after s u f f i c i e n t t r a n s i t i o n s have occured from one state to another. By steady state, we imply that the probability of finding the process in state j , afte r a large but f i n i t e number of tr a n s i t i o n s , tends to a fixed p r o b a b i l i t y , independent of the beginning state. It i s important to note that steady state p r o b a b i l i t y does not mean that the process " s e t t l e s " down into 84 one state. The process continues to make tr a n s i t i o n s from one state to another, and at any step n, the t r a n s i t i o n probability from state i to state j i s s t i l l p t^ . Mathematically, P i j = n j ( 9 ) n 0 0 with the following observations H. > 0 (9a) J m H. S I n. p. . , for j • 0 , 1 . . . . .M (9b) m I n. = 1 (9c) j=0 J where "FTj h are the steady state p r o b a b i l i t i e s of the Markov chains. Their reciprocals provide an estimate of the expected recurrence time for the states, that i s , the expected time which the stochastic process takes to return to i t s starting state on a long-term basis. Knowing the t r a n s i t i o n matrix, the steady state (long-run) properties of the Markov chain can be determined. 85 3.7.2 Stochastic model for drainage design Crop y i e l d s are i n d i r e c t l y related to water table fluctuations. High water tables lead to oxygen deficiency in the active root zone of the plants thereby hindering plant growth in a number of ways. S i m i l a r l y , low water tables lead to less moisture in the active root zone causing a decline in crop y i e l d s . This amounts to saying that although not a perfect "indicator" of the drainage system performance from a t h e o r e t i c a l point of view, use of water table fluctuations in p r a c t i c a l designs i s reasonable (see section 3.6.4). Water tables fluctuate in response to the weather. Both p r e c i p i t a t i o n and evapotranspiration and groundwater recharge flow rates at the bottom of the system (in the absence of any other l o c a l source or sink) govern i t s movement. Water tables r i s e or f a l l due to the net effect of these hydraulic forcing functions. The phenomenon of r i s i n g and f a l l i n g water tables also occurs in drained f i e l d s . Therefore, the design of a drainage system i s intended to minimize those sequences of fluctuations that are harmful to plants. Dependence upon weather being evident, the risk associated with weather becomes a risk in the drainage system design. A f i n i t e state Markov chain model for water table fluctuations i s proposed in t h i s study. The stochastic model is described below and i t s application in p r a c t i c a l drainage system designs in the Lower Fraser Valley i s i l l u s t r a t e d with the help of examples. 86 3.7.2.1 Transition matrix The essence of simulating a stochastic process by Markov chains l i e s in the construction of a t r a n s i t i o n matrix. The matrix contains p r o b a b i l i t i e s of going from one state, which in th i s case i s defined as an elevation, to another in a given time period. The period of the t r a n s i t i o n matrix selected in the study was one day. It was arrived at by considering p r e c i p i t a t i o n patterns in the Valley for spring months. The autocorrelation c o e f f i c i e n t for the net recharge decays very fast for t h i s period and i s v i r t u a l l y zero at lags higher than two. At lag of unity, i t . i s positive and f i n i t e but for almost a l l locations in the Valley, i t is less than 0.30. A computer program was written to solve the Boussinesq equation numerically. The output from the program i s the net recharge that i s necessary to move the water table from a given state to a l l others for a given drainage system and for a given s i t e . A sample output from t h i s program is shown in Table 4. Another computer program was written to analyze net p r e c i p i t a t i o n data, that i s p r e c i p i t a t i o n minus evapotranspiration, in winter and spring in the Valley. Output from t h i s program gives a cumulative probability d i s t r i b u t i o n function of net recharge. The function i s obtained by doing a frequency analysis of d a i l y net recharge values over the whole length of the h i s t o r i c a l record available for the s i t e . The pro b a b i l i t y of occurrence of any event that i s less than or equal to a stated value can now be obtained from i t . Figure 5 87 TABLE A OUTPUT FROM THE DRAINAGE MODEL BASED ON THE BOUSSINESQ EO.UATION Net Rai nfal1 (mm) Midspan water t a b l e height above the impermeable layer a f t e r one day (mm) 21.112 215.0 15.9^0 205.0 10.810 195.0 5.750 185.0 0.620 175.0 - *».370 165.0 - 9.370 155.0 - Ht.310 11*5.0 - 19.250 135.0 - 2k.]20 125.0 - 29.000 115.0 - 33.750 105.0 - 36.190 100.0 NOTE: I n i t i a l midspan water t a b l e height above the impermeable layer = 215.0 cm 88 89 shows these functions for Abbotsford Airport in B.C. for winter and spring periods. Net recharge causing the water table to move from one elevation to another can now be obtained from Table 4. The probability of the net recharge i s thus obtained from Figure 5. At t h i s stage, the t r a n s i t i o n matrix can be completed from the above analysis for both periods. Transition matrices thus constructed for winter and spring periods are shown in Tables 5 and 6. There was some concern with regard to the number of states to be used in the construction of the t r a n s i t i o n matrix. The greater the number of states to be considered, the more accurate the analysis. However, the computational burden also increases exponentially. For a given s i t e , the underlying p r i n c i p l e i s that the states should not be so "coarse" that i t i s p r a c t i c a l l y impossible for the water table to move from one state to another in the period of the t r a n s i t i o n matrix. There may be weather constraints for the movement such as lower p r e c i p i t a t i o n or evapotranspiration rates or there may be l i m i t a t i o n s with respect to the conducting capacity of the s o i l . The number of states in the t r a n s i t i o n matrix should be such that the states can "communicate" with adjacent states for that s i t e . I n i t i a l l y , a fine sandy s o i l with a state i n t e r v a l of 30 cm was used to check the hypothesis that i f adjacent states can communicate with each other, then any int e r v a l length less than or equal to the maximum int e r v a l can be used. Simulations were repeated using a state i n t e r v a l of 10 cm. There was v i r t u a l l y Midspan w a t e r t a b l e h e i g h t above the d r a i n c e n t r e a t b e g i n n i n g o f the p e r i o d (cm) 06 TABLE 6 SPRING TRANSITION MATRIX (March and A p r i I ) 0) *-* 4-1 ro <D i _ *-* c 0) u c •— e (D u l_ TJ TJ 0) O .c 4-» 1-8. > o <u nj 4-J -C o cn •— O) <u c .c c 0) c -D cn (TJ <u •M U 0) 4-1 (0 3 c (TJ Q. tn XJ . 0 0 5 . 0 0 5 . 0 1 0 . 0 0 5 . 0 1 0 . 0 1 0 . 0 1 0 . 0 2 5 . 0 3 0 . 0 6 0 . 1 3 0 . 7 0 0 . 0 0 5 . 0 1 0 . 0 1 0 . 0 0 5 . 0 1 0 . 0 1 5 . 0 2 5 . 0 2 0 . 0 7 0 . 1 1 0 . 5 0 0 . 2 2 0 . 0 1 0 . 0 0 5 . 0 1 0 . 0 1 0 . 0 1 5 . 0 1 5 . 0 3 0 . - 3 5 . 0 o 0 . 2 3 0 . 4 8 0 .070 . 0 1 5 . 0 1 0 . 0 1 0 . 0 1 0 . 0 2 0 . 0 2 0 . 0 4 0 . 0 7 5 . 1 7 0 . 5 3 0 . 1 0 0 . 0 0 0 . 0 2 0 . 0 1 0 . 0 1 0 . 0 1 5 . 0 2 0 . 0 3 5 . 0 5 0 . 1 2 0 . 4 1 0 . 3 1 0 . 0 0 0 . 0 0 0 . 0 3 0 . 0 1 0 . 0 1 0 . 0 1 5 . 0 2 5 . 0 4 0 .080 . 2 4 0 . 4 8 0 . 0 7 0 . 0 0 0 . 0 0 0 . 0 3 5 . 0 1 0 . 0 1 5 . 0 3 0 . 0 2 0 . 0 8 0 . 1 5 0 . 5 4 0 . 1 2 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 4 0 . 0 1 0 . 0 3 0 . 0 2 0 . 0 6 0 . 1 2 0 . 4 0 0 . 3 2 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 4 5 . 0 2 0 . 0 2 0 .090 . 0 2 5 . 1 8 0 . 5 1 0 . 1 1 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 5 5 . 0 2 5 . 0 3 0 . 0 6 0 . 1 2 0 . 5 0 0 . 2 1 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 6 5 . 0 2 5 . 0 8 5 . 0 4 5 . 1 8 0 . 5 2 0 . 0 8 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 8 0 . 0 2 5 . 0 6 0 . 1 2 5 . 4 9 0 . 2 2 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 . 0 0 0 2:120 110 100 90 80 70 60 50 40 30 20 10 M i d s p a n w a t e r t a b l e h e i g h t above t h e d r a i n c e n t r e a t t h e end o f the p e r i o d (cm) 92 no difference in results from the two simulations. However, for standard purposes, an interval of 10 cm i s recommended. For the Valley, i t encompasses almost a l l possible combinations of s o i l properties and net recharge. As stated e a r l i e r , any other physically possible combination can be used as may be required for any p a r t i c u l a r s i t e . It must be noted that these combinations are a function of the period of the t r a n s i t i o n matrix. If the period needs a change due to any pa r t i c u l a r reason, the number of states w i l l be changed. For i l l u s t r a t i o n purposes, a medium textured s o i l (fine sand) was used in the study. The re s u l t i n g t r a n s i t i o n matrices for winter and spring are shown in Tables 5 and 6 respectively. It i s noted that according to the requirements of a f i n i t e state t r a n s i t i o n matrix, the t r a n s i t i o n p r o b a b i l i t i e s in any row add up to unity. Physically, i t implies that the water table can move only to one of the various states considered herein from any s t a r t i n g state. Having considered a l l possible states, the p r o b a b i l i t i e s in a row should add to unity. A zero entry in the matrix implies that the water table can not move from another state to t h i s state due to physical constraints (low evapotranspiration rates in early spring e t c . ) . An entry of one in the matrix s i g n i f i e s a "trapping" state, that i s , once the system moves into this state from any other state, i t can not move out . Any other entries correspond to the t r a n s i t i o n p r o b a b i l i t i e s for moving from one state to another in the period of the stochastic matrix (1 day). A check was applied to determine i f the t r a n s i t i o n matrix 93 exhibits s t a t i o n a r i t y . The steady state property can be checked by repeatedly multiplying the one-step t r a n s i t i o n matrix by i t s e l f . The entries in the resulting matrices should cease to change after few ite r a t i o n s indicating thereby that steady state has been reached. A further check is usually made to see i f the steady state values are independent of the beginning state. A l t e r n a t i v e l y , the steady state p r o b a b i l i t i e s can be computed by using equations ( 9 ) , ( 9 a ) , (9b) and ( 9 c ) . Both checks were made on the t r a n s i t i o n matrices for winter and spring periods. The steady state was reached rapidly (25 days on an average) indicating that the assumption of steady state i s v a l i d for such cases. A variety of beginning states were employed and yet steady state was reached very soon and the steady state values were the same irrespective of the starting state. Two t r a n s i t i o n matrices were constructed: one each for winter (November to February) and spring (March and Ap r i l ) periods. Both of them exhibited s t a t i o n a r i t y and passed the tests for steady state. For the winter t r a n s i t i o n matrix, steady state was reached after about 25 days from the beginning in November. It implies that the assumption of steady state at the end of winter period i s more than s a t i s f i e d . The resulting steady state p r o b a b i l i t i e s for the winter period are shown in Table 7. As stated e a r l i e r , these are the long-term p r o b a b i l i t i e s that the water table w i l l be in one of the various states for a given drainage system and for a given s i t e . These long-term p r o b a b i l i t i e s can be read d i r e c t l y from the steady TABLE 7 STEADY STATE PROBABILITIES FOR THE WINTER TRANSITION MATRIX 01 JZ * J 4-J <0 V L-4-1 c Ol u c .— <o u i_ T 3 -a 01 0 JZ * J i_ 01 01 a. > o Ol .o x: u-JZ o O l .— 0) c x: c 01 c - O O ! (0 0> +J u 01 +J 10 3 c 10 a ul '£ 10 20 30 140 50 60 70 80 90 100 no >120 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.000 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0.120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.000 0.130 0.220 0.170 0.120 0.0*10 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0. 120 0.0*40 0.030 0.0*10 0.060 0.020 0.030 0.050 0.090 0.130 0.220 0.170 0. 120 0.0*10 0.030 I-l 1 0.0*10 Midspan water table height above the drain centre at the end of the period (cm) 95 state p r o b a b i l i t y matrix. For example, there i s 22% pr o b a b i l i t y that the water table w i l l be located in a depth i n t e r v a l of 50 cm to 60 cm from the s o i l surface on a long-term basis. Similar deductions may be made for other states. It should be kept in mind that t h i s p robability i s r e l a t i v e to a p a r t i c u l a r drainage system in a certain location. Results are not transferable and a new pr o b a b i l i t y w i l l result i f either the drainage system design or the location i s changed. Steady state p r o b a b i l i t i e s were assumed to occur at the end of the winter period (that i s , end of February) because they were achieved within 25 days of the beginning of the winter period. These served as i n i t i a l conditions for t r a n s i t i o n s in the water table to various other states in the spring as governed by the spring t r a n s i t i o n matrix. Ideally, the pr o b a b i l i t y of obtaining the f i r s t ten workable days in the month of March can be estimated from simple p r o b a b i l i s t i c manipulations. That i s , given that at the beginning of March, the water table is not in an undesirable state, the probability of i t not going to t h i s state after the end of one day can be estimated using the spring t r a n s i t i o n matrix. The repeated use of t h i s process can y i e l d a pro b a b i l i t y for any desired number of consececutive days that the water table i s not in that undesirable state. However, farmers may not be interested in the f i r s t ten days of March. What they may require i s any ten consecutive days of workable period in March. There are twenty one ways in which t h i s i s possible, that i s , any of the periods 1-10, 2-11, 3-12 etc. i s 96 of interest to them. They may or may not have any special regard for a s p e c i f i c period. In most cases, what they might want i s one workable period in early spring. Therefore, the pr o b a b i l i t y of getting at least one workable period i s a l l that may be required. It i s assumed that although more than one workable period i s a p o s i t i v e event for the farming enterprise, the design may be based upon obtaining at least one such period. This problem can be solved by making use of the union pro b a b i l i t y formula from s t a t i s t i c s . It can y i e l d a pr o b a b i l i t y of occurrence of at least one event when there are many events that can occur either at the same time or any one of them occurs. The formula can be used only i f the number of events under consideration i s small because the analysis becomes very complicated and perhaps intractable for events greater than ten; For example, l e t us look at the p r o b a b i l i t y of occurrence of at least one event given that the t o t a l number of such events i s two, that i s , P[A U B ] = P ( A ) + P ( B ) - P ( A f \ B) Where P [ A U B ] i s the p r o b a b i l i t y that either event A occurs or event B occurs or both occur; P ( A ) i s the probability that event A occurs; P ( B ) i s the pr o b a b i l i t y that event B occurs and P ( A I \ B ) i s the j o i n t p r o b a b i l i t y of events A and B, that i s , the p r o b a b i l i t y that both A and B occur. 9 7 The analysis i s quite tractable as long as the number of events i s small. The generalized form of the above formula i s p {AjU A 2U . . . . uAn> = z p(Ak) - i H^r\\) +1 P ^ f t A . . ^ ) +(-!)" P C A ^ f X AA n) where £.P(An) is the sum of p r o b a b i l i t i e s of from 1 to n, $ P ( A j O A^) i s the sum of p r o b a b i l i t i e s of A^ f\ A ^ with j and k from 1 to n and k > j , and so on. The above formula, although i t looks quite simple, i s d i f f i c u l t to apply when the number of events i s large. An attempt was made to write a computer program to use the above formula but d i f f i c u l t y soon became obvious. The analysis became very complicated and keeping track of the various terms became almost impossible. On the advice of the sta f f of the Computing Centre at the University of B r i t i s h Columbia, a simulation approach was used. To estimate the required p r o b a b i l i t i e s a computer program was written to perform the s t a t i s t i c a l simulation. Random numbers were generated to id e n t i f y the st a r t i n g condition for the simulation. Once i t was known, another set of random numbers were generated to keep track of the tr a n s i t i o n s from the i n i t i a l 98 c o n d i t i o n and so on. Some s t a t e s were tagged as u n d e s i r a b l e s t a t e s from the w o r k a b i l i t y p o i n t of vi e w . An account was kept w i t h r e g a r d t o t r a n s i t i o n s from a d e s i r a b l e i n i t i a l c o n d i t i o n t o o t h e r d e s i r a b l e s t a t e s w i t h t i m e . The pro c e d u r e was abandoned i f a t r a n s i t i o n t o an u n d e s i r a b l e s t a t e from an i n i t i a l l y d e s i r a b l e s t a t e took p l a c e o r i f the i n i t i a l c o n d i t i o n was found t o be u n d e s i r a b l e . The whole p r o c e d u r e was r e p e a t e d and by c o u n t i n g t he number of s u c c e s s e s i n a g i v e n number of t r i a l s , the p r o b a b i l i t y of o b t a i n i n g any r e q u i r e d number of c o n s e c u t i v e workable days was o b t a i n e d . The above procedure assumes o n l y two outcomes f o r the p r o c e s s , t h a t i s , s u c c e s s o r f a i l u r e c o r r e s p o n d i n g t o d e s i r a b l e and u n d e s i r a b l e s t a t e s r e s p e c t i v e l y from the s o i l w o r k a b i l i t y p o i n t of view. Such a p r o c e s s can be assumed t o f o l l o w a b i n o m i a l d i s t r i b u t i o n . Then the number of t r i a l s r e q u i r e d t o o b t a i n a g i v e n a c c u r a c y i n the e s t i m a t i o n of p r o b a b i l i t y can be det e r m i n e d as f o l l o w s : L e t us say t h a t we want t o be 95% c e r t a i n t h a t the number of t r i a l s s h o u l d be s u f f i c i e n t such t h a t t he v a l u e of the p r o b a b i l i t y t h a t we get from t h e t r i a l s i s w i t h i n 3% of the a c t u a l p r o b a b i l i t y . Assuming t h a t the p r o c e s s can be r e p r e s e n t e d by a b i n o m i a l d i s t r i b u t i o n , we have P{|(X/N-p) / y^q7N"| < 1.96} « 0.95 where 1.96 corresponds to 95% certainty in binomial distr ibution or -1.96 < T < 1.96 where T stands for the term on the l e f t side of the above equation. 99 or -1.96 v/pq7N < X/N - p < 1.96 /pajti Therefore, the spread (s) around the desired value is s = 2 * 1.96 * /pq/rT = 2 * 0.03 But f 1/2 (q - 1-p) It is possible i f we have number of t r i a l s such that 2 * 1.96 * /pq7N < 2 * 1.96 * 1/2 * 1 / » C < 2 * 0.03 or AT > 1.96/0.06 or N > 1068 It implies that i f the number of t r i a l s i s greater than or equal to 1068, the prob a b i l i t y obtained from the program w i l l be within 3% of the actual p r o b a b i l i t y as would have been obtained i f the number of t r i a l s were very large. Similar deductions can be made for other tolerance l e v e l s . The methodology can be summed up in the following steps: (1) Decide the number of states that are feasible for the given location (2) Use an appropriate drainage model to estimate the net recharge which would cause the water table to move from one state to another for a given drainage system at a 100 given s i t e (3) Estimate cumulative pr o b a b i l i t y d i s t r i b u t i o n functions of net recharge for the given s i t e and for periods of interest by performing frequency analysis on precipitation/evapotranspiration records (4) Construct t r a n s i t i o n matrices for each period. Entries for t r a n s i t i o n matrices are made by using information obtained in steps (1), (2) and (3) (5) Check the resulting matrices for s t a t i o n a r i t y and steady state (6) If they are found stationary and exhibit steady state, proceed with step (7); otherwise the analysis must be abandoned (7) Using steady state p r o b a b i l i t i e s from the winter t r a n s i t i o n matrix and the spring t r a n s i t i o n matrix, perform the s t a t i s t i c a l simulation to estimate the p r o b a b i l i t y of obtaining at least one N consecutive day workable period in March (8) Repeat th i s process from steps (1) to (7) for d i f f e r e n t drainage system designs and obtain corresponding p r o b a b i l i t i e s of obtaining at least one workable period (9) Choose the preferred design on the basis of i n i t i a l cost of the system and the r i s k 101 (that i s , the pr o b a b i l i t y of not obtaining a workable period in March) a farmer i s w i l l i n g to accept in the farm enterprise. To validate the stochastic model, a comparison was made between the new model and the "long-term deterministic" model of Chieng(1975). Two s o i l types were considered for t h i s purpose, namely a clay s o i l and a loam soil.. Input data used for model valid a t i o n i s shown in Table 8. Chieng's model was run for both cases using the d a i l y weather data obtained from 20 years of h i s t o r i c a l record for Abbotsford Airport. In both cases, Chieng's model predicted that there w i l l be at least one ten-day workable period (when the water table stays below 60cm from the s o i l surface). The Markov chain model gave p r o b a b i l i t i e s of 89% and 92% for clay and loam s o i l s respectively. These results are not too d i f f e r e n t from those obtained using Chieng's model. Therefore, i t i s believed that the proposed model is capable of y i e l d i n g equally good or better solutions to the drainage problem .than the present design tools. The beauty of t h i s approach l i e s in the fact that repeated long-term simulations that would be required for d i f f e r e n t drainage designs are avoided. Every time a new drainage system design i s to be performed, the procedure given e a r l i e r is followed. Gain in computational time with t h i s approach w i l l vary with d i f f e r e n t drainage models used to simulate the drainage process. For example, the proposed model requires about one-third or less computational time for a single run as compared to Chieng's model. There i s a d e f i n i t e reduction in the computational time 102 TABLE 8 INPUT DATA USED IN THE COMPARISON OF CHIENG'S MODEL AND THE PROPOSED MODEL Parameter Va1ue fo r Clay Soi1 Loam Soi1 Drainage C o e f f i c i e n t Depth o f the impermeable l a y e r below the d r a i n cent Drain depth D r a i n a b l e p o r o s i t y Drainage f a c t o r (C = i» K/S 2) re 12 mm/day 1000 mm ]k00 mm 0.0k 0.000004 mm day 12 mm/day 1000 mm \k00 mm 0 . 0 7 0.000004 mm day 1 03 with the model proposed herein because we do not have to run a drainage model for the f u l l record length r e p e t i t i v e l y as i s required in almost a l l p r o b a b i l i s t i c models of drainage design. Design graphs are shown in Figure 6 for some l o c a l s o i l s . The p r o b a b i l i t y of obtaining a workable period in March can be d i r e c t l y read from the graph corresponding to a given drainage system and s o i l type. The importance of good s o i l and water management i s also apparent from the.figure. The p r o b a b i l i t y of obtaining a workable period in March increases s i g n i f i c a n t l y with saturated hydraulic conductivity of the s o i l for a given drainage design. In other words, improvement in s o i l structure can now be expressed in quantifiable terms and underlines the importance of good s o i l and water management on the farm. An observation that has a strong bearing on p r a c t i c a l designs was also made from Figure 6. For fine textured s o i l s , the drainage system w i l l ensure at least one workable period in March i f the drain spacing obtained from conventional methods using steady state theory i s made smaller by 25%. This finding should be of interest to designers who do not want to indulge in complicated design procedures due to data l i m i t a t i o n s but, at the same time, want to upgrade design of drainage systems. The d i s t r i b u t i o n of workable periods of.varying length and the corresponding p r o b a b i l i t i e s of their occurrence for a given drainage system are shown in Figure 7. This can aid designers to develop drainage design procedures that are relevant to the requirements of individual farmers. They can make recommendations on a drainage design depending upon the number 104 k f (cm/day) 10 0.03 30 0.04 50 0.045 75 0.05 100 0.06 6 8 10 12 14 16 18 20 22 24 26 28 30 Drain spacing (m ) Probability of obtaining at least one workable P e ^ j " M a r c h Jj!.* function of drain spacing for some lowland soils of the Fraser Valley. 105 k = 10 cm/day f = 0.03 1 2 3 4 5 6 7 8 9 10 11 12 Consecutive workable days re 7. Relationship between the number of consecutive workable days and their probability of occurrence in March for a fine textured soil. Drain spacing 5.0 m 6.0 m 7.5 m — — 10.0 m 13.0 m 106 of consecutive workable days required for a given crop and the nature of farm implements. A family of such graphs should also enable them to recommend a drainage system depending upon the risk a farmer i s ready to accept in his farming enterprise. Output from the repeated use of the stochastic model for d i f f e r e n t designs can provide a useful quantitative index of the drainage system performance. If some monetary value can be assigned to the timeliness of farming operations, a benefit-cost analysis could be done to a r r i v e at an optimal spacing. However, e x p l i c i t estimation of t h i s loss function for l o c a l conditions i s d i f f i c u l t to obtain. In the meantime the model provides a meaningful index of the performance of possible drainage designs for a given s i t u a t i o n . It can be used in a tradeoff against cost to decide which design is most suitable under the circumstances in question. 1 07 CHAPTER IV PARAMETER UNCERTAINTY IN DRAINAGE DESIGN A major dilemma facing p r a c t i c i n g engineers these days i s how to incorporate v a r i a b i l i t y of s o i l properties in p r a c t i c a l designs. Among the various s o i l parameters in drainage design, characterization of saturated hydraulic conductivity i s indispensable. However, i t i s well known fact that t h i s property exhibits a large degree of s p a t i a l and perhaps temporal variation (Freeze, 1975; de Vries, 1982). Spacing of drains is highly dependent on s o i l hydraulic conductivity (proportional to the square root of hydraulic conductivity - Hooghoudt's formula). U n t i l now, drainage designs have been performed mostly on a deterministic framework by using equivalent s o i l parameters. There has been an on-going discussion among a g r i c u l t u r a l engineers, hydrogeologists and s o i l p h y s i c i s t s with regard to the meaning of th i s equivalence in groundwater flow through porous media. The concept holds a lot of promise in actual .designs because r e l a t i v e l y fewer data is required in such analyses. The characterization of the v a r i a b i l i t y of s o i l parameters is very important in drainage designs. The v a r i a b i l i t y of s o i l parameters, among others, led the International Drainage Workshop held at Wageningen, The Netherlands in 1978 to conclude that the exis t i n g drainage theory can be expected to give an order of magnitude difference in the required drainage system 108 design for a location. Methods of analysing uncertainty in drainage design are discussed below for both steady and unsteady saturated flow analyses. 4.1 Methods of Analysing Uncertainty Many e f f o r t s have been made in the l a s t decade to incorporate v a r i a b i l i t y of s o i l properties in deterministic analyses (Cornell, 1972; Freeze, 1975; Peck et a l . ,1973; Warrick et a l . ,1977; Bakr et a l ,1978; Delhomme, 1978; Sagar, 1978a and 1978b; Dagan, 1979, Smith and Freeze, 1979; P h i l i p , 1980; Tung and Mays, 1980; Cordova and Bras, 1981; Russo and Bresler, 1981). Most of the above works seem to have more academic than immediate p r a c t i c a l value. This i s due to either large data requirements which are seldom, i f ever, available in practice or excessive computational e f f o r t s required as prerequisites to the solution of the problem. Some of the methods require knowledge of advanced mathematics as an essential prerequisite and therefore, they may not be p r a c t i c a l to employ. There are others that can be applied to p r a c t i c a l designs quite readily without most of the above constraints, but only at the cost of some accuracy. However, such a trade-off c e r t a i n l y allows stochastic modeling to be used in p r a c t i c a l designs where i t was not feasible before; or a complete analysis of the problem can now be done, with obvious merits, instead of a p a r t i a l analysis. There are very few stochastic models for designing drainage systems that can cope with the v a r i a b i l i t y of s o i l properties. 109 Most of these models have either p r a c t i c a l constraints with respect t o data requirements and computational e f f o r t s or they can perform only a p a r t i a l analysis of the problem (that i s , only few parameter uncertainties can be considered simultaneously). This often leads to a model which can either simulate the response to weather only (Chieng, 1975) or one which can deal with parameter uncertainties only (Sagar, 1978b). The composite response of the system can not be obtained when both the above uncertainties are present at the same time. More often, only a few parameters are treated as random variables at one time. This is partly due to the n o n - f e a s i b i l i t y of an a l y t i c a l solutions to those cases where many parameters are random simultaneously. In general, an approximate solution of the whole problem i s better than an exact solution of only a portion of i t . Further, i t seems that models that can incorporate v a r i a b i l i t y of equivalent parameters have far greater immediate p r a c t i c a l u t i l i t y as compared to models that tend to incorporate micro-level heterogeneities. To thi s end, the use of f i r s t and second order analyses to incorporate parameter uncertainty in drainage design based upon both steady saturated and nonsteady saturated analyses w i l l be demonstrated. 4.1.1 F i r s t and second order analysis F i r s t and second order analyses require information about the f i r s t two moments of the random variable(s) only (mean and variance-covariance) and can provide information about the same moments for the solution to the problem. It i s often very 1 10 d i f f i c u l t to obtain f u l l p r o b a b i l i t y d i s t r i b u t i o n s of the random variables as i s required in more exact models, and i t i s comparatively easier to estimate their f i r s t two moments. However, i t may be noted that only a normal d i s t r i b u t i o n can be completely defined by these moments. If a random variable i s not normally d i s t r i b u t e d , the method w i l l result in approximate solutions only. Cornell (1972) observed that i f the c o e f f i c i e n t of v a r i a t i o n remains small the method w i l l provide accurate r e s u l t s . It must be noted that f i r s t order analysis i s based on the mean and variance-covariance of a random function as calculated from the f i r s t order Taylor series expansion. Second order analysis employs the mean based upon second order Taylor series expansion but the variance-covariance i s s t i l l obtained from the f i r s t order Taylor series expansion. Therefore, the mean derived in these analyses w i l l be d i f f e r e n t but the variance-covariance w i l l not. Let us start with a random function Y functionally related to another random function X, that i s , Y = f(X) If the above function behaves s u f f i c i e n t l y well, and provided the c o e f f i c i e n t of variation i s not large, the function can be expanded by Taylor series about the expected value of X as 111 Y = f(X() + f'(X-X|) + 1/2 f ' ( X - X | ) 2 + ... where f and f ' are the f i r s t and second derivatives of f with respect to X, evaluted at X|, the expected value of X. Neglecting second and higher order terms as being small compared to the f i r s t two terms and taking expected values on both sides, we get E[Y] =1 E[f(X|) + f'(X-X|)] where =1 stands for the f i r s t order approximation. E[ ] denotes the expected value operator. Since expected value i s a linear operator, we have E[Y] =1 f(X|) + f'(E(X) - X|} = ' f (X | ) (since E[X] =' X | ) The second moment can be estimated as Var(Y) = or Var(Y) = E[{f(X) - f(X|)} 2] E[{f(X|) + f'(X-X|) E[{f'(X-X|)} 2] [ f ' ] 2 E[(X-X|) 2] [ f ' ] 2 Var(X) - f(x|)} 2] S i m i l a r l y , a second order analysis of the mean can be obtained using a second order Taylor series expansion, that i s 1 12 Y = 2 f(X|) + f'(X-X|) +l/2f " (X-X|) 2 where =2 stands for second order approximation. Again taking expectations on both sides, we get E[Y] =2 f(X|) + l / 2 f " E[(X-X|) 2] = 2 f(X|) + l/2f' ' Var(X) The above estimate of the mean i s more accurate beacuse i t i s conditional on the mean-and variance of X (Dettinger and Wilson, 1981). Similar approximations can be made in multivariate sit u a t i o n s . A multidimensional Taylor series expansion i s used in this case. Considering a function Y = f(U,V) Expanding i t by Taylor series Y = f(U|,V|) + [{f(U).(U-U|)} + {£*(V).(V-V|)}] + l/2[{f " (U).(U-U| ) 2} + {f' ' (V).(V-V| ) 2}] + Where f'(U) and f'(V) are the f i r s t derivatives of f with respect to U and V evaluated at the mean value of U and V respectively. S i m i l a r l y , f''(U) and f''(V) are the second 1 13 derivatives of f with respect to U and V evaluated at the mean values of U and V respectively. Taking expectations on both sides E[Y] = 2 f(U|,V|) + l/2[f " (U).Var(U) + f " (V).Var(V)] S i m i l a r l y , variance can be calculated as Var(Y) = 2 [f'(U)] 2.Var(U) + [f'(V)] 2.Var(V) If the variables U and V are correlated, we have E[Y] =2 f(U|,V|) •+ l / 2 [ f ' ' (U,U) .Cov(U,U) + 2f''(U,V).Cov(U,V) + f'*(V,V).Cov(V,V)] and Var(Y) =2 [f' (U)] 2.Cov(U,U) + 2 . [f* (U).f' (V)].Cov(U,V) +[f'(V)] 2.Cov(V,V) Where f'(U,U), f'(V,V) and f'(U,V) are the second derivatives with respect to U only, V only and U and V respectively evaluated at the mean values of the variables. S i m i l a r l y Cov(U,U) and Cov(V,V) are the variances of variables U and V respectively. Cov(U,V) i s the covariance between the random, variables U and V. Generalizations can be made ea s i l y for multivariate 1 1 4 situat i o n s . Therefore, using the above simple relationships, the f i r s t or second order estimates of the moments of the dependent variable can be made ea s i l y once these moments are known for independent random variables. In the rest of the chapter, f i r s t and second order analysis w i l l be applied to drainage theory. Examples w i l l be given for both steady and nonsteady state formulas. 4.1.1.1 Steady state theory Most of the present-day drainage designs are based on steady state theory (Hooghoudt's method). It i s generally written as S 2 = [4 K (2 d h + h 2) ] / R Where S i s the drain spacing, K i s the saturated hydraulic conductivity, R i s the drainage c o e f f i c i e n t , h i s the permissible midspan water table height above the drains and d is the equivalent depth to the impermeable layer. For s i m p l i c i t y only, assuming that most of the flow takes place below the drains, we have S 2 = (8 K h d) / R or h = (S 2 R) / (8 K d) In the above equation, h is the dependent variable and we can investigate the uncertainty in h due to uncertainty in other parameters for a given drain spacing. The above transformation 1 1 5 has been performed only for the sake of convenience. A change in the dependent variable does not a l t e r the analysis. The parameter h i s usually determined from both the workability requirements of the s o i l and the crop water requirements. Other parameters such as K, d, and R are assigned as random variables and any uncertainty in characterizing them w i l l result in an uncertain drainage system performance measured by the parameter h for a given drain spacing. Parameters K and R may exhibit weak co r r e l a t i o n between them because in most cases, a high K value i s matched with a high R value and vice-versa. F i r s t and second order analyses w i l l be applied to the above equation for the following cases : (1) K as random variable only; (2) R as random variable only; (3) both K and R as random variables and; (4) d as random variable only. Case (3) w i l l be treated for both uncorrelated and correlated K and R. 1. Random hydraulic conductivity In Hooghoudt's formula, drain spacing varies d i r e c t l y as the square root of hydraulic conductivity.- It i s noted here that the hydraulic conductivity that we are interested in i s some equivalent K of the s o i l which i s not e x p l i c i t l y known for obvious reasons. A p r a c t i c a l way of dealing with the problem may be to try to estimate i t s f i r s t two moments by whatever means that are available. It may include both the auger hole tests and subjective judgements. The estimation of the moment should be easier than the estimation of the complete pr o b a b i l i t y 1 16 d i s t r i b u t i o n . It i s suggested that the mean and the variance of hydraulic conductivity can be estimated from a low (L), probable (P) and high (H) estimates of the s o i l property. These estimates can be obtained either by quizzing an experienced designer in the area or from few auger hole tests or both. A low value may stand for a 90% pr o b a b i l i t y that the actual value of the parameter can not be less than a prescribed value (a general convention). A high value may represent a 90% pr o b a b i l i t y that the actual value can not be higher than the given value, again a general convention. A probable value is the value seems most l i k e l y to represent the actual value. Knowing L,P and H estimates of the property, i t s mean and variance can be estimated using formulas given by Benjamin and Cornell (1970). However, in estimating these moments from L,P and H estimates, some assumption needs to be made with regard to the prob a b i l i t y d i s t r i b u t i o n of the L,P and H estimates. Selection of the d i s t r i b u t i o n to be used to obtain the mean and the variance of the property i s not f u l l y understood at this time and needs some future research. The normal or i t s s i s t e r d i s t r i b u t i o n s may be used but a better i d e n t i f i c a t i o n of the underlying d i s t r i b u t i o n when people make guesses from few known facts and i n t u i t i o n should be studied. There are some guidelines which can be observed while obtaining these estimates from f i e l d personnel (Russell, 1981). However, i t may be noted that the f i r s t and second order analyses as such do not depend on the d i s t r i b u t i o n used in a r r i v i n g at the estimates of moments. What i s required i s a good estimate of 1 1 7 these moments. Of course, better estimates w i l l lead to more r e l i a b l e r e s u l t s . Rewriting Hooghoudt's formula for the case when most of the flow takes place below the drains, h = (S 2 R) / ( 8 K d) and r e c a l l i n g that the expected value of the dependent variable in the f i r s t order analysis is E[Y] =2 f(X|) + 1/2 f''.Var(X) and variance i s Var(Y) =1 V a r ( X ) . [ f ' ] 2 Therefore, we need to calculate the f i r s t and the second derivatives of the function with respect to the random variable (K, in our case), that i s h' = (S 2 R) / (8 d).(-1).K- 2 h'' = (S 2 R) / (8d).(-1).(-2).K" 3 Let us assume a drainage system design problem where h=60 cm, d=200 cm, R=0.903 cm/day and K estimates are L=20, P=40 and H=100 cm/day. Assuming reciprocal normal d i s t r i b u t i o n (that i s the reciprocals are normally distributed) for K, we obtain 118 estimates of mean and varianace as 49.16 cm/day and 549.32 respectively. It should be noted that the writer has no special reason for using the inverse normal d i s t r i b u t i o n for the K estimates. Some assumption i s needed with regard to the d i s t r i b u t i o n in order to estimate the f i r s t two moments. It i s i hoped that future research in thi s area w i l l be able to provide some ra t i o n a l answer to t h i s problem. Using Hooghoudt's formula (with equivalent K of 49.16 cm/day), a drain spacing of 22 m w i l l be required for the drainage system. Using i t as the design drain spacing, we w i l l investigate the uncertainty in h due to uncertain equivalent K in the drainage system design. The expected value of the permissible midspan water table height can be obtained as E[h] = [(S 2 R) / (8 K d)]| +1/2[(S 2 R) / (8 d).(-1).(-2).K- 3]|.Var(K) or =55.56 + 12.67 =68.19 cm Had we used a single equivalent hydraulic conductivity value of 49.16 cm/day, h would be h = 55.56 cm Therefore, in an uncertain situation where equivalent K i s not completely known, the expected value of the midspan water table height i s 68.19 cm. The variance in h can be estimated as 1 19 Var(h) = [(S 2 R) / (8 d).(-1).K~ 2] 2|.Var(K) = 701 .8 or the standard deviation in h i s S.D. = 26.49 cm The variance or standard deviation may be regarded as the risk in drainage system design due to uncertainty about hydraulic conductivity. For a normal d i s t r i b u t i o n , i t implies that there i s a 90% pro b a b i l i t y that h would not exceed 102.1 cm (mean + 1.28*S.D.). . S i m i l a r l y , there i s a 90% pr o b a b i l i t y that h would exceed 34.3 cm (mean - 1.28*S.D.). 2. Random drainage c o e f f i c i e n t The drainage c o e f f i c i e n t (R) i s the amount of water that must be removed by a drainage system in one day. It depends upon the climate of the region, s o i l type, crop water requirements and f i e l d machinery operations (Chieng, 1978). In short, R i s a measure of the hydraulic loading on the drainage system. Although advanced methods of estimating R for a given location are available (Chieng et a l . , 1978) i t i s l i k e l y that the estimates w i l l be uncertain due to p r a c t i c a l l i m i t a t i o n s . F i r s t and second order analyses can be performed to estimate uncertainty in h due to uncertainty in R. 1 20 Rewriting the reduced form of Hooghoudt's equation h = (S 2 R) / (8 K d) The p a r t i a l derivatives required for the uncertainty analysis are h' = S 2 / (8 K d) h 1 1 = 0 Estimating a L,P and H values of R as 0.5, 0.8 and 1.5 cm/day and again assuming reciprocal normal d i s t r i b u t i o n for R estimates, we obtain a mean value of 0.903 cm/day and variance of 0.108. Uncertainty analysis can now be applied as E[h] = {(S 2 R) / (8 K d)}| + 1/2(0}.Var(R) =55.56 cm and Var(h) = {S2 / (8 K d)} 2|.Var(R) = 408.9 or S.D. (h) = 20.22 cm Therefore, for a normal d i s t r i b u t i o n the 90% prob a b i l i t y bounds for ' the high and low values of h are 75.78 and 35.34 cm 121 respectively. 3. Random drainage c o e f f i c i e n t and hydraulic conductivity Both R and K are known to vary in p r a c t i c a l situations. Therefore, such a case w i l l be studied in t h i s section. The analysis can be performed in two ways (a) R and K are uncorrelated and (b) R and K are correlated. Each case w i l l be dealt with separately in the following sections. Rewriting Hooghoudt's formula as h = (S 2 R) / (8 K d) We w i l l need the following derivatives for the analysis h'(R) = S 2 / (8 K d) h''(R) =0 h' (K) = {(S 2 R / (8 d).(-1).K- 2} h' '(K) = {S2 R) / (8 d).(-1).(-2).K~ 3} h"(R,K) = h"(K,R) = [{S 2 / (8 d)}.(-l).K" 2] Now, we w i l l deal with each case separately. (a) Independent R and K or E[h] = {(S 2 R) / (8 K d)}| + 1/2 [0 + {(S 2 R / (8 d)}.(-1).(-2).K- 3]|.Var(K) E[h] = 68.20 cm 1 22 and Var(h) [S 2 / (8 K d)} 2|.Var(R) + {(S 2 R) / (8 d).(-1).K- 2} 2|.Var(K)] or Var(h) 1111 or S.D.(h) 33.33 cm Therefore, for a normal d i s t r i b u t i o n the 90% confidence bounds for high and low values of h are 110.9 and 25.54 cm respectively. (b) Correlated R and K In nature, R and K are correlated to some extent because in most locations where R i s high, K i s high and vice-versa. Although no data i s available to quantify the above relationship, a correlation c o e f f i c i e n t of 0.3 is assumed in the present example to i l l u s t r a t e the use of uncertainty analysis in such cases. E[h] = [ (S 2 R) / (8 K d)]| + 1/2[{(S 2 R) / (8 d).(-1).(-2).K- 3}|.Var(K) + 2{(S 2 / (8 d).(-1).K" 2}|.Var(K,R) + 0] 71.08 cm and Var(h) = [{S 2 / (8 K d)} 2|.Var(K) + 2{S 2 / (8 K d).(S 2R / 8 d).(-1).K- 2}| 1 23 .Cov(R,K) + {(S 2 R) / (8 d).(-1).K" 2}|.Var(K)] = 789.6 or S.D.(h) = 28.10 cm Again, had the d i s t r i b u t i o n been normal, the 90% confidence l i m i t s for high and low values are 107.05 and 42.98 cm respectively. It should be noted that there i s a reduction in the variance in permissible midspan water table height due to consideration of correlation between R and K as compared to the case when R and K were uncorrelated. 4. Random depth to the impermeable layer The depth to the impermeable layer i s seldom known in most p r a c t i c a l designs. It i s estimated either from experience in adjacent areas or by digging a hole up to a certain depth at few locations in the area. Assuming a L,P and H estimates of 150, 200 and 300 cm for d and using reciprocal normal d i s t r i b u t i o n , we obtain a mean value of 213.56 cm and a variance of 2712.67 cm. Rewriting Hooghoudt's formula as h = (S 2 R) / (8 K d) Various derivatives required for the analysis are h' = (S 2 R) / (8 K).(-1).d- 2 h " = (S 2 R) / (8 K ) . ( - l ) . ( - 2 ) . d - 3 1 24 Applying the uncertainty analysis E[h] = (S 2 R) / (8 K d) | + 1/2{(S2 R) / (8 K).(-1).(-2).d- 3}|.Var(d) or = 55.14 cm . Var(h) = {(S 2 R) / (8 K).(-1).d~ 2} 2|.Var(d) = 197.5 or S.D.(h) = 14.05 cm For a normal d i s t r i b u t i o n , the 90% confidence bounds for high and low values of h are 73.14 and 37.15 cm respectively. It appears from the above results that d does not affe c t h s i g n i f i c a n t l y . It is because the effect of d on the drainage system design decreases as i t s depth below the drains increases. The greatest effect of d on the design w i l l be observed when i t i s closer to the drains. But in that case, i t may not be a random variable because i t s location w i l l be known. Results from the above uncertainty analysis are shown in Table 9. The respective c o e f f i c i e n t s of va r i a t i o n of the uncertain parameters as well as those of results are also shown. Though these c o e f f i c i e n t s appear to be somewhat higher (and that may lead to some error) the importance of taking into 125 TABLE 9 RESULTS FROM THE PARAMETER UNCERTAINTY ANALYSIS OF HOOGHOUDT'S EQUATION INPUT OUTPUT (he ight o f the water t a b l e above the d r a i n cent re ) (cm) Random V a r i a b l e Mean Value Coef f i c i ent o f V a r i a t i o n Mean Value C o e f f i c i ent o f V a r i a t i o n H y d r a u l i c c o n d u c t i v i t y 49.16 cm/day 0.48 68.19 0.39 Drainage c o e f f i c i e n t 0.90 cm/day 0.36 55.56 0.36 H y d r a u l i c c o n d u c t i v i t y 6 d ra inage c o e f f i c i e n t a) u n c o r r e l a t e d I I 11 68.20 0.49 b) c o r r e l a t e d I I 11 71.08 0.40 Depth o f the impermeable l a y e r below the d r a i n cen t re 213.56 cm 0.24 55.14 0.25 Determi ni s t i c S o l u t ion 49.16 cm/day 0.90 cm/day - 55.56 -126 consideration the stochastic nature of input data in drainage design i s apparent. 4.1.1.2 Nonsteady state drainage design In f i e l d s i t u a t ions, steady state flows are seldom encountered. Most of the time, we are dealing with fluctuating water tables. Therefore, in t h i s section we w i l l apply the f i r s t and second order analyses of analysing input uncertainty to such designs. For saturated flows and for the case where the D-F assumptions are applicable, the Boussinesq equation i s generally solved for the drainage initial-boundary value problem. It can be written as f 3h t K _ L _ (h 111) 3t dx dx where f i s the drainable porosity, h i s the height of the water table above an impermeable layer, K i s the saturated hydraulic conductivity, X i s the horizontal distance and t i s time. This i s a nonlinear equation and i s quite d i f f i c u l t to solve in closed form. Therefore, i t has been l i n e a r i z e d quite frequently (Polubarinova-Kochina, 1962) to bring h outside the parenthesis on the right side of the equation. It i s noted here that there are a n a l y t i c a l solutions available to the nonlinear Boussinesq equation for the case when drains l i e on the 1 27 impermeable layer. However, there i s only one such solution available when drains do not rest on the impermeable layer (Van Schilfgaarde, 1974). This method of solution i s applicable only i f the water table recession can be considered as a sequence of small steps. Therefore, numerical solutions have been obtained recently by solving the Boussinesq equation on the computer (Moody, 1966; Skaggs, 1975). These solutions have very few constraints in comparison to the "analytical solutions. The f i r s t and second order analyses w i l l be applied to both solution methodologies in the following pages. 1. A n a l y t i c a l solution to the drainage design problem The a n a l y t i c a l solution to the Boussinesq equation when drains do not rest on the impermeable layer and when water table recessions can be considered as a sequence of small steps is given by Van Schilfgaarde (1964) as t = (f S 2) / (9.K d) ln [{ml (2 d + m2)} / {m2 (2 d + ml)}] Where t i s the time taken by the water table to recede from an i n i t i a l height above drains of ml to a new height m2, f i s the drainable porosity, K i s the saturated hydraulic conductivity, d is the equivalent depth to the impermeable layer and S i s the drain spacing. Let us solve the example that was discussed in the last section, that i s , S = 22 m, d = 200 cm, mean K = 49.16 cm/day 1 28 and var(K) = 549.32, f = .04, ml = 120 cm and m2 = 80 cm. The following derivatives w i l l be required for the uncertainty analysis t' = ( - D . ( f S 2) / (9 d).ln [{ml (2 d + m2)} / {m2 (2 d + ml)}]. K'2 t " = (-1 ) . (-2) . (f S 2) / (9 d).ln [{ml (2 d + m2)} / {m2 (2 d + ml)}]. K-3 Applying the uncertainty analysis E f t ] = (f S 2) / (9 K d).ln [{ml(2 d + m2)} / {m2 (2 d + ml)}3| + 1/2[(-1).(-2).(f S 2) / (9 d).ln {ml (2 d + m2)} / {m2 (2 d + ml)} . K- 3]|. Var(K) = 0.87 days and Var(t) = Var(K).[(-1 ) . ( f S 2) / (9 d) •In {ml (2 d + m2)} / {m2 (2d + ml)}. K" 2] = 0. 1153 or S.D.(t) = 0.34 day Therefore, for normal d i s t r i b u t i o n the 90% confidence l i m i t s for high and low values of t are 1.31 days and 0.44 day respectively. 129 2. Numerical solution to the drainage design problem F i r s t and second order analyses can be incorporated in numerical schemes provided the f i r s t and second derivatives of the dependent variable can be expressed in terms of the random independent variables. Dettinger and Wilson (1981) have proposed a procedure that can incorporate f i r s t and second order uncertainty analyses in groundwater flow problems by employing the v e c t o r i a l expansion of Taylor series as suggested by Vetter(l970, 1971 and 1973). However, the analysis gets very complicated and may not be p r a c t i c a l for design purposes. The method of analysing uncertainty can be applied very e a s i l y to the numerical models of drainage system design (or any other design) by estimating the required derivatives numerically. The model can be run i t e r a t i v e l y around the mean value of the random variables and thus the derivatives can be approximated. Of course, i t w i l l add to the computational burden in solving the model, but the increased costs may be j u s t i f i e d considering the costs of running other available stochastic models (Freeze, 1975) or the knowledge of advanced mathematics required in certain other models (Bakr, 1976; Sagar, 1978a and 1978b). The Boussinesq equation was solved numerically in the present work using an alternate predictor-corrector method which is a very e f f i c i e n t method (no i t e r a t i o n s are required to obtain the solution at the next time step). In the equation, h i s the dependent variable and only K i s assumed to be a random variable in the present analysis. However, any other input variable or 1 30 variables can be considered simultaneously without a f f e c t i n g the methodology. The various steps involved in the uncertainty analysis are as follows: 1. The governing p a r t i a l d i f f e r e n t i a l equation (Boussinesq equation) is solved using the mean value of saturated hydraulic conductivity. 2. The equation i s solved again using a hydraulic conductivity value of (K + DK), where DK is a very small deviation in K around the mean of K. 3. The equation i s resolved using a hydraulic conductivity value of (K - DK). 4. From the above, the required derivatives can be estimated numerically from the d e f i n i t i o n of a derivative, that is h' = {h(K + DK) - h(K - DK)} / 2DK h'' = {h(K + DK) - 2 h(K) + h(K - DK)} / DK2 5. The expected value as well as the variance in h can now be estimated using the f i r s t and second order analyses. The example given in the last section was solved using the numerical techniques. The results are shown in Figures 8, 9 and 10. In Figure 8, the expected value of the midspan water table expected value T 1 r 1 r 5 6 7 8 9 t, in days Figure 8. Water table drawdown at the mid-spacing computed from the Boussinesq equation. 100 -90-80-70 60 H 50 40-30-20 j 10 J G t, in days 6 7 - T 8 T 9 —t 10 Figure 9. Growth of uncertainty in the output as a function of time. Soil surface 120 y xxxx 100 H 80 H 60-^ 40 20 ^ Expected value Expected value ± one standard deviation Expected value ± two standard deviation XXXXs / / o / / y _ o \ 0.1 o — — o -"ST 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Nondimensional horizontal distance (x/s) Figure 10. Standard deviation of water table height above the drain centre line (two days after the beginning of drawdown). 134 height above drains is plotted as a function of time. Results from a deterministic run of the model using the mean saturated hydraulic conductivity are also shown in the figure. Both results are i d e n t i c a l at t=0 because the i n i t i a l condition i s assumed to be completely known. As time progresses, the expected values diverge from the deterministic values and a maximum difference of 7.5 cm i s observed after 4 days. Results are also shown for +1 and -1 standard deviation in the figure. These may be regarded as approximately 90% confidence l i m i t s . Figure 9 shows a plot of standard deviation of midspan water table height as a percentage of i t s expected value versus time. The curve starts with zero standard deviation at zero time because the i n i t i a l condition i s assumed to be completely known. At t=4 days, i t i s found to be 66% of the expected value. Thereafter, the standard deviation increases very slowly while the expected value decreases. Since as time increases, the midspan water table height above the drains tends to approach the drain l e v e l , which i s completely deterministic, both the expected value and standard deviation decrease. The standard deviation i s not only a function of time but is also a function of distance from the drain axis. This i s shown in Figure 10. The standard deviation is zero at the drain axis because of the boundary condition h=d for a l l times greater than or equal to zero. It increases as we move toward the midspan. Again, curves are drawn in Figure 10 for +1, +2, -1 and -2 standard deviations in the figure. The two-standard deviation p r o f i l e s may be taken as approximately 95% confidence 1 35 l i m i t s . Therefore, i f the design c r i t e r i o n i s based on the f a l l i n g water table concept, that i s for example, the water table should drop from the s o i l surface at least 60 cm in two days, the assumed spacing of 22 m w i l l be inadequate. Although the expected value of the midspan water table height meets t h i s requirement, the prob a b i l i t y of i t being exceeded i s not negl i g i b l e and therefore there i s some ris k involved in such a drainage system design. The procedure may be repeated to arr i v e at a spacing that meets the design c r i t e r i o n with an allowable r i s k . It may be noted that the standard deviation remains f i n i t e at a l l times. In drainage design, i t s magnitude may be regarded as a measure of the risk involved in selecting that design. Knowledge of the risk due to uncertain input data should be of considerable help in a r r i v i n g at a design that meets the requirements of a pa r t i c u l a r farmer. In the above sections, hydraulic conductivity has been treated as the single uncertain parameter in the drainage design. It was • selected because i t i s the most important parameter in the drainage design. However, i t could be any other input parameter such as r a i n f a l l , drainable porosity, depth to the impermeable layer. Any dependencies among input parameters can also be considered. 1 36 CHAPTER V LINKING PARAMETER AND CLIMATIC UNCERTAINTIES IN DRAINAGE SYSTEM DESIGN It should be possible to incorporate uncertainty due to both s o i l parameters and climate in a r a t i o n a l drainage design procedure. However, the writer i s not aware of any ex i s t i n g models that consider these uncertainties in drainage system designs. In this chapter, an attempt i s made to link the climatic uncertainty model developed in chapter III with the parameter uncertainty model developed in chapter IV and t h i s i s then used to estimate confidence l i m i t s for the pr o b a b i l i t y of obtaining at least one workable period in the month of March. The confidence l i m i t s are measures of the uncertainty in the solution due to uncertain s o i l parameters. For s i m p l i c i t y , only saturated hydraulic conductivity w i l l be considered as the random s o i l property. However, the procedure i s amenable to incorporate uncertainty in any other s o i l parameters or combinations thereof. 5.1 The Stochastic Drainage Model In t h i s section, an attempt i s made to link the uncertainty in s o i l parameters with the uncertainty in climate in a rat i o n a l way. The procedure w i l l be explained with the help of an example. The following drainage system design i s considered: d = 100 cm D = 120 cm 1 37 L P H K = 6.25 8.33 10.42 cm/hr f = 0.07 S = 24 m ( f i r s t t r i a l ) Various steps in the design are as follows: 1. From the L, P, and- H estimates of the saturated hydraulic conductivity, i t s f i r s t two moments are estimated, namely the expected value and the variance. As discussed e a r l i e r , some suitable assumption is required with regard to the probability d i s t r i b u t i o n of the estimates. For the above example, i t was assumed that the reciprocals of the estimates are normally d i s t r i b u t e d . Any other d i s t r i b u t i o n that may seem suitable can be used. 2. Using the mean and the variance of saturated hydraulic conductivity, the drainage parameter uncertainty model as described in chapter IV i s run to get estimates of the mean and the variance of the d a i l y r a i n f a l l amount that can cause the water table to move from one state to another. The state i n t e r v a l used in the example i s 30 cm. It i s somewhat on the high side of the recommended int e r v a l size but considering the high saturated hydraulic 138 conductivity of the s o i l , i t appears to be v a l i d . However, a smaller interval could have been used without a f f e c t i n g the methodology. The causes of uncertainty in net r a i n f a l l are a " r e f l e c t i o n " of the uncertainty in the s o i l parameter(s). 3. Matrices can be constructed now that can provide information with respect to the net r a i n f a l l required to move the water table from one state to another. Three such matrices were constructed, one -each for the expected value of the net r a i n f a l l , the expected value plus two standard deviations and the expected value minus two standard deviations. The plus/minus two standard deviations around the expected value can be considered approximately as 95% confidence l i m i t s . The results are shown in Table 10. 4. With the help of the cumulative probability d i s t r i b u t i o n functions of net r a i n f a l l for winter and spring, six tr a n s i t i o n matrices were constructed ( one each for spring and winter for the three net r a i n f a l l matrices obtained from step 3). These matrices are shown in Table 11. 5. After checking the steady state features 139 TABLE 10 NET RAINFALL MATRICES a) Deterministic Solution Ol t l c 0 c £ c Ol V x XI Ol *—-V E x u X t l n 0 x> IB V u. *-> t l a L. C t l 41 ti *-i O X IB *J X c u -e ID 0 10 i . o . - D t l X 30 60 90 2:120 B0.0k 53.33 27.89 3.63 66.58 39-33 13.31 -11.55 53.81 26.03 -0.53 -25-94 41.67 13.39 -13.67 -39.62 »120 90 60 30 0 M i d s p a n water table height above the drain centre at the end of the period (cm) b) Expected Value Solution Ol tl c > — O C xt c IS — Oi tl JZ XI Ol — ti E ti x y x *« t l T3 — n o x — IB tl k. M L V fc» Q. L. C tl tl tl 4J U X c o o &-o in V tl — X X w 30 60 90 79.03 52.34 26.89 3.63 67.58 39.33 12.32 -10.56 52.82 26.05 -0.53 -25.9* 40.67 14.38 -14.67 -39-61 £120 90 60 30 Midspan water table height •bove the drain centre at the end of the period (cm) 140 cn II c > — £ i IC — cn - M V £ U X *> i) 4-* ^ — It O £ — IB tl 1-«J I- tl u c l i II V n c n o A -•o tl — X TAR1F 10 (cont'd) MET RAINFALL MATRICES c) Expected Value plus Two Standard Deviations 89.<t1 59-30 30.92 1 5.58 80.64 48.61 18.79 -6.53 ~~ 68.57 38.01 8.27 -19-60 ) — 59. n 28.9' -3.31 Ir. -30.82 n Midspan water table height above the drain centre at the end of the period (cm) d) Expected Value minus Two Standard Deviations cn v c > — i I IB — en +> ti X JD Ol - ~ — H E tl X o £. *< ~— II ^ — IB 0 -O — IB II L <v >- t> *> a k c ti ti ti w O X IB « C IB O 68.65 45.38 22.86 1.68 5«t.52 30.05 5.85 -14.59 ) • — 37.07 14.08 -9.32 -32.29 J • 22.24 -0.15 -26.02 -48.40 in 1 „,d»p.n water table heloht .bove the drain centre at the and of the period (cm) 141 TABLE 11 TRANSITION MATRICES - DETERMINISTIC SOLUTION a) Winter 0.000 0.040 0.115 0.845 0.015 0.060 0.525 0.400 0.040 0.140 0.820 0.000 0.065 0.858 0.350 0.000 Midspan water table height above, the drain centre at the end of the period (cm) b) Spring £120 0.000 0.015 0.075 0.910 0.005 0.035 0.360 0.600 0.015 0.080 0.905 0.000 0.035 0.415 0.550 0.000 £120 90 60 30 Midspan water table height above the drain centre at the end of the period (cm) 142 TABLE H (cont'd) TRANSITION MATRICES - EXPECTED VALUE SOLUTION c) Winter V c > — & i m — x J! 11 x u — • O X — W tl k - i 8. k c t l tl II U X • M B c te-e n o a k a -o •o ti — X 30 60 90 £120 £ 1 2 0 0.000 0.015 0.075 0.910 0.005 0.035 0.360 0.600 0.015 0.080 0.905 0.000 0.035 0.515 0.450 0.000 90 60 30 Midspan water table height above the drain centre at the end of the period (cm) d) Spring 30 60 90 £ 1 2 0 0.000 0.015 0.075 0.910 0.005 0.035 0.360 0.600 0.015 0.080 0.905 0.000 0.035 0.515 3.450 0.000 60 30 £ 1 2 0 90 Midspan water table height above the drain centre at the end of the period (era) 143 TABLE 11 (cont'd) TRANSITION MATRICES - EXPECTED VALUE PLUS TWO STANDARD DEVIATIONS e) Winter U c > —£ c IB • -O l ** 4> X x O l — V E V X u X tl fc> •o — IB 0 X — IB tl k w k tl +* a <- c tl tl ti 4J U X IB j c C "is o IB 1-a TJ l/l T l tl — X X *> 30 60 90 2120 0.000 0.030 0.005 0.885 0.005 0.045 0.300 0.650 0.015 0.070 0.915 0.000 0.035 0.155 0.810 0.000 Midspan water table height above the drain centre at the end of the period (cm) f) Spring O l tl c > —o c X C O l *J tl X X O l — ti E tl X o X tl " •o — IB o X IB II k <U k II a k C II tl u *> o X IB J C c ~n 0 VI •o t> — X X f 30 60 90 2120 0.000 0.010 0.060 0.930 0.000 0.025 0.185 0.790 0.005 0.045 0.950 0.000 0.010 0.090 0.900 j 0.000 60 30 2120 90 Midspan water table height above the drain centre at the end of the period (cm) 144 TABLE 1) (cont'd) TRANSITION. MATRICES - EXPECTED VALUE MINUS TWO STANDARD DEVIATIONS SOLUTION g) Winter 30 60 90 2120 0.010 0.040 0.150 0.800 0.035 0.080 0.885 0.000 0.075 0.405 0.520 0.000 0.220 0.780 0.000 0.000 2120 90 60 30 M i d s p a n w a t e r t a b l e h e i g h t a b o v e t h e d r a i n c e n t r e a t t h e e n d o f t h e p e r i o d ( c m ) h) Spring cn V c o c £ c 10 cn t l X -o cn _ t i E t i X Ll X ' t l *t n o X IS t i k. L. t l a. C V t l t i 4V u X n X c u-e te O • L. CL T> m T3 t i X X 30 60 90 2120 0.000 0.025 0.085 0.890 0.010 0.060 0.928 0.002 0.040 0.280 0.680 0.000 0.115 0.885 0.000 0.000 2120 90 60 30 Midspan water table getgh above the drain cenre at the end of the period (cm) 145 of the expected-value t r a n s i t i o n matrices, the procedure of Chapter III is employed to estimate the expected pr o b a b i l i t y of obtaining at least one workable period in March. 6. The procedure of step 5 i s repeated for the t r a n s i t i o n matrices corresponding to the expected value plus two standard deviations and the expected value minus two standard deviations. The p r o b a b i l i t i e s of obtaining at least one workable period thus obtained form the approximate 95% confidence l i m i t s for the solution. 7. The procedure of steps 1 to 6 i s repeated for other drain spacings. For the present example, d i f f e r e n t drain spacings used were 14 m, 24m, 28 m and 35 m. The results are shown in Table 12. For a drainage system design with drain spacing of 24 m, there i s a 72.5% pro b a b i l i t y that the design w i l l provide at least one workable period in March. However, due to our i n a b i l i t y to f u l l y map the uncertain nature of the saturated hydraulic conductivity of the s o i l , the approximate 95% confidence l i m i t s are 94% and 69% respectively. The deterministic solution to the problem i s also shown in the table. The expected value solution and the deterministic solution to the problem are nearly the same. However, benefits 146 TABLE 12 PERFORMANCE OF VARIOUS DRAINAGE SYSTEM DESIGNS IN THE LOWER FRASER VALLEY FROM A WORKABILITY POINT OF VIEW Drain Probability of obtaining at least on workable period in March {%) (m) Deterministic Expected Value Expected Value +2 S.D. Expected Value -2 S.D. ]k 100.00 100.00 100.00 97.50 2k 73.90 72.50 9^.00 69.00 28 30.^0 20.15 75.00 0.00 35 0.30 0.05 11.50 0.00 1 47 of the uncertainty analysis l i e in the estimation of the variance in the required p r o b a b i l i t y which i s not feasible in a deterministic analysis. The variance may be regarded as the risk in a given design and an awareness of i t s extent should help designers. It also leads us to a point where an indi v i d u a l i z e d drainage system design would be possible depending upon the amount of ri s k a p a r t i c u l a r farmer i s prepared to accept in the farming enterprise. The above results are also shown in Figure 11. The expected value solution to the problem along with the 95% confidence l i m i t s are shown in the figure. It appears that for the example under consideration, large benefits can be derived by decreasing the drain spacing from 30 m to 24 m( net gain in the pr o b a b i l i t y of obtaining at least one workable period in March i s 0.6). Thereafter, the increase in pr o b a b i l i t y per unit decrease in drain spacing diminishes considerably. It seems from the figure that there w i l l be no benefits to the farming enterprise i f drain spacing i s further reduced from 11 m. However, a guiding factor in the selection of a drainage system design i s i t s i n i t i a l i n s t a l l a t i o n cost. Selection of a system that f u l l y meets the design requirements such as a drain spacing of 11 m may not be economically f e a s i b l e . A trade-off should be made between what i s f u l l y e f f i c i e n t from the performance point of view and what is feasible economically. It i s suggested that tables l i k e Table 12 or figures l i k e Figure 11 can aid prac t i s i n g engineers to do designs better suited to p a r t i c u l a r situations than are possible at present. A benefit-cost 1 4 8 Expected value E[ J-2S .D. E[ ] + 2 S.D. Drain spacing (m) Figure 11. Probability of obtaining at least one workable period in March as a function of drain spacing. 1 49 analysis can be done to select an optimal drainage system design for a given location. 1 50 CHAPTER VI SUMMARY, CONCLUSIONS AND RECOMMENDATIONS In t h i s chapter, the proposed procedure for the design of a drainage/subirrigation system for the Lower, Fraser Valley i s summarized. Conclusions drawn from the study are l i s t e d and recommendations are given for future work. The present study f i r s t attempts to formally define the problem of drainage/subirrigation system design. On recognizing that i t was almost impossible to pursue t h i s l i n e of approach, i t was decided to follow a "branch and bound" approach where a sequence of small studies i s undertaken rather than a complete study and an attempt is made to simplify the design procedures by eliminating or regrouping those alternatives that are of secondary importance. A study was undertaken to check whether drainage or subirr i g a t i o n i s the l i m i t i n g factor in a drainage/subirrigation system design. Because of the inherent numerical complexities associated with the transient saturated/unsaturated flow models, i t was decided to apply a steady saturated/unsaturated flow analysis for the study. However, the unsaturated flow properties of the s o i l were required in the study and since they were not available for the l o c a l s o i l s , i t was decided to use the properties of some similar s o i l s obtained from the l i t e r a t u r e for the analysis. Both subi r r i g a t i o n systems and drainage systems were designed for a range of s o i l types on the basis of steady state flow theory using generally accepted c r i t e r i a . It was found that the design drainage spacing was 151 smaller than the design subirrigation spacing in a l l cases with medium to fine textured s o i l s except that the two designs approached each other as the s o i l texture became finer and the subi r r i g a t i o n requirements seemed to override drainage requirements for very fine textured s o i l s . Subirrigation, however, would not be used with very fine s o i l s . For most s o i l s , the height required to be maintained in the outlet to get the s u b i r r i g a t i o n system working was found to be quite feasible for use of the drainage systems for s u b i r r i g a t i o n . Thus, i t was established from the above simple study that for almost a l l p r a c t i c a l situations in the Valley, a design developed to meet drainage requirements should be more than s u f f i c i e n t to s a t i s f y s ubirrigation requirements. Therefore, i t is not necessary to consider sub i r r i g a t i o n when designing a dual-purpose drainage/subirrigation system in the Lower Fraser Valley. However, the l a t e r a l s should not be made too long to function properly for s u b i r r i g a t i o n . This w i l l be guided by the height to which water table should be raised in the outlet for s u b i r r i g a t i o n . Also, the location for water l e v e l control chambers and water f i l t e r s must be considered while designing a dual-purpose system to s a t i s f y drainage requirements only. Another problem that needs more attention from a g r i c u l t u r a l engineers i s " W i l l the fact that subirrigation w i l l keep s o i l around the drain pipes saturated for the hot part of the year, and therefore, w i l l prevent the s o i l from drying out and developing shrinkage cracks and improved s o i l structure have any 1 52 detrimental effects on the 'entry resistance' or 'exit resiatance' near the drain tubes?". Next a study was undertaken to formulate a p r a c t i c a l design c r i t e r i o n for designing drainage systems. The design requirements to be met by a drainage system during d i f f e r e n t seasons were investigated for l o c a l conditions. By deductive reasoning, i t was suggested that early spring requirements (workable conditions in March) override a l l other requirements. Winter requirements were not considered. Therefore, a drainage system that can meet early spring requirements should be more than s u f f i c i e n t to .provide an optimal growth environment in the summer and t r a f f i c a b l e conditions in the summer and in the f a l l . An advantage of the above study i s that we do not have to know the nature of any crop loss function, the single most d i f f i c u l t and hence usually unknown parameter in drainage design. The next step was to quantify the early spring requirements for the Valley. A d i s t i n c t i o n was made between a t r a f f i c a b l e and a workable s o i l . It was suggested that a drainage system should ensure workable conditions of the s o i l during seed-bed preparations and other t i l l a g e operations in March. An obvious question was "When does a s o i l become workable?". From the l i t e r a t u r e , i t was found that a saturated/unsaturated flow analysis could best provide the above information. However, due to the inherent d i f f i c u l t i e s in such analyses, i t was decided to use depth to the water table as an indicator of workability, following the work of Paul and deVries (1983) who established from a saturated/unsaturated flow analysis that the depth to the 153 water table from the s o i l surface may be used as an indicator of traf f i c a b i l i t y . The next step was to determine how many consecutive workable days would be required to perform the necessary f i e l d work in March. After consultations with various l o c a l experts, including farmers, i t was decided to use ten days as a s t a r t i n g value. However, recognizing that depth to the water table might not by i t s e l f be a f u l l y r e l i a b l e indicator of workability, a safety margin of two days was added to the ten days requirement. Therefore, the objective of a good drainage system was defined as providing at least one workable period of twelve consecutive days with water table more than 60 cm from the s o i l surface in March in most years of i t s operation (that i s , 80% or more pro b a b i l i t y of obtaining a workable period). The writer has not considered whether two work periods of six days in March could be just as s a t i s f a c t o r y as one work period of twelve days. Following the same argument, three work periods of fiv e days instead .of one work period of twelve days were also not considered. However, the proposed methodology can be e a s i l y modified to simulate such strategies. The next step was to estimate the p r o b a b i l i t i e s of obtaining at least one workable period in March from a given drainage system and for a given s i t e . For t h i s , a f i r s t order Markov chain model was proposed that uses t r a n s i t i o n p r o b a b i l i t i e s of water table moving from one state to another in order to a r r i v e at the p r o b a b i l i t y of obtaining at least one workable period in March. A s t a t i s t i c a l simulation was used to 1 54 get the required probability due to the complexity of d i r e c t l y computing pr o b a b i l i t y when the number of events i s greater than twelve. Another study was undertaken to study the e f f e c t s of parameter uncertainty in drainage design. After establishing that f u l l d i s t r i b u t i o n analyses such as using Monte Carlo methods are impractical for p r a c t i c a l design, f i r s t and second order methods of analyzing uncertainty were used in the study. Examples were given for both steady state and nonsteady state methods of designing drainage systems. A simple approach was suggested for performing such analyses on numerical models of groundwater flow by using numerical d i f f e r e n t i a t i o n . The approach was applied to the numerical model of drainage design using the Boussinesq equation and the nature of uncertainty due to uncertain s o i l parameters was studied. The next step was to formulate the design c r i t e r i o n for spring in the presence of s o i l parameter uncertainty. It was shown how the climatic uncertainty model can be linked with the s o i l parameter uncertainty model to refine the drainage system design. Confidence l i m i t s on the results are made possible by such a linkage, and these can help a designer to develop designs that w i l l meet the requirements of individual farmers. During the course of the study, a novel way of designing drainage systems on a steady state basis was discovered when both the s o i l hydraulic conductivity and the location of the impermeable layer are unknown. The method uses a p i l o t drainage system on which measurements are taken. Also, the writer became 1 55 aware of the importance of studying some of the possible enviornmental impacts of a g r i c u l t u r a l drainage such as peak flows, t o t a l runoff volume and water quality on the watershed hydrology. However, since the main thrust of the thesis was on the transient drainage/subirrigation system design under uncertainties due to climate and s o i l properties, the above topics are covered in the Appendices. 6.1 Conclusions 1. For the Lower Fraser Valley, a drainage/subirrigation system may be designed to meet drainage requirements only; subirrigation requirements w i l l be i m p l i c i t l y met in such a design. 2. For the Lower Fraser Valley, a drainage system designed to meet early spring requirements should be more than s u f f i c i e n t to meet other crop water management requirements during the crop growth and harvesting periods. Winter requirements were not considered. 3. The early spring requirements for the Lower Fraser Valley can be quantified by designing a system that can provide at least one workable period of any number of consecutive days (twelve days in the present study) in March in most years. 4. A f i r s t order Markov chain model that incorporates 156 t r a n s i t i o n p r o b a b i l i t i e s of the water table moving from one state to another can be used to estimate long-term p r o b a b i l i t i e s of the water table being in one of the various states. 5. For most s o i l s in the Valley, a state i n t e r v a l of 10 cm may be used while constructing t r a n s i t i o n matrices. 6. S t a t i s t i c a l simulation can be used to estimate the p r o b a b i l i t y of obtaining any number of consecutive days during which the water table stays in any desired range. 7. A comprehensive computer model was developed that can be used with any s o i l properties and weather data to estimate the p r o b a b i l i t y of obtaining at least one workable period of any number of consecutive days for a given drainage system. 8. Design graphs were obtained for some l o c a l s o i l s that can aid designers to develop drainage designs depending upon the risk a farmer i s ready to accept in the farming enterprise. 9. The above graphs underline the importance of a good s o i l and water management on the f i e l d . There i s a s i g n i f i c a n t gain in the probability of obtaining a workable period when the saturated hydraulic conductivity of the s o i l i s high. 157 10. The methodology can be used to obtain p r o b a b i l i t i e s for any number of consecutive workable days for a given drainage system. It can be used by designers to develop a drainage design that w i l l ensure a required number of workable days as guided by the crop requirements or the farm implements or both. 11. It appears that for fine textured s o i l s , the probability of obtaining a workable period increases s i g n i f i c a n t l y i f the drain spacing obtained from the steady state theory i s decreased by about 25%. Similar deductions may be made with respect to the drainage c o e f f i c i e n t or the depth to the water table at midspacing instead of drain spacing. 12. F i r s t and second order methods of analyzing uncertainty can be used to study the e f f e c t s of uncertainties in s o i l parameters on the drainage system design. 13. The method of analyzing uncertainty can be used for both a n a l y t i c a l formulations and numerical models of groundwater flow. 14. A design procedure for a drainage system can be formulated that incorporates uncertainties due to both weather and s o i l properties. 15. Confidence bounds can be established on the results due to 158 uncertain s o i l properties. These can aid designers develop drainage-system designs that w i l l meet the requirements of individual farmers (depending upon the ri s k they are ready to accept in the farm enterprise). A new s i m p l i f i e d methodology for designing drainage systems on the basis of steady state theory that can be used without any pre-knowledge about the saturated hydraulic conductivity of the s o i l and the location of the impermeable layer i s suggested in Appendix B. A g r i c u l t u r a l drainage can aff e c t watershed hydrology. There may be a positive c o r r e l a t i o n between the peak flows and the surface drainage a c t i v i t y in a basin. The nature of such a correlation between the peak flows and the subsurface drainage on a basin-wide scale i s not yet understood. This question i s examined in Appendix C. 6.2 Recommendations In order to translate the results of t h i s study into some monetary equivalents to which farmers can e a s i l y relate q u a n t i f i c a t i o n of the various benefits of a g r i c u l t u r a l drainage to the farm enterprise should be studied in d e t a i l . 159 2. The computer model developed in the study should be tested on d i f f e r e n t s o i l s and at di f f e r e n t locations to gain confidence in i t s predictions. Design curves given in the thesis should be checked for p r a c t i c a l applications. 3. The unsaturated flow properties for the l o c a l s o i l s should be estimated. They w i l l a i d greatly in future drainage/subirrigation designs. 4. The new methodology for designing drainage systems on the basis of steady state flow theory should be tested for lo c a l s o i l s . The method holds promise for f i e l d designs where the designer i s always confronted with the problem of quantifying the s o i l properties and the boundaries of the flow system. 5. E f f o r t s should be made to incorporate uncertainty due to s o i l parameters in design methods on a "macro basis"; methods that require a f u l l p r o b a b i l i t y d i s t r i b u t i o n of the random s o i l properties should be discouraged for p r a c t i c a l designs, because such data w i l l probably not be available in the near future for most s o i l s . Instead, e f f o r t s should be directed to search for some equivalent s o i l properties that are not known to the designer. An approach i s discussed in the present study that uses low, probable and high values of these properties. E f f o r t s should be made to refine the above analysis by estimating the prob a b i l i t y 1 60 d i s t r i b u t i o n s of equivalent s o i l properties. The environmental impacts of a g r i c u l t u r a l drainage should not be ruled out from future studies in s o i l and water engineering. It appears that watershed hydrology may be affected by these micro-level manipulations and steps should be taken by s o i l and water engineers to aid the personnel in environmental d i s c i p l i n e s to detect the nature of these e f f e c t s . 161 BIBLIOGRAPHY Anat, A., H. R. Duke and A. T. Corey. 1965. Steady upward flow from water tables. Colarado State University, Fort C o l l i n s , Colarado. Hydrology Paper No. 7:34 pp. Anon. 1972. B r i t i s h Columbia Drainage Guide. Anon., 1980. Fraser Valley I r r i g a t i o n Design C r i t e r i a Project. Soilcon Laboratories Ltd., Vancouver, B.C., Canada. Baehr, B. E. 1980. Personal Communication. A g r i c u l t u r a l Engineering Branch of the B.C. Ministry of Agriculture and Food, Abbotsford, B.C., Canada. Baier, W., G. W. Robertson and M. F. Clarke 1969. A Climatological Analysis of I r r i g a t i o n Requirements in the Lower Fraser Valley, B.C.. Publication 1318. Research Branch, Agriculture Canada, Ottawa, Canada. Bakr, A. A. 1976. Stochastic analysis of the effects of s p a t i a l variations in hydraulic conductivity on groundwater flow. Thesis submitted to New Mexico Inst i t u t e of Mining and Technology, Socorro, New Mexico, in p a r t i a l fulfilment of the requirements for the degree of Doctor of Philosophy. Bakr, A. A., L. W. Gelhar, A. L. Gutjahr and J. R. Macmillan 1978. Stochastic analysis of s p a t i a l v a r i a b i l i t y in subsurface flows, 1, Comparison of one-dimensional flows. Water Resour. Res. 14(2):263-271 . Basak, P. 1979. An a n a l y t i c a l solution for the transient d i t c h drainage problem. J. Hydrol. 41:377-382. Bear, J. 1979. Hydraulics of Groundwater. McGraw-Hill, New york. Benjamin, J. R. and C. A. Cornell 1970. Probability, S t a t i s t i c s and Decision Analysis for C i v i l Engineers. Mcgraw-Hill, New York. Bhattacharya, A. K. and R. S. Broughton 1979. Spacing drains to maximize revenue increases. Transcations of the ASAE 22(2) :351-356. Bhattacharya, A. K. and R. S. Broughton 1979. Variable drainable porosity in drainage design. Journal of the I r r i g a t i o n and Drainage D i v i s i o n , ASCE 105(IR1):71-85. Bierhuizen, J . F. and R. 0. Slatyer 1965. Effects of atmospheric concentration of water vapour and CO in determining transpiration-photosynthesisrelationships of cotton leaves. A g r i c u l t u r a l Meteorology, Vol. 2:259-270. Black, T. A. and M. J. Goldstein 1977. E f f e c t of ditch water l e v e l on root zone water regime and crop water use. Report to B.C. Department of Agriculture and Food, Department of S o i l Science, University of B r i t i s h Columbia, Vancouver, B.C., Canada as c i t e d by Oke, T. R. 1978. Boundary layer climate. Halstead Press Book, pp 372. Boersma, L. 1965. F i e l d measurement of hydraulic conductivity above a water table. In: C. A. Black (Ed.) Methods of S o i l Analysis, Part I. Agronomy 9:234-252. American Society of Agronomy, Madison, Wisconsin. Bouma, J., W. A. Z i e b e l l , W. G. Walker, P. G. Olcott, E. 162 McCoy and F. D. Hole 1972. S o i l Absorption of Septic • Tank Effluent - A F i e l d Study of Some Major S o i l s in Wisconsin. Information Circular No. 20. Geological and Natural History Survey, S o i l Survey Divi s i o n , Madison, Wisconsin, pp 235. Boumans, J . H. 1979. Drainage calculations in s t r a t i f i e d s o i l s using the anisotropic s o i l model to simulate conductivity conditions. Proceedings of the International Drainage Workshop, ILRI, Wageningen, The Netherlands. Boussinesq, J . 1903. Sur le debit. En temps de secheresse, d'une source alimentee par une nappe d'eaux d ' i n f i l t r a t i o n . Compt. Rend. Sceances Acad. S c i . 136:1511-1517. Bouwer, H. 1961. A double tube method for measuring hydraulic conductivity above a water table. S o i l S c i . Soc. Amer. Proc. 25:334-339. Bouwer, H. 1962. F i e l d determination of hydraulic conductivity above a water table with the double tube method. S o i l S c i . Soc. Amer. Proc. 26:330-335. Bouwer, H. 1964a. Unsaturated flow in groundwater hydraulics. ASCE Proc. 90(HY5):121-144. Bouwer, H. 1964b. Measuring horizontal and v e r t i c a l hydraulic conductivity of s o i l with the double-tube method. S o i l S c i . Soc. Amer. Proc. 28:19-23. Bouwer, H. 1966. Rapid f i e l d measurement of air-entry value and hydraulic conductivity of s o i l as s i g n i f i c a n t parameters in flow system analysis. Water Resour. Res. 2:729-738. Bouwer, H. 1974. Developing drainage design c r i t e r i a . In: Jan van Schilfgaarde (Ed.) Drainage for Agriculture. Agronomy Monograph 17, American society of Agronomy, Madison, Wisconsin. pp 67-79. Bouwer, H. and R. D. Jackson 1974. Determining s o i l properties. In: Jan van Schilfgaarde (Ed.) Drainage for Agriculture. Agronomy Monograph 17:611-672. American Society of Agronomy, Madison, Wisconsin. Broughton, R. S. 1976. Land drainage in Canada. Paper No. 76-2512 presented at the December 14-16 ASAE Meeting held at Palmer House, Chicago, 111. Bulavko, A. G. 1971. The hydrology of marshes anr marsh-ridden lands. UNESCO. Nature and Resources, Vol. 111(1). Chieng, S. T. 1975. The effects of subsurface drain depths and drainage rates on water table l e v e l s . Thesis presented to the McGill University, Montreal, Canada, in p a r t i a l fulfilment of the requirements for the degree of Master of Sc ience. Chieng, S. T. 1980. Computer-aided subsurface drainage system design and dra f t i n g . Thesis submitted to McGill University, Montreal, Quebec, Canada, in p a r t i a l fulfilment of the requirements for the degree of Doctor of Philosophy. Chieng, S. T., R. S. Broughton and N. Foroud 1978. Drainage rates and water table depths. Journal of the I r r i g a t i o n and Drainage D i v i s i o n , ASCE 104(IR4):413-433. Childs, E. C. 1952. The measurement of the hydraulic 163 conductivity of saturated s o i l in s i t u , 1, P r i n i c i p l e s of a proposed method. Roy. Soc. (London), P r o c , A. 215:525-535. Childs, E. C. 1960a. The nonsteady state of the water table in drained land. J. Geophys. Res. 65(2):780-782. Childs, E. C. 1960b. A treatment of c a p i l l a r y fringe in the theory of drainage. J . S o i l S c i . 1 9(1 ) :83-1 00,. Childs, E. C. 1960c. A treatment of the c a p i l l a r y fringe in the theory of drainage I I . J . S o i l S c i . 11(2):293-304. Childs, E. C , A. H. Cole and D. H. Edwards 1953. The measurement of the hydraulic permeability of saturated s o i l in s i t u , I I . Roy. Soc. (London), P r o c , 216:7.2-89. C i s l e r , J. 1969. The solution for maximum ve l o c i t y of isothermal steady flow of water upward from water table to s o i l surface. S o i l S c i . 108:148. Cordova, J. R. and R. L. Bras 1981. Physically based p r o b a b i l i s t i c models of i n f i l t r a t i o n , s o i l moisture, and actual evapotranspiration. Water Resour. Res. 17(1):93-106. Cornell, C. A. 1972. F i r s t order analysis of model and parameter uncertainty. In: International Symposium on Uncertainties in Hydrologic and Water Resources Systems, University of Arizona, Tucson, Arizona. Dagan, G. 1979. Models of groundwater flow in s t a t i s t i c a l l y homogeneous porous formulations. Water Resour. Res. l5(2):47-63. Davis, S. N. and R. J . M. De Wiest 1966. Hydrogeology. John Wiley & Sons, New York, 463 pp. De Vries, J. 1981 and 1982. Personal Communications. S o i l Science Department, University of B r i t i s h Columbia, B.C., Canada. Delhomme, J. P. 1978. Kriging in hydrosciences. Advances in Water Resources 1 (5):251-266. Dettinger, M. D. and J. L. Wilson 1 9 8 1 . F i r s t order analysis of uncertainty in numerical models of groundwater flow, 1, Mathematical development. Water Resour. Res. 17(1):149-161. Douglas, J. J r . and B. F. Jones 1963. On predictor-corrector methods for nonlinear parabolic d i f f e r e n t i a l equations. Jour. SIAM 11:195-204. Driehuyzen, M. G. 1982. Personal communications. B. C. Ministry of Agriculture and Food, S o i l s Branch, Cloverdale, B. C. Duke, H. R. 1972. Cap i l l a r y properties of s o i l s - Influence upon s p e c i f i c yoeld. Transactions of the ASAE 1 5(4)-.688-691 . Dumm, L. D. 1954. New formula for determining depth and spacing of subsurface drains in i r r i g a t e d lands. Agr. Eng. 35:726-730. Dumm, L. D. 1964. Transient-flow concept in subsurface drainage: i t s v a l i d i t y and use. Transcations of the ASAE 17(2):142-145, 151. Dumm, L. D. 1968. Subsurface drainage by transient flow theory. Journal of the I r r i g a t i o n and Drainage Di v i s i o n , 164 ASCE 94(IR4):505-519. Eagleson, P. S. 1978. Climate, s o i l and vegetation, 5, A derived d i s t r i b u t i o n of storm surface runoff. Water Resour. Res. 14(5):741-748. Ernst, L. F. 1950. A new formula for the ca l c u l a t i o n of the permeability factor with the auger hole method. T.N.O. Groningen 1950. Translated from the Dutch by H. Bouwer, Cornell University, Ithaca, New York, 1955. Feddes, R. A. and A. L. M. Van Wijk 1977. An integrated model-approach to the eff e c t of water management on crop y i e l d . Agriculture Water Management 1:3-20. Federer, C. A. 1979. A soil-plant-atmosphere model for transpiration and a v a i l a b i l i t y of s o i l water. Water Resour. Res 15(3):555-562. Fogel, M. M., L. Duckstein and W. T. O'Brien 1979. Optimal design and operation of drainage systems. Paper No. 79-2071 Presented at the 1979 Summer Meeting of the ASAE and CSAE at the University of Manitoba, Winnipeg, Manitoba, Canada. Freeze, R. A. 1975. A stochastic-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media. Water Resour. Res. 11(5):725-741. Freeze, R. A. and J . A. Cherry 1979. Groundwater. Prentice-Hall Inc., Englewood C l i f f s , New Jersey. 604 pp. Frevert, R. K. and D. Kirkham 1948. A f i e l d method for measuring the permeability of s o i l below the water table. Highw. Res. Bd., Proc. 28:450-452. Gardener, W. R. 1958. Some steady solutions of the unsaturated flow equation with application to evaporation from a water table. S o i l S c i . 85:228-233. Gardener, W. R. and M. Fireman 1958. Labortary studies of evaporation from s o i l columns in the presence of a water table. S o i l S c i . 85:244-249. Goldstein, M. J. 1978. Evaluation and prediction of subi r r i g a t i o n using a subsurface drainage system. Thesis submitted to the University of B r i t i s h Columbia, Vancouver,B.C., Canada, in p a r t i a l fulfilment of the requirements for the degree of Bachelor of Science in Ag r i c u l t u r a l Sciences, in the Department of S o i l Science. Gutjahr, A. L. and L. W. Gelhar 1981. Stochastic models of subaurface flow: I n f i n i t e versus f i n i t e domains and st a t i o n a r i t y . Water Resour. Res. 17(2):337-350. H a l l , A. E. and M. R. Kaukmann 1975. Regulation of water transport in the soil-plant-atmosphere continum, Ecological studies 12:187-210. Hare, F. K. And M. K. Thomas 1974. Climate Canada. Wiley Publishers of Canada Limited, Toronto, Ontario. Hayhoe, H. N. 1980. Calculation of workday p r o b a b i l i t i e s by accumulation over sub-periods. Canadian A g r i c u l t u r a l Engineering 22(0:71-75. H i l e r , E. A. 1969. Quantitative evaluation of crop drainage requirements. Transactions of the ASAE 12(4):499-505. H i l e r , E. A. and R. N. Clark 1971. Stress day index to characterize effects of water stress on crop y i e l d s . 165 Transactions of the ASAE 14(4):757-761. H i l l , A. R. 1976. The environmental impacts of a g r i c u l t u r a l land drainage. J. Environ. Mgmt. 4:251-274. H i l l i e r , F. S. and G. J . Lieberman 1981. Introduction to Operations Research. Holden-Day, Inc., San Francisco, C a l i f o r n i a . Hoffman, G. J . and G. O. Schwab 1964. T i l e drain spacing prediction based on drain spacing. Transactions of the ASAE 13(4):444-447. Hooghoudt, S. B. 1938. Bijdrage tot de kennis van enige natuurkundige grootheden van den grond. V e r s l . Landb. Onderz., 43(13)B, 215 pp. Hooghoudt, S. B. 1940. Bijdragen tot de kennis van enige natuurkundige grootheden van de grond nr. 7. V e r s l . Landbouwk. Onderz. 46:515-707. Hubbert, M. K. 1956. Darcy's law and f i e l d equations of the flow of underground f l u i d s . J. Petr o l . Technol. 8:222-239. Kirkham, D. 1954. Measurement of the hydraulic conductivity in place. Symposium on the Permeability of S o i l . Amer. Soc. Test Mater Spec. Tech. Pub. 163:80 97. Kirkham, D. 1958. Seepage of steady r a i n f a l l through s o i l s into drains. Trans. Am. Geophys. Union 39:829-908. Kirkham, D. 1967. Explanation of paradoxes in Dupuit-Forchheimer seepage theory. Water Resour. Res. 3:609-622. Kirkham, D. and C. H. M. Van Bavel 1948. Theory of seepage into auger holes. S o i l S c i . Soc. Amer. Proc. 13:75-81. Kirkham, D., S. Toksoz and R. R. Van Der Ploeg 1974. Steady flow to drains and wells. In: Jan Van Sch.ilfgaarde (ed.) Drainage for Agriculture. Agronomy Monograph 17:203-244. American Society of Agronomy, Madison, Wisconsin. Krayenhoff van de Leur, D. A. 1958. A study of nonsteady groundwater flow with special reference to a reservoir c o e f f i c i e n t . De Ingenieur 70 B:87-94. Luthin, J. N. 1980. Drainage Engineering. Robert E. Krieger Publishing Company, Huntington, New York. Luthin, J. N. and D. Kirkham 1949. A piezometer method for measuring permeability in s i t u below a water table. S o i l S c i . 68:349-358. Massland, M. 1959. Water table fluctuations induced by intermittent recharge. J. Geophys. Res. 64:549-559. M i l l e r , J. E. 1981. Personal Communications. U. S. Geological survey, Water Resources Division, Bismark, North Dakota. Moody, W. T. 1966. Nonlinear d i f f e r e n t i a l equation of drain spacing. ASCE Proc. 92(IR2):1-9. Moore, I. D. and C. L. Larson 1979. Effect of drainage projects on surface runoff from small depressional watersheds in the North Central Region. B u l l e t i n No. 99. Water Resources Research Center, University Of Minnesota, St. Paul, Minnesota. Moore, I. D. and C. L. Larson 1980. Hydrologic impact od draining small depressional watersheds. Journal of the 1 66 I r r i g a t i o n and Drainage D i v i s i o n , ASCE 106(IR4):345-363. Morgan, T. H., A. W. Biere and E. T. Kanemasu 1980. A dynamic model of corn y i e l d response to water. Water Resour. Res. 16(1):59-64. Morton, F. I. 1975. Fundamentals for a national basin program, Unpublished paper of the National Hydrology Research Institute, Environment Canada, Ottawa, Ontario, Canada. Morton, F. I. 1979. The hydrological effects of a g r i c u l t u r a l and r u r a l land drainage, Memorandum to Mr. D. H. Lennox, Director, National Hydrology Research I n s t i t u t e , Environment Canada, Ottawa, Ontario, Canada. Musy, A. and L. Duckstein 1976. Bayseian approach to t i l e drain design. Journal of the I r r i g a t i o n and Drainage D i v i s i o n , ASCE 102(IR3):317-334. N a i r i z i , S. and J. R. Rydzewski 1977. E f f e c t s of dated s o i l moisture on crop y i e l d s . Experimental Agriculture 13:51-59. Oosterbaan, R. J. and G. P. Wind 1979. Design and research. In: Jan Wesseling (Ed.) Proceedings of the International Drainage Workshop. ILRI, Wageningen, The Netherlands. Paul, C. L. and J. DeVries 1979. Effect of S o i l Water Status and strength on T r a f f i c a b i l i t y . Can J. S o i l S c i . 59: 313-324. Paul, C. L. and J. DeVries 1983. S o i l t r a f f i c a b i l i t y in spring, 2. Prediction and the ef f e c t of subsurface drainage. Accepted for publication in Can. J . S o i l S c i . Peck, A. J . , R. J. Luxmoore and J. L. Stolzy 1977. Effects of s p a t i a l v a r i a b i l i t y of s o i l hydraulic properties in water budget modeling. Water Resour. Res. 13(2):348-354. Penkava, F. 1974. Drainage. In: The a l l o c a t i v e c o n f l i c t s in water-resource management, University of Manitoba, Winnipeg, Manitoba, Canada. Agassiz Center for Water Studies. pp 121-135. Penman, G. L. 1967. Discussion on the groundwater, recharge study in Notheastern Jordon, Proc. Inst. C i v i l Engrs. 37:705-706. P h i l i p , J. R. 1980. F i e l d heterogeneity: Some basic issues. Water Resour. Res. 16(2):443-448. Pikul, M. F., R. L. Street and I. Remson 1974. A numerical model based on coupled one-dimensional Richards and Boussinesq equations. Water Resour. Res. 10(2):295-302. Polubarinova-Kochina, P. Ya 1962. Theory of Groundwater Movement. Princeton University Press, Princeton, New Jersey. 613 pp. Prasher, S. 0. 1981. Exploratory study of the ef f e c t of extensive a g r i c u l t u r a l drainage on the outflow hydrograph of a basin. Report submitted to National Hydrology Research Institute, Environment Canada, Ottawa, Ontario, Canada. Prych, E. A. 1980. Perturbation solutions for water table aquifers. Journal of Hydraulics D i v i s i o n , ASCE 106(HY4):603-608. Raats, P. A. C. 1971. Some properties of flows in 1 67 unsaturated s o i l s with an exponential dependence of the hydraulic conductivityupon the pressure head. J . Hydol. 14:129-138. Ravelo, C. J . 1978. A r a t i o n a l approach for incorporating crop needs into drainage system design. Thesis submitted to Texas A & M University, College Station, Texas, in p a r t i a l fulfilment of the requirements for the degree of Doctor of Philosophy. Reeve, R. C. and N. R. Fausey 1974. Drainage and timeliness of farming operations. In: Jan Van schilfgaarde (Ed.) Drainage for A g r i c u l t u r a l . Agronomy Monograph 17, Americam society of Agronomy, Madison, Wisconsin. PP 55-66. Remson, I., G. M. Hornberger and F. J. Molz 1971. Numerical Methods in Subsurface Hydrology. John Wiley & Sons, New York. Russell, S. 0. 1981. Class notes for the course- Decision Analysis for C i v i l Engineers (CIVL 495), C i v i l Engineering Department, University of B r i t i s h Columbia, Vancouver, B.C., Canada. Russo, D. and E. Bresler 1981. S o i l hydraulic properties as stochastic processes: I. An analysis of f i e l d s p a t i a l v a r i a b i l i t y . S o i l S c i . Soc. Amer. J. 45:682-687. Sagar, B. 1978a. Analysis of dynamic aquifers with stochastic forcing function. Water Resour. Res. 14(2):207-216. Sagar, B. 1978b. Galerkin f i n i t e element procedure for analyzing flow through random media. Water Resour. Res. 14(6):1035-1044. Sagar, B. 1979. Solution of l i n e a r i z e d Boussinesq equation with stochastic boundaries and recharge. Water Resour. Res. 15(3) :6l8-624. Sagar, B. and C. C. K i s i e l 1972. Limits of deterministic p r e d i c t a b i l i t y of saturated flow equations. In: Proceedings of the Second Symposium on Fundamentals of Transport Phenomenain Porous Media, Vol. 1:194-205. International Association of Hydraulic Research, Guelph, Ontario, Canada. Sagar, B. and C. P r e l l e r 1980. Stochastic model of water table under drainage. Journal of the I r r i g a t i o n and Drainage Division, ASCE 106(IR3):189-202. Sakkas, J. G. 1975. Generalized nomographic solution of Hooghoudt's equation. Journal of the I r r i g a t i o n and Drainage Division, ASCE 101(IR1):21-39. Sakkas, J. G. and V. Z. Antonopoulos 1981. Simple equations for Hooghoudt equivalent depth. Journal of the I r r i g a t i o n and Drainage Division, ASCE 107(IR4):411-414. Shady, A. M., R. S. Broughton, Y. Bernier, U. D e l i s l e and G. Sylvestre 1976. Experience with production scale f i e l d measurements of hydraulic conductivity. Paper Presented at the 3rd National Drainage Symposium, ASAE, Palmer House, Chicago, I l l i n o i s . Shaykewich, C. F. and L. Stroosnijder 1977. The concept of matric flux potential applied to simulation of evaporation from s o i l . Neth. J . ag r i c . S c i . 25:63-82. Sieben, W. H. 1964. Het Verband tussen onwatering en 1 68 ophbrengst b i j de jonge zavelgronden in de Noordoostpolder. Van Zee tot Land. 40, Tjeenk Willink V, Zwolle, The Netherlands. Singh, R. N. and S. N. Rai 1980. On subsurface drainage of transient recharge. J. Hydrol. 48:303-311. Singh, S. D. 1981. Moisture-sensitive growth stages of dwarf wheat and optimal sequencing of evapotranspiration d e f i c i t s . Agronomy Journal 73(3):387-391. Singh, S. R. and C. M. Jacob 1976. Numerical solution of Boussinesq equation. J. Eng. Mech. Div., ASCE 102(EM5):807-824. Skaggs, R. W. 1975. Drawdown solutions for simultaneous drainage and ET. Journal of the I r r i g a t i o n and Drainage D i v i s i o n , ASCE 91(IR3):9-22. Skaggs, R. W. 1976. Determination of the hydraulic conductivity-drainable porosity r a t i o from water table drawdown measurements. Transactions of the ASAE 19(1):73-80,84. Skaggs, R. W. 1978. A water management model for shallow watertable s o i l s . Technical Report No. 134, Water Resources Research Institute of the University of North Carolina, N.C. State University, Raleigh, N.C. Skaggs, R. W. 1979. Factors a f f e c t i n g hydraulic conductivity determinations from drawdown measurements. Paper presented at the 1979 Summer Meeting of the ASAE and CSAE, University of Manitoba, Winnipeg, Manitoba, Canada. Skaggs, R. W. and Y. K. Tang 1976. Saturated and unsaturated flow to p a r a l l e l drains. Journal of the I r r i g a t i o n and Drainage Divi s i o n , ASCE 102(IR2):221-238. Smedema, L. K. 1979. Drainage c r i t e r i a for s o i l workability. Neth. J. agric. S c i . 27:27-35. Smiles, D. E. and E. G. Young 1963. A multiple-well method for determining the hydraulic conductivity of a saturated s o i l in s i t u . J. Hydrol. 1:279-287. Smith, L. and R. A. Freeze 1979. Stochastic analysis of steady state flow in a bounded domain, 1, One-dimensional simulations. Water Resour. Res. 15(3):521-528. Sn e l l , A. W. and J. Van schilfgaarde 1964. Four-well method of measuring hydraulic conductivity in saturated s o i l s . Transactions of the ASAE 7:83-87,91. S o i l Conservation Service 1973. Drainage of A g r i c u l t u r a l Land. U. S. Department of Agriculture. Tang, Y. K. and R. W. Skaggs 1977. Experimental evaluation of t h e o r e t i c a l solutions for subsurface drainage and i r r i g a t i o n . Water Resour. Res. 13(6):957-965. Tanner, C. B. 1981. Transpiration e f f i c i e n c y of potato. Agronomy Journal, Vol. 73:59-64. Taylor, G. S. 1960. Drainable porosity evaluation from outflow measurements and i t s use in drawdown equations. S o i l S c i . 90(6):338-343. Taylor, S. A. and G. L. Ashkroft 1972. Physical Edaphology. The Physics of Irrigated and nonirrigated s o i l s . W. H. Freeman and Company, San Francisco, pp 533. T e r z i d i s , G. 1968. Discussion: F a l l i n g water table between 169 t i l e drains. ASCE 94(IR1):159. Tung, Y. and L. W. Mays 1980. Risk analysis for hydraulic design. Journal of the I r r i g a t i o n and Drainage D i v i s i o n , ASCE 106(HY5)-.893-913. UN Water Conference, 1977. Mar del Plata Action Plan. Report of the United Nations Water Conference, Extract 1. Van Beers, W. F. J. 1965. Some nomographs for the cal c u l a t i o n of drain spacing. B u l l . No. 8, ILRI, P. 0. Box 45, Wageningen, The Netherlands. Van Beers, W. F. J. 1970. The auger-hole method. B u l l . No. 1, ILRI, Wageningen, the Netherlands. Van Schilfgaarde, J. 1963. Design of t i l e drainage for f a l l i n g water tables. ASCE Proc. 89(IR2):1 -11. Van Schilfgaarde, J. 1964. Closing discussion: Design of t i l e drainage for f a l l i n g water tables. ASCE Proc. 90(IR3):71-73. Van Schilfgaarde, J. 1965. Transient design of drainage systems. ASCE Proc. 91(IR3):9-22. Van Schilfgaarde, J. 1970. Theory of flow to drains. Advances in Hydroscience 6:43-106. Acedemic Press, Inc., New York. Van Schilfgaarde, J. 1974. Nonsteady Flow to Drains. In: Jan Van Schilfgaarde (Ed.) Drainage for Agriculture, Agronomy Monograph 17:245-270. American Society of Agronomy, Madison, Wisconsin. Van Schilfgaarde, J. 1979. Progress and problems in drainage design. In: Jan Wesseling (Ed.) Proceedings of the International Drainage Workshop. ILRI, Wageningen, The Netherlands, pp 633-644. Vauclin, M., D. Khanji and G. Vachaud 1979. Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem. Water Resour. Res. 1089-1101. Vetter, W. J. 1970. Derivative operations on matrices, IEEE Trans. Autom. Control, AC-15:241-244. Vetter, W. J . 1971. Corrections to: Derivative operations on ' matrices. IEEE Trans. Autom. Conteol, AC-16. Vetter, W. J. 1973. Matrix calculus operations and Taylor expansions. SIAM Rev. 15(2) : 352-369. Wallis, C. H. 1981. Modeling forage growth in Peace River region. M.Sc. Thesis submitted to the S o i l Science Department of the University of B r i t i s h Columbia, Vancouver, B.C., Canada. Warren, J . E. and H. S. Price 1961. Flow in heterogeneous porous media. Soc. Petrol . Eng. J. 1:153-169. Warrick, A. W., G. J. Muller and D. R. Nielson 1977. Scaling field-measured s o i l hydraulic properties using a similar media concept. Water Resour. Res. 13(2):355-362. Werner, P. W. 1957. Some problems in non-artesian groundwater flow. Amer. Geophys. Union, Trans. 38(4):511-518. Wesseling, J . 1964. A comparison of the steady state drain spacing formulas of Hooghoudt and Kirkham in connection with design practice. J. Hydrol. 2:25-32. Wesseling, J. 1974. Crop growth and wet s o i l s . In Jan van Schilfgaarde (ed.) Drainage for Agriculture. Agronomy 1 70 Monograph 17:7-37. American society of Agronomy, Madison, Wisconsin. Williamson, R. E. and J. R. Kriz 1970. Response of a g r i c u l t u r a l crops to flooding, depth of the water table, and s o i l gaseous composition. Transactions of the ASAE 13:216-220. Wind, G. P. 1976. Application analog and numerical models to investigate the influence of drainage on workability in spring. Neth. J. agr i c . S c i . 24:155-172. Wind, G. P. and W. Van Doorne 1975. A numerical model for the simulation of unsaturated v e r t i c a l flow of moisture in s o i l s . J . Hydrol. 24:1-20. Winger, R. J., J r . 1960. In-place permeability tests and thei r use in subsurface drainage. Int. Congr. Comm. I r r i g . and Drainage, 4th, Madrid, Spain, 48 pp. Wiser, E. H., R. C. Ward and D. A. Link 1974. Optimized design of a subsurface drainage system. Transactions of the ASAE 17(1 ) : 1 75-178. Wisler, C. 0. And E. F. Brater 1965. Hydrology. John Wiley and Sons, N. Y., PP 55-56. Young, E. G. 1968. Shape factors for Kirkham's piezometer methos for determining the hydraulic conductivity of s o i l in s i t u for s o i l s overlying an impermeable floor or i n f i n i t e l y permeable membrane. S o i l S c i . 106:235-237. Young, T. C. and J. T. Ligon 1972. Water table and s o i l moisture p r o b a b i l i t i e s with t i l e drainage. Transactions of the ASAE 15(3):448~451. 171 APPENDIX A SOIL PROFILE DESCRIPTION Ladner S i l t y Clay Loam Horizon Depth, cm Description Ap 0-15 S i l t loam; primary structure - moderate grade and medium subangular blocky; secondary structure - moderate grade and subangular blocky; l i g h t gray (5Y 7.5/2 moist); mottles - common, fine size and d i s t i n c t contrast; contains fragments of Aegj. Aegj 15-25 S i l t loam; primary structure - moderate grade and coarse subangular blocky; secondary structure - moderate grade and fine subangular blocky; mottles - many, coarse size and prominent contrast. Btg1 25-40 S i l t y clay loam; primary structure moderate to strong grade and medium prismatic; secondary structure - weak to moderate grade and fine subangular blocky; l i g h t o l i v e gray (5Y 6/2 moist); mottles -many, coarse size and prominent contrast. Btg2 40-55 S i l t y clay loam; primary structure moderate grade and coarse angular blocky; secondary structure - weak grade and fine 1 72 subangular blocky; mottles - many, medium and prominent contrast. BCg 55-115 S i l t loam; primary structure - moderate grade and coarse subangular blocky; secondary structure - moderate grade and fine subangular blocky; mottles - many, medium size and prominent contrast; concretions, nodules and casts - few, medium size and irregular shape. Cg1 115-120 Loamy sand; massive single grained structure; mottles - few, fine size and fai n t contrast; concretions, nodules and casts - common, medium size and irregular shape. Cg2 120-130 Sandy loam; massive single grained structure; mottles - few, fine size and faint contrast; concretions, nodules and casts - common, medium size and irregular shape. 173 APPENDIX B A NEW METHOD FOR DESIGNING DRAINAGE SYSTEMS In practice, subsurface drainage systems are designed on the basis of either steady saturated flow or nonsteady saturated flow theories. Common inputs to these design tools include both s o i l properties and boundaries of the flow zone. These are, among others, the saturated hydraulic conductivity of the s o i l and the depth of "the impermeable layer below the drains. However, both of these properties are seldom e x p l i c i t l y known in most p r a c t i c a l designs. The saturated hydraulic conductivity i s a random property of the s o i l and i s known to vary both in space and with time by several orders of magnitude (Freeze and Cherry, 1979; DeVries, 1982). The location of the impermeable layer i s known e x p l i c i t l y only i f i t l i e s within few meters from the s o i l surface. It i s usually guessed in most p r a c t i c a l designs. Both of these design parameters are so undeterministic that they seem to consitute weak lin k s in the drainage theory. Some researchers have attempted to consider them in p r a c t i c a l designs in a somewhat nndirect way (Skaggs, 1979). Few drain l a t e r a l s or sometimes a single main l i n e are l a i d f i r s t and the i r effect on the surrounding water table i s observed. The water table recession measurements are used with any relevant drainage theory to obtain a few e f f e c t i v e design parameters that would explain the observed water table fluctuations. Generally, estimates are made about the ef f e c t i v e saturated hydraulic conductivity in t h i s way and are used in the 174 subsequent analyses to upgrade drainage system designs. It i s contended by some researchers (Skaggs, 1979) that errors due to the imperfect knowledge about the location of the impermeable layer below the s o i l are somewhat reduced in such analyses. In t h i s section, we w i l l discuss an indirect method of designing drainage systems by steady state theory which assumes no pre-knowledge about the saturated hydraulic conductivity and the depth to the impermeable layer below the drains. The method is based upon the Hooghoudt's formula for designing drainage systems on the basis of steady state theory. Under steady state flow conditions, the t o t a l inflow into the system equals the t o t a l outflow from the system, but the flux density into and out of the system are not the same. Under steady state flow conditions, Hooghoudt's formula may be written as S 2 = 4 K h (2d + h) / R where S i s the drain spacing, K i s the saturated hydraulic conductivity of the s o i l , h is the permissible height of the water table at the mid-spacing, d i s the depth of the impermeable layer below the drains and R i s the drainage c o e f f i c i e n t . To account for the convergence of flow l i n e s near the drains, the depth of the impermeable layer below the drains is replaced by some equivalent depth in p r a c t i c a l applications 175 of t h i s formula. The method proposed in t h i s section i s intended for use when the saturated hydraulic conductivity of the s o i l and/or the location of impermeable layer are unknown. The l a t t e r w i l l be the case when the impermeable layer is located at greater distances from the s o i l surface such that i t s detection i s not possible with the usual tools at the disposal of designers. In such cases, i t can be assumed that a major proportion of the flow takes place below the drains, which are generally 1.2 m below the s o i l surface. For such cases, the Hooghoudt's formula can be written as 2 _ = 8 K h d / R or ' 2 _ = C h / R (1 ) where C i s a c o e f f i c i e n t . The Hooghoudt formula as such assumes no entry resistance to groundwater flow entering the drains. The entry resistance to flow i s due to the number and shape of drain openings and the nature of the drain envelope. Hooghoudt (1938) suggested that the entry resistance to flow can be considered by re-defining the parameter K as an e f f e c t i v e K and that i t may be estimated from water table observations at the mid-spacing and volume outflow rate after i n s t a l l i n g a few drains in the f i e l d . It i s also imperative that effects of 176 convergence of the flow lines near the drains are i m p l i c i t l y considered in such methods. In the proposed method, we also define K as an e f f e c t i v e K and, therefore, the c o e f f i c i e n t C not only includes the saturated hydraulic conductivity of the sub-soil and the depth of the impermeable layer but also the entry resistance to flow and the r a d i a l flow near the drains. The method discussed in t h i s section hinges upon the determination of the c o e f f i c i e n t C for a given location. Once C i s known, di f f e r e n t drainage systems designs can be tested quickly to check their effectiveness in c o n t r o l l i n g the water table on a steady state basis. The following procedure may be used to estimate the c o e f f i c i e n t : 1. At least three or at the most five drains may be l a i d at a drain spacing and drain depth that has been found to be appropriate for similar s o i l s . The i n i t i a l guess about the spacing requires pertinent f i e l d experience in that location because i t is made without any s p e c i f i c information about the design parameters. It w i l l be shown later that even i f the i n i t i a l guess is off by up to one hundred percent, the proposed method w i l l lead to a correct design for the rest of t h e ' f i e l d . However, the proposed method w i l l only lead to correct designs for the rest of the f i e l d i f the s o i l s are resonably uniform in the 1 77 f i e l d . Some measurements of the hydraulic conductivity may be made in other parts of the f i e l d to roughly test the assumption of uniform s o i l s . 2. Daily observations should be made on the midspan water table height and the corresponding drain outflow for a few days. The c o e f f i c i e n t C can be estimated from equation (1) by substituting the observed drainage rate for R and the corresponding midspan water table height for h. 3. The representative value of the c o e f f i c i e n t can be estimated by averaging the C values obtained from step 2. 4. Once a representative C value i s known, the design drain spacing can be calculated from equation (1) that w i l l result in a desired midspan water table height above the drains on a steady state basis. Ideally, the average value of the c o e f f i c i e n t C would be determined from f i e l d observations of water table fluctuations after i n s t a l l i n g a few drains. However, i t was not possible for most of the examples that were considered. Therefore, Kirkham's method (Kirkham et a l . , 1974) was used to obtain the steady state midspan water table height above the drains that w i l l result from the i n i t i a l design. 178 The above method has been tested on a variety of hypothetical but relevant drainage system designs. It appears to provide correct estimates of the drainage system design in a l l cases. In the following, a number of design examples are discussed that were used in the validation of the method: Example 1. We need to design a drainage system for an area with a drainage c o e f f i c i e n t of 1.2 cm/day. We assume that the saturated hydraulic conductivity of the s o i l i s 40 cm/day, the depth to the impermeable layer below the drains i s 200 cm and the diameter of the drain t i l e (2r) is 10 cm. The above data w i l l be required in the model vali d a t i o n phase where the model results w i l l be compared with those of the Kirkham method. It w i l l be c a l l e d the Kirkham method in the rest of the chapter. However, i t should be kept in mind that the designer has no information about the above s o i l properties or the boundaries of the flow zone as w i l l be the common case in most designs in pr a c t i c e . Let us assume that the i n i t i a l guess about the possible drain spacing i s 20 m for the problem. The steady state height of the water table above the drains w i l l then be 126.8 cm (from Kirkham's method). Now, the c o e f f i c i e n t C can be estimated from equation (1) as 37855 cm2/day. For spring t r a f f i c a b i l i t y , the permissible water table depth below the surface i s approximately 50 cm. For the example, i t implies a permissible midspan water table height above the drains of 70 cm. Therefore, the drain spacing that 179 w i l l meet th i s design requirement can be obtained from equation (1) by substituting a value of 70 cm for h. A drain spacing of 15m i s obtained. Had the designer been provided with the information with respect to the s o i l properties and the boundaries of the flow zone, a drain spacing of 14m would have been required from the Kirkham method. The difference between the two designs i s about 7%. Had the i n i t i a l guess about the drainage system design been 25 m (approximately double the actual spacing), a drain spacing of 16m would have resulted. Again the difference between the two designs i s about 14%. As a matter of interest, had we assumed a drain spacing of 40 m, a spacing of 17 m would have been obtained (a difference of 21%). Example 2. The data considered in t h i s design example are as follows: K = 400 cm/day d = 300 cm D = 120 cm We obtain an average value of C as 664599 cm2/day for an i n i t i a l drain spacing of 50 m (permissible midspan water table height above drains under steady state flow conditions i s obtained as 45.14 cm for the spacing from Kirkham's formula). Using the above value of C and for an h of 70 cm, we obtain a drain spacing of about 62 m. R = 1.2 cm/day h = 70 cm 2r = 10 cm 180 From the Kirkham's method, a drain spacing of 66 m w i l l result by employing the above input data. The difference between the two designs i s about 6%. Example 3. The input data i s the same as in example 2 with the exception that the impermeable layer i s 50 cm below the drains (d = 50 cm). Using an i n i t i a l guess of 20 m as the drain spacing, we obtain the design spacing of 29 m. From the Kirkham's method, a drain spacing of 30 m w i l l be required to meet the design requirements. The difference between the two methods i s about 3%. Had our i n i t i a l guess about the drain spacing been 40 m, the design spacing would be 32 m. Again, the difference between the two methods i s about 7%. Example 4. The input data i s the same as in the last example but the saturated hydraulic conductivity i s 50 cm/day. With an i n i t i a l guess of 6 m as the drain spacing, we obtain a design spacing of 9 m. From the Kirkham's method, a drain spacing of 10m i s obtained to meet the design requirements. Had we used a drain spacing of 20 m as our f i r s t guess, the correct spacing of 11 m would have been obtained. Example 5. In th i s section, we w i l l assume a transient case when drains have just been l a i d and the observations are made on the midspan water table height above ,the drains for few days. 181 These observations can be used to estimate the average value of the c o e f f i c i e n t C which, in turn, w i l l be used in the estimation of correct drainage system design. However, such an exercise w i l l require a relationship between the drain discharge per unit f i e l d area and the midspan water table height above the drains. Such data was obtained from the l i t e r a t u r e (Chieng, 1975; Figures 5 and 12). Assuming a drain spacing of 36.58 m (120 f t ) , C values of 32636, 48658 and 64747 cm2/day are obtained for drainage c o e f f i c i e n t s of 0.2, 0.4 and 0.6 cm/day corresponding to h values of 82, 110 and 124 cm. Therefore, the average C value for the location is 48680 cm2/day. The other input data for the s i t e i s assumed as follows: d = 200 cm K = 22 cm/day We obtain a drain spacing of about 37 m using the Hooghoudt's formula. Using an average C value of 48680 cm2/day, a drain spacing of 39 m i s obtained. It i s very close to the spacing as would have been resulted provided there was an e x p l i c i t knowledge about the various design parameters. It i s assumed in the analysis that a major proportion of the groundwater flow takes place below the drains. However, the assumption may be invalidated i f the i n i t i a l drain spacing i s very large as compared to the correct drain spacing. A large h = 80 cm R = 0.253 cm/day 1 82 drain spacing w i l l result in very high water tables above the drains and therefore, the flow above the drains can not be neglected. Also, the entry resistance to flow w i l l be quite d i f f e r e n t (from the case when the correct drain spacing i s used) because of the increase in the hydraulic gradient responsible for the flow of groundwater towards the drains. On the other hand, i f the i n i t i a l drain spacing i s very small as compared to the correct drain spacing the entry resistance to groundwater flow w i l l again be very d i f f e r e n t due to the decrease in hydraulic gradients above the drains. However, in practice l o c a l designers should be able to ident i f y the i n i t i a l drain spacing close enough to the correct spacing for most cases. The analysis appears to be applicable as long as the i n i t i a l drain spacing i s not d i f f e r e n t from the correct drain spacing by more than one hundred percent. The l i m i t a t i o n s of this method are mostly the same as for the Hooghoudt's formula except that the proposed method may not be applicable when the s o i l heterogeneity is such that a fine textured s o i l l i e s under a coarse textured s o i l . Such a case w i l l invalidate the major assumption in the method that most of the groundwater flow takes place below the drains. The method is applicable to those cases where a coarse textured s o i l l i e s under a fine textured s o i l . In the l a t t e r case, the method i s applicable even though the drains are located in the lower coarse textured layer rather than at the boundary between the two layers. However, the method i s not applicable when drains are located in the upper layer (then the Hooghoudt's formula 183 i t s e l f i s not applicable - Van Beers, 1965). The above examples indicate that the proposed method can be used to design drainage systems in the absence of any d e f i n i t e knowledge about the s o i l parameters or the flow boundaries. The method compares favourably with Kirkham's method for a variety of drainage system designs. It provides a p r a c t i c a l way of designing drainage systems. Moreover, the method makes use of the phenomenon that i s to be predicted ultimately and therefore, should be more accurate than d i r e c t methods that use point values of the design parameters. The prodedure, though somewhat cumbersome, should be amenable to day-to-day p r a c t i c a l designs. 184 APPENDIX C ENVIRONMENTAL IMPACTS OF AGRICULTURAL DRAINAGE Ag r i c u l t u r a l drainage has been practised for centuries but i t s hydrologic effects have generally been studied in only an a g r i c u l t u r a l context, i . e . lim i t e d to the top s o i l and i t s a g r i c u l t u r a l productivity. The lack of understanding of i t s environmental effects can be ascribed very d e f i n i t e l y to the lack of understanding of the hydrology of the unsaturated zone and the soil-vegetation-atmosphere interface. In Canada, only 4.3% of the t o t a l land area i s suitable for c u l t i v a t i o n (Broughton, 1976). About half of the cu l t i v a t a b l e land in the humid areas of eastern Canada and some areas of western Canada need drainage improvement. Subsurface drains are being l a i d at the rate of 62,000 ha/year or about 1% of the cu l t i v a t a b l e land needing drainage improvement (Broughton, 1976). Although large areas of land have been drained and there i s an ever increasing trend among the farmers to drain t h e i r lands, l i t t l e emphasis i s being given to the possible environmental impact of such a drastic change in the hydrologic c h a r a c t e r i s t i c s of the land surface. Most texts on drainage (e.g. Luthin, 1980; Van Schilfgaarde, 1974; S o i l Conservation Service, 1973) provide l i t t l e or no reference to such e f f e c t s . To my knowledge, the only text book which addresses th i s problem in more than a cursory way is that of Wisler and Brater (1965) which discusses the effects of such change studies on the watershed. They also try to es t a b l i s h the causative factors and 185 a n a l y t i c a l problems in these studies. Drainage improvement changes the s o i l moisture regime and thereby influences natural vegetation, w i l d l i f e , regional water tables, both quantity and quality of runoff ( H i l l , 1976). In th i s chapter, the writer w i l l be addressing the role that drainage plays in the generation of runoff in a watershed. S p e c i f i c a l l y , a possible c o r r e l a t i o n between a g r i c u l t u r a l drainage and flood peaks on a basin-wide scale w i l l be addressed. It must be noted that we are not too concerned about major flood events l i k e 25-, 50- or 100-year floods. In such major events the role that a g r i c u l t u r a l drainage plays i s probably i n s i g n i f i c a n t . Our concern l i e s with the moderate floods of 10-or 15- year recurrence interval or less which may be aggravated by the increased a g r i c u l t u r a l drainage a c t i v i t y in the basin. The nature of such a c o r r e l a t i o n i s a bone of contention among experts in engineering and related d i s c i p l i n e s (Morton, 1979). There are two main schools of thought regarding these change effects on a watershed. The f i r s t school argues in favour of drainage and i s widely supported by personnel in a g r i c u l t u r a l d i s c i p l i n e s . Their oulook i s that drainage plays a corrective role in the hydrologic cycle by decreasing flood peaks. They j u s t i f y their claim on the basis of the presence of an unsaturated zone above the water table in drained lands. The rain water or snowmelt must f i r s t f i l l the available pore space before being transmitted to the water table or becoming available for runoff. The second school argues against drainage 186 and i s widely supported by personnel in various environmental d i s c i p l i n e s . They agree with farmers l i v i n g downstream in a basin that drainage may be one of the factors causing higher flood peaks in recent years (South-Nation basin, Ontario; Red River basin, Manitoba). They argue that drainage provides new and faster flow paths for the storm water, thereby delivering that water to the downstream end which.otherwise would have stayed on the upstream 'side or slowly drained toward the downstream end. In other words, they f e e l that a g r i c u l t u r a l drainage increases the e f f e c t i v e size of the basin by providing a link between areas which previously flowed to "blind channels" and the downstream channel. H i l l (1976) presents the following reasons as to why a single consensus has not been achieved: "The view i s often expressed that land drainage inevitably results in increases in ri v e r flood peaks as drainage reduces the o r i g i n a l surface retention capacity of the land and accelerates the flow of water into natural water courses. However, there are a variety of reasons why increases in flood peaks may not occur following land drainage. Drainage-induced changes in stream hydrographs are clo s e l y related to the amount of land drained in proportion to the size of the rive r basin. Small drainage schemes are not l i k e l y to produce major changes in flow at the downstream s i t e s . The location of drained land within the r i v e r basin i s also a s i g n i f i c a n t factor. The drainage of considerable areas in the lower portion of a rive r basin may actually reduce flood peaks at the basin outlet because water from t h i s area w i l l be removed 187 before the a r r i v a l of flood peaks from the headwater areas (Wisler and Brater, 1965). Furthermore, land drainage can have a p o s i t i v e effect by increasing the i n f i l t r a t i o n rates and the water retention capacity of s o i l s which may compensate in part for the losses in surface retention (Penkava, 1974)." A g r i c u l t u r a l drainage a f f e c t s streamflows due to a variety of reasons. It changes runoff factors of the basin by a f f e c t i n g the flow paths of storm water. Also, such a change brings about a change in land use where a land formerly under hay i s now put under more valuable crops l i k e corn. This, in turn, has the following e f f e c t s on the hydrologic cycle: (a) a change in interception pattern; (b) a change in evapotranspiration; and (c) a change in root zone storage due to extraction of s o i l moisture from increased depths. The cumulative e f f e c t of the above changes brought about by a g r i c u l t u r a l drainage on a watershed w i l l make the usefulness of h i s t o r i c streamflow records questionable. In most hydrologic studies, we assume st a t i o n a r i t y in h i s t o r i c a l records. Had drainage been a f f e c t i n g the streamflows (increased runoff, change in land use) then we would have every reason to suspect non-stationarity in streamflow records of the l a s t two or three deacdes. This non-stationarity in h i s t o r i c streamflow records should be detected and taken into account in analysis, before such records are used in the design of various hydrologic structures. The above effects on basin hydrology are v i s i b l e conceptually only in the case of surface or d i t c h drainage. There i s l i t t l e c o n f l i c t among hydrologists about i t s e f f e c t s on 188 the time of concentration of the storm water in a basin. But in the case of subsurface drainage, t h i s relationship i s not f u l l y understood yet and needs more research before statements can be made about i t s environmental impacts on the basin. The environmental impacts of a g r i c u l t u r a l drainage may have gone undetected for such a long time due to the fact that generally a drainage system becomes i n e f f i c i e n t in two or three years after i t s i n s t a l l a t i o n because of poor maintainance and a corresponding reduction in carrying capacity of the water courses. Weed growth and slumping are the main causes of the deterioration of these channels. In properly maintained drainage projects, when ditches are dredged or vegetation k i l l e d off at regular interva l s , impact of a g r i c u l t u r a l drainage on streamflows w i l l become more permanent. Other reasons for inconclusive or contradictory results are due to our i n a b i l i t y to map the e f f e c t s of differences in drain configuration, s o i l s , topography and climate. Above a l l , inadequate theory and i n s u f f i c i e n t data l i m i t our predictive a b i l i t y . From the above, i t i s clear that our knowledge i s scanty as to the side-effects of a g r i c u l t u r a l drainage. With the current rate of land going under a g r i c u l t u r a l drainage, i t is surprising that we have been ignoring i t s possible effects on the hydrology of the watersheds for so long. We must remember that although these e f f e c t s are not very v i s i b l e now, there i s enough reason to suspect that ignoring these e f f e c t s now may prove very costly l a t e r . By studying t h i s phenonmenon now, we might end up with a 189 strategy for developing watersheds in such a manner that undesirable effects on streamflows are minimized. 1. Extent of the Problem Both planned and unplanned manipulations of watersheds affe c t our water resources. Those effects which are due to a change in land use or c u l t u r a l practices may prove b e n e f i c i a l or non-beneficial to water resources. In either case, quantitative predictions of these effects are mandatory before s i g n i f i c a n t land use modifications are undertaken or allowed to occur. In the context of a g r i c u l t u r a l drainage, these effects have been recognized very recently in North America and abroad. Certainly, t h i s problem does ex i s t , to a varying degree, in many parts of the world wherever extensive a g r i c u l t u r a l drainage i s prevalent or planned.- In th i s section, the problem w i l l be reviewed for Canadian conditions. For a somewhat global view of the problem, the reader is referred to a report submitted by the writer to Environment Canada in 1981 t i t l e d "Exploratory Study of the Ef f e c t of Extensive A g r i c u l t u r a l Drainage on the Outflow Hydrograph of a Basin". In Canada, the problem has recently caught the attention of personnel involved in land and water conservation, environmental and natural resources management. It i s of a varying importance in d i f f e r e n t provinces depending upon the extent of a g r i c u l t u r a l drainage and the susceptability to flooding. Evidently, i t is a more serious problem in the central provinces than elsewhere mainly due to the large areas which are under a g r i c u l t u r a l 1 90 drainage and an ever-increasing tendency among farmers in these provinces to drain their lands for economic reasons. There are four main types of drainage works that can be distinguished for Canadian conditions: a) Drainage works to minimize the adverse ef f e c t s of flooding of a g r i c u l t u r a l land by improving drainage channels and structures for rapid carry-off of peaks flows (Lower Fraser Valley, B.C.). This is a t y p i c a l engineering problem and as such f a l l s outside the scope of t h i s chapter. There are no s o c i a l c o n f l i c t s associated with t h i s type of drainage improvement. b) Drainage works that provide outlets for wetlands and potholes in areas of internal drainage, thereby increasing the e f f e c t i v e drainage area of a water course ( P r a i r i e Provinces). This type gives r i s e to s o c i a l c o n f l i c t s between landowners in the drained upland areas and the property owners in the flood plains, adding the burden of regulating and assesing impacts of drainage improvements on the water resource manager and the public adminstrators. 191 c) Intensive drainage of r e l a t i v e l y f l a t farmlands for improved s o i l management and e a r l i e r access to f i e l d s for farm machinery. Drainage may be underground by means of subsurface drains (Ontario, Quebec, B r i t i s h Columbia) or by f i e l d ditches (Manitoba), both types feeding into open ditches and large municipal drains. The choice between subsurface and ditch drainage depends largely on the permeability of the s o i l , and economic reasons. There i s a f a i r l y strong but unverified consensus that ditch drainage increases flood peaks, and no consensus at a l l on the effe c t of subsurface drainage; the difference however may be largely due to the higher v i s i b i l i t y of ditch drainage; and even in the case of ditch drainage, qu a n t i f i c a t i o n remains d i f f i c u l t and inaccurate. Social problems associated with this type of drainage improvement are similar to those described under b), perhaps with the note that perception and public awareness of the problem are higher in the case of ditch drainage because of i t s higher v i s i b i l i t y . 192 d) Mining of peat, thereby draining marshes. Although not s t r i c t l y belonging to a g r i c u l t u r a l drainage, i t i s mentioned here because of i t s possible future importance in the Maritimes, where mining of peat i s increasing and i t s potential as a fuel for power generation is under study. Social problems may arise as under b) and c ) . The bulk of the research reported in the l i t e r a t u r e that deals s p e c i f i c a l l y with the topic of this report i s based on comparative basin studies and trend analyses. In the comparative analysis, a comparison i s made between simultaneous runoff observations in adjacent drained and undrained areas (Domain Crop Area Demonstration Project, Manitoba). In runoff trend analysis, multiple regression models are used to investigate a change in the response function of a watershed (Boyne River at Carman and Rat River near Otterburne in Manitoba). These studies may or. may not give an answer to the question: "Does drainage improvement increase t o t a l discharge or peak flows in t h i s p a r t i c u l a r basin?" However, their findings are not transposable to other basins, and do not add s i g n i f i c a n t l y to our understanding of the physical processes involved. Other research has t r i e d to use available hydrologic models to simulate the. effects of changing land use. However, such models depend heavily on careful c a l i b r a t i o n against h i s t o r i c records and may or may not have predictive c a p a b i l i t y in the 193 basin for which they are cal i b r a t e d . U n t i l there i s a clear understanding of which transfer parameters of the model are affected by t i l e or ditch drainage, and in which way and how much, their predictive c a p a b i l i t y in a changing watershed i s suspect. B a s i c a l l y , the hydrologic response of a basin i s the resu l t , by integration over time and space, of a complex set of routing or p a r t i t i o n i n g decisions - governed by the laws of physics - moving water from i t s points of entry in the basin, through a complex network, to the outlet or the main water course. Time and place of drainage improvement are therefore of great importance in the generation of peak discharge; with t h i s concept in mind some q u a l i t a t i v e statements can be made: a) Surface drainage works in the lower reaches of a basin tend to carry off pr e c i p i t a t i o n before water from the upper reaches can be routed through the network, thereby tending to "smooth" out discharge peaks; on the other hand, drainage works in the upper reaches tend to make discharge "catch up" with the discharge from the lower reaches, thereby increasing peaks. b) In large basins, drainage works, at any given time, w i l l be i n s t a l l e d , more or less randomly in the component sub-basins; 1 94 l o c a l l y , individual, small drainage basins may suffer increased peakflows in the lower order watercourse, but in the higher order watercourse some of these l o c a l e f f e c t s w i l l tend to cancel because of their d i f f e r e n t timing. Thus, the effect of drainage improvement w i l l be r e l a t i v e l y less severe in large basins; t h i s tendency has been observed in the Minnesota study (Moore and Larson, 1979). c) E f f e c t s of a g r i c u l t u r a l drainage can be minimized by promoting a balanced development of drainage improvement, in which drainage improvement in the upper reaches i s balanced by simultaneous drainage improvement in the lower reaches. In th i s connection, consideration should be given to the p o s s i b i l i t y of providing special catchment areas in the upper reaches of the basin, thereby lessening the impact of drainage on the peak flows as well as providing measure of drought proofing at the same time. At present, such basin management i s more of an art than a science. It should be noted that although these statements may seem plausible they are not supported by any r e l i a b l e s c i e n t i f i c 195 evidence. Therefore, there is a chance that they could be invalidated by the kind of unknown interactions or paradoxical cause-and-effect relationships described in the next section. 2. Modeling Approach In t h i s section, we w i l l discuss possible answers to the question "How can we predict the hydrologic e f f e c t s of changing land use?" We w i l l also look at the types of models at our disposal to make such predictions. A f i r s t step towards modeling a change study i s to d i f f e r e n t i a t e between input and watershed nonstationarities. We must convince ourselves that the observed change in watershed response i s not due to input disturbance alone. A careful evaluation of changes in input i s prerequisite to meaningful detection of any watershed changes; otherwise nonstationarity in input may get f a l s e l y assigned to the watershed. Having established that a change study i s required to aid in the futher watershed management, hydrologists seem to have followed mainly two routes, namely, s t a t i s t i c a l approach and mathematical modeling. The s t a t i s t i c a l approach i s simple and involves manipulating the h i s t o r i c data to detect new trends. It is valuable in i d e n t i f y i n g the nature of the change by helping to answer questions about s i g n i f i c a n t differences between two mean values, variances, regression slopes, and other parameters in either watershed input or output. We must have a very long record for t h i s type of analysis and i t s usefulness probably 196 l i e s in the detection of a change, not in i t s prediction. The approach has i t s l i m i t a t i o n s when the effects being investigated are r e l a t i v e l y small and d i f f i c u l t to separate due to other factors (Moore and Larson, 1980). Mathematical models require a sound mathematical basis and a detailed knowledge as well as mathematical description of a l l hydrologic processes involved. For the problem under investigation, i t implies a thorough knowledge about processes l i k e i n f i l t r a t i o n , runoff generation, overland flow and streamflow routing, snowmelt runoff, unsaturated flow in the s o i l and evapotranspiration. Once, a mathematical model i s formulated that gives due consideration to these processes there should not be any problem in using i t to make predictions i f one of the processes such as runoff generation i s affected by a change in land use. However, most of the physical changes are complex in nature and their multifarious e f f e c t s cannot be represented adequately by changing one or two parameters in the model. For example, a g r i c u l t u r a l drainage may not only a f f e c t the flow paths of the excess water in the basin but may also lead to a change in crop use practice which, in turn, may a f f e c t evapotranspiration. To represent a l l these effects in mathematical models has always been and perhaps always w i l l be a d i f f i c u l t job for hydrologists. To circumvent these problems, we simplify our models by invoking certain assumptions which may or may not be v a l i d . As a result, numerous hydrologic models have come to the surface and more are coming every year. It has become d i f f i c u l t to 1 97 ide n t i f y which model i s applicable to a problem at hand. Often, indiscriminate use of . the models i s made without ascertaining beforehand whether a given model i s appropriate for the given problem or not. The choice of an appropriate model must be dictated by the manner in which i t considers physical processes and their i n t e r r e l a t i o n s h i p s , by data requirements and by the desired l e v e l of accuracy. M i l l e r (1981) studied the effects of a g r i c u l t u r a l drainage on the discharge of the Red River in N. Dakota. He evaluated many watershed models on the basis of how they take into account relevant physical processes. He concluded, to no one's surprise, that none of the models was suitable for t h i s study. He recommended the formulation of a new model which addresses each involved processe accurately. From my personal experience, I f u l l y agree with M i l l e r ' s suggestion that a new model is needed for t h i s type of study. None of the watershed models considers various physical processes such as i n f i l t r a t i o n , evapotranspiration according to the present state of knowledge. A watershed model may be formulated by making f u l l use of the present state of the art yet i t s refined data requirements w i l l rule out i t s use as a p r a c t i c a l management t o o l . The increase in computational time w i l l be another constraint on i t s a p p l i c a b i l i t y . Moreover, in such studies, there are so many parameters and variables which can never be mapped for a watershed that a refined model may not be j u s t i f i a b l e . However, i f such a model i s formulated and the above constraints are 198 overcome by spending considerable amounts of money, time and e f f o r t , i t w i l l enable us to see more c l e a r l y the e f f e c t s of man-made physical changes on our watersheds in a short time. Of course, i t w i l l be an engineering study which w i l l y i e l d better solutions to the problem in a short time. There i s l i t l e or no room for research on various involved processes in such an approach. Now, a question may ari s e about the v a l i d i t y of our current knowledge "Can we characterize the various physical processes accurately enough without current knowledge?" This question is very d i f f i c u l t to answer and hydrologists may have d i f f e r e n t opinions about i t . I think that although a l o t of research has been done in the l a s t two decades on basin 'hydrology, very l i t t l e progress has been made in acquiring fundamental knowledge about the cause and e f f e c t relationships in the soil-vegetation-atmosphere system. Morton (1975) gives an interesting summary of the present state of the art in hydrology and current modeling approaches. "The sophisticated models attempt to synthesize a l l aspects of the hydrologic cycle. However, with the present state of knowledge about hydrologic processes and their interactions, such models are merely conglomerations of plausible but untested hypotheses. Therefore, confidence in the results requires confidence in the conventional wisdom. The problem is that the conventional wisdom cannot be trusted because of paradoxical interactions in the soil-vegetation-atmosphere system. This can be demonstrated by 199 known interactions associated with evapotranspiration - a process that normally consumes 50% to 95% of the p r e c i p i t a t i o n input to a basin. It i s well established that a decrease in atmospheric humidity w i l l increase the evaporation from a continuosly moist surface ( i . e . the potential evaporation). Therefore, i t has seemed reasonable to assume that such a decrease would also increase the flux of moisture from the s o i l through the vegetation to the atmosphere. In watershed models thi s i s expressed by a relationship which has evapotranspiration proportional to potential evaporation and some function of s o i l moisure, with the vegetation acting as a passive wick. However, the recent research reviewed by H a l l and Kaufmann (1975) demonstrates that decrease in humidity cause the stomatal apertures of leaves to contract. The amount of contraction depends on the kind of leaves and their temperature but can be enough to decrease transpiration. If t h i s should happen over a large enough area, the reduction in vapour flux to the atmosphere w i l l cause . a further decrease in humidity with resultant positive feedback. Therefore, low atmospheric humidity may be a cause of high potential evaporation, a cause of high resistance to transpiration and an effect of low evapotranspiration. Although the l a t t e r two functions are by far the most important in terms of the hydrologic cycle, they are so paradoxical, so confusing and so d i f f i c u l t to explain or assess that they are p r a c t i c a l l y unknown to hydrologists. This is unfortunate since they indicate that the relationship between evapotranspiration and potential evaporation should be negative 200 and may even be complimentary, a concept that is contrary to the assumptions used in watershed models. The existence of a p o t e n t i a l l y large number of feedback interactions in the soil-vegetation-atmosphere system would not be too serious i f they could be detected by testing the models. However, th i s i s precluded by the need for a large number of insubstantial c o e f f i c i e n t s (or parameters) that must be derived from h i s t o r i c a l data for each basin by mathematical or trial - a n d - e r r o r optimization. - For example, the optimization process and the pronounced seasonal patterns of radiation, temperature, vapour pressure d e f i c i t and vegetation development for any one basin can combine to mask the inadequacy of the conventional assumption that evapotanspirat ion i s proportional to potential evaporation. This i s evident in a discussion by Penman (1967) where i t is stated that "in terms of hydrology, one could be very wrong in the estimates of potential evaporation and s t i l l get the right answer i f one had the right value of D, in which D i s a s o i l moisture parameter that cannot be measured. With l i t t l e chance that unknown interactions or paradoxical cause-and-effect relationships can be detected in the v e r i f i c a t i o n process, the sophisticated watershed models have doubtful s c i e n t i f i c value." The above discourse points towards the p o s s i b i l i t y of flaws in our conceptualizations of the physical processes. To detect these flaws, we need to encourage " o r i g i n a l " research on hydrologic processes and ' t h e i r interactions. Therefore, the second approach to solve the problem i s a long-term research 201 project which w i l l have enough room for doing fundamental research on characterizing hydrologic processes and their interactions. Two d i f f e r e n t approaches to solve the problem have been discussed. The f i r s t is an engineering approach wherein we can formulate a model based on the present state of the art and use i t to predict the watershed response to man-made changes (a g r i c u l t u r a l drainage). This approach must y i e l d better solutions than those possible with the existing models because our characterization of the various hydrologic processes w i l l be more precise. The second approach i s more research-oriented. F i r s t we should increase understanding of the physical processes and the paradoxical nature of the interactions among them and then formulate a model. It involves re-evaluation and an in-depth study of each of the processes involved and a modification, i f necessary, of the assumptions i m p l i c i t in their mathematical charaterizations. Such an approach i s very demanding in terms of personnel, money, time and e f f o r t s in data c o l l e c t i o n but w i l l provide physical basis for operational and water management decisions. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo 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-0076973/manifest

Comment

Related Items