.5. 317...! talc 3ng; 52.": .. :L HAT. . . EL: 0.. 3 5 .¥ . ‘ . v. wank. ‘ . . . zihxfiw”. I, <‘_. q . . swdflrr £33.. .21 .. H. “as . x 3...: . m... aaefiwfiixz. .dI- 4.)‘l|. {’3' s £32. dflhrfit "dun. . 51%;“. a3» 33.. . II to: ‘36:: .5) «I $33.“? .. «.3 .15}. «a “H i $4... ha»: sari-sawinvx .02. 3...... :3 :...1t.; : tints}. .. i at italVOn i. 1.23.1...1. Illllllflilillillllfil This is to certify that the dissertation entitled DIGITAL TERRAIN ANALYSIS AND SIMULATION MODELING TO ASSESS SPATIAL VARIABILITY OF SOIL WATER BALANCE AND CROP PRODUCTION presented by Bruno Basso has been accepted towards fulfillment of the requirements for Ph.D Crop and Soil Sciences degree in 0 Major professor Date 5 I77“? 1—?” MS U is an Affirmative Action/Equal Opportunity Institution 0-12771 LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE ‘ DATE DUE | APR $45 [2me 3 m’sb’ififfi Wm FEB 03? 20077 $32—04: 23192: W SEW ZIP! moo chIRCIDateOuopfiS—p.“ DIGITAL TERR AIS . Elm VARLABHJT Ir | DIGITAL TERRAIN ANALYSIS AND SIMULATION MODELING TO ASSESS SPATIAL VARIABILITY OF SOIL WATER BALANCE AND CROP PRODUCTION By Bruno Basso A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Crop and Soil Sciences 2000 DIGHAL TERRAIN Silt U- HRHBIUT lmm charactcn ~. 733 ohm modify cm IIUI‘ station. a well as C I Learning saturated n 3::tttergence of both s 323133 a contentional. ( 292110 mluate the hxd ABSTRACT DIGITAL TERRAIN ANALYSIS AND SIMULATION MODELING TO ASSESS SPATIAL VARIABILIT Y OF SOIL WATER BALANCE AND CROP PRODUCTION By Bruno Basso Terrain characteristics and landscape position control soil physical properties. They often modify environmentally sensitive processes such as leaching, erosion and sedimentation, as well as dynamic factors affecting crop production. The likelihood of soils becoming saturated increases at the base of slopes and in the depression where there is a convergence of both surface and subsurface flow. The objective of this study was to combine a conventional, one-dimensional soil water balance model with a terrain analysis model to evaluate the hydrological and agricultural processes occurring on sloping land surfaces. A new digital terrain model, TERRAE-SALUS was developed to study and model how the terrain affects the vertical and lateral movement of water occurring on the land surface and in the shallow, subsurface regimes. This study evaluated the capability of TERRAE-SALUS applied at a field scale with rolling terrain where the soil water content was measured. The model was able to partition the landscape into an interconnected series of element network from a grid DEM. TERRAE-SALUS was evaluated using three different scenarios to gain a better understanding of the factors affecting the runoff-runon processes. The high elevation point consistently showed lower water contents compared to the upper and lower saddles and depressions. The subsurface lateral flow was highest on the saddles between two peaks, indicating the :rtst. perfarmzmce of if. ”m located on the p am: tamed from 0.22 i strap srmuiatton mode .. the ability of thc 522$} variable landscaj Elam modeling and 31531 inability. Mitch L1? correct performance of the model in predicting the contribution of water from the elements located on the peaks. The RMSE between measured and simulated soil water content varied from 0.22 cm to 0.68 cm. A second experiment was carried out applying the crop simulation model CROPGRO in combination with remote sensing data to evaluate the ability of the model to identify factors responsible for the yield variation in a spatially variable landscape. Results from this study showed that the combination of crop simulation modeling and remote sensing can identify management zones and causes for yield variability, which are prerequisites for zone-specific management prescription. To my family, my friends, and my love Valentina! Irish to express m) smc iia'temet Joe tnelxe )e. Sreemett. we became \ feyears I “as in East L." . | 2:: are no words able II c for many years. trm Cl muEedge in various SCI Feels-an let me make ti Inert to thank Dr. Fran . aserttttttee member an I 3132! “ant to extend a sri Em, . , "ment of Agricultu . v 1“‘:. ct m 1d lCe d Un ACKNOWLEDGMENTS I wish to express my sincere gratitude and appreciation to my advisor Dr. Joe T. Ritchie. I have met Joe twelve years ago when I came to MSU for the first time to learn English. Since then, we became very good friends and a special relationship developed throughout the years I was in East Lansing. Joe and Ann basically became my family and I think there are no words able to express my gratitude to them. I had the fortune to be next to Joe for many years, travel around the world with him, and learned not only the knowledge in various scientific fields but also their applications to real world problems. He always let me make the decisions and face the consequences, trusting in my skills. I want to thank Dr. Fran J. Pierce for his guidance, suggestions for my research, serving as committee member and for his friendship that developed throughout the years. I also want to extend a special word of appreciation and gratitude to Dr. Jim Jones of the Department of Agricultural and Biological Engineering at University of Florida, for accepting to be on my graduate committee, and for his valuable advises on crop modeling. Thanks to Dr. Jeff Andresen and Dr. Dave P. Lusch, members of my graduate committee for their tutoring in the field 'of Agricultural Climatology and Remote Sensing. They were excellent teachers and I will benefit from their lectures throughout my future career. A special thanks goes to Dr. John C. Gallant of CSIRO-Land and Water in Canberra, Australia for the development of the digital terrain model TERRAE. I have visited John in Canberra twice during my Ph.D. project and his time, help and hospitality were inestimable. in? mm ll'IEI‘tCIS dUT‘ Water. Bflfln BC‘L’J I H 3-...“ ed lelp, litanlgs Brian? “rip help and su-p’ rt cappuccmo to; lilac-unit to. extent m) 33m and lCRlS A7. CereofCDNlT and 'l 72:13; 10 All my N0“ 11? little of such profe lli'i \rlnocur and A)“. limb“ 0f pCOple u e: Wally like to thank l~ Enigma} Thanks to B7 mg me making l‘ll .. Lit bur net least. my ch. 53m . .‘Ltel'S, GlOVam-u‘ 21' I 'valf‘j’l‘. the fix her [Dye a“ I I made many friends during my stay at MSU, and the friendship with some of them will last forever. Brian Bear is one of them; I cannot ever thank him enough for his friendship and help. Thanks Brian! A very special friend that I would like to thank for his friendship, help and support is EdMartin. . Ed and I had many memorable conversations enjoying cappuccino together. Grazie fratello! I also wish to extent my gratitude for their support to the International Research Center of CINIMYT and ICRISAT. A special thanks goes to my dear friends Tim Reeves and Peter Grace of CIMMYT and to Prabhakar Patak, Piara Singh and Dr.S. Virmani of ICRISAT. Thanks to all my Nowlin Chair colleagues for their friendship. Sharing the office with colleagues of such professional and human dimensions as Carlos Paglis, Scott Piggott, Marta Vinocur and Ayman Suleiman was a privilege. A number of people were involved in the research project at the Durand site. I would especially like to thank Ricardo Braga for the memorable and fun days we had during the field study. Thanks to Brian Long and Cal Bricker for the field harvest and John Anibal for letting me making his field look like a Swiss cheese with my neutron probe access tubes. Last, but not least, my deepest gratitude goes to my parents, Francesco and Angela, and my brothers, Giovanni and Claudio, for their love and support. I gratefully thank Valentina for her love and my friends in Napoli for their inestimable friendship. vi Uni OF TABLES ........ Ill OFHGIRES ...... SST OF ABBREVIAII llltl’lIRl [\TROD7 ll Background an; 12 litpotheSts ........ 1.3 Objective .......... ~IAP’IIRII LlTERA' 2.1 Sntlwatcrbalar 22 Dtgtttu‘ Terrain 3.2.1 DigitalE. 2.22 Datasour 2.2.3 Spatial reI . I 32.4 Dlgllillli 2.2.5 landscap-I 13- ltttatn-basedh pa L» J 5—3 5/) l C (/2 ( '7‘ TABLE OF CONTENTS LIST OF TABLES ............................................................................................................ ix LIST OF FIGURES .......................................................................................................... x LIST OF ABBREVIATIONS ........................................................................................... xv CHAPTER] INTRODUCTION ..................................................................................... l 1.1 Background and Rationale ................................................................................... l 1.2 Hypothesis ............................................................................................................ 8 1.3 Objective .............................................................................................................. 10 CHAPTER II LITERATURE REVIEW ......................................................................... 12 2.1 Soil water balance modeling ............................................................................... 12 2.2 Digital Terrain Analysis ...................................................................................... 16 2.2.1 Digital Elevation Model ............................................................................ 17 2.2.2 Data source of Digital Elevation Model ................................................... 21 2.2.3 Spatial resolution and accuracy of Digital Elevation Model .................... 22 2.2.4 Digital Terrain Modeling .......................................................................... 24 2.2.5 Landscape topographic attributes .............................................................. 25 2.3 Terrain-based hydrological modeling ................................................................ 28 CHAPTER III TERRAE—SALUS ................................................................................... 33 3.1 TERRAE: A new method for element network .................................................. 33 3.2 SALUS and Spatial soil water balance model .................................................... 39 3.2.1 SALUS vertical soil water balance model ............................................... 40 3.2.2 Spatial soil water balance model ............................................................... 44 CHAPTER IV ASSESSING AND MODELING SOIL WATER BALANCE IN A SPACIALLY VARIABLE TERRAIN USING TERRAE-SALUS ........ 48 4.1 Introduction ........................................................................................................... 48 4.2 Methodology ......... . ............................................................................................... 54 4.2.1 Models description .................................................................................... 54 4.2.2 Model simulation ...................................................................................... 59 4.3 Results and discussions ......................................................................................... 62 4.3.1 Topographic attributes .............................................................................. 62 4.3.2 Model simulations ..................................................................................... 64 4.3.3 Validation Study—Scenario 3 .................................................................. 67 4.4 Conclusions ........................................................................................................... 7O vii LETTER V’ UNDER CROP .\ 5.1 Introduction ....... 52 Materials and in; 5.2.1 Site descr 5.2.2 Remote s. 5.2.3 Crop gI'O'. 5.2.4 Simulatt-t- 53 Results and dis. 5.3.1 Field me. 5.3.2 SimUlLiIlt 5.3.2.1 Seen,“ 5.3.2.2 .\'D\'j 5.4 Conclusion ....... Cam \I CONCI APPBDDI AWN... satin B .............. ' 3"“!- ‘ 0d CHAPTER V UNDERSTANDING SOYBEAN YIELD VARIABILITY USING CROP MODELS AND REMOTE SENSING ........................................ 99 5 - 1 Introduction ........................................................................................................... 99 5 -2 Materials and methods .......................................................................................... 102 5.2.1 Site description and field measurements ................................................. 102 5.2.2 Remote sensing data .................................................................................. 104 5.2.3 Crop growth model ................................................................................... 105 5.2.4 Simulation experiments ............................................................................. 106 5 -3 Results and discussion ........................................................................................ 106 5.3.1 Field measurements ................................................................................... 106 5.3.2 Simulation experiments ............................................................................. 109 5.3.2.1 Scenarios 1-4 ....................................................................................... 109 5.3.2.2 NDVI Classes ...................................................................................... 109 5.4 Conclusion .......................................................................................................... 110 CHAPTER VI CONCLUSIONS AND RECOMMENDATIONS ................................ 121 APPENDIX A ................................................................................................................. 126 APPENDDI B ................................................................................................................. 130 REFERENCES .................................................................................................................. 162 viii Tail: ”.1. Incal topogra Au 5 71:52:22. Resional top. Q -23. Catchment or 1- 5. ..\lodel inputs -n‘2a Sentiranogr. T21652b. Correlation T eixntint ICIIIIOIOIII- 5:53. Variables me. 1., 1‘54. Summary tab area weighted aver; LIST OF TABLES Table 2.1. Local topographic attributes .................................................... 26 Table 2.2. Regional topographic attributes ................................................ 27 Table 2.3. Catchment oriented topographic attributes .................................... 28 Table 5.1. Model inputs, number of model runs and RMSE for each simulation scenarios ....................................................................................... 112 Table 5.2a. Semivariograms parameters for the variables measured in the field. . . .. 113 Table 5.2b. Correlation matrix among variables measured in the field at the 52 grid-point ...................................................................................... 113 Table 5.3. Variables measured in each NDVI class and average for the 52 grid points ........................................................................................... 117 Table 5.4. Summary table for measured and simulated yield for each NDVI Class and area weighted average across the field. 118 ix E33241. Elevation at Figure 4.2. Slope for the 13334.3. Profile curt .: 51:34.4. Plan curvatu- E::e 4.5. Maximum r‘ Egtre 4.6. Location in .zitRiE for the Stud} 4 are 4.71 Soil water . :3: rairuall. no restricti :gre 4.7b. Soil water .‘ 1.21. High rainfall. no re fire 3.8a Cumulatm .. nitorm soil type. it Est-e ~ 117.2481). CUIDUlEIIthl we2193.33 surface . . rating soil laver i :17: 4. 10 11":11 -5Ut‘face RSmCting S 55:: 4. 11 S ‘3‘ finial] n HIT-21C C‘ IESIT‘ICtt cc 0 g 3011 lal'erl 55‘4".13& . “S- O >- 31713111. 1:218 : Cb 5&4 ”e ~13b' Sol] in ”all, no It: ate LIST OF FIGURES Figure 4.1. Elevation at l m grid resolution for the study area, Durand, MI ................ Figure 4.2. Slope for the study area, Durand, MI. .............................................. Figure 4.3. Profile curvature of the study area, Durand, MI ................................... Figure 4.4. Plan curvature of the study area, Durand, MI ...................................... Figure 4.5. Maximum ponding capacity for the study area in Durand, MI .................. Figure 4.6. Location in the landscape of the element network created by TERRAE for the study area in Durand, MI ...................................................... Figure 4.7a. Soil water content (0-26 cm) on day-1 for scenario I (l uniform soil type, high rainfall, no restricting soil layer) ............................................................ Figure 4.7b. Soil water content (26-77 cm) on day-1 for scenario I (1 uniform soil type, high rainfall, no restricting soil layer) ...................................................... Figure 4.8.a. Cumulative surface water flow out of the elements on day-l for scenario 1 ( l uniform soil type, high rainfall, no restricting soil layer) ................................. Figure 4.8.b. Cumulative surface water flow onto the elements on day-l for scenario 1 ( l uniform soil type, high rainfall, no restricting soil layer) ................................... Figure 4. 9. Net surface flow on day-1 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer) ............................................................................. Figure 4. 10. Surface ponding for day-1 for scenario I (l uniform soil type, high rainfall, no restricting soil layer) .................................................................. Figure 4. 11. Subsurface lateral flow on day-1 for scenario I (1 uniform soil type, high rainfall, no restricting soil layer) ............................................................ Figure 4. 12. Drainage on day-1 for scenario 1(1 uniform soil type, high rainfall, no restricting soil layer) ............................................................................. Figure 4. 13a. Soil water content (0-26) on day-2 for scenario I (1 uniform soil type, hlgh rainfall, no restricting soil layer) ........................................................... Figure _4. 13b. Soil water content (26—77) on day-2 for scenario 1 (1 uniform soil type. h‘gh 1' m“fall, no restricting soil layer) ............................................................ 71 71 72 72 73 73 74 74 74 74 76 76 77 77 78 78 Ezra 4.14. Subsurface 11:2 runfall. no ICSIIICI‘ Ear: 4.15. Drainage 0' reaming 5011 later E;tre 4.161 Soil IA ate: agranfall. no restnct' Egr: 4.161). Soil u ate: 1;; rainfall. no restrict Legs: 4.12. Sum of the s " . I xnigh rainfall. no re 334.18. Drainage 0? assisting soil l3\€l I E;.‘rei.l9a. Soil is ate: tees. high rainfall. u 11.” :53: 41%. Soil \\ ater -§t>.higlt rainfall. uttl’ Fire I 3;; 4~03 Runoff or. 'Eimc’dng soil late :33 ‘- 5'- :flhRunon on ....attting soil late Figure 4.14. Subsurface lateral flow on day-2 for scenario I (l uniform soil type, high rainfall, no restricting soil layer) ............................................................ Figure 4.15. Drainage on day-2 for scenario I (1 uniform soil type. high rainfall, no restricting soil layer). ........................................................................... Figure 4.16a. Soil water content (0-26) on day-7 for scenario 1 (l uniform soil type, high rainfall, no restricting soil layer) ............................................................ Figure 4.16b. Soil water content (26-77) on day-7 for scenario I (1 uniform soil type, high rainfall, no restricting soil layer) ............................................................ Figure 4.17. Sum of the subsurface lateral flow on day-7 for scenario I (1 uniform soil type, high rainfall, no restricting soil layer) ...................................................... Figure 4.18. Drainage on day-7 for scenario I (1 uniform soil type, high rainfall, no restricting soil layer) ............................................................................. Figure 4.19a. Soil water content (0-26 cm) on day-l for scenario 2 (3 different soil «types, high rainfall, with restricting soil layer at 120 cm) ...................................... Figure 4.19b. Soil water content (26-77 cm) on day-l for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ...................................... Figure 4.20a. Runoff on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) .............................................................. Figure 4.20b. Runon on day-l for scenario 2 (3 different soil types, high rainfall, with restricting soil layer) ........................................................................... Figure 4.21. Surface ponding on day-l for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ................................................... Figure 4.22. Net surface flow (cm) on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ................................................... Figure 4.23. Subsurface lateral flow on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ............................................. Figure 4. 24. Drainage on day-l for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) .............................................................. figure 4. 25a. Soil water content (0-26 cm) on day-2 for scenario 2 (3 different soil types high rainfall, with restricting soil layer at 120 cm) ...................................... xi 79 79 80 80 81 81 82 82 83 83 84 84 85 85 86 Eats 4.2513. 5011 \4 ate: “65.111231F31nl113l. \ttt" 'e 4.26. Subsurface :3 axial]. with restn. E33427. Drainage or it: resisting 5011 la} 1‘ Fig324.25a. Soil u ate.’ geshzgh rainfall. u 112’ E :e 4.2831 Soil u at." res high rainfall. in it? 19:3: 4.29.1 Suhsurfa.‘ .51. mull. truth restr. ..;:=4.29b. Sum of sul pea. ugh rainfall. Wll.’ '- “fin-a“ .5... 4.3021. Simulatec .anerent soil npes. ltl L“: 43-01). Simulatet: . l. ulitrent soil ttpes .5432 ' 2:41;] DTRInage 0 1.33.13 . Bit“? 3 Slmulatec Figure 4.25b. Soil water content (26-77 cm) on day-2 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ...................................... Figure 4.26. Subsurface lateral flow on day-2 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ............................................. Figure 4.27. Drainage on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) .............................................................. Figure 4.28a. Soil water content (0-26 cm) on day-7 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ...................................... Figure 4.28b. Soil water content (26-77 cm) on day-7 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ...................................... Figure 4.2%. Subsurface lateral flow on day-7 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ............................................. Figure 42%. Sum of subsurface lateral flow on day-7 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm) ...................................... Figure 4.30a. Simulated soil water content (0-26 cm) on day 184 (July 3) for scenario 3 ‘ (3 different soil types, low rainfall, with restricting soil layer at 120 cm) ................... Figure 43%. Simulated soil water content (26-77 cm) on day 184 (July 3) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm) ................ Figure 4.31. Subsurface lateral flow on day 184 (July 3) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm) ....................................... Figure 4.32. Drainage on day 184 (July 3) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm) ................................................... Figure 4.33a. Simulated soil water content (0-26 cm) on day 188 (July 7) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm) ................... Figure 4.33b. Simulated soil water content (26-77 cm) on day 188 (July 7) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm) ................. Figure 4.34. Subsurface lateral flow on day 188 (July 7) for scenario 3 (3 different soil type, 10W rainfall, with restricting soil layer at 120 cm) ........................................ Figure 4.35. Drainage on day 188 (July 7) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm) ................................................... xii 86 87 87 88 88 89 89 9O 91 91 92 92 93 93 We 4 36. Measured s i.;.i. .. R3437. Measured s 5?:rire438a. Error map taunt? 3 13 different 5 Ear: 4.38hMeasured season on the high elex; "3 4.393. Measurec 36—77 am for the hi gh Fi:;ie 4.3%. Measured iIiWiiii for the med; Er 4.403. Measured 257cm for the med; E53 44%. Measured ..'t~7‘.7cmgi for the med: 3 3‘ .QECd m, . rigid w Of Potentia » . | «a s i— .i we e alse‘COlOr a are . . _". IS NDVLi Fee: Figure 4.36. Measured soil water content (0-26 cm) on day 184 (July 3) .................... Figure 4.37. Measured soil water content (26-77 cm) on day 184 (July 3) .................. Figure 4.38a. Error map of soil water content (0-26 cm) on day 184 (July 3) for scenario 3 (3 different soil type, low rainfall, with restricting soil layer at 120 cm) ....... Figure 4.38b.Measured and simulated soil water content (26-77 cm) for the entire season on the high elevation point (peak) ........................................................ Figure 4.39a. Measured and simulated water content for the soil profile (0-26 cm) and (26—77 cm) for the high elevation zone (peak) for the entire season .......................... Figure 4.3%. Measured and simulated water content for the soil profile (0-26 cm) and (26—77 cm) for the medium elevation zone (upper saddle) for the entire season ............ Figure 4.40a. Measured and simulated water content for the soil profile (0-26 cm) and (26-77 cm) for the medium elevation zone (lower Saddle) for the entire season ........... Figure 4.40b. Measured and simulated water content for the soil profile (0-26 cm) and (26—77 cm) for the medium elevation zone (depression) for the entire season .............. Figure 4.41a. Water content change as function of elevation for 0-26 cm soil profile for day l 84 (July 3), day 224 (August 13) , and day 240 (August 28) ........................... Figure 4.41b. Water content change as function of elevation for 26-77 cm soil profile for day 184, day 219 (August 7), and day 248 (September 4) ................................. Figure 5.1. Krigged map of measured soybean yield (a); krigged map of LAI measured on August 8 (b); krigged map of plant population measured on July 4 (c) .................. Figure 5.2. Kn'gged map of soil depth (cm) (a); krigged map of clay content (%) (b); krigged map of potential extractable soil water (PESW-mm) (c). ........................... Figure 5.3a. False-color composite image taken on July 18, 1997 (original image deli vered by Emerge) ................................................................................ Figure 5.3b. NDVI map from false-composite image taken on July 18, 1997 ............... Figure 5.3c. Reclassified NDVI map from false-composite image taken on July 18, 1997. The areas in white is NDVI-Class l, the areas in grey is NDVI-Class 2 and the areas in black is NDVI-Class 3 ..................................................................... Figure 5 .4. Kriged map of simulated yield using measured plant population and measured soil type for the 42 grid point (a); measured yield for the three NDVI classes (b); SImUIated yield for the three NDVI classes using average measured inputs ............ xiii 94 94 95 95 96 96 97 97 98 98 114 115 116 116 116 119 .3 55 Correlatior - LU . F. F' . . igure 5.5. Correlation between NDVI (image taken on July 18) and yield... 120 xiv LIST OF ABBREVATIONS cm Centimeters GIS Geographic Information System km Kilometers MI Michigan M)( Mexico USDA United States Department of Agriculture USGS United States Geological Survey XV 1.1. Background 8 Most demographefb’ re chronically mal r challenges; it occurs m‘tdses the demar population pressure. CHAPTER I INTRODUCTION 1.1. Background and Rationale Most demographers agree that almost 800 million people of the world’s population today are chronically malnourished. In the developing nations, poverty presents increasing challenges; it occurs in more advanced countries as well. As the global population increases, the demand for food expands, with a growing diversity in diets. As a result of population pressure, the world’s finite resources are taxed to the limits by those same peOple whose existence depends on them. The World Resources Institute (1992) estimates that there are about a billion people, roughly one in every five, who must survive on less than the equivalent of one US dollar per day. Equally alarming is the deteriorating condition of the natural resources that underpin our current agricultural Pdeuction systems. We are now witnessing a never-before-seen rate of increase in the World population; nearly 200 new residents are added to this crowded planet every minute. Moreover, global food stocks, as a percent of utilization, are at their lowest level since we began keeping such records. A 13130“ b} the “ world is heading '- degmdation. and ‘ years alone. an CL} bzihon heetaresl h impaired through from soil loss b) e overgrazin g. de for a‘tn'itres that ha\ e heshu ater resourc: resources is in jeor hmares of land 311(- 5011.111ng POVCTI} CO’u‘ mes requires r. According to the U 7rd uu the key to (then mast have the incoh 3nd 'LhCn m\erSed I A report by the World Resources Institute (1992) presents data that clearly show the world is heading for an uncertain future if present trends in population growth, land degradation, and water and atmospheric pollution continues unabated. In the past 50 years alone, an equivalent to almost 11 percent of the Earth’s vegetated surface (1.2 billion hectares) has had its capacity to support productive forest and agricultural systems impaired through human-induced soil degradation. This degradation has resulted mainly from soil loss by erosion and from chemical and physical deterioration caused by overgrazing, deforestation, and inappropriate agricultural practices. Many of the same activities that have degraded soil resources have diminished the quality and availability of freshwater resources as well. Clearly, the long-term productivity of soil and water resources is in jeopardy. According to Rosensweigh and Hillel (1998) almost five billion hectares of land globally have been degraded in the past half century. Solving the poverty problem and, thus the food insecurity problem in developing countries requires rapid increases in food production, income, and employment. According to the United Nations, agricultural production-the engine of the devel0pment and the key to alleviating poverty-must be tripled within the next 50 years, and the people must have the income to buy it. The erosion of the natural resource base must be halted and then reversed. If the environmental degradation is to be curtailed, and if the food demands of a growing human population are to be met, agricultural land management PraCtices that sustain and enhance the long-term productivity of the natural resource base must be implemented. Successful implementation of these practices, referred to as sustainable land management (IBSRAM, 1991), will require quantitative evaluation of the factors that determine whether an agricultural system is sustainable or unsustainable. Only by idenufyin 3' long-term perform. task. The issue of \- transeends concern cultural. and politi. found to be sustain Ill-58081 of agrieu; h} increasing CrOp mmnml more C} from scientific TESe mndenr‘ and {her km We 31? inunfi that affeq them. C1 and soil nutrient pr disciplines to trans Only by identifying and measuring these factors will it become possible to evaluate the long-term performance of a given management practice. This, however, is not an easy task. The issue of what constitutes sustainable land management is complex and transcends concerns of a physical-chemical-biological nature to include socio-economic, cultural, and political concerns. Because of this complexity, a land management practice found to be sustainable at one site might not be equally sustainable at another site. The goal of agricultural research is to improve agricultural productivity and sustainability by increasing crop yields and by using inputs (water, fertilizer, labor, and farm machinery) more efficiently and at less cost to the farmer and the environment. Results from scientific research based on a single discipline are site specific and time and weather dcpendent, and they often require many experiments at different locations over several Years. We are inundated with information regarding site-specific crop yield and factors that affect them. Climatically driven crop growth simulation models which quantify water and 80i1 nutrient processes, integrate experimental results and knowledge from various disciNines to transfer the scientific knowledge to other sites and other years. These mOdCIS offer predictions and recommendations independent of location, season, crop, cultivar’ and management. Crop growth simulation models are increasingly used to SUPPOI't field research, extension, and teaching. The number of costly, multi -treatrnent, fime‘Consuming field trails can be substantially reduced by crop simulation as crop models can quantify the magnitude and variability in response to treatments (Ritchie et al,, 1989; Ritchie 1991; Jones and Ritchie, 1991). From 1970, when computers became easily available to help us deal with the complexity of crops, crop simulation modeling dei-eloped rapidly. ' and combined simr The other was seen ph}siological proce 19%). These No a; functional and meet- llechanrstic models- rro'olied as current. Darcy'slaw, and st Because of the large assumptions. mecha aei‘eropers. For this except for academic Fm‘nonal models 3313- ..3 to model then developed rapidly. Two distinct types of models emerged: one was essentially practical, and combined simple relationships and rules of thumb to predict the behavior of crops. The other was seemingly scientific in spirit, and sought to represent the biological and physiological processes thought to occur in plants and their environments (Passioura, 1996). These two approaches correspond to what Addiscott and Wagenet (1985) termed functional and mechanistic in their analysis of leaching models. Mechanistic models incorporate most fundamental mechanisms of the processes that are involved as currently understood. For example, soil water flow would be modeled using Darcy’s Law, and solute transport would involve mass flow and diffusion-dispersion. Because of the large amount of input information and the uncertainty of some assumptions, mechanistic models are usually not used by those other than their developers. For this reason, mechanistic models are seldom used for problem solving except for academic purposes. Functional models, on the other hand, represent the same processes but use simplified ways to model them. These types of models may be able to express a process as accurately as mechanistic models, although they use less input data and require much less cakilllation. Because of that, they can usually be used by those others than the developers WithOUt much difficulty. The best functional models might be thought of as containing rational empiricism to express rather complex relationships. The well-established CERES family of crop models (Ritchie et al., 1985) is predominantly functional rather than mechanistic, as a matter of fact they are built with the minimum data set concept (Nix, 1983). This minimum data set consists of information on weather, soil, crop genetic and m. agement. To model for several strategic planning analysis. and der’ir panicular cultix ar and development 4 pl ’ll grouth and t or’the [oral produ; dfficits on gm“ th 330W) and lOngqc consequences of Q The soil Water bala l . ,3... . «thing to the gm! wfibhshlng how rr it)" and how rr management. The CERES models were developed to provide users with an operational model for several purposes: assistance with farm decision making, risk analysis for strategic planning, within-year management decisions, large area yield forecasting, policy analysis, and definition of research needs. Models can predict the performance of a particular cultivar sown at any time, on any soil, in any climate. They estimate the growth and development of the crop, including the duration of each stage of growth, the rate of plant growth and they calculate the part of growth that becomes the economic component of the total production (Ritchie and NeSmith, 1991), the impact of soil water and nitrogen deficits on growth and development. Simulation models that accurately describe crop growth and long-term soil processes are also increasingly valuable tools for studying the consequences of global climate change (Adams et al., 1990). The soil water balance is calculated to evaluate the possible yield reduction caused by soil and plant water deficits, and the environmental impact of nutrients and pesticides leaching to the groundwater. The soil water balance model requires inputs for establishing how much water the soil will hold by capillarity, how much will drain out by gfaVitY. and how much is available for root uptake (Ritchie, 1985). Water content in any 3°“ laYer can decrease by soil evaporation, root absorption, or flow to adjacent layers. The 1iInits to which water can increase or decrease are input (volumetric fraction) for ”Ch Soil layer as the lower limit (LL) of plant water availability, the limit where capillary forces are greater than gravity forces, the drained upper limit (DUL), and for field saturation (SAT). The values used for these limits must be appropriate to the soil in the fireld, and accurate values are quite important in situations where the water input suoph is marginal hate been describe (1%. The CERES familj sails when the dra: such models is on? andrime. To use sr min. the spatial {11115.1 be addTCSSed There are three ma s 3 . Hortom'an overlant fateofthe soil; am onduetive 8011 Pre- acontinuum of Prt hm 1973). 1n . Precipitation fallin. supply is marginal (Ritchie, 1981). Determinations of these soil water extraction limits have been described in Ritchie et al., (1986), Ritchie and Crum (1989), Ritchie et al., (1999). The CERES family models have proven to be effective in simulating the water balance of soils when the drainage is vertical, often an unrealistic assumption. Runoff produced by such models is only from a point in space and there is no account for the water over space and time. To use such models for erosion estimates and for poorly drained, sloping terrain, the spatial and temporal relationship between various hydrological processes must be addressed. There are three main mechanisms that produce storm runoff: saturation overland flow that occurs when a rising water table intersects the soil surface, generating exfiltration; Hortoni an overland flow that occurs when the rainfall intensity exceeds the infiltration rate 0f the soil; and subsurface flow in which water flows laterally through a highly condUCtive soil profile (Horton, 1933; Dunne, 1970, 1983). These mechanisms are part of 3 continuum of processes and may operate singularly, but more often in combination (Freeze, 1972). In the case of saturation overland flow and Hortonian overland flow, Precipitation falling directly on the saturated zone at the soil surface produces surface runoff or overland flow. These saturated areas may occupy only a portion of a catchment and may vary in size depending upon soil properties such as saturated hydraulic conductivity, organic content, depth to restricting soil layers, antecedent soil water conttint, and topography. The Hortonian mechanism of runoff generation is most important in ser OC‘Clll'S. Marry hydrologr; nature of natural hydrologie proce hydrologic and e: power Of these m tithe catchment l Variation of inflot lhmtapMecon\3 ”all” impact 0n nnmnmnohenr and Sodimentation accounting for the tumble. Terrain ; Efitmomhological benefit proCeSScS l}, Etta] Elevathn \ t 3‘ ‘ J:‘pbgrapth ann- b I Pallal dlsmbu' important in semiarid and desert areas, and on agricultural land when surface sealin g OCC UI'S . Many hydrologic and water quality models crudely represent the three-dimensional nature of natural landscapes and therefore crudely represent spatially distributed hydrologic processes. As transport modeling becomes increasingly important in hydrologic and environmental assessment, this becomes a limiting factor in the predictive power of these models. Not only do we need to know the temporal variation in discharge at the catchment outlet, we also need to be able to accurately predict the temporal Variati on of inflow depths and flow velocities throughout the catchment. The effect of topographic convergence and divergence on flow characteristics in natural landscapes has a major impact on the values of hydrologic variables (Moore and Grayson, 1991). The tCrrain often modifies environmental sensitive processes such as leaching, erosion and sedimentation as well as dynamics factors affecting crop production. Without accounting for the terrain characteristics, accurate prediction of water quality is not POSSible. Terrain analysis is becoming increasingly important in the hydrological, geomorphological and ecological sciences for examining the spatial relationships bctween processes occuning on the land surface and the shallow subsurface regime. Digital Elevation Models (DEMs) provide the basic data for characterizing the topographic attributes of landscapes. A Digital Terrain Model (DTM), instead, includes the 8partial distribution of terrain attributes. A DTM is a topographic map in digital fon'nat, consisting not only of a DEM, but also the type of land use, settlements, types of didnage lines a processing. DT? Algorithms halt canilure. aspe proximity l'ana'f For closer COrTes ”55 0f panltlo christenzed by; 1~3~Hll>othesis The Ol'cra“ hYPOl EUHUOI $0” phVSlt Mon also deter M . “11‘1“ how lor mimics or ms 4‘33] SUbSUrfECe drainage lines and so on. They are a major constituents of geographical information processing. DTMs help to model, analyze and display phenomena related to topography. Algorithms have been developed to estimate slope gradient, plan curvature, profile curvature, aspect direction, specific catchment area, and a variety of surface drainage proximity variables from the DTMs (Moore et al., 1993; Bell et al., 1994). For closer correspondence in scale between model predictions and measurements, there is a need of partitioning landscapes into small areas where the hydrological processes, and the soil and topographic characteristics can be considered uniform or at least can be characterized by simple relationships. 13- HyPothesis The overall hypothesis of this study is that terrain characteristics and landscape position control soil physical properties through organic matter accumulation and formation of 50“ hOlizon and soil structure, which highly influence the soil water balance. Landscape POSitiOn also determines how much precipitation infiltrates into the soil profile because it reg‘llates how long water can pond on the surface, as well as how much can pond before it infiltrates or runs off to other areas in the landscape. Hortonian overland flow and latei‘al subsurface flow of water are most likely to be significant on the backslope, adding Water to soils in toeslope position. Thus, soil water is influenced by the terrain characteristics due to the effect of runoff-runon processes. In this stud.“ ‘ mmClTlClll at 2 mid consequt‘r water from 501 balance. Accur optimum mans not like a buck: influence of gr: the plants as my tones in soil. Water infiltratl n4 lateral movemen msenL This nest mienslly exceeds consequently Cau In this study, it is also hypothesized that the partitioning between vertical and lateral movement at a catchment level will help us to better predict the complete water balance and consequently the available water for the plants over space and time. Plants derive water from soil through their roots in an attempt to maintain a favorable hydraulic balance. Accurate evaluation of available soil water reservoir is vital to developing optimum management for rainfed crop production. However, the soil water reservoir is not like a bucket. Some water may percolate down out of the root zone under the influence of gravity. All water remaining in the root zone reservoir cannot be taken up by the plants as rapidly as needed because it is held too tightly by cohesive and adsoprtive forces in soil. Water infiltrating into the soil profile moves vertically and laterally. The subsurface lateral movement occurs when a low conductivity soil profile or shallow water table are present. This restriction in the soil forces the water to move laterally and if the rainfall intensity exceeds the infiltration rate an overland flow occurs, increasing runoff and consequently causing erosion problems. The horizontal movement of water varies with the soil properties and with the terrain attributes. Accurate prediction of the terrain characteristic will lead us to a better prediction of the water balance and water quality. Factors affecting both crop production and environmental sensitivity vary in both space and time. Topographic convergence and divergence of water flow characteristics in natural landscapes have a major impact on the values of hydraulic variables. movement 0 understandtn affecting crop 1.3. Objective The ObjCCll\‘C c balance model ' agricultural prot Amwmgmhe m5 I“Tam affect surface and m [h mducn‘lt)’ 501’] CH“ l 3 dptf 4 Comprfi wdopfllem Of d aeration models Chapter The combination of crop growth models and DTM will enable us to predict the movement of water over space and time at a catchment level and will provide a greater understanding of the environmental sensitive processes as well as dynamic factors affecting crop production. 1.3. Objective The objective of this study is to combine a conventional one-dimensional soil water balance model with more realistic terrain analysis to evaluate the hydrological and agricultural processes occurring on a sloping land surface. A new digital terrain model, TERRAE-SALUS, was developed to study and model how the terrain affects the vertical and lateral movement of water occurring on the land surface and in the shallow, subsurface regime, where a shallow water table or a low conductivity soil layer may exist. The thesis is constructed as follows: Chapter 2 comprehensively reviews studies pertaining to soil water balance modeling, the development of digital terrain analysis, topographic attributes derived from digital elevation models and terrain-based hydrological modeling; Chapter 3 presents TERRAE-SALUS: the new element network for deriving flow lines and for modeling the spatial variability of the soil spatial water balance; 10 Chapter 4 desen Durand. Mlchrg; Chapter 5 descn' lifflill) the taste agriculture conte Chapter6 glles t‘ Agricultural CCOS amplify the com] relationships occr other sites where adlanees in men lllfOpardy. Chapter 4 describes the application of TERRAE-SALUS at agricultural field scale in Durand, Michigan, USA to test the hypothesis of this research; Chapter 5 describes and evaluate the application of crop models and remote sensing to identify the factors responsible for the spatial variability of crop production in a precision agriculture context; Chapter 6 gives conclusions and implications of results from this dissertation. Agricultural ecosystems are very complex entities. The scope of this thesis is not to simplify the complexity of nature but to be able to explore and explain how complex relationships occurring on the dynamic environment can be described and transferred to other sites where the need for such understanding is highly beneficial. Technology and advances in research should be applied to properly manage the natural resources that are in jeopardy. ll soil uater be assessing app lhelHerature {OT Soil pater nSapphcation. D3933] Elel'am related applicat. CHAPTER II LITERATURE REVIEW This chapter reviews the finding of other scientists that have studied problems related to soil water balance, hydrological modeling and the use of digital terrain analysis for assessing appropriate management strategies of natural resources in agricultural sciences. The literature review is divided in four sections. Section 2.1 discusses previous attempts for soil water balance modeling. Section 2.2. describes Digital Terrain Analysis (DTA), its applications in environmental modeling and the topographic attributes derived from Digital Elevation Models (DBMS). Studies pertaining to hydrological modeling and GIS- related application are presented in section 2.3. 2.1. Soil Water Balance Modeling More than 99% of the water on earth is salty and the remaining fresh water is unevenly distributed. Humid regions have an abundance of it, so that frequently the problem is how to dispose of excess water. Arid and semiarid regions, on the other hand, are afflicted with a chronic shortage of it; and in some areas of arid regions, fresh is water is so limited that life is unbearable. Indeed, of all the major physico-chemical resources needed 12 by plants. water rater from prev resenorr affect r Because of the c models of the sc to interpret long Pastors at‘fectln; Ritchie ((1972) c in his model for incomplete cox-t mtitration and .931“ SOTglmm l mm at 3 “ate soil “are, COnte Ute e‘l'tlatl'On Wit by plants, water is the most limiting factor. The temporal variation in the supply of soil water from precipitation, and the spatial variability in the soil factors affecting soil water reservoir affect crop productivity and create a risky environment for growing crops. Because of the difficulty in obtaining long-term variability in yield, computer simulation models of the soil water balance and crop growth are necessary for agricultural scientists to interpret long-term productivity. Factors affecting the soil water balance have been exploited by several researchers. Ritchie (1972) developed a program to separate soil evaporation from plant evaporation in his model for calculating the daily evaporation rate from a crop surface with incomplete cover. The model provided close agreements between the simulated evaporation and the evaporation measured in the field using a weighing lysimeter in a grain sorghum plot. Richardson and Ritchie (1973) evaluated the soil water balance model at a watershed scale to study the effect of soil water content on runoff. The total soil water content was predicted accurately by substituting the simulated runoff needed in the equation with the runoff measurements available in the study. In 1985, Ritchie presented a model of the soil water balance that estimates the daily change in the storage capacity (AS) of the profile. This soil water balance model (SWBM), commonly called the Ritchie model, is a major component of the CERES (Crop Estimation through Resource and Environmental Synthesis) family of crop-soil- atmosphere models which are used in the Decision Support System for Agrotechnology Transfer (IBSNAT, 1986, Hoogenboom, 1998). Ritchie (1987) and Kovac’s et al. (1995) 13 used the Ritci applicable for (1993) used it able to elosel; cm i. Yates (1996) l tfl973) equalth the potential it saturation usi: study area. Ph based modelln; arid l’eé'lorls, ‘ ll. lung \r aided 0“ field It Gill} fol’ the V63] Soil water balan C ii“Ear models am 5&1 re ‘ ‘ l 5 EX {6] used the Ritchie soil water balance to simulate nitrate leaching. The model proved to be applicable for simulating the downward flow to the groundwater. Gerakis and Ritchie (1998) used the SWBM in the simulation of atrazine leaching and concluded that it was able to closely simulate the soil water content at three depths (13 cm., 26 cm., and 67 cm). Yates (1996) presented a water balance model (W atBal) that used the Priestly-Taylor (1972) equation to estimate potential evaporation. The WatBal model was used to assess the potential impact of climate change on a river basin runoff. The model required intense calibration using test data sets, a major limitation for interpreting the results outside the study area. Physically based models, like Ritchie (1972) are more adaptable to a GIS- based modeling environment because there is little needed calibration for individual sites and regions. Shanhoultz and Younos (1994) used a water balance approach to study the influence of tillage practices on soil water. In this study, evaporation was estimated with a model based on field measurements of pan evaporation. Results of this model were reported only for the years when the measurements were available, making the model’s usefulness limited. Soil water balance models can be stand-alone models as well as components of other larger models and their validation is often done through their evaluation. As matter of fact, there is extensive literature on the application and validation at different scales of 14 Ills? “filer b3 models. Diet Fed'des et al centenl and C use-ti Ritchie i sell sater and hater balance various pedOCl and a more aCC aehteted. This from bale or m: soils and comp: this study conch Stlll. uhile the ac the experimental H998) tested the the water balance components of the DSSAT models as well as other crop growth models. Diercks et al. (1988) used a soil water balance model (SWATRE) developed by Feddes et al. (1978) in conjunction with SUCROS crop model to estimate soil water content and crop yields under different irrigation strategies. This water balance approach used Ritchie (1972) model parameters. The model provided satisfactory results for both soil water and yield data. Gabrielle et al. (1995) evaluated the components of Ritchie’s water balance at a field scale. The model was tested against field data collected from various pedoclimatic regimes in France. The authors modified the drainage coefficient and a more accurate prediction of the soil water storage and surface water content was achieved. This was confirmed by comparing the model output against independent data from bare or maize (Zea mays L.) cropped conditions and for silt loam or sandy loam soils and compared the modified version of the model with the original one. Results from this study concluded that the original water balance model preformed accurately in sandy soil, while the accuracy of the simulations performed with the modified model fell within the experimental error in the measurements for silty-loam soils. Maraux and Lafolie (1998) tested the ability of a model, mechanistic with respect to soil-water flow and empirical for soil-plant and plant-atmosphere interactions, to predict soil-water balance components for long periods of time when input parameters are measured or estimated independently. A data set gathered in Nicaragua during several months was used for this purpose. Soil hydraulic properties were measured independently and parameters taken from the literature were used for plant processes modeling. The model predicted reasonably well the soil-water balance for maize (Zea mays L.) sorghum [Sorghum bicolor (L.) Merr.] sequence and for a grass sod. Singh et al., (1999) used Ritchie’s 15 model for ion Verttc Inceptl ohsened paint The soil water nater balance I use such mode necessar} to ad hydrological pr mt‘reasrngly imp Clarmning the s 12' Digital Te” ’ll kt . m3} k appI-Oac h! ”031mm ts at timer] , 001‘ atli WW. 1991 . l h‘anbeaPOu 0 \ model for long-term simulation of the water balance in a soybean-chickpea rotation on a Vertic Inceptisols. The authors observed that the simulated results were fairly close to the observed patterns. The soil water balance models discussed have proven to be effective in simulating the water balance of soils when the drainage is vertical, often an unrealistic assumption. To use such models for erosion estimates, for a poorly drained soil in a sloping terrain, it is necessary to address the spatial and temporal relationships between the various hydrological processes occurring on the landscape. Terrain analysis is becoming increasingly important in the hydrological, geomorphological and ecological sciences for examining the spatial relationship between processes occurring on the land surface and the shallow subsurface regime. 2.2. Digital Terrain Analysis Terrain analysis is the quantitative analysis of topographic surfaces. The purpose of a digital terrain system is the digital representation of terrain so that "real world" problems may be approached accurately and efficiently through automated means. Most attempts at modeling landscapes have been unsuccessful because the landscape was either looked at in little detail or the landscape was considered in two dimensions (Hall and Olson, 1991). Three-dimensional data patterns have a very high information content and can be a powerful vehicle for conveying essential landscape surface information. 16 lomnanhie at tunaiure can h erosion (Beret. al.1993d; Mo. al.. 1993; Wrist . properties of th tegetation (Mo. , I Basically. dlglld environmental r management dep 2.2-.1. Digital E Topographic attributes, including specific catchment area, slope, aspect, and plan curvature can be calculated and used to predict spatial patterns of soil water content and erosion (Beven and Kirkby, 1979; Moore et al., 1991; Moore and Wilson, 1992; Moore et al., l993d; Moore 1995; Wilson and Gallant, 1996); solar radiation estimation (Moore et al., 1993; Wilson and Gallant, 1998); spatial distribution of physical and chemical properties of the soil (Moore et al., 1993; Gessler, et al., 1995); spatial distribution of vegetation (Moore et al., 1993) and prediction of vegetation types (Brown, 1994). Basically, digital terrain analysis provides the basis for a wide range of landscape-scale environmental models, which are used for solving research-related problems as well as management decisions. 2.2.]. Digital Elevation Model There is a long history of studying landscape surfaces and an abundant knowledge and technology to measure topographic attributes have been developed. Digital Elevation Models are the source of the primary data used as a source of topographic surfaces information alone (Pike, 1988), for landscape modeling (Moore et al., 1991, 1993) as data layers in a GIS (Wiebel and Heller, 1991) and as ancillary data in remote sensing image analysis (Franklin, 1991). In principle, a digital elevation model (DEM) describes the elevation of any point in a given area in digital format. A discrete representation of a spatially continuous surface is merely a sample of values from the continuous surface. The sample is a finite set of spatial points with definite value (x, y, z) in a given 17 coordinate S) s“ sampled to pr:- suriaee is imp representatls'e surface. A dis.- is generated fr ll993bl stated t ' accurately re ' minimize d3: ' maximize dd coordinate system. A continuous surface has an infinite number of points that could be sampled to precisely represent the surface. Sampling the infinite points of the continuous surface is impractical and unnecessary; indeed a sampling method is used to extract representative points to build a surface model that approximate the actual continuous surface. A discrete sampling set of a continuous surface can still retain the continuity if it is generated from the original surface by following certain sampling procedure. ESRI (1993b) stated that a discrete surface model should: 0 accurately represent the surface; 0 minimize data storage requirements; 0 maximize data handling efficiency; The type of spatial surface dictates the representation and sampling method of the surface. No matter how smooth the landscape surface appears, it is not a mathematical surface, and cannot be represented using a single mathematical function. A landscape surface is a very particular continuous surface but no single mathematical function can be used to describe it. It is a product of the composition of many geological processes (faulting, erosion, and sedimentation). Geological young terrains typically have sharp ridges and valleys, in contrast to older terrains which have been smoothed by prolonged exposure to erosional forces (ERSI, l993b). There are three principal ways used to represent a surface in digital form: contour lines, arrays of equally spaced sample points, and irregularly spaced sample points (ESRI, 18 l993b). The {stored as C specified ele storage. the photographs manipulation is recorded 0: The most pop The surface i elevation moi- Single point a] Carter (1988) Surfaces, T0p( Imill“ is nor F has to be adju: 1min, I! is a] and efflClenu}. 1993b). The Vector or Contour line model describes the elevation of terrain by contours (stored as Digital Line Graphs, DGLs), the x,y coordinate pairs along each contour of specified elevation. Vector DEMs are based on the most common form of elevation data storage, the topographic map. Topographic maps are prepared directly from aerial photographs or field surveys so the information has undergone the minimum of manipulation, therefore minimizing errors. In the digital contour structure the elevation is recorded only once per contour string. The most popular way to represent a surface is an array of equally spaced sample points. The surface is represented by a "regular grid", or matrix, of elevation values. Gridded elevation models can be distributed as simple matrices of elevation, with the location of a single point and the grid spacing, implying the horizontal locations of all other points. Carter (1988) describes the methodologies for the digital representation of topographic surfaces. Topographic surfaces are non-stationary, more specifically, the roughness of the terrain is not periodic but changes from one land type to another. A regular grid therefore has to be adjusted to the roughest terrain in the model and be highly redundant in smooth terrain. It is apparent that, if one has to model these non-stationary surfaces accurately and efficiently, one must use a method that adapts to this variation. In response to this problem the Triangulated Irregular Network (TINs) was created. Tle are based on "coordinate random" but "surface specific" sample points. The location of these models would be dictated solely by the surface being modeled. By "surface specific" it is meant that they would be clustered in those regions of maximum roughness. In its most common form, the TIN is a set of irregularly spaced points connected into a network of 19 gdgesthatfo'." according to . I Thiessen pol}. planar. The tf pnll'lllTV ad\ Lil surface is rug;p where the sun include import along ridges 2 storage). The: l'nlilie a regt of the matrix Bill. 1975'). l dimenSion triangles). aI’Plications Wartime Emmorph 0. a1“ (1991). edges that form space-filling, non-overlapping triangles. The points are usually connected according to a Delaunay triangulation, a procedure that joins the centers neighboring Thiessen polygons (Delaunay, 1934; Thiesse 1911). The facets are usually assumed to be planar. The irregular nature of the TIN has many advantages and disadvantages. The primary advantage is variable resolution: a TIN can include many points where the surface is rugged and changing rapidly, but at the same time, only a few points in areas where the surface is relatively uniform. Another significant advantage is the ability to include important surface points (peaks, pits, passes, road and stream intersections, points along ridges and drains) at their exact locations (due to the precision of the coordinate storage). These advantages are countered by complexities in storage and manipulation. Unlike a regular grid, which provides an implicit neighborhood through the mechanism of the matrix, a TIN system would have to include this neighborhood explicitly (Peucker et al., 1975). Indeed, the location of every point in a TIN must be specified in the x,y, and 2 dimensions, as well as the topology of the points (the edges and adjacencies of the triangles). An intensive comparison between these three structures, together with applications of terrain analysis methods based on these structures for calculating topographic attributes and terrain-based indices of a variety of hydrological, geomorphological and biological processes is discussed by Kumler (1990) and Moore et al., (1991). 20 7 Dilltl St {5 ‘1 5.3. In principle. 2 a DEM data elelation moi frame points stream chann or stereo S] Complementa data in the fc “my al’alla poolished ml United States DEM data: 2.2.2. Data Source of Digital Elevation Models In principle, any data that contains the elevation information with location context can be a DEM data source. Practically, the main source of data for producing the digital elevation model are topographic contour lines, randomly distributed elevation points, the frame points of land surface such as peak, sinks, passes, points of change in slope, ridges, stream channels and shorelines, as well as stereoplotter data (e.g. stereo aerial—photo pair or stereo SPOT image pair) etc. Stereocorrelated DEMs are created from two complementary images, aerial photographs, or satellite images (Schenk, AF, 1989). Raw data in the form of stereo photographs or field survey (the accurate data source) are not readily available to potential end users of a DEM. Therefore, most users must rely on published topographic maps or DEMs produced by government agencies such as the United States Geological Survey (USGS). USGS produces several standard types of DEM data: 0 75-minute DEMs have a 30-by-3O meter point spacing in x and y; o 30-minute DEMs have 2-by-2 arc second point spacing, approximately 60-by-60 meter point spacing in x and y; o 1-degree DEMs have 3-by—3 are second point spacing, approximately lOO-by-IOO meter spacing in x and y. 21 3.2.3. Spatial The distance and }' horizon spatial resolu another 15. Sp increased. GI DEM. The n increasing th “N“ the pn processmg ti “"6 Optima broad €305 2.2.3. Spatial resolution and accuracy of digital elevation model The distance between two adjacent cells, or the geometric size of a cell or pixel in the x and y horizontal directions is called the spatial resolution of the DEM (or "grain"). The spatial resolution of a DEM is higher than another if its cell size is smaller than the another is. Spatial resolution is refined if cell size is decreased, or coarsened if cell size is increased. Generally, the finer the spatial resolution is the higher the accuracy of the DEM. The number of cells of a DEM covering a certain area will be increased when increasing the spatial resolution, and vice versa. The spatial resolution is very dependent upon the primary data used to produce the DEM, and the cost of computer storage and processing time. The optimum spatial resolution of a DEM is closely related to the spatial scale of the landscape pattern analysis and goo-modeling. For example, when soil properties with broad geographic extent are required, then a DEM with relatively coarse spatial resolution is indicated. To model detailed spatial distribution of soil properties, instead, a fine spatial resolution DEM will be needed. The topographic attributes computed from DEMs are dependent on the resolution of the elevation data from which they are computed. A regular grid is not an ideal representation of topographic surfaces for the study of scale effects. Gallant and Hutchison (1997) pointed out that when we subsample an elevation grid to obtain another grid at coarser resolution, beside the intended change in losing fine scale features of the surface, we also change the number of square cells into which the surface is divided. This is of particular importance when studying a "specific 22 catchment a: is important properties oil appropnatels' shape of feati The accurac source data 1 have been e- llte contour high. Howe could be it Sampling p Generally catchment area" that is computed by accumulating cell areas from adjacent cells. Thus, it is important not to confuse scale effects with grid effect if the objective is to study scale properties of a topographic surface. Gallant and Hutchinson (1997) suggested that to appropriately represent a topographic surface for the analysis of scale effects, the size and shape of features should be assessed at different scales. The accuracy of DEMs in representing the land surface is mainly dependent upon its source data spatial resolution (USGS, 1987). If we build the DEM from contour data that have been captured directly from aerial photographs as primary data using a stereoplotter, the contours are highly accurate (ESRI, l993b) and the accuracy of the DEM could be high. However, when the contours have been generated from point data, the accuracy could be lower because contours must be interpolated. A DEM usually uses discrete sampling points with raster structure to represent the relief of the landscape surface. Generally, it is difficult using discrete sampling points to represent every detailed feature and anomaly such as streams, ridges, peaks, and pits. Consequently, the higher the spatial resolution, the more detailed information content the DEM could represent and therefore the more accurate the DEM is. Conversely, a DEM with lower spatial resolution will miss more detailed information of the land surface. This generalization reduces the ability to recover position of specific features less than the interval spacing. Theoretically, for a given source data set, the only way to enhance the representation of detailed information of the landscape surface is to refine the spatial resolution of the DEM; as the spatial resolution is refined, there is an increasing likelihood that significant features of land surfaces will be represented. Nevertheless, it is not possible for a DEM to obtain more 23 detailed inft hon DEM the spatial rt surface \dl‘i surface usu: resolution \\ appropnate : POllllS. Ther on this stag SUIlEtCC that that are not Only becomi DEM. The E the PTOducti detailed information than that contained in the source data. Hutchinson (1996) shows how DEM resolution can be matched to information content of source data. Moreover, the spatial resolution of a DEM required to contain detailed information of a landscape surface varies with roughness characteristics of natural landscape surface. A rough surface usually needs a DEM with relatively fine resolution, while a coarse spatial resolution will be required by a smooth surface. After selecting the source data at the appropriate scale, the final stage is to interpolate the source data to a grid of elevation points. There are many choices here, and the quality of the DEM is critically dependent on this stage. General-purpose interpolation methods such as Kriging will produce a surface that is reasonably consistent with the data but may contain features such as sinks that are not really present in the real topography. They may also introduce biases that only become apparent when deriving terrain attributes such as slope and aspect for the DEM. The attention to shape and the drainage characteristics of the surface are critical to the production of a high quality DEM. 2.2.4. Digital Terrain Modeling Digital Terrain Models (DTM) have been used in geoscience application since the 1950s (Miller and Laflamme 1958). Since then, they have become a major constituent of geographical information processing. They provide a basis for a great number of applications in the earth and the engineering sciences. In 018, DTMs provide an opportunity to model, analyze and display phenomena related to topography. Indeed, DTMs include the spatial distribution of terrain attributes. The spatial distribution of 24 I topographt. vanahiltt} t 7.3.5. (and Landscape I mptesent [ht technology 1 SPCight (19 ‘31 al.. ( I99 primal}.- am tram Clevat Secondaryl “@355 that 0n the land; Th C ”lather can be fOUn Topographi EPOgraphiC ‘c topographic attributes can thus be used as a direct or indirect measure of spatial variability of these processes. 2.2.5. Landscape topographic attributes Landscape topographic attributes are spatial variables that are used to describe and represent the shape and pattern of the landscape surface. Digital terrain analysis and GIS technology provide tools to quantitatively define landscape attributes. Speight (1974) described over 20 attributes that can be used to depict landforms. Moore et al., (1991, 1993) also described terrain attributes and divided them into categories: primary and secondary or compound attributes. Primary attributes are directly calculated from elevation data and include variables such as elevation, slope, aspect, curvature etc. Secondary or compound attributes involve combinations of the primary attributes and are indices that describe or characterize the spatial variability of specific processes occurring on the landscape such as soil water content or the potential for sheet erosion. The mathematical representation of most attributes and the methods for calculating them can be found in Moore (1991, 1993), ESRI (1993), Gallant and Wilson (1996, 2000). Topographic attributes can also be divided in local, regional and catchment. Local topographic attributes are those that can be calculated from a small neighboring area surrounding the DEM cell using certain algorithm. The neighboring area is usually 3x3 25 cells. Table attributes is Regional to; lmcr geon: I topographic IOpOgtaphic Catchment 0 I0 ll'lC “hc characten stic Wain com't 10mg aphic z cells. Table 2.1 gives most of these attributes. The accuracy of the local topographic attributes is closely related to the spatial resolution of the DEM. Regional topographic attributes are those attributes that are calculated using considerably larger geometric area than the local topographic attributes (Table 2.2). The regional topographic attributes are less sensitive to the spatial resolution of the DEM than local topographic. Catchment oriented topographic attributes (Table 2.3.) are those attributes that are related to the whole catchment area, and are the measurement of certain catchment characteristics. The output value of the attribute at each DEM cell is calculated from certain combinations of all of DEM cells in the catchment. The catchment oriented topographic attributes have very little sensitivity to the spatial resolution of the DEM. Table 2.1. Local topographic attributes Attribute Definition Altitude Elevation above sea level Slope Maximum rate of change in elevation from each DEM cell Aspect Direction of the maximum rate of change in elevation from each cell DEM Surface curvature Measure of the surface convexity or concavity Profile curvature Curvature of a surface in the direction of steepest slope Plan curvature Curvature of a surface perpendicular to the direction of steepest slope Tangent curvature Plan curvature multiplied by the slope 26 Table 2.2 Regional topographic attributes Attribute Definition Upslope area Catchment area above a short length of contour Upslope slope Mean slope of upslope area Upslope height Mean height of upslope area Upslope length Mean length of flow paths to a point in the catchment Dispersal area Area downslope from a short length of contour Dispersal slope Mean slope of dispersal area Dispersal length Distance from a point in the catchment to the outlet Flow path length Maximum distance of water flow to a point in the catchment Specific catchment area Upslope area per unit width of contour Elevation percentile Ranking of the central point elevation compared to all the points in the surrounding region with a given area radius Elevation difference Difference between the central point elevation compared to all the points in the surrounding region with a given area radius Elevation deviation Elevation difference scaled by the standard deviation of elevation of the surrounding region with a given area radius Elevation standard deviation Standard deviation of the surrounding region with a given area radius Elevation semi—variance Two-dimensional semi-variogram of the surrounding region with a given area radius. It is an appropriate measure of the two-dimensional fractal dimension of the region 27 lahle 1.3 C selfml reg:- soil mm, Sand Slli a; Table 2.3 Catchment oriented topographic attributes Attribute g Definition Catchment area Area draining to catchment outlet Catchment slope Average slope over the catchment Catchment length Distance from highest point to catchment outlet Several researchers (Bell et al., 1992; Moore et al., 1993; Gessler et al.; 1995; Xu, 1999) have found high correlation between changes in these terrain variables and changes in soil drainage characteristics, A horizon depth, organic matter content, extractable-P, pH, sand, silt and soil taxonomic classes. 2.3. Terrain-based hydrological modeling In recent years, considerable progress has been made in the development of computer- based mathematical and computational techniques to model hydrological processes for various scales of analysis. GIS technology has become widely used in hydrological and water quality modeling. Hillslope hydrologists have long assumed that the downslope movement of water can be described by surface topography since gravitational potential largely dominates hydraulic gradients in steep terrains. Hence with the increased availability of DTMs, surface topography is driving many popular hydrological models (Moore et al., 1993; Vertessy et al., 1993; Gallant and Wilson, 1996;). Since the first 28 h computer . htdrologts: lcatchmcn: equations c for the ph computatior concept “a. {imposed a (lines of stet bl this nettt contour line: all elem-em C tumour line Elfimcms for lllfir CalChn- Elll'ernmen‘ ltginmng at SutCESSivelv has been USE mmensmnal hZ‘idlologic n he firm is nepressionle: computer-based model hydrologic models were developed in the early 1960’s, hydrologists have been attempting to use micro-scale process descriptions in meso-scale (catchment scale) hydrology. The massive computational effort required to solve equations describing processes in three dimensions and the intensive inputs requirement for the physically based model has limited the success of such models. However, computations may be reduced if the dimensions can be reduced from three to two. This concept was first applied by Onstad and Brakensiek (1968) and Onstad (1973). The proposed a flow net of gravitational potential between contours and their orthogonals (lines of steepest slope). Water was routed laterally down strips of land elements defined by this network and they termed this approach "stream path" or "stream tube". Adjacent contour lines and streamlines define irregularly shaped elements. Surface runoff enters an element orthogonal to the upslope contour line and exists orthogonal to the downslope contour line. Flow from one element can then be successively routed to downslope elements formed by the same stream tube. Moore et al., 1993 adopted this approach in their catchment partition model: TAPES-C (T opographic Analysis Programs for the Environmental Sciences-Contour. TAPES-C performs the partitioning of the catchment beginning at the contour line of lowest elevation and ending at the highest contour line, successively determining the elements for each adjacent pair of contour lines. TAPES-C has been used for distributed hydrological modeling that accounts for the effect of three- dimensional terrain on storm runoff generation. THALES (Grayson et al., 1992) is the hydrologic model that is coupled with TAPES-C. This DTM has two major limitations: the first is that it cannot handle depression for the flow network, thus requires a depressionless DEM, which is not a reality in many agricultural fields. The second 29 limitation available. t'anablcs a lhc TAPE TAPES-G : l SlaIlC mode gcncration ent'ironmen' 3Vailable in limitation is that the model is mechanistic, it requires several inputs that are often not available. Also, there is inconsistency in scale between the measurements of field variables and the way they are applied in the models. The TAPES model has also a grid version, TAPES-G (Gallant and Wilson, 1996). TAPES-G generates primary and secondary attributes from a DEM. It is considered a static model since it does not contains a dynamic water balance model. Through the generation of topographic attributes, TAPES-G has been applied in a variety of environmental modeling applications. For hydrological modeling, flow routing is available in TAPES-G with four different algorithms. Flow is routed from one cell to one and only one of its eight neighbor cells is based on the deepest descent. This algorithm called D8 produces parallel lines of flow along preferred directions. A second algorithm for flow directions (Rh08) aims to break up the parallel flow lines by introducing a random disturbance to the flow direction. The Rh08 algorithm is stochastic, indeed it produces a different flow network each time it is run. Flow dispersion is introduced in FD8 and FRh08, where the fractional amount of flow dispersed to each of the neighbors depends on the slope from the center cells to the neighbor. TAPES-G also computes the terrain wetness index (TWI), helpful in identifying areas of divergence and convergence based on the slope gradient. Where the slope gradient are low, the soil becomes wetter because the water is not removed to other downslope elements. Moore et al., (1988) found a strong correlation between this index and the distribution of surface soil water content. Gessler et al. (1995) found that the index, along with plan curvature, is a fairly good predictor of soil properties (A horizon depth, solum depth). 30 With a s dct'elopeti id salt t4 framework uscthchar lateral satti Ct'apOtranSp ”much a di‘ “Ill the Un Costly. lime guessed (Rt PTEdtct “'31 COmpared v mllOff b} Is With a similar approach of TAPES-C, TOPOG, an ecohydrological model, was developed by CSIRO in Australia to predict plant growth and the three-dimensional water and salt balance of heterogeneous catchment. Vertessy et al., 1993 describe the framework of this physically based, distributed parameter catchment model. The model uses Richard’s equations for vertical moisture flow, in multilayered soils, Darcy’s Law for lateral saturated flow, the convection-dispersion equation for solute transport, and evapotranspiration based on the Penman-Monteith model. Soil water extraction is through a distributed root system from the multilayered soil, and there is water exchange with the underlying aquifer system. The model demands significant input data that are costly, time consuming and difficult to measure, so most of the model inputs have to be guessed (Refsgaard et al., 1992). Vertessy et al., (1993, 1996) have used TOPOG to predict water yield from a mountain ash forest. Modeled and observed daily runoff compared well. Over the full period of simulation (12 years) the model overpredicted runoff by less 5%. Beven and Kirky, (1979) developed an hydrological model called TOPMODEL with the general thinking that variable source areas could be identified and the process of modeling basin hydrology be simplified, by summarizing the saturation potential, based on topographic position. Several other terrain-based overland flow, runoff and non-point source pollution model have been reported in the literature, including the TIN-based models of Jones et al. 31 (1990); g? H996). \V H998 l. The h}dro'. approach h] Nathan. 191 the concern ofthcm are SOll CVapom Present (Ritl dfsmbed inl (1990); grid-based models such as SHE (Abbott el al., 1986), MEDRUSH, Kirky et al., (1996), WEPP, Laflen et al., 1997, Cochrane and Flanagan, (1999), Wang and Hjelmfelt (1998). The hydrological models examined in this review were all physically based and such approach has come under scrutiny in recent years (Grayson et al., 1992 a, b, Grayson and Nathan, 1993). There is a considerable skepticism about their use in hydrology, because the concerns related to the scarcity of appropriate input and validation datasets. Also most of them are based on Richards equations for water flow, that can produce good results for soil evaporation, but it cannot predict plant evaporation as well when the root system is present (Ritchie and Johnson, 1990). An alternative to the models described above is described in Chapter 111. 32 This chapttv new methot elements ( C TERRAE-S CHAPTER III TERRAE-SALUS This chapter contains two sections. Section 3.1. presents the principles of TERRAE: the new methods for deriving flow lines and constructing a network of interconnected elements (Gallant, 1999). Section 3.2. discusses the spatial soil water balance of TERRAE-SALUS. 3.1. TERRAE: A new method for element network TERRAE (Gallant, 1999) constructs a network of elements by creating flow lines and contours from a grid DEM, that is the only required input. A flow line is a line of steepest descent down the surface that represents the flow of water. TERRAE can create contours at any elevation in the grid and does not rely on pre-defined contours. Each element created by TERRAE is an irregular polygon with contours as the upper and lower edges and flow lines as the left and right edges. The elements are connected so that the flow out of one element flows into the adjacent downslope element. A regular grid digital elevation model (DEM) provides the elevation data for TERRAE. Currently TERRAE reads floating point binary data exported from ARC/INFO using the GRIDFLOAT command. ANUDEM (Hutchinson, 1989), also available as TOPOGRID 33 ‘1‘”— within Al 01 COlllOU. TERRAE channel he uhtch of T To constru This is an e saddles are maximum. on Each 355 Peak TWO 1 mm at this the Saddle. 5 within ARC/INFO, is the recommended method for creating the DEMs from spot height or contour data. TERRAE is controlled by a parameter file that specifies the names of the DEM file, channel head file and boundary file. This parameter file also contains settings that control which of TERRAE’s functions are active and thresholds for sink clearance. To construct the element network, TERRAE first identifies all flat points in the surface. This is an essential feature for modeling water routing across the landscape. Peaks and saddles are recognized as critical points for computing flow network. A peak is a local maximum. A saddle is mixed extreme with a minimum along a ridgeline and a maximum on each associated drainage line. Topography is more complex near a saddle than near a peak. Two regions of convergent topography and two regions of divergent topography meet at this point. Topography is divergent in the upper part and convergent after turning the saddle. Slope lines turn sharply if they pass close to the saddle. TERRAE applies a user-specified boundary polygon, creates streamlines down from saddles and channel heads then adds ridge lines up from stream junctions and saddles. The resulting network of lines defines polygons that are then subdivided into elements to form the final element network. 34 Surface TERRAE. contours. ' results in . clctatton . All l‘ltll p0; are critical identified b steepest (165 m} be deft [NGEXER HUD) eaCh 5 Surface TERRAE uses a continuous surface derived from the DEM to construct its flow lines and contours. This surface is a quadratic B-spline with continuous first derivatives, which results in smoothly curved flow lines and contours. The surface exactly matches the elevation at every grid point in the DEM. All flat points in the spline surface are identified before creating any flow lines, as these are critical points defining the topology of the surface. Peaks, sinks and saddles are identified by locating exactly all points in the surface where the slope is zero. The lines of steepest descent and ascent from saddles are also determined. A boundary for the analysis may be defined using the supplied boundary file. This file is in ARC/INFO UNGENERATE format. TERRAE builds a polygon from the supplied lines. Sinks and depressions From each saddle, TERRAE follows a flow line down the surface either side of the saddle until it terminates at the edge of the DEM or at a sink point. Once all the lines connecting saddles to sinks are known the lowest draining saddle for each sink is determined by working upwards from the lowest connected saddle to the highest. A saddle is considered connected when it drains to the edge of the DEM. The sinks that are drained by connected saddles are then marked and the other saddles flowing into these sinks are also marked as connected. This process is repeated until every sink has a lowest 35 training 5 real. \l‘hc the flou l: side of the constructei determine the param 5:7?de Stream lt file. and 1 '1 T mIOLisl ,. ”“3 Crear Junctions draining saddle. Depressions around sink points are considered to be either spurious or real. When constructing further flow lines, spurious depressions are cleared by following the flow line from the sink up to its draining saddle and continuing down from the other side of the saddle. Real depressions remain as features in the surface, and elements can be constructed within them. The depth and distance from the sink to the draining saddle determines whether a depression is classified as real or spurious using values specified in the parameter file. Streams Stream lines are constructed down from the channel heads defined in the channel head file, and from real depressions out over their draining saddles. When a line approaches a previously created line (including the user-specified boundary line) it can connect to that line creating a stream junction. Streamlines can also connect at sink points. These stream junctions become part of the set of critical points defining the surface topology. For catchment modeling, the polygons defined from channel heads and stream junctions combine to form the modeling area. Catchment outlet points may also be provided; these are treated like stream junctions, so divide lines are constructed upwards either side of the outlet point. 36 the OUilli'l' l Ridge llllt' Ridgelines filial; 0r th bOUlldar}; ll ”Wards fro Cross, 30 the Smituilines. C tOllOVt‘g [he 1 CalChmentc LEnChanneled In agricultural applications channel heads will frequently not be used because the areas are smaller than natural first order catchment. In this case the boundary polygon provides the outline for the modeling area. Ridge lines Ridgelines are created from either side of each saddle and followed until they reach either a peak or the edge of the DEM. When a line approaches an existing ridgeline or the boundary, it can connect to it, creating a ridge junction. Ridgelines are also constructed upwards from stream junctions to form catchment divides. It is important that lines do not cross, so these catchment divide lines are started by interpolating between the incoming streamlines. When the interpolated line is sufficiently far from the streamlines, TERRAE follows the surface upslope as for other ridgelines. Catchment divide lines may also be created for channel heads to delineate the unchanneled contributing area for the channel head. Elements The network of lines created at this point forms a series of polygons covering the surface. TERRAE builds the polygons by tracing these lines from each stream junction, channel head and real depression. At each critical point (sink, saddle, peak, stream junction or 37 an C0} Tr. COl ridge junction) there may be several connected lines. TERRAE finds the next line in an anti-clockwise direction and follows each successive line until it returns to the starting point, forming a closed polygon. If a line terminates at the DEM edge, the polygon is considered to be open and is ignored. The polygons created by this procedure are converted to elements by locating the lowest and highest point along each polygon boundary, which then form the lower and upper boundary of the element. These points are considered to be contours with zero length. These elements can then be subdivided into smaller elements for hydrologic modeling by constructing contour lines and flow lines within the initial elements. (Figure). The subdivision of large elements into smaller elements can be done manually or automatically. The automatic method subdivides elements until they are smaller than a threshold area specified in the parameter file. The manual method displays the initial polygons and allows the user to interactively place contour lines and flow lines that subdivide elements. At each subdivision, the connections to adjacent elements must be determined. When an element is split by a contour line, the entire contribution of the upper element is directed to the lower element. When an element is split by a flow line, both the inwards contribution from the element above and the flow out to the elements below must be connected correctly. The proportion of flow is determined by the relative lengths of contour between the two elements. 38 At the end of the processing, a number of files are written. The most important of these are the files that describe the geometry, attributes and connectivity of each element. In each of these files, the elements are identified by a unique element number. The element geometry file (.elemgeom) contains the coordinates of the element boundaries. The element attribute file (.elemattr) specifies the slope, aspect, area, coordinates (x, y and z) of the centroid, upper and lower boundary lengths and other properties of each element. The element connectivity file (.elemconnect) contains the element numbers of the downslope neighbors of every element with the corresponding proportions of flow, and the sink number for elements that drain into real depressions. Depressions also need to be. described to permit modeling of ponding. The sink properties file (.sinks) contains the location and elevation of the sink, along with the relationship between depth and ponded volume and the number of the element it flows into when it drains over its draining saddle. 3.2. SALUS and Spatial Soil Water Balance Model Reference pertaining to the development and validation of the soil water balance model were cited in chapter 11 (Ritchie, 1972; Ritchie, 1985 and Ritchie, 1998). This section contains two main parts, the first part describes the principle of soil water balance as 39 llil sul Sit. described by Ritchie 1998 and some of the revisions made recently in the SALUS soil water balance; the second parts discusses the spatial components for the surface and subsurface lateral movement of water: the main core of this research. 3.2.1. Soil water balance model The one-dimensional (vertical) soil water balance model is calculated to account for soil and plant water stress at each point. The model calculates the profile water balance on a daily basis using the equation: dS/dt=P+I-R-Es-Ep-D where dS/dt = the change in water storage (S) in time period t P = precipitation I = irrigation R = runoff Es = evaporation from bare soil surface Ep = transpiration by plants D = drainage from the profile The soil water is distributed in several layers with fixed depths of : 40 ht Th: dec sat: The WI] The the surf lDa The 0-2 cm, 2-7 cm, 7-15 cm, 15-26 cm, 26-40 cm, 40-57 cm, 57-77 cm, 77-100 cm, 100-125 cm, 125-150 cm, 150-175 cm, 175-200 cm. Water content in any soil layer can decrease by root absorption, flow to an adjacent layer, or by soil evaporation in the case of layer 1. The input required by the model are the soil water limits to which water can increase or decrease (saturation, SAT; the drained upper limit, DUL, and lower limit, LL) and the saturated hydraulic conductivity (KSAT) for each layer. Definition and determination of these soil water extraction limits are described in Ritchie, 1998 and Ritchie et al., 1999. The use of KSAT has recently been introduced to calculate runoff based on the time-to- ponding approach instead of the Soil Conservation Service Curve Number (CN) method. The CN method was proven to be inadequate in representing variation in infiltration characteristics associated with differences in tillage and residue management. The soil surface KSAT varies as function of tillage, soil compaction, surface residue amounts (Dadoun, 1993) and it is the main parameter controlling the time-to-ponding curve. The time-to-ponding approach was first described by White et al., 1989. The application of the semi-empirical version of White equation used in the soil water balance model is discussed in chapter 4 of this dissertation. The time-to-ponding (TP) curves relate rainfall intensity to infiltration rate and define the point at which cumulative rainfall intensity exceeds the infiltration capacity of the soil (White et al., 1989, 1990; Chou, 1990), at which time water ponding in micro-depressions in the soil surface occurs. After ponding begins, infiltration is equal to the amount predicted by the TP curve as long as rainfall rate exceeds the infiltration capacity. When the rainfall rate becomes less than the 41 ah. infiltration capacity, rainfall plus the surface ponded water are infiltrated until the amount ponded is depleted. Runoff is predicted using a new methodology (Ritchie and Gerakis, personal communication). In this methodology a new parameter ("a") is introduced that varies with the time of the year. The a coefficient can vary for each month of the year and is obtained by calculating the slope for every month of the curve of cumulative rainfall (cm) on the Y axis, and the cumulative rain hours on the X axis. The a coefficient is then calculated using the following equation: a = EXP (-l.3*ln (S)-5.9 where S = slope of the curve of the cumulative rain (cm) vs cumulative rain time (hours) The slope of the curve of the daily runoff (cm hr'l) vs daily rain is described with the model: RS: EXP (a* KSAT) Finally the runoff (cm) can be estimated by the following equation: R: RS * (P-PM) 42 where R = Runoff (cm) P = Daily Precipitation (cm) PM = Ponding capacity (cm) Water is moved downward from the top soil layer to lower layers using a cascading approach. Water entering a layer in excess of the holding capacity of that layer (SATi- DUL i) is passed directly to the layer below by saturated flow. The drainage coefficient (K) is also been recently modified. The calculation for the Es and Ep are taken primarily from the work of Ritchie (1972) and by using Priestly-Taylor type equation (1972) instead of the Penman equation to calculate the potential evaporation (E0). E0 is calculated as function of the air temperature and solar radiation levels. Potential soil evaporation is a function of the potential evaporation E0 and the current leaf area index (LAD. LAI is the ratio of leaf area to ground area. As LAI increases, potential soil evaporation is decreased because of the protection from the wind and the shading from the leaf cover. The root water uptake routine has also been modified, but it is not described here. Ritchie, 2000 (unpublished data), discusses details of those modifications. 43 a. 4... ox 3.2.2. Spatial Soil Water Balance Model The element network created by executing TERRAE allows for the lateral movement of water across the landscape. Surface runoff and subsurface lateral movement is routed from one element to the next starting from the top element and moving downward. The spatial soil water balance model allows the presence of different soil types for each of the elements created if needed. The spatial routines initialized by reading information from the file "filenameelements" produced by TERRAE containing the element attributes. The element attributes are: the element number, the area of the element, the slope of the element, the X,Y, and Z coordinate of the center of the element and the topology (the connections of the elements). Surface Runoff The daily time loop is initiated by reading the weather file and by calculating the soil water balance for the downward flow for each of the element. The surface runoff produced by each element will move laterally to the next downslope element. The amount of surface runoff is calculated by multiplying the surface runoff of the upslope element by the area of the element. This amount of water will be added onto the next downslope elements as additional precipitation. If there is no downslope element, the surface water runs off to field outlet. Subsurface lateral flow The downward flow is calculated by introducing a correction factor to account for the slower flow that occurs at the deeper layers. The correction factor consists in separating the KSAT variable into a KSAT for the effective vertical flow (KSAT-Vert) and a KSAT for the saturated flow (KSAT-Macro). The correction factors is calculated as follow for the various depth: KSAT-Vert at 0 em = KSAT-Macro * 1 + 0*SIN (slope) KSAT-Vert at 50 cm = KSAT-Macro * 0.75 + 0.25*SIN (slope) KSAT-Vert at 100 em = KSAT-Macro * 0.50 + 050me (SIOpe) KSAT-Vert at 150 cm = KSAT-Macro * 0.25 + O.75*SIN (slope) KSAT-Vert at 200 em = KSAT-Macro * 0 + 1*SIN (slope) At the soil surface there is no difference between KSAT-Vert and KSAT-Macro, thus there is not need for a correction factor. At 200 cm the correction factor will be the SIN of slope, creating a lower conductivity. The subsurface lateral flow is computed using the following equation: SLF= Kef (dH/dx) * (Aup/Adw) Where 45 SFL = Subsurface lateral flow (cm day") Kef = Saturated hydraulic conductivity calculated as harmonic mean between Ksat of the upslope element and the downs10pe element (cm day") dH = distance between the saturated layer and the soil surface dx = distance between the center of the upslope element and the downslope element Aup = area of the upslope element (m2) Adw = area of the downslope element (m2) The hydraulic head (dB) is calculated by the soil water balance model, while dx is calculated by TERRAE. The Kef is calculated as harmonic mean as follows: Kef = 2/ (l/KSATup + l/KSAwa) where KSATup = KSAT of the upslope element KSAwa = KSAT of the downslope element The subsurface lateral flow occurs only when the SIN of the slope is greater than 0. The subsurface lateral flow is added to the next downslope elements into the layer that has the capacity to take it in, starting from the bottom layer and moving upwards. 46 If the SIN of the s10pe is zero and the amount of water is greater than the KSAT-Macro, then the water backs up within the same element. Subsurface lateral flow will occur again if a hydraulic head is created; dH/dx is then equal to the distance between the saturated layer of the upslope element and the saturated layer of the downslope element. Applications of the spatial soil water balance described above are presented in Chapter 4. The routines and the codes for the spatial soil water balance simulation are given in Appendix B. 47 CHAPTER IV ASSESSING AND MODELING SOIL WATER BALANCE IN A SPATIALLY VARIABLE TERRAIN USING TERRAE-SALUS 4.1. Introduction Water has been long known to be essential for plants’ life and at the same is one of the most limiting factors for their growth. In many agricultural regions of the world, the supply of water is highly variable due to variations that occur spatially and from year to year. One of the most important properties of the soil is that it is a reservoir for water. Without access to such a reservoir, most plants would not survive periods between rains. The factors that affect the soil water content include (1) soil characteristics: soil water limits (saturation-SAT, drained upper limit-DUL, lower limit-LL), saturated hydraulic conductivity (KSAT), thickness of the hydrologically active zone; (2) topography: local slope (a measure of the hydraulic gradient), specific catchment area (a measure of the potential maximum water flux), plan curvature (a measure of the rate of flow convergence and divergence), profile curvature (a surrogate for the rate of change of hydraulic gradient), and aspect and topographic shading (together with slope these influence the amount of solar radiation and in turn, evapotranspiration); (3) vegetation: variation in surface cover and water use characteristics; and (4) weather: net rainfall, net radiation, wind, and temperature (Moore, et al., 1991; Barling et al., 1994). 48 unr and E511 ant The 311:1 lher L’tar Hon, rateC (Grim 5 The ability to characterize the spatial variability of soil water content is of major importance. Models that consider the dynamics of soil water balance and crop growth have been extensively used to quantify the risk related to the uncertainty in water supply (Ritchie 1994, Jones and Ritchie, 1996). The CERES family models have proven to be effective in simulating the water balance of soils when the drainage is vertical, often an unrealistic assumption. Runoff produced by such models is only from a point in space and there is no account for the water over space and time. To use such models for erosion estimates and for poorly drained, sloping terrains, the spatial and temporal relationship between various hydrological processes must be addressed. The water infiltrating into the soil profile moves vertically and laterally. The lateral movement occurs when a low conductivity soil profile or shallow water table are present. This restriction in the soil, forces the water to move laterally and if the rainfall intensity exceeds the infiltration rate an overland flow occurs, increasing runoff and consequently causing erosion problems. The horizontal movement of water varies with the soil properties and with the terrain’s attributes. There are three main mechanisms that produce storm runoff: saturation overland flow that occurs when rising water tables intersect the soil surface, generating exfiltration; Hortonian overland flow that occurs when the rainfall intensity exceeds the infiltration rate of the soil; and subsurface flow in which water flows laterally through a highly conductive soil profile (Horton, 1933; Dunne, 1970, 1983). These mechanisms are part of a continuum of processes and may operate singularly, but more often in combination 49 the pit; lUlli COll ‘L i. llCCl SlE’ll \I‘OH uni Sun the = T. (Freeze, 1972). In the case of saturation overland flow and Hortonian overland flow, precipitation falling directly on the saturated zone at the soil surface produces surface runoff or overland flow. These saturated areas may occupy only a portion of a catchment and may vary in size depending upon soil properties such as saturated hydraulic conductivity, organic content, depth to restricting soil layers, antecedent soil water content, and topography. The Hortonian mechanism of runoff generation is most important in semiarid and desert areas, and on agricultural land when surface sealing occurs. Hortonian overland flow and lateral subsurface flow of water is most likely to be significant on the backslope adding water to soils in toeslope position. Thus, soil water is influenced by the terrain characteristics due to the effect of runoff-runon processes. Subsurface storm flow is generally considered to occur as lateral movement of water in the upper soil layers. Van de Griend and Engman (1985) reviewed the reasons for this and reported the influence of hard pans (plow pans) and impeding layers. When subsurface flow converges, the capacity of the soil to transmit the flow is exceeded and saturated areas are formed. These saturated areas are impermeable so in addition to exfiltration, all rainfall on them becomes runoff. Many hydrologic and water quality models crudely represent the three-dimensional nature of natural landscapes and therefore crudely represent spatially distributed hydrologic processes. As transport modeling becomes increasingly important in hydrologic and environmental assessment, this becomes a limiting factor in the predictive power of these models. Not only do we need to know the temporal variation in discharge 50 at the catchment outlet, we also need to be able to accurately predict the temporal variation of in flow depths and flow velocities throughout the catchment. The effect of topographic convergence and divergence and divergence on flow characteristics in natural landscapes has a major impact on the values of these hydrologic variables (Moore and Grayson, 1991). Topography can also affect the location of zones of surface saturation and the distribution of soil water content (i.e., the soil water content overlying an impermeable or semi permeable soil horizon at shallow depth). The likelihood of soils becoming saturated increases at the base of slopes and in depressions where there is a convergence of both surface and subsurface flow (Kirkby and Chorley, 1967; Moore et al., 1988a). Hall and Olson (1991) determined the effects of landscape morphology on soil physical and chemical properties and soil water movement across the landscape. Without accounting for the terrain characteristics, accurate prediction of the soil water balance was not possible. The automation of terrain analysis and the use of Digital Elevation Models (DEMs) has made it possible to easily quantify the topographic attributes of the landscape and to use topography as one of the major driving variables for many hydrological models. These topographic models, commonly called Digital Terrain Models (D'I‘Ms), partition the landscape into a series of interconnected elements based on the topographic characteristics of the landscape and they are usually coupled to mechanistic soil water balance models (Moore et al., 1993 Grayson, etc al., 1993, Kirkby et a1, 1979; Vertessy et al.,l996). These DTMs have two major limitations: the first is that they cannot handle depression for the flow network, thus requiring a depressionless DEM, which is not a 51 reality in many agricultural fields. These DTMs were designed for large-scale applications and for quantifying water quality running into streams, thus the sinks and depressions are filled to have a continuos flow of water down to the streams. The second limitation is that the mechanistic soil water balance models require several inputs that often are either not available, costly, time consuming and difficult to measure, so in most of the cases they have to be estimated (Refsgaard et al., 1992). Also, there is inconsistency in scale between the measurements of field variables and the way they are applied in the models. There is a considerable skepticism about their use in hydrology, because the concerns related to the scarcity of appropriate input and validation datasets. Also, most of them are based on Richards equations for water flow, that can produce good results for soil evaporation, but it cannot predict plant evaporation as well as water extraction from the root system (Ritchie and Johnson, 1990). The idea of creating a DTM that would include the topographic effect on the soil water balance and would be coupled with a functional soil water balance to spatially simulate the soil water balance became clear from the reasons mentioned above. This lead to the development of TERRAE-SALUS, a DTM for predicting the spatial and temporal variability of soil water balance (Chapter 3). The model requires a DEM for the creation of the element network for landscape partitioning, weather and soil information for the soil water balance simulation. 52 Soil information required by the model include (SAT, DUL, LL, KSAT). These parameters can be obtained through measurements or estimated using empirical equations (Ritchie et al., 1999). For an accurate soil water balance simulation, the depth of the lowest KSAT should also be included. Indeed, significant contributions of water from subsurface lateral flow occur in saturated conditions. This information is difficult to obtain, but it can be estimated based on topography, or using historical information on the site (knowing the dry and wet areas across the field). Soil surveys could also be helpful providing information on the drainage characteristics of the soil (i.e. soil poorly drained indicates the presence of shallow water table or low conductivity layer). If no information is available on the soil, crop data could be used as indicators of stresses through remote sensing, where imagery interpretation can help identify areas to be sampled and determine soil information necessary for the model. The overall hypothesis of this study is that the terrain characteristics and landscape positions control soil physical properties through organic matter accumulation, formation of soil horizons and soil structure that highly influence the soil water balance. Landscape position also determines how much precipitation infiltrates into the soil profile and for how long water can pond on the surface, as well as how much water can pond before it infiltrates or runs off to other areas in the landscape. In this study, it is also hypothesized that the partitioning between vertical and lateral movement at a field-scale level will help us to better predict the complete soil water balance and consequently the available water for the plants over space and time. Accurate 53 predictions of the terrain characteristics will lead to better predictions of the soil water balance. The objective of this study is to combine a conventional one-dimensional soil water balance model with a more realistic terrain analysis to evaluate the hydrological and agricultural processes occurring on a sloping land surface. A new digital terrain model, TERRAE-SALUS, was developed to study and model how terrain affects the vertical and lateral movement of water occuning on the land surface and in the shallow, subsurface regimes, where a shallow water table or a low conductivity soil layer exist. This study evaluates the capability of TERRAE-SALUS applied at field scale with rolling terrain where the soil water content was extensively measured. The model was evaluated using three different scenarios (Scenario 1, Scenario 2, Scenario 3) to gain a better understanding of the factors affecting the runoff-runon processes occurring on the landscape. 4.2. Methodology 4.2.1. Models description Digital Terrgin Model: mAE TERRAE is a new method for creating element networks where landscape depressions are included. TERRAE constructs a network of elements by creating flow lines and 54 contours from a grid DEM. A flow line is a line of steepest descent down the surface that represents the flow of water. TERRAE can create contours at any elevation in the grid and does not rely on pre-defined contours. Each element created by TERRAE is an irregular polygon with contours as the upper and lower edges and flow lines as the left and right edges. The elements are connected so that the flow out of one element flows into the adjacent downslope element. A regular grid digital elevation model (DEM) provides the elevation data for TERRAE. To construct the element network, TERRAE first identifies all flat points in the surface. This is an essential feature for modeling water routing across the landscape. Peaks and saddles are recognized as critical points for computing flow network. TERRAE then applies a user-specified boundary polygon, creates streamlines down from saddles and channel heads then adds ridge lines up from stream junctions and saddles. The resulting network of lines defines polygons that are then subdivided into elements to form the final element network. These elements can then be subdivided into smaller elements for hydrologic modeling by constructing contour lines and flow lines within the initial elements. The subdivision of large elements into smaller elements can be done manually or automatically. The automatic method subdivides elements until they are smaller than a threshold area specified in the parameter file. The manual method displays the initial polygons and allows the user to interactively place contour lines and flow lines that subdivide elements. 55 At the end of the processing, a number of files are written. The most important of these are the files that describe the geometry, attributes and connectivity of each element. In each of these files, a unique element number identifies the elements. The element geometry file (.elemgeom) contains the coordinates of the element boundaries. The element attribute file (.elemattr) specifies the slope, aspect, area, coordinates (x, y and z) of the centroid, upper and lower boundary lengths and other properties of each element. The element connectivity file (.elemconnect) contains the element numbers of the downslope neighbors of every element with the corresponding proportions of flow, and the sink number for elements that drain into real depressions. Depressions also need to be described to permit modeling of ponding. The sink properties file (.sinks) contains the location and elevation of the sink, along with the relationship between depth and ponded volume and the number of the element it flows into when it drains over its draining saddle. Spatial Soil Water Balance Model The element network created by executing TERRAE allows for the lateral movement of water across the landscape. Surface runoff and subsurface lateral movement is routed from one element to the next starting from the top element and moving downward. The 56 spatial soil water balance model allows the presence of different soil types to a maximum equal to the number of the elements created. Basically, each element created could have different soil characteristics if necessary. The spatial routine is initialized by reading information from the file "filenarne.elements" produced by TERRAE containing the element attributes. The element attributes are: the element number, the area of the element, the slope of the element and the x, y and z coordinates of the center of the element and the topology (the connections of the elements). The daily loop is initiated by reading the weather file and by calculating the soil water balance for the downward flow for each of the element. The surface runoff produced by each element is moved laterally to the next downslope element. The amount of surface runoff is calculated by multiplying the surface runoff of the upslope element by the area of the element. This amount of water is added onto the next downslope elements as additional precipitation. If there is not a downslope element, the surface water runs off to the field outlet. The downward flow is calculated by introducing a correction factor to account for the slower flow that occurs at the deeper layers. The correction factor consists in separating the KSAT variable into a KSAT for the effective vertical flow (KSAT-Vert) and a KSAT for the saturated flow (KSAT-Macro). The correction factors are discussed shown in Chapter 3. 57 The subsurface lateral flow is computed using the following equation: SLF= Kef (dH/dx) * (Aup/Adw) where SFL = Subsurface lateral flow (cm day'l) Kef = Effective saturated hydraulic conductivity calculated as harmonic mean between Ksat of the upslope element and the downslope element (cm day") dH = distance between the saturated layer and the soil surface dx = distance between the center of the upslope element and the downslope element Aup = area of the upslope element (m2) Adw = area of the downslope element (m2) The hydraulic head (dH) is calculated by the soil water balance model, while dx, the distance, is calculated by TERRAE. The effective saturated hydraulic conductivity (Kef) is calculated as a harmonic mean as follows: Kef = where 2/ (l/KSATup + 1/KSAwa) KSATup = KSAT of the upslope element ( cm day") KSAwa = KSAT of the downslope element (cm day") 58 The subsurface lateral flow occurs only when the sine of the slope is greater than zero. The subsurface lateral flow is added to the next downslope elements into the layer that has the capacity to take itin, starting from the bottom layer and moving upwards. If the sine of the slope is zero and the amount of water is greater than the KSAT-Macro, then the water backs up within the same element. Subsurface lateral flow will occur again if an hydraulic head is created; dH/dx is then equal to the distance between the saturated layer of the upslope element and the saturated layer of the downslope element. 4.2.2. Model simulation The first simulation run of TERRAE-SALUS was done using a single, uniform soil type with no restricting soil layer for the entire area with a high rainfall (76 mm) occuring on the first day. This simulation done was chosen to demonstrate the ability of the model to partition the vertical and horizontal subsurface flow. The second simulation run of TERRAE-SALUS used three different soil types with a low conductivity layer (KSAT=0.01 cm hr-l) at 120 cm. The soil types were a shallow sandy soil for the high elevation zones and peaks; a medium sandy-loam for the medium elevation zones and saddles areas; and a loamy soil for the low elevation areas and depressions. The rainfall was the same as scenario 1 (76 mm). This scenario was selected to have a direct comparison with scenario 1 but with a restricting layer at 120 cm that altered the vertical flow. 59 The final simulation run of TERRAE-SALUS was done to perform a model validation at field scale. Similar to scenario 2, the model was set up using three different soil types with low conductivity layer (KSAT=0.01 cm hr-l) at 120 cm. The soil types were a shallow sandy soil for the high elevation zones and peaks; a medium sandy-loam for the medium elevation zones and saddles areas; and a loamy soil for the low elevation areas and depressions. Field measurements of profile soil water content were taken on a three ha portion of a field located 10 km south of Durand, M1, to compare with model values. The field was planted with soybeans on May 5, 1997. A digital elevation model (DEM) was created for the site using a high accuracy differential global positioning system (DGPS) at 1 m grid resolution (F.J Pierce and T.G.Mueller, personal communication, 1997). Using the DEM, the following topographic attributes were determined for the site: elevation, slope, plan curvature and profile curvature. A regular grid consisting of 28 grid locations spaced 30.5 m apart was imposed on the experimental area. Latitude, longitude and elevation of each grid points were determined with DGPS. Neutron probe access tubes were installed at each of the 28 grid locations. A neutron moisture gauge was used to measure the spatial variability of soil water content at 15-cm increments to the depth of the C horizon or a maximum of 150 cm depth, which ever occurred first. Neutron probe calibration was obtained by filling four large metal cylinders with soil collected from two different location in field. The locations were 60 selected based on the soil type. Two cylinders were filled with the predominant soil type present in the field (sandy loam), and the other two were filled with the sandy soils that characterized the high elevations points and the peak. Each cylinder was carefully filled reproducing the field bulk density. For each soil type, one cylinder was filled with air-dry soil and the other was brought to saturation. This methodology is the most appropriate for neutron probe calibration. It allows the correct determination of the slope of the lines that joins the the driest and wettest point for that soil. It also decreases the errors obtained by fitting a curve through a clouds of points as observed from the traditional field methods measurements. Measurements on soil water were taken on a weekly basis throughout the season. During the installation of the neutron probe access tubes, soil samples were taken at the intersection of each of the 28 grid points in 25-cm increments and stored for analysis. Soil samples were air-dried and passed through a 2 mm sieve. Particle size was determined for each segment of each soil profile using the hydrometer method (Gee and Bauder, 1986). The upper and lower limit of soil water availability was determined using soil water measurements taken in the field, and from empirical equations based on soil texture (Ritchie et al., 1999). A datalogger (Licor 1000) was installed to collect weather data on solar radiation, minimum and maximum temperature, and precipitation, which are required as model inputs. Precipitation was measured with an electronic tipping bucket rain gauge every five minutes to record rainfall intensities as well as daily total amounts. 61 The spatial structure for each parameter was assessed using a semivariance analysis. Soil water measurements taken on each grid point were interpolated using punctual kriging technique available in GS+ Version 3.1a (Gamma Design Software, 1999). The first day selected for the model validation was July 3, 1997, day of the year (DOY) 184. The rainfall that occurred on this day was 7.5 mm. This day was chosen to test the performance of the model under a low rainfall amount. The simulations started with the soil profile at DUL for the first 100 cm of the soil profile and at saturation for 100 cm to 200 cm for all scenarios. The performance of the model was evaluated by using the RMSE between the predicted and observed values. 4.3. Results and Discussions 4.3.1. Topographic attributes The topographic attributes are shown in Figure 4.1 through 4.6. “Images and figures of this dissertation are presented in color”. The elevation map (Figure 4.1) shows that the field had an elevation relief of 3.6 m. The north part of the field is the highest point but two other small areas (peaks) also have high relief. These two peaks can be observed in the slope map (Figure 4.2). The slope of the field varies from zero in the flat areas to 3.4 62 % on the backSIOpes of the peaks. Surface curvatures are shown in Figure 4.3 (profile curvature) and Figure 4.4 (plan curvature). They can be thought of as the curvature of a line formed by the intersection of a plane and the topographic surface. This intersection is in the direction of the maximum slope for the profile curvature and transverse to the slope for the plan curvature. Profile curvature is negative for slope increasing downhill (convex flow profile, typically on upper slopes) and positive for slope decreasing downhill (concave flow profile, typically on lower slopes). Plan curvature or contour curvature measures the topographic convergence and divergence and hence, the concentration of water in the landscape. The plan curvature is negative for diverging flow (on ridges) and positive for converging flow (in valleys). The maximum ponding capacity is purely a function of s10pe. Indeed, this can be observed in the maximum ponding capacity map (Figure 4.5) that shows an opposite trend from the slope map. The maximum ponding capacity varied from zero observed on the backslope of the peaks to 3.0 on the flat areas. Figure 4.6 depicts the location and the number of the elements created by executing TERRAE study area. The highest point in the landscape is represent by one element while the bottom of landscape is represented by several elements that describe the lowest elevation point of the field (depressions). 63 4.3.2. Model Simulations Scenario 1-- Day 1. The model results for scenario 1 (uniform soil type across the landscape, no restricting layer, 76 mm rainfall on the first day) are shown in Figures 4.7 through 4.17. The units used in the outputs for all the variables are in cm (height of water). The soil water content for the 0-26 cm (Figure 4.7a) is quite uniform across the field, except for the low elevation areas, which are higher due to accumulation of surface flow onto the elements. These areas also showed higher water content for the profile at the 26-77 cm depth (Figure 4.7b). The cumulative surface leaving the element was high, as expected, due to the quantity of rainfall (Figure 4.8a). However, the surface water balances out as can be seen from the map of the cumulative surface flow out of the each element (Figure 4.8b). The model predicted that water not infiltrated on the element located on top of the landscape runs off to the next element downslope as runon. This explains the balance observed between flow out and flow in. Both maps clearly show the effect of the landscape in the surface water routing. The highest amount of water leaving and entering each element is 65 cm and it is observed in the depression areas due to the contributions from the upperslope elements. The net surface flow (Figure 4.9) is calculated by subtracting the amount of water coming into the element from the one leaving the element. The highest value (-5) is observed on top of the landscape since those elements do not have water running onto them. Figure 4.10 shows the surface ponding. The model was able to correctly determine that the depression areas have higher surface ponding capacities. The subsurface lateral flow is shown in figure 4.11. The highest amount of horizontal flow is observed in the depressions due to high soil water content present at these locations. The vertical drainage is depicted in figure 4.12. The drainage amount predicted is quite small throughout the landscape. This may be due to the rapid occurrence of saturation in each soil layer due to the high rainfall amount. If the elements have a slope greater than zero, the model allows the water to flow horizontal through the KSAT corrections factors that decreases the vertical flow. _D_ay_; The soil water content for the second day of simulation is shown in figures 4.13a and 4.13b. The 0-26 cm soil water content (Figure 4.13a) decreased from the previous day on the high elevation zones and peaks, but did not greatly change in the saddles and depressions due to the higher water flow coming onto the elements and ponding conditions occurring at these zones. Similar phenomenon was also observed for the 26-77 cm soil profile (figure 4.13b). Ponding conditions disappeared on this day, as well as the flow in and out of the elements. The subsurface lateral flow (figure 4.14) also decrease for the second day of simulation but the highest amount (0.3 cm) was still observed in the low elevation areas. The drainage (Figure 4.15) did not vary significantly from the previous day, both for the amounts and locations of occurrence. MThe soil water content at 0-26 cm (Figure 4.16a) indicates that soil surface dried out quite uniformly across the field. The 26-77 cm soil water content (Figure 4.16b) shows a high water content (14.10 cm) at the lowest point on the landscape as result of the higher water content, ponding, surface water flow onto the element and subsurface lateral flow affecting this area. Figure 4.17 shows the sum of the subsurface lateral flow for the seven days of simulation. The highest value (3.6 cm) was observed in the lower 65 Fig COE SCCr lheh elevation areas and depressions as expected. The drainage for day 7 is shown in Figure 4.18. The drainage consistently decreased from the previous days but was still present. Subsurface flow, however, had terminated by day 4. Scenario 2--Day 1. The simulated results for scenario 2 are shown in Figures 4.19 through 30. Figure 4.19a and 4.1% show the soil water content for the 0-26 cm and 26- 77 cm. The model predicted a higher amount of soil water on depression areas and sinks as compared to scenario 1 for the 0-26 cm and 26-77 cm soil profiles. Since the rainfall amount used for scenario 2 was the same for scenario 1, the surface water flow in and out of the elements Figure 4.20a and 4.20b showed similar values and locations to the maps of flow in and out produced for scenario 1 on day 1. The surface ponding (Figure 4.21) is also similar to the previous scenario. The net surface flow is shown in figure 4.22 and showed the same values predicted for scenario 1, since the values of runoff and runon were also similar with the ones observed in scenario 1. The lowest negative value was observed on the high elevation areas since there was no water contribution from the upslope elements. The subsurface lateral flow, shown in figure 4.23, was higher on the depressions due to the higher water content consistently present on the low elevation areas. The value for the subsurface lateral flow was lower compared to the ones from scenario 1 since most of the water remained on the surface as ponding for scenario 2. Figures 4.24 depicted the drainage that occurred on first day. The values were lower compared to scenario 1 due to the restricting soil layer present at 120 cm depths. In scenario 2 the effects of the three different soil types were not visible. This may be due to the high rainfall amount that minimized the soil type influence. 66 Bag. The 0-26 cm soil water content (Figure 4.253) did not change from day 1. Similar results were also found for the 26-77 cm soil depth (Figure 4.25b). Runoff did not occur on day 2. The subsurface lateral flow was highest on the saddles between two peaks, indicating the contribution of water from the elements located on the peaks (Figure 4.26). Drainage values for day 2 (Figure 4.27) were low due to the low conductivity layer. Also on day 2, the effect of the different soil types was not present. ng_7_. The soil water content for the 0-26 cm depth (Figure 4.28a) decreased on the high elevation zones, but it remained the same for the low elevation zones and depressions. The soil water content for the 26-77 cm depth (Figure 4.28b) did not change significantly across the field. The vertical drainage terminated on day 4, while the subsurface later flow was still occurring (Figure 4.29a). The sum of the subsurface flow (Figure 4.2%) is higher (4.6 cm) compared to the one for scenario 1 (3.6), showing the influence from the low conductivity layer. 4.3.3. Validation study-- Scenario 3 Scenario 3 includes a validation with field measured data on soil water content. The simulated soil water content for the 0-26 cm depth is shown in Figure 4.30a. The model predicted higher water content for the saddles and depression areas. The simulated soil water content varied from 1.8 cm to 6.2 cm. Figure 4.3% shows the simulated soil water content for the 26-77 cm soil depth. The highest values of soil water were also observed at the saddle point and depressions and varied across the field from 3.2 to 13.4. The 67 subsurface lateral flow was higher in the saddles and depression area (Figure 31) due to the contribution from runon that increased the amount of potential infiltration on those areas. The values predicted by the model for the vertical drainage (Figure 32) were lower than those predicted for the lateral flow. The locations in the field that had the highest amount of water draining vertically were the areas located on the high and medium elevation. TERRAE-SALUS correctly simulated higher amount of vertical drainage in the areas occupied by the sandy soil. The soil water content slightly changed on DOY 188 (July 7) both for 0-26 cm and 26-77 cm (Figure 4.33a and b). The simulated subsurface lateral flow for day 188 is shown in Figure 4.34. The saddles and depressions showed the highest amount of lateral flow (0.0045 cm). The lowest value of flow was a result of the low rainfall occurring on day 184. Figure 4.35 shows the drainage that occurred on day 188. Its values are also small, almost identical and at the opposite locations where the subsurface lateral flow took place, indicating correct partitioning being employed by the model. Figure 4.36 depicts the error maps from the measured and observed soil water content for day 184 and for the 0-26 cm depth. Figure 4.37 shows the predicted versus the observed measurement compared to a 1:1 line for the soil water content at 26-77 cm depth and for the entire season on the highest elevation point. The model provides accurate results for the entire season but it slightly underestimates the soil water content measured. Figures 4.38 through 4.41 show the measured and simulated results for the soil water content for 0-26 and 26-77 cm soil depth for the entire season using four points along a streamline (from the top-peak, to the bottom of landscape-depression). The model performance was compared using the root mean square error (RMSE). Figures 4.38a shows the comparison between the measured 68 and simulated results for the 0-26 cm and 26-77 cm for the point located on the highest elevation point of the field (264 m). The RMSE observed was 0.51 cm, for the 0-26 cm depth and 0.62 cm for the 26-77 depth. The simulated soil water content for the point located in the upper saddle (263 m) are compared with soil water measurements in Figure 4.38 b. The RMSE observed for this comparison were 0.39 cm for the 0-26 cm depth and 0.52 cm for the 26-77 cm depth. Figure 4.39 shows the comparison between simulated and simulated soil water content for the lower saddle point (262 m). A RMSE of 0.46 cm and 0.49 was observed for this comparison for the 0-26 cm and 26-77 cm depth. An evaluation of the model performance was also done for the depression area of the streamline selected (260 m) The RMSE observed for this evaluation were 0.47 cm for the 0—26 cm and 0.59 cm for the 26-77 cm depth. T116 soil water content simulations and field observations were also compared as function of elevation. The RMSE for this evaluations are shown in Figure 40 a and b. The high elevation point consistently showed a lower water content compared to the upper and lower saddles and for the depressions. The days used in this final evaluation of the model, the soi ] water content did not change significantly for the saddle points and for the depressions. This is due to the contributions of water running downhill fiem the peaks. 69 4.4. Conclusions This chapter discussed the application of TERRAE-SALUS, a digital terrain model with a functional spatial soil water balance model, at a field scale to simulate the spatial soil water balance and how the terrain affects the water routing across the landscape. The first part of the chapter described the principles for the two models. The second parts discussed the capability of the model tested under three different scenarios. The scenarios were selected to evaluate the model sensitivity with a soil having no vertical drainage restrictions and for one that had basically almost no vertical flow. The model was able to partition the subsurface lateral flow and the vertical drainage differently for the two scenarios, but the high amount of rainfall seemed to have an higher effect through the amount of runoff and ponding that occurred on the first day, making the rest of the days qui (:6 similar between the scenarios. The model provided excellent results when compared to the field measured soil water content. The RMSE between measured and simulated results varied from 0.22 cm to 0.68 cm. The performance of TERRAE-SALUS is very promising and its benefits can be quite substantial for the appropriate management of water resources as well as for identifying the areas across the landscape that are more susceptible for erosion. It is necessary to further validate the model at different sites with different soils and weather characteristics. 70 Figure 4.1. Elevation at 1 m grid resolution for the study area, Durand, MI. Figure 4.2. Slope for the study area, Durand, MI. 71 “tin » Figure 4.3. Profile curvature of the study area, Durand, MI. Figure 4.4. Plan curvature of the study area, Durand, MI. 72 Figure 4.5. Maximum ponding capacity for the study area in Durand, MI. Figure 4.6. Location in the landscape of the element network created by TERRAE (Gallant, 1999). 73 Figure 4.7.a. Soil water content (0-26 cm) on day-l for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). (cm) l4_l0 Figure 4.7.b. Soil water content (26-77 cm) on day—l for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). 74 Figure 4.8.a. Cumulative surface water flow out of the elements on day—1 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). Figure 4.8b. Cumulative surface water flow onto the elements on day-l for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). 75 Figure 4.9. Net surface flow on day-1 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). Figure 4.10. Surface ponding for day-1 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). 76 Figure 4.11. Subsurface lateral flow on day-1 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). Figure 4.12. Drainage on day-l for scenario [(1 uniform soil type, high rainfall, no restricting soil layer). 77 (cm) (cm) l4.l0 13.70 49‘ ‘ Figure 4.13.a. Soil water content (0-26) on day-2 for scenario 1(1 uniform soil type, high rainfall, no restricting soil layer). Figure 4.13.b. Soil water content (26-77) on day-2 for scenario 1(1 uniform soil type, high rainfall, no restricting soil layer). 78 Figure 4.14. Subsurface lateral flow on day-2 for scenario I (1 uniform soil type, high rainfall, no restricting soil layer). Figure 4.15. Drainage on day-2 for scenario 1(1 uniform soil type, high rainfall, no restricting soil layer). 79 Figure 4.16 a. Soil water content (0-26) on day-7 for scenario 1(1 uniform soil type, high rainfall, no restricting soil layer). Figure 4.16 b. Soil water content (26-77) on day-7 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). 80 Figure 4.17. Sum of the subsurface lateral flow on day-7 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). Figure 4.18. Drainage on day-7 for scenario 1 (1 uniform soil type, high rainfall, no restricting soil layer). 81 (cm) I 1.50 10.50 Figure 4.19.a. Soil water content (0-26 cm) on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). Figure 4.19.b. Soil water content (26-77 cm) on day-l for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). 82 70.“) 65.“) 60.11) 55.“) 50.00 45.“) 40.00 35.00 30.00 25.11) 20.“) l 5 .11) 10.11) 5.00 Figure 4.20.21. Cumulative surface water flow out of the elements on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). Figure 4.20. b. Cumulative surface water flow onto the elements on day-1 for scenario 2 (3 different soil types, high rainfall,with restricting soil layer at 120 cm). 83 Figure 4.21. Surface ponding on day-l for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). Figure 4.22. Net surface flow (cm) on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). 84 \ ., ’5' *o O O jg§fia e o:¢’. ’. 9W ’o‘o too». MW ' 000% fi,‘ ' a 9 s .o..:r»..:.. I, . .Q . Q - tea-aw. . e 0‘. o‘OQ"... $.Q . €43§3£39§eb ‘3%§§§€59¥° Tfi:§§5§3§€§*&efifi§fis§§§9€;. b - Q Q‘ ‘ a.o.:.;:::. ‘_ ‘ 3338‘ - - o -; - \\ .°'°.::’.£“;;;“ 'bfiqgfihflhsfitks":‘ ‘s?: a 0’09.“\“: '0..:.‘:.“\ o ow‘53‘ 0*. 9 Figure 4.23. Subsurface lateral flow on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). Figure 4.24. Drainage on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). 85 Figure 4.25.a. Soil water content (0-26 cm) on day-2 for scenario 2 (3 different soil types,high rainfall, with restricting soil layer at 120 cm). Figure 4.25 .b. Soil water content (26-77 cm) on day-2 for scenario 2 (3 different soil tYDeshigh rainfall, with restricting soil layer at 120 cm). 86 Figure 4.26. Subsurface lateral flow on day-2 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). Figure 4.27. Drainage on day-1 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). 87 (cm) (cm) Figure 4.28. a. Soil water content (0-26 cm) on day-7 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). Figure 4.28. b. Soil water content (26-77 cm) on day-7 for scenario 2 (3 different Soil types, high rainfall, with restricting soil layer at 120 cm). 88 (cm) Figure 4.29. a. Subsurface lateral flow on day-7 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). Figure 4.29. b. Sum of subsurface lateral flow on day—7 for scenario 2 (3 different soil types, high rainfall, with restricting soil layer at 120 cm). 89 49‘ Figure 4.30a. Simulated soil water content (0-26 cm) on day 184 (July 3) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). Figure 4.30.b. Simulated soil water content (26—77 cm) on day 184 (July 3) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). 90 0.0 l 20 0.01 10 0.0100 0.11190 0.“)80 011170 0.(X)60 01050 0.0040 0.0030 0.0020 0.(X) 10 0.0000 Figure 4.31. Subsurface lateral flow on day 184 (July 3) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). Figure 4.32. Drainage on day 184 (July 3) for scenario 3(3 different soil types, low rainfall, with restricting soil layer at 120 cm). 91 Figure 4.33 a. Simulated soil water content (0-26 cm) on day 188 (July 7) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). Figure 4.33 b. Simulated soil water content (26—77 cm) on day 188 (July 7) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). 92 Q‘ Q s" s s .‘20‘:“:o 9“.‘.~‘.““ . . 9...... 9" \‘ s‘ .‘ 0 <3 0 ,Q ‘Q." ‘ ‘ ‘."s‘..‘ 0‘ Q _¢ o" ,» ‘e c‘ "..: \‘¢‘ \‘ o 9%“stsfi — \ ‘ . a.“ e \ é‘q’ hex‘o; Q 9 r » ‘\ o .e‘ ‘ Q Q v s‘ s‘o‘o $ \s‘.$..p‘. . c e;.~t.¢‘ \s‘:¢ ~‘ _.s;.o“s - .\.. .§ ~s ¢3¢¢3¢3 3‘ Q‘ 0 ‘ vk‘ ’19 “0 o..¢.s..0,» ‘ . O. .‘. 5‘ o‘ ‘3 . v.46»?- / .‘0‘35‘3’5 e»: .ekgtasgrofi- \ 9‘9. ““‘“ ‘\\ 33353;,“ - V, ‘7’ dogs} -3 .. >‘ 1‘ ‘e’.%‘ $3: . a A“ \¢ 1 ¢~ 4e \ -r'z’ - ' ‘.n,”" '. _ Q’HV " a $33.. $30: \ .s 7.. '~°.<~:~;‘1~‘1‘:e3 g. ,‘ Figure 4.34. Subsurface lateral flow on day 188 (July 7) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). Figure 4.35. Drainage on day 188 (July 7) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). 93 §V“:\“\ -.- Q ‘ \s‘.o‘-¢®‘. . c .; \Q‘.‘ . .9 \s‘ - - ...:20;‘s‘2:::.‘::.::::‘,¢::‘:.\ . .o .3 ..§ ,s ‘. ‘.s .s sios‘ -. .‘.«‘s‘c‘o Q 5 ‘ ¢ ‘ o 9 - s . ‘¥.°.$.9¢ ““.s‘v a “«’3‘3°$‘:3‘$°‘<‘~ “ » g Q * ¢ 0‘s: 0‘ ’o' i "0.4 0“, 9“ a ‘9 g :a (cm) Figure 4.36. Measured soil water content (0—26 cm) on day 184 (July 3). Figure 4.37. Measured soil water content (26-77 cm) on day 184 (July 3). 94 15.“) [4.50 ”.00 l3.50 13.00 12.50 Mel-ma Sol Wu" Content (cm) 0 I I 0 2 4 6 8 l0 l2 Sllllmd Sol Woul- Come"! (cm) Figure 4.38.21. Error map of soil water content (0-26 cm) on day 184 (July 3) for scenario 3 (3 different soil types, low rainfall, with restricting soil layer at 120 cm). Figure 4.38.b. Measured and simulated soil water content (26—77 cm) for the entire season on the high elevation point (peak). 95 264 murtl’eak) l2.” ~ , , » —D— Mcasued 0-26 cm 10( RMSE=0.51 --|-- smmozo cm ' ) —A— Mcasued 26-77 cm mA- Simhtcd 26-77 cm .. 8.0 ~ E .2 I g 6.0 ~ In a 3 4.0 4 2.0 4 0.0 . . . f 7 s . 1 l 80 I90 200 2 H) 220 230 240 250 260 Day of (In Year 263 natal-(Upperaadtle) 12.0 _._,___. —O— Measued 0—26 cm ‘ - - In Situated 0-26 cm 10.0 4 —a— Mcasued 26-77 cm “L MA 5mm 26-77cm 8.0 4 'O-. .... ...... ----- u .......... ‘...-~---‘ l RMSE=O.39 ‘ Water Council! (an) o o 4.0 ~ 2.0 ~ RMSE=O.52 0.0 . T . . . . . fl I80 190 200 210 220 230 240 250 260 Day of!!! Year Figure 4.39 a. Measured and simulated water content for the soil profile (0-26 cm) and (26-77 cm) for the high elevation zone (peak) for the entire season. Figure 4.39 b. Measured and simulated water content for the soil profile (0-26 cm) and (26—77 cm) for the medium elevation zone (upper saddle) for the entire season. 96 262 murthower sadde) 16.0 ~ ~. . ,1 14.0 4 RMSE=O.46 +Mwm0 Mn 1 A.. ---I-- 5m tod0-26cm | . —o— Measu'cd 26-77 cm l 12'" i ---A- snuh' ted 2077 cm II E 100 ~ 2 e 3 8.0 J L RMSE=0.49 s 6.0 « 4.0 1 2.0 ~ 0.0 . Y . 7 . . . . 180 190 200 210 220 230 240 250 260 Day ofthe Year 260 utter (Demsion) 160 -——————— — — - ——~———— - — — - 4 14.0 i RMSE=0.47 —0— Measured 020 cm 5" ---I-~ smmnzscm 120 —A- Measured 26-77 cm -----A Situated 2077 cm E 10.0 « 3 e g .0. I- i- l 3 6.0 4.0 4 2.0 4 0.0 . . . , fl - - . I80 190 200 210 220 230 240 250 260 DayonheYear Figure 4.40 a. Measured and simulated water content for the soil profile (0-26 cm) and (26-77 cm) for the medium elevation zone (lower Saddle) for the entire season. Figure 4.40 b. Measured and simulated water content for the soil profile (0-26 cm) and (26-77 cm) for the medium elevation zone (depression) for the entire season. 97 Water Comm Chm vs Elevation 7.5 4 7.0 4 6.5 4 ------ .. RMSE=O.34 6.04 5.5 4 —D—— Meas. 026cml'XW-lu - - -I - ~ Snail Ohm-DOWN —a— Men 01’ cmDOY 2'15 ----h Saml OhmDOYZE —O— Mus 026 cmDOY 240 - - 4. -- Saml (mev-DOYZ-IO Water Content (cm) P J RMSE=0.57 2.5 “ . 2.0 4 1,5 - . . . . . 259 260 26l 262 263 264 265 Elevation Water Code! Cline vs Elevation mm 1 ____._-_.__.__..____..-_- _ -.-__--..,,.. —~ -—- ,_._. .__ u- , , ,, 14.00 4 12.“) 4 g —-o— Meas. 26774;!!ka nu E m... smutzs-ncwoovm " -—o—Mm. 2671ch 219 5 10.00 ‘ RMSE=O-5l ----ASmn1.2677cu»-00Y2|9 § k .. . A —o— Meas, 2677cm-DOY m i --v.--Simil.2677cm-DOY248 8.00 4 6.00 ~ ‘ 4,“) v r f 1 g 1 259 260 261 262 263 264 265 Elevation (In) Figure 4.41 a. Water content change as function of elevation for 0-26 cm soil profile for day 184 (July 3), day 225 (August 13) , and day 240 (August 28). Figure 4.41 b. Water content change as function of elevation for 26-77 cm soil profile for day 184, day 219 (August 7), and day 248 (September 5). 98 CHAPTER V UNDERSTANDING SOYBEAN YIELD VARIABILITY USING CROP MODELS AND REMOTE SENSING 5.1. Introduction Agricultural production systems are inherently variable due to spatial variation in soil properties, topography, and climate. To achieve the ultimate goal of sustainable cropping systems, variability must be considered both in space and time because the factors influencing crop yield have different spatial and temporal behavior (Pierce and Nowak, 1999). Advances in technologies such as Global Positioning Systems (GPS), Geographic Information Systems (GIS) and remote sensing have created the possibility to assess the spatial variability present in the field and manage it with appropriate site-specific practices (Verhagen et al., 1995). Site-specific management (SSM) strategies may be able to optimize production, but their potential benefits are highly dependent on the accuracy of the assessment of such variability (Pierce and Nowak, 1999). Traditional analytical techniques, such as regression of static measurements against yield, have failed to explain the reasons for yield variability because the dynamic, thus temporal, interactions of stresses on the crop growth and development cannot be accounted for (Jones et al., 1989; Cambardella et al., 1996; Sudduth et al., 1996). Process 99 oriented crop simulation models, such as the CERES and CROPGRO models (Ritchie eta], 1998; Boote et al., 1998), have the capability to integrate the effects of temporal and multiple stress interactions on crop growth processes under different environmental and management conditions. Even though crop models have shown high potential for optimizing production and minimizing environmental impact, their application for SSM has been limited thus far. Crop models can be used for understanding yield variability, leading to a more sustainable environment (Sadler and Russell, 1997, Cora et al., 1998). Batchelor et al. (1998) and Paz et al., (1999), used CERES-Maize and CROPGRO- Soybean simulation model to determine the effect of soil water variation throughout the season on yield spatial variability optimizing for soil water limit parameters. The differences between measured and simulated yield for 224 grid points over a 3 year period ranged from i 10% for 70% of the grids and i 20% for 96% of the grids in maize and for soybean from i 10% for 84% of the grids and i- 20% for 92% of the grids. Recent advances on the resolution and availability of remote sensing imagery, coupled with a decrease in its associated costs, have allowed the collection of timely information on soil and crop variability by examining spatial and temporal patterns of vegetation indices (Moran et al., 1997). Such information can be useful to derive inputs for crop models in a GIS environment (Moran et al., 1997; Barnes et al., 1998, Johannsen et al., 1999). Vegetation indices involvemathematical relationship between near infrared and red reflectances and they have been extensively used with the goal of estimating vegetation amount (Wiegand et al., 1990; Jackson and Huete, 1991; Price, 1992). Among vegetation indices, the normalized difference vegetative index (NDVI) is the one most 100 commonly used to quantify canopy vigor and density (Price, 1992; Carlson and Ripley, 1997). NDVI is defined as: _ NIR — RED NIR + RED ND VI where NIR and RED represent the surface reflectance averaged over ranges of wavelengths in the near infrared (A ~ 0.8 pm) and the visible (A ~ 0.6 pm) regions of the spectrum, respectively. NDVI increases almost linearly with increasing leaf area index (LAI, leaf area per unit land area) until LAI exceeds values of 3-4, above which NDVI rapidly approaches an asymptotic limit (Liu and Huete 1995; Jasinski 1996, Carlson and Ripley, 1997). NDVI analysis performed on images taken at critical times during a growing season can help characterize spatial variability in crop performance. Clearly, the goal of crop simulation in precision agriculture is to explain the spatial variability of crop performance mapped with grain yield monitoring systems and to help guide in management decisions related to the site-specific management of crop inputs. It is also clear that crop simulations cannot be performed everywhere given that the cost and the availability of detailed inputs would be prohibitive. A more balanced approached to the application of crop simulation models to precision agriculture would be to delineate zones within the field of similar crop performance. One approach may be to obtain vegetation indexes derived from remotely sensed imagery during critical times during the growing season to select spatial patterns to sample and to use the results of the target sampling as inputs for the models. Model validation can be performed at selected sites within these delineated management zones. 101 The objective of this study was to examine a new procedure for spatial validation of crop models for use in precision agriculture that uses the CROPGRO-Soybean model to simulate soybean performance using a progressive increase of spatial inputs. The procedure also uses the crop model to validate management zones across the field delineated using a NDVI classification procedure. 5.2. Materials and Methods 5.2.1. Site Description and Field Measurements The study area consisted of a seven ha portion of a field located 10 km south of Durand, MI. The field has been cropped to a com-soybean rotation for more than 10 years. Soils are variable containing five soil map units and considerable spatial variability in soil fertility (Pierce et al., 1995; Pierce and Warncke, 2000) with the major soil type in the experimental area mapped as Capac loam (Udollic Ochraqualf fine, loamy, mixed, mesic). Soybean was grown in 1997. The field was planted on May 5 by direct drilling soybean (Variety Asgrow 1901, a Roundup Ready variety) in 37-cm rows at a seeding rate of 494,000 seeds ha". A regular grid consisting of 52 grid locations spaced 30.5 m apart was imposed on the 7 ha experimental area after planting. Latitude, longitude and elevation of each grid point were determined with a high-resolution differential GPS. Neutron probe access tubes 102 were installed at each of the 52 grid locations. A neutron moisture gauge was used to measure the spatial variability of soil water content at lS-cm increments to the depth of the C horizon or a maximum of 150 cm depth, which ever occurred first. Measurements were taken on a weekly basis throughout the season. During the installation of the neutron probe access tubes, soil samples were taken at the intersection of each of the 52 grid points in 25-cm increments and stored for analysis Soil samples were air-dried and passed through a 2 mm sieve. Particle size was determined for each segment of each soil profile using the hydrometer method. Soil organic matter was determined on the surface 25 cm of each soil profile by dry combustion using a CHN analyzer (Carlo Erba Instruments, Italy). The upper and lower limits of soil water availability were determined using soil water measurements taken in the field during the season, and also from using Ritchie’s simple model to estimate soil water limits using soil texture data (Ritchie et al., 1999). Soil depth for each grid point was determined using the deepest depth observed during the installation of neutron probe access tubes. Potential extractable soil water (PESW) was determined by subtracting the lower limit of plant water availability from the upper limit for each soil layer and integrated for the entire profile. A 5 m2 area was delineated at each grid location for selected plant measurements. Plant population and the distance between plants were measured at emergence (May 20) and at the 3rd leaf stage of development (June 15). A non-destructive optical device with a fish- eye sensor (LAI- 2000; LI-COR) was used to quantify the LA] at the 52 grid points on July 15 and August 10. Soybean yield was obtained by harvesting four rows along a 20 m 103 length centered on each grid point using a plot combine. Grain moisture was obtained after harvest on a subsample from each harvested area. A datalogger (Licor 1000) was installed to collect weather data on solar radiation, minimum and maximum temperature and precipitation, which are required as model input. Precipitation was measured with an electronic tipping-bucket rain gauge every five minutes to record rainfall intensities as well as daily total amounts. Standard statistical analyses were done for the variables measured in the field. The spatial structure for each parameter was assessed thorough a semivariance analysis. Measurements taken on each grid point were interpolated using the punctual kriging technique available in GS+ Version 3.1a (Gamma Design Software, 1999). Correlation matrices were developed to determine the relationships among variables for each single class and for the 52 grid points. 5.2.2. Remote Sensing Data The airborne false color composite images in the green, red and NIR portion of the spectrum were taken on June 1, June 28, July 18, July 29, August 13, August 29, September 15 at 1 meter pixel resolution. The images provided spatial information about the condition of the crop throughout the season. Each image was used to generate NDVI maps of the field and to identify spatial patterns across the field . The false color composite image taken on July 18 (figure 5.1a) was selected for quantifying areas with 104 similar reflectance by grouping areas into classes of similar NDVI values (Figure 5.2) using supervised classification technique available in Idn'si v32 software (Clark Labs, 1999). Pixels of similar reflectance were queried across the field after trying various ranges of values able to reproduce the spatial patterns visible in the original false—color composite image (Figure 5.3). 5.2.3. Crop Growth Model CROPGRO-Soybean v.3.5 is a process oriented model that simulates plant responses to environmental conditions (soil and weather), genetics and management strategies. A detailed description of the inputs required for the model is described in Ritchie and Dent (1990). This model is part of the Decision Support System for Agrotechnolgy Transfer (DSSAT 3.5, Hoogenboom et al., 1998) that provides several tools for model application. The soil water limits used to run the simulation experiments varied spatially and according to the observed data of soil texture and soil water content at the 52 grid points. The model was evaluated using the Root Mean Square Error (RMSE): RMSE =[}1;2(Y.- _ y )1“: i=1 where y, are the measurements, y ,. the predictions, and n is the number of comparisons. 105 5.2.4. Simulation Experiments The soybean model was used to simulate yields in the field using progressive increase of spatial inputs. Yields simulations were made using five different scenarios. These scenarios varied from one that assumed uniform soil and management conditions across the field to one that used field-measured, spatially variable inputs for the soil water balance parameters (LL, DUL, SAT and soil depth) and plant populations, to one that simulated three areas identified by the NDVI analysis. The five cases are described in table 5.1. 5.3. Results and Discussion 5. 3.1 . Field Measurements The yield was spatially variable across the field (Figure 5.23), ranging from 1900 to 3600 kg ha’1 with a mean value of 2953 kg ha‘1 and a standard deviation of 433 kg ha]. The spatial distribution of yield was consistent with other field measurements (LAI and plant population) and by the remote sensing image that showed high reflectance in the high yielding areas (Figure 5.1 and 5.2). Field measurements of LAI on August 8, (Figure 5.2b), reflected the different soil water regimes present across the field. The highest value of LAI was 4.6 and it was observed in 106 areas of high plant population, deeper soil and high PESW. The areas of the field with rocky soil and highest elevation had the lowest LAI value of 1.7. The mean and the standard deviation for LA] were 3.6 and 0.6 respectively. The areas showing high LAI corresponded with the remote sensing image areas with high reflectance as depicted in figures 5.1 and 5.2. Record cool weather in May delayed soybean emergence and resulted in variable population densities across the field (Figure 5.2c). Plant populations varied from 22 to 60 plants in2 with a mean value for the 52 grid-points of 47 plants In2 and a standard deviation of 8 plants m'z. Plant stand is highly influenced by the environmental condition (soil and weather) at planting time. Soil textural analysis from the 52 grid-points showed high spatial variability for sand and clay particles. The clay content varied across the field from 8% in soil of the high relative elevation areas to 25% in soil in low relative elevation areas. Sand percentage varied from a minimum of 40% to a maximum of 82% and logically had an opposite spatial distribution from the clay content. Based on textural analysis results and their spatial distribution, three main soil types were detected across the field. A deep-dark sandy-loam soil located in lower elevation areas of the field, a sandy loam characterizing the flat areas and a sandy-rocky soil present in the higher elevation areas. Soil depth measurements also showed the presence of high spatial structure across the field (Figure 5.3a). Spatial variations in soil depth had similar trends as clay (Figure 107 5.3b). Peaks had lower soil depth due to erosion phenomena. Soil depths ranged from 95 cm to 150 cm with a mean value for the 52 grid-points of 130 cm and a standard deviation of l4-cm. Potential extractable soil water (PESW) values are shown in Figure 5.3c. PESW is a function of soil depth and soil texture, thus the spatial distribution of these variables were similar. A maximum PESW of 140 mm was observed in the low elevation areas, while the lowest value of 70 mm was found at high elevations. The mean PESW value for the field was 111 mm with a standard deviation of 19 mm. The field is characterized by a rolling terrain that caused high spatial variability of soil properties. Landscape position and relative elevation highly influenced soil physical properties thorough erosion processes that occurred over the years (Mueller, 1998). The spatial dependence was determined for each soil and crop variable measured in the field. Geostatistical analysis revealed spatial structure for all the variables giving ranges of distance that varied from 60 m for the plant population to 150 m for the yield (Table 5.2a). The spherical model fitted the semivariance well. Table 5.2b reveals that the yield measured at the 52 grid points is highly correlated to LAI, PESW and NDVI as shown by the correlation coefficients of 0.86, 0.87 and 0.97. 108 5.3.2. Simulation Experiments 5.3.2.1. Scenarios 1-4 Error in yield prediction decreased as the level of input detail increased for the simulation scenarios tested (Table 5.1). The field average of 2995 kg ha'1 was underestimated under scenario 1 which predicted a soybean yield of 2530 kg ha]. The RMSE for scenario 1 was 465 kg ha". Under scenario 2, adding site-specific plant population data as model input improved model performance by decreasing the RMSE to 296 kg ha", a reduction of 36 % over scenario 1. Using site-specific soils data at constant plant populations in scenario 3 improved yield prediction 18 % more than scenario 2 as evidenced by an RMSE of 245 kg ha’1 and 47% over scenario 1. Using both site-specific soils and plant population input further reduced RMSE improving the prediction of soybean yield over scenario 1 by 58% (Figure 5.4a). 5.3.2_..2. NDVI Classes The 18 July composite image and the corresponding NDVI image clearly show spatial variability in soybean performance (Figure 5.1b). The correlation between NDVI and crop yield for the 52 grid points was very strong (Figure 5.5). Classification of the NDVI image indicated three classes of importance in this field. Note that the areas of different classes are not contiguous. Table 5.3 summarizes the areal distribution and properties of these classes as well as the data for soils and plant populations used as input in the crop 109 simulations under scenario 5. Yield predictions for the three NDVI classes were very good as evidenced by an RMSE values for the three classes (Table 5.4, Figure 5.1c). These RMSE were the lowest of all crop simulation scenarios evaluated showing that the NDVI reclassification procedure adopted in the study was appropriate and proven to be a useful way of creating zones in the field. Multi-year simulation with this approach would allow for the characterization of the field. The progressive use of site specific model inputs combined with the NDVI- reclassification has a major advantage over an issue that thus far has limited the power and application of simulation model in precision agriculture: scale! The site-specific input approach is scale-independent because the scale is controlled by the observed variation in the field and that is the scale at which the model will be applied. 5.4. Conclusion Analysis of remote sensing imagery processed into the NDVI identified the spatial patterns of crop growth variability. The variability in soybean populations within the field provided validation of the plant population and soil type effects on soybean yields predicted by the CROPGRO- Soybean model. The model gave accurate predictions of the yield across the field when the correct inputs were used showing great potential for use in yield map prediction and interpretation in the context of site-specific management. This 110 study showed that the soil types present in the field could be managed differently using different cultivars or plant density to achieve higher yield while minimizing the costs. It is clear that zone-specific management to optimize production can be developed through a combination of remote sensing and simulation models. This is a more affordable alternative to the use of traditional soil sampling and micro-scale sensing. It also answers questions related to the scale issue by applying the model at the scale of variability observed through remote sensing and NDVI image reclassification. lll Scenario Model Input Model RMSE No. Runs Kg ha" 1 Average of grids 1 465 (Average soil type and target plant population) 2 Average soil type and grid plant population 52 296 3 Grid soil type and target plant population 52 245 4 Grid soil type and grid plant population 52 198 5 Average soil and plant population for 3 101 3 NDVI Classes Table 5.1. Model inputs, number of model runs and RMSE for each simulation scenarios. 112 458-2% Nm 2: 8 22.. 2: E 9:38.: 83%? muons SEES 5:20:00 dad 038. _ wwd w _ d- odd and Ed. 3d nod 5d 5G2 _ de- @— d mmd cvd- 5d Ed mad 3%.: _ .wd- mad- 3 .d m _ .d- 2d- w. d- AEumhdv team as _ C d m _ .d- mdd 5d nod EB mhdv .35 a.» _ M: .o- .3 63 mac 5%; =8 _ ovd- 3d- Ned- gage—m— _ mod owd 23> _ and cede—=9..— =35 _ :4— ;GZ 3mm.— AEomhdv deem 6x6 :5 mhdd >20 6R. 53: aem 5532m— 209 cote—:3.— .car— :3 EBBSS .23 05 E 35.82: moist: 05 a8 B20888 mEEwotmEEom .m~.m £an 113 8. 8... sec 3823 3mm: 2 _ 83 23 .85.} Eu 8_ a; ._ 83 .822} 58a =8 8 48.. o. S .853 8:22....— .5: c: on... 2.3 .853 :3 o2 32 33 36.22% 2...» owes— Em «alum—.2 .352 836.23» Figure 5.1 a. False-color composite image taken on July 18, 1997 (original image delivered by Emerge). Figure 5.1 b. NDVI Image. Figure 5.1 c. Reclassified NDVI Image. The areas in white are NDVI-Class 1, the areas in grey are NDVI-Class 2 and the areas in black are NDVI-Class 3. 114 257350 257400 257450 257500 257550 257600 257650 474785 257350 25 7400 257450 257500 257550 257600 257650 474785 257350 257400 257450 257500 257550 257600 257650 Figure 5 .2. Krigged map of measured soybean yield (Kg ha'l)(a); LAI measured on August 8 (b); plant population measured on June 5 (plant m'2) (c). 115 257350 257400 257450 257500 257550 257600 257650 474785 257350 257400 25 7450 257500 2575 50 257600 257650 257350 257400 257450 257500 257550 257600 257650 Figure 5.3. Krigged map of soil depth (cm) (a); clay content (%) (b); and potential extractable soil water (mm) (c). 116 £an 2% mm 2: cod owns; 28 mam—o 592 some 5 ©2335 333:5 .m.m 03E. 858 Ed Ed o. .2 3 dm. N.d Omen mmv mmdm x N.v dd On I Macao 25 mm mdd omd ~_ on. 2 Nam. 5d vNcm Q: mm; N. dm wmd mod. mdA 053 m dd 3.0 0 ea N. mm. vd Numom N._m 0.6.5 N. av 3d Worm md-m~d wmmmv N mdd 2d 0 MN. mm d: wd ,6..ch omm mmd. m. cm 3d 2N mmdd mew—m _ NEE :55 2.3 3.3 :5 :5 N..£_mv_:_.2wv: $5.2: $5.3 500.6 :32 .>oQ..m :82 sonam =82 sun—gm :32 500.5 :82 460.5 :32 3.0.5 :82 owes— . 68.0 35.8 E ~>GZ =00 awe—U 53a— zeta—ago.— N>GZ 3mm.— =om gouge—m 20; :3..— _<‘_ 117 .20: on. £88 owfiozw BEEP.» we.“ was 820 N >0 Z some a8 22% dos—:86 new @2838 cod 038 bee—cam .vd 2an .o. Sam 9% 8.6 mm 836 :52 ..a sea omens... 82.363 6. . 82 2.9. 62 an m 2.5.392 55.3 22.... 6.3.62 N. 83 £8 a... m. ~ 25.5% 55.; 22.... 69262 N. o. N. 2... . N3 5 . §.u._>nz 55.3 33:. omega. / .22. no: :2. we: .22. we: .2: .56.. 3.5 ”.52: 26; 33.2.5 26; 3.3.8:. 8:. .38. cc L.2...=.z 6.55.282... 5.3.2.5 / 18 1 4747850.“ 4747800." . 1 ’ ' a I l l 257350.00 257400.00 257450.00 257500.00 257550.00 257600.00 257650.00 ': g. E3 if 5‘? b a? 1 l l l I 257350 257400 257450 257500 257550 257600 257650 474 C 474 25760 25765 25735 25740 25745 25750 25755 Figure 5.4. Kriged map of simulated yield using measured plant population and measured soil type for the 52 grid point (a); measured yield for the three NDVI classes (b); Simulated yield for the three NDVI classes using average measured inputs. 119 3yn< 3an< Yield (kg ha") 2am~ UG)‘ 1am Figure 5.5. Correlation between NDVI (image taken on July 18) and yield. 2flD* R2 = 0.9353 0| 02 03 04 NDVI 120 05 06 07 08 CHAPTER VI CONCLUSIONS AND RECOMMENDATIONS This chapter summarizes the importance of using simulation modeling and digital terrain analysis to evaluate the effect of topography and soil physical properties on spatial soil water balance. Chapter one presented the rationale and background of this dissertation. A detailed literature review (chapter two) was done for the terrain analysis, reviewing existing digital terrain models and soil water balance models. From this chapter, it was shown that a new hybrid model that combined a digital terrain model with a spatially sensitive soil water balance was needed in order to better simulate water movement over the terrain. The terrain model and the spatial component of the soil water balance models were discussed in chapter three. The model was able to partition the landscape into an interconnected series of element network from a grid DEM. The soil spatial water balance of SALUS was able to partition the downward water flow in the vertical and horizontal dimensions. The model was rapid in its simulations and output outcome. 121 Chapter four was the primary focus of this dissertation. In this chapter, data were presented on the newly developed TERRAE—SALUS model and its ability to simulate the spatial soil water balance as affected by landscape topography. The model was applied at a field scale in Durand, Michigan, were an extensive data set on soil water measurements was available. The RMSE between model and observed results varied from 0.22 cm to 0.68 cm of water. In addition, two scenarios were presented illustrating the model capabilities to partition the subsurface lateral flow and the vertical drainage as well as the surface water runoff-runon as affected by landscape positions and by rainfall amount. TERRAE-SALUS was able to simulate satisfactorily the soil water content. The biggest advantage of this model appears to be its simplicity and at the same time accuracy. Due to the functional approach of the soil water balance, the data inputs requirements are also minimum and easy to obtain. Chapter five described the integration of the current technology available in agriculture to predict the spatial variability of soybean yield. The crop simulation model CROPGRO was applied in combination with remote sensing data to evaluate the capability of the model to identify factors responsible for the yield variation in a spatially variable landscape. Results from this study showed that a combination of crop simulation model and remote sensing can identify management zones and causes for yield variability, which are prerequisites for zone-specific management prescription. The processes of modeling the soil water balance over space and time is of crucial importance for the appropriate management of the water resources, especially in areas 122 where they are in continuous jeopardy. Although, this is not an easy task due to the difficulties related to the complexity of the soil-water-atmosphere systems, uncertainty of weather, and lacking of good quality data to be used for the model evaluation. The scale issue is also a restriction on the power of the existing simulation models if applied at the inappropriate scale to simulate a process characterized by different scale characteristics. Moreover, it must not be forgotten that any model is a simplification of reality. If the biological processes modeled were not dynamic, the modeling of such system would have been much easier. In the case of this study, TERRAE-SALUS has been shown to produce satisfactory results in its field application. The accuracy of the model results is highly dependent on the quality of inputs used, especially pertaining to soil characteristics, elevation data and weather information. Further investigations are recommended to evaluate approaches to the problems of up scaling model simulations. The up-scaling is not solely achieved by running the model patch scale models for larger areas consisting of many patches, but that different processes and connectivities emerge as dominant as we move from the plot scale to catchment scale. The promising results showed by TERRAE-SALUS demonstrated that its application can be beneficial in water resources management. My vision for the future is that digital terrain modeling will be become increasingly important in simulating the most sustainable soil and water resources management practice. The application of DTM can help in identifying areas across the landscape more susceptible to erosion and with the highest surface runoff. Moreover, the model does not need to be applied on the entire landscape, requiring large amount of inputs but it can be executed 123 on small areas and then extend the output to areas that are alike across the landscape, avoiding in this way repetitive data collection, not necessarily needed. 124 APPENDICES 125 APPENDIX A USING TERRAE TERRAE requires a grid DEM. TERRAE is run with a single argument specifying the parameter file: Unix> TERRAE dem.params The parameter file contains all the information required to run the program, including the names of the DEM, channel head and boundary files as well as numeric parameters such as the maximum sink depth and logical parameters specifying which components of TERRAE are active. The most commonly used parameters are: 0 Grid file = dem.flt The name of the DEM file exported from ARC/INFO using GRIDFLOAT. TERRAE expects to find the header file (eg dem.hdr) in the same directory as the binary DEM file. A full path name may be used. This parameter must be set (there is no default). 0 Channel head file = channels.txt The channel heads, as ASCII x y values one per line. This is an optional parameter; if it is not set or the file cannot be opened, channel heads will not be used. 0 Diversion line file = boundary.ung The boundary polygon in ARC/INFO ungenerate format. This is an optional 126 parameter; if it is not set or the file cannot be opened, no boundary polygon will be used. Remove sinks = yes or no If yes, sinks will be cleared if they are within the depth and/or distance thresholds. Default is yes. Sink depth threshold = number The maximum depth of a depression that will be considered a spurious sink. There is no default value, so if “Remove sinks” is set, a value must be provided. Sink draining distance threshold = number The maximum distance from sink to saddle for a spurious sink. The default value is twice the grid spacing. Sink threshold combine logic = AND or OR Use AND if both thresholds must be satisfied to drain the sink; use OR if either can be satisfied. Default is OR. Construct polygons = yes or no Construct elements = yes or no These two settings must be yes for elements to be created. The default for both is yes. Subdivision size = number This sets the target area for automatically subdividing elements. The initial set of elements will be automatically subdivided until each element is smaller than the specified area. 127 Create flow lines out from depressions = yes or no If this is set to yes, stream lines will be created from the depression over its draining saddle as if the depression was a channel head. Default is yes. Spaced flowline spacing = number Apart from creating element networks, TERRAE can operate in a simpler mode by creating a set of flowlines at fixed spacing across the landscape. In this mode, the lines do not connect to form ridge and stream junctions. Set this parameter to a number and set “Construct polygons” to no to use this mode. After its processing, TERRAE writes the following files: Dem.sinks, dem.saddles, dem.peaks The locations and properties of each flat point. Dem.ridgelines, dem.valleylines The set of lines generated by TERRAE Dem.streamj, dem.ridgej The stream and ridge junction points Dem.polygons The outlines of each polygon Dem.elemgeom, dem.elemattr and dem.elemconnect The element description files (described above) Dem.topology Description of the line network topology - each peak, sink, saddle, stream junction 128 and ridge junction and the lines they connect to; sinks and saddles have additional information. These files are all in simple x y z format (apart from the .topology file) that can be plotted directly using GNUPLOT or converted to ungenerate format to read into ARC/INFO. 129 APPENDIX B Program Spatial Soil Water Balance C (SALUS-TERRAE) C DEFINITIONS OF VARIABLES ARE LISTED AT THE END OF TI-fl3 PROGARM USE elemwatbal IMPLICIT NONE CHARACTER I-IEAD*32,FNAl\/IE*8,doystring*3 type(soil_type) 2: Soils(MaxSoils) type(element_type) :: Elements(MaxElems) REAL LAT,LONG,ALT,Press REAL CROPEC,VPDMIN,RGRATE,RTLMR REAL MEANDEP, CUMDEP, Depth REAL DTI‘,LAI,RLD(MaxLayer), RootDep, RootDepIncr REAL SOLAR,TMAX,TMIN,RAIN REAL OutFlow REAL SWSTR REAL SWTot,SWO_26,SW26_77,SW77_Bot INTEGER YR,DOY,L,I,EMDAY,ELEM,NUMELEMS, l elemnum,elemdown,downnum,numsoils,soil REAL FINDY,LAI_TAB(2,IOO),REAL_DOY integer month !AG real a(12), b(12), runoff_slope(maxsoils) !AG !*************************READIN’G INPUT DATA*************************** WRITE(*,*)’ ENTER WEATHER (_.WTH) FILE NAME WITHOUT EXTENSION’ READ(*,’(A8)’)FNAME OPEN(UNTT =1 l,F1LE=trim(FNAME)//’.WTH’, STATUS=’OLD’) OPEN(UNIT=12,FILE=’SOIL.DAT’,STATUS=UNKNOWN’) OPEN(UNIT =14,FILE=’CROP.DAT’,STATUS=’UNKNOWN’) OPEN(UNIT=13,FILE=trim(FNAME)//’.PRN’, STATUS=UNKNOWN’) OPEN(UNIT =1 5,FILE=trim(FNAME)/l’.PRF’, STATUS=UNKNOWN') OPEN(UNIT =17,FILE=trim(FN AME)” ’.PRL’, STATUS=UNKNOWN’) OPENCUNIT=18,FILE=’RUNOFF.DAT’, STATUS=UNKNOWN’) !SPATIAL: Need files to read elements information from OPENCUNIT= l 9,FILE=trim(fname)//’.elemparams’,STATUS=’OLD’) OPEN (UNIT=1 lO,F1LE=trim(fname)// ’.elemsoils’,STATUS=’OLD’) 130 ! ***** Read values for LAI and root depth for given days of the year*** OPEN (UNIT:120,FILE=’LAI.DAT’,STATUS=’OLD’) CALL READXY(120,LAI_TAB) CLOSE(120) READ(l l,’(A32)’)I—IEAD ! Read till the line with Latitude, Longitude and elevation is found 700 READ(1 l,5,ERR=700,END=900)LAT,LONG,ALT IF(LAT .EQ. 0.0 .AND. LONG .EQ. 0.0 .AND. ALT .EQ. 0.0) GOTO 700 5 FORMAT(8X,2F9.3,F7.1) ! Read the runoff parameters: !AG read (UNIT = 18, fmt = 300) (month, a(i), b(i), i = l, 12) !AG 300 format (I, 12(1x, i2, 2(1x, f7 .4), I» !AG READ(14,*) READ(14,*) READ(l4,*)CROPEC,VPDMIN,RGRATE,RTLMR READ(12,*) soil=l DO WHILE(.TRUE.) READ(12,*,END=790) IF(soil.GT.maxsoils) THEN WRITE(*,*) Too many soils: maximum is ’, maxsoils STOP ENDIF READ(12,*) Soils(soil).NLAYR,Soils(soil).COEFW, l Soils(soil).KSMICRO,Soils(soil).PONDMAX READ(12,*) CUIVfl)EP=O. !SPATIAL read element data !*******************mmALIZING AND READING SOILS m*************** WRITE(*,*) ’Soil number ’, soil DO I;l,Soils(soil).NLAYR READ(12,*)Soils(soil).DLAYR(L),Soils(soil).SWLL(L), l Soils(soil).SWDUL(L),Soils(soil).SWSAT(L), 2 Soils(soil).KSMACRO(L),Soils(soil).RI-IF(L), 3 Soils(soil).Il\lITSW(L),Soils(soil).BD(L) CUMDEP=CUMDEP+Soils(soil).DLAYR(L) Soils(soil).ZLayr(L) = CumDep Soils(soil).SWAD(L)=O.44*Soils(soil).SWDUL(L)**2 MEANDEP = (CUMDEP + (CUMDEP - Soils(soil).DLAYR(L)))l2 131 Soils(soil).FLOWUCO(L)=O.63/MEANDEP* *2* l (O.495+EXP(2.804- l0.76*Soils(soil).SWDUL(L))) IF(L.LT.5)THEN Soils(soil).RUCO(L)=Soils(soil).RHF(L)*0.15 ELSE Soils(soil).RUCO(L)=Soils(soil).RHF(L)*O.10 ENDIF RLD(L)=0.0 WRITE(*,10)CUMDEP,Soils(soil).DLAYR(L),Soils(soil).SWLL(L), l Soils(soil).SWDUL(L),Soils(soil).SWSAT (L), 2 Soils(soil).KSMACRO(L),Soils(soil).RHF(L), 3 Soils(soil).INITSW(L) 10 FORMAT(2F6.0,3F6.3,F6. 1,F6.3,F6.3) ENDDO Soils(soil).FLOWUCO(l)=Soils(soil).FLOWUCO( l )* 1 (0.82-4.7*(O.45-Soils(soil).SWDUL(l))**2) WRITE(l3,’(A32)’)I-IEAD WRITE(IS,’(A32)’)HEAD WRITE(17,’(A32)’)HEAD WRIT'E(13,18)LAT,LONG,ALT,CROPEC,VPDMIN,Soils(soil).PONDMAX, l Soils(soil).KSMICRO 18 FORMAT(’LAT’,F5. 1,’ LONG’,F6. 1,’ ALT’,F5.0,’ CROPEC’,F5. 1,’ VPD 1M1N’,F6.2,’ PONDMAX’,F5. l,’ KSMICRO’,F5. l) WRIT‘E( 13,30) WRITE(15,31) WRITE(17,32) WRITE(17,33) WRITE(13, l 5)(Soils(soil).II\lITSW(I),I=l ,12) 15 FORMAT(8X,’- ----- menu/day ---------- cm’,2X,12F5.3) ! Soils(soil).PONDY=0.0 soil=soil+l ENDDO !SPATIAL - Read information from TERRAE about element attributes(area, slope, etc.) ! and connection (topology, which is the downslope elements, etc) from the ’filename.Elements’. 790 CONTINUE numelems = 0 ! The Element array is sparse. The USED variable will tell the program ! which elements have valid data Elements.Used = .FALSE. DO WHILE(.TRUE.) 132 READ(19,* ,end=79 l) elemnum, Elements(elemnum).centreX, Elements(elemnum).centreY,Elements(elemnum).centrez, Elements(elemnum).Area, 3 Elements(elemnum).ndownslope,Elements(elemnum).slope Elements(elemnum).Used = .TRUE. WRITE(*,*) ’Elem ’,elemnum,’has ’, l Elements(elemnum).ndownslope,’connections’ if(Elements(elemnum).ndownslope .gt. maxdown) then WRITE(*,*) Too many connections: maximum is ’, maxdown STOP endif do downnum=l, Elements(elemnum).ndownslope READ(19,*) Elements(elemnum).downslope(downnum), 1 Elements(elemnum).downfrac(downnum) WRITE(*,*) downnum, Elements(elemnum).downslope(downnum), 1 Elements(elemnum).downfrac(downnum) enddo READ(110,*) elemnum,soil Elements(elemnum).Soil = soil do i=l,Soils(Soil).nlayr Elements(elemnum).sw(i) = Soils(Soil).initsw(i) enddo Elements(elemnum).pond = 0.0 Elements(elemnum).pondMax = Soils(Soil).PondMax numelems = max(numelems, elemnum) enddo Nu— ! initialise depth to saturated layer in all elements do elem=l ,numelems soil = Elements(elem).Soil Elements(elem).SatDepth = 9999 Elements(elem).SatLayer = Soils(soil).nlayr + 1 do L=Soils(soil).nlayr,1,-l if (Elements(elem).sw(L) .lt. Soils(soil).swdul(L)) then exit else Elements(elem).SatDepth = Soils(soil).ZlayrG.) - l Soils(soil).Dlayr(L) Elements(elem).SatLayer = L endif enddo ! convert depth from cm from the top of the soil to meters above some ! base level Elements(elem).SatDepth = Elements(elem).CentreZ - 133 l (Elements(elem).SatDepth/100.0) enddo 791 CONTINUE PRESS=101.-0.0107*ALT ROOTDEP=0.0 EMDAY=0 YR=l ! write out Arc LUT for element parameters OPEN(UNIT =1 1 l,FILE=trim(FNAIVIE)//’.area’, STATUS=UNKNOWN’) do elemnum = l,numelems if (Elements(Elemnum).Area .gt. 0) then WRITE(111,*) elemnum, Elements(Elemnum).Area, 1 Elements(Elemnum).centrex,Elements(Elemnum).centrey, 2 Elements(Elemnum).centrez, Elements(Elemnum).Slope endif enddo close(UNIT=l l l) ! Open a Sufer DAT file to store some non-volitle data Open(20,file=trim(fname)//’_fi xed.dat’,status= ’unknown ’) write(20’9(7A l 3),)I'X" :IOY" 9, ”'Z" :mSIOW" 9, ”'Area" 9, l ”'PondMax" ’,’"ElemNum"’ do elemnum = l,numelems if (Elements(elemnum).Used) then soil = Elements(elemnum).Soil ! Initialize variable needed for the water balance CALL ElemCSWB(INIT, DOY, YR, SOLAR, TMAX, TMIN, RAIN, LAT, LONG, ALT, Press, dTT, LAI, RLD, SWSTR, Elemnum, Soil,Soils(soil),Elements(elemnum),CROPEC,RTLMR,RGRATE, 3 Elements,Soils) write(20,’(6fl3.3,ll3)') Elements(Elemnum).centreX, Elements(Elemnum).centrey,Elements(Elemnum).centreZ, Elements(Elemnum).Slope,Elements(Elemnum).Area, Elements(Elemnum).PondMax,ElemNum endif enddo close (20) Nu— WNW 134 ! ***************DAILY LOOP (READ WEATHER AND CALCULATE WATER BALANCE******* DO WHILE(YRGTO) 800 READ(l l,20,ERR=800,END=900)YR,DOY,SOLAR,TMAX,TMIN,RAIN IF (DOY .EQ. O) GOTO 800 20 FORMAT(12,I3 ,4F6. l) REAL_DOY = REAL(DOY) dTT = amaxl(0.,(Tmax + Tmin)/2 -lO.) LAI = FINDY(REAL_DOY,LAI_TAB) do elem=l ,numelems Elements(elem).surfinflow = O Elements(elem).subsurfinflow = O Elements(elem).inflowdepth = 9999 enddo outflow = O WRITE(doystring,’(I3.3)’) doy ! OPEN (UNIT=1 l l ,FILE=trim(FNAME)// ’.wat’//trim(doystrin g), ! l STATUS=UNKNOWN’) ! Open a Sufer DAT file to store some daily data Open(20,file=trim(fname)//trim(doystring)//’.dat’, 1 status=’unknown’) write(20,’(25Al3))"X"’,"'Y"’,’"ElemNum"’,’"Runoff"’,"'RunOn"’, 1 "NetSurfFl"’,"Ponding"’, 2 ’"SW 0-26"’,’"SW 26-77"’,"'SW 77-Bot"’,"'Tot. SW'", 3 ’"SubFlow"’,’"SumSubFlow"’,’"Drainage"’ DO ELEM=1,NUMELEMS if (.not. Elements(ELEM).Used) then ! no element specified, go on to next one cycle endif soil = Elements(ELEM).Soil WRITE(13,*) ’Water balance for element ’,elem I*********************************************************************** ! RATE calculations l***************************III11"."!**************************************** CALL ElemCSWB(RATE, DOY, YR, SOLAR, TMAX, TMIN, RAIN, 1 LAT, LONG, ALT, Press, d'IT, LAI, RLD, SWSTR, Elem, 2 Soil,Soils(soil),Elements(elem),CROPEC,RTLMR,RGRATE, 3 Elements,Soils) I****33*!!!*#************************##1##********************************** 135 INTEGERATION calculations 1*********************************************************************** 1 2 3 CALL ElemCSWB(INTEG, DOY, YR, SOLAR, TMAX, TMIN, RAIN, LAT, LONG, ALT, Press, d'IT, LAI, RLD, SWSTR, Elem, Soil,Soils(soil),Elements(elem),CROPEC,RTLMR,RGRATE, Elements,Soils) ! SPATIAL if there is no downslope element, we let the surface 45me Chm-huts.)— water run off and the subsurface water drain using the slope of the element. Outflow is reported as a volume WRITE(13,*) ’Runoff = ’, Elements(elem).runoff if(Elements(elem).ndownslope .eq. 0) then outflow = outflow + Elements(elem).runoff * Elements(elem).Area else do downnumzl, Elements(elem).ndownslope elemdown = Elements(elem).downslope(downnum) Elements(elemdown).surfinflow = Elements(elemdown).surfinflow + Elements(elem).runoff * Elements(elem).downfrac(downnum)* Elements(elem).Area/Elements(elemdown).Area enddo endif SWTot = sum(elements(Elem).SW(l:Soils(Soil).NLayr) * Soils(soil).DLayr(l :Soils(Soil).NLayr)) SWO_26 = sum(elements(Elem).SW(l:4) * Soils(soil).DLayr(l :4)) SW26_77 = sum(elements(Elem).SW(5:7) * Soils(soil).DLayr(5:7)) SW77_Bot = sum(elements(Elem).SW(8:Soils(Soil).NLayr) * Soils(soil).DLayr(8:Soils(Soil).NLayr)) if (swtot < 0.) then write(*,*) ’Something wrong’ endif write(20,’(2f13.3,113,20F13.3)) Elements(Elem).centrex, Elements(Elem).centrey,ElemElements(Elem).runoff, Elements(Elem).surfinflow, Elements(Elem).surfinflow - Elements(Elem).runoff, Elements(Elem).Pond,SWO_26,SW26_77,SW77_Bot,SWTot, Elements(elem).SumSubInFlow - Elements(elem).SumSubInFlowYest, Elements(elem).SumSubInFlow, 136 7 Elements(elem).FlowD(Soils(soil).NLayr) ENDDO Close(20) ! close(UNIT=l 1 1) ENDDO 900 CONTINUE STOP 30 FORMATC YRDOY RAIN ROFF DRAN ES EP LAIRTDEP SW1 SW2 lSW3 SW4 SW5 SW6 SW7 SW8 SW9 SW10 SW11 SW12) 31 FORMAT(’ YR DOY RAIN RNOFS RNOFD POND FLOW1 FLOW2 FLOW3 F 1LOW4 FLOWS FLOW6 FLOW7 FLOW8 FLOW9 FLOW10 FLOW11 FLOW12’) 32 FORMATC ROOT LENGTH DENSITY FILE «units cm/cm"3’) 33 FORMAT(’ YR DOY DEPl DEP2 DEP3 DEP4 DEP5 DEP6 DEP7 DEP8 D 1EP9 DEPlO DEPll DEPlZ’) END SUBROUTINE READXY (INUNIT,POINTS) Subroutine to read an arbitrary number of xy points and put them into the array POINTS 0000 INTEGER I,J,INUNTT,NUMPOINTS REAL POINT S(2,*) C Read number of points for this equation from the file I = 0 READ(INUNIT ,*,ERR=100,END=100) NUMPOINTS Repeat until error (ie. a line without two real numbers on it) 000 DO I = l,NUMPOlNTS READ(INUNIT‘,*,ERR=IOO,END=IOO) (POIN'I‘S(I,J),I=1,2) END DO 100 CONTINUE C C Place end of points marker C IF (J .EQ. 0) THEN 137 OOOO GOO 000 C POINT S(1,1) = -9999.0 ELSE POINTS(1,J) = -9999.0 ENDIF IF (J .LE. 1) THEN WRITE(*,*) ’ No points read!’ ENDIF RETURN END REAL FUNCTION FINDY (X,POINT S) From the arbitrary X,Y values in the POINTS array find a Y value for a given X value. INTEGER I REAL POINT S(2,*),WEIGHT ,X !variableXadded If input X is less than minimum X set Y to value of first X IF (X .LE. POINTS(1,1)) THEN FINDY = POINTS (2,1) ELSE I: 2 DO WHILE ((X .GT. PONTS(1,I)).AND.(POII\ITS(1,I) .GE. -9998.0)) I=I+l END DO If input X is greater than maximum X set Y to value of last Y IF (POINTS(IJ) .LT. -9998.0) THEN FINDY = POINTS (2,1-1) ELSE . WEIGHT = (X-POINTS(1,(I—l)))/(POINTS(1,I)-POINTS(1,(I-l))) FINDY = POINTS(2,I-l) + (POINT S(2,I) - POINTS(2,I-1))*WEIGHT ENDIF ENDIF RETURN END C C NAILUJ, Subroutine 138 OOOOOOOOOOOOOOOOO0000000000000 Determines Julian date Revision history 1. Written 2 Modified by 3. Header revision and minor changes P.W.W. 2-8-93 INPUT :JULD,NYRCHK LOCAL :RNAME(),NSUM,JCOUNT,NDIF,MON() OUTPUT : NDAY RMON Called : SEHARV SENS SEPLT SETIIVIE OPDAY Calls : None DEFINITIONS RMON : RNAMEO: NSUM : JCOUNT: NDIF : NYRCHK : NDAY : MON () : SUBROUTINE NAILUJ (JULD,NYRCHK,month) IMPLICIT NONE CHARACTER*3 RMON,RNAME(12) INTEGER NSUM,JCOUNT,NDIF,JULD,NYRCHK,NDAY,MON(l2), month DATA MON /3 1 ,28,3 1,30,31,30,3 1,31,30,3l,30,31/ DATA RN AME I’J AN ’,’FEB’,’MAR’,’APR’,MAY’,’JUN’, & ’JUL’,’AUG’,’SEP’,’OCT’,’NOV’,'DEC’/ IF (MOD(NYRCHK,4) .EQ. 0) THEN MON(2) = 29 ELSE 139 MON(2) = 28 ENDIF NSUM = 0 DO JCOUNT = 1, 12 NDIF = JULD - NSUM IF (NDIF .LE. MON (JCOUNT)) THEN NDAY = NDIF RMON = RNAME(J COUNT) month = jcount RETURN ENDIF NSUM = NSUM + MON (JCOUNT) END DO RETURN END C DEFINITIONS C EO Potential evapotranspiration (cm/day) C ESO Potential soil evaporation (cm/day) C ES Actual soil evaporation (cm/day) C EPO Potential plant evaporation (cm/day) C EP Actual plant evaporation (cm/day) C VPD Estimated vapor pressure deficit for a day (kPa) C VPDAY Estimated vapor pressure for a day (kPa) C VP Estimated saturation vapor pressure for a day (kPa) C VPDMIN A crop specific threshold vapor pressure deficit below which the aero- C dynamic component in the potential evaporation equation is zero (kPa) C CROPEC A crop specific constant--slope of the aerodynamic component relation- C ship in the potential evaporation equation (MJ/kPa) C RADCOM The radiation component in the potential evaporation equation (mm/day) C AEROCOMS The aerodynamic component in the potential evaporation equation for C bare soil surfaces (mm/day) C AEROCOMC The aerodynamic component in the potential evaporation equation for C crops surfaces (mm/day) C GODPG Gamma divided by delta plus gamma in the potential evap. equation C PRES The estimated atmospheric pressure based on elevation (kPa) C RTLMR The root length to mass ratio for the empirical root growth routines (cm/g) C RUNOFFS Daily runoff related to lack of surface infiltration (cm/day) C RUNOFFD Daily runoff related to lack of drainage at depth (cm/day) C LAI Leaf area index (cm2/cm2) C LAI_TAB Table of LAI for given dates (input data) C ROOTDEP,TAB Table of root depth for given dates (input data) C RLD(L) Root length density (cm/cm3) 140 C COEFD Empirical coefficient used to calculate new water content after C infiltration (unitless) C CUMDEP Cumulative soil depth at bottom of layer (cm) C MEANDEP Soil depth at center of layer (cm) C SW(L) Soil water content (cm3/cm3) C SWT Temporary variable for a new soil water content (cm3/cm3) C SWTD Temporary variable for calculating water content during drainage (cm3/cm3) C SWTI Temporary variable for calculating water content during infiltration (cm3/cm3) C SWDUL(L) Soil water content at the drained upper limit (cm3/cm3) C SWSAT(L) Soil water content at saturation (cm3/cm3) C RHF(L) Root hospitality factor input from soils file unitless C RUCO(L) Root uptake coefficient unitless C SWDIFU(L) Soil water content difference resulting from upward flow (vol frac) C SWD[FD(L) Soil water content difference resulting from downward flow (vol frac) C SWDIFR(L) Soil water content difference resulting from root uptake (vol frac) C SWDIF(L) Net soil water content difference from all factors (vol frac) C FLOWU(L) Flow rate upward at the layer being evaluated (cm/day) C FLOW(L) Net flow rate at the layer being evaluated (cm/day) C FLOWEXC Excess water that cannot flow downward (cm/day) C FLOWUCO(L) Coefficient for each soil layer used to calculate evaporation rate - unitless C RUCO(L) Coefficient for each soil layer used to calculate root uptake - unitless C REDCOU Reduction coefficient to reduce upward flow to potential soil evap. C REDCOR Reduction coefficient to reduce root uptake to potential plant evap. C ROOTDEP Depth of rooting cm C INFILT Amount of infiltration during a given period (mm). Temp.var. C KSMACRO(L) Saturated hydraulic conductivity with macropores (cm/day). C KSMICRO Saturated hydraulic conductivity at surface without macropores (cm/day). C (Should be lower than KSMACRO) C MAXRAIN Maximum rate of rainfall (cm/day) C NUMSTEPS Number of Steps in a rainstorm C PINF Daily precipitation infiltrated (cm) C PIP Precipitaion infiltrated from ponding (cm) C POND Amount/depth of water held in ponding on the surface (cm) C PONDMAX Maximum amount of water that can be stored in ponding (cm) C RAIN Daily measured rain from weather input file (mm). C is assumed to infiltrate and the ponding routine is bypassed (mm). C PPRECIP Amount of precipitation for a given period (cm). C PRECIP Sum of RAIN and irrigation (cm) convert from m. C RAINDUR Length of rain storm (hrs) C RUNOFFS Daily runoff (cm). C TIME Current time into the storm (days) C TIMESTEP Amount of time per step in a storm (days) C TPAMT For a given period amount of water that could infltrate (cm) C TRUNOFF Amount of water runoff during a given period (cm). 141 C LR Depth increment in depth counter (L) where root depth is located Element Spatial Soil Water Balance module elemwatbal implicit none !**** CONSTANTS !* MAX_LAYER: Maximum number of soil layers INTEGER,PARANIETER INTEGER,PARAMETER INTEGER,PARAMETER INTEGER,PARAMETER INTEGER,PARAMETER INTEGER,PARAMETER INTEGER,PARAMETER REAL,PARAMETER Type soil_type :: MaxLayer = 20 :: MAXELEMS= 100 :: maxdown=20 :: maxsoils=100 ::INlT=1 ::RATE= 2 ::INTEG= 3 ::=RockBD 2.65 !Bulk density of rock material 1n soil integer :: NLAYR ! Number of layers Real :: COEFW Real :: KSMicro Real :: PondMax Real :: DLAYR(MaxLayer),ZLayr(MaxLayer) Real :: BD(MaxLayer) Real :: SWDUL(MaxLayer) Real :: SWLL(MaxLayer) Real :: SWSAT(MaxLayer) Real Real Real Real Real Real :: KSMACRO(MaxLayer) :: SWAD(MaxLayer) :: FLOWUCO(MaxLayer) :: RHF(MaxLayer) :: RUCO(MaxLayer) ' ..InitSW(MaxLayer) end type soil _type type element_type logical :: Integer :: Integer :: Integer :: Used NDownSlope DownSlope(MaxDown) Soil 142 Integer :: EMDAY Integer :: LR Integer :: OutFlowLayer Integer :: SatLayer !Highest (non-perched) saturated layer Real :: DownFrac(MaxDown) Real :: Pond Real :: PondMax Real :: Pondy !Used for Error checking Real :: CentreX, CentreY, CentreZ Real :: Area Real :: Slope Real :: SatDepth !Depth in meters above a base elevation Real :: surfinflow Real :: subsurfinflow,SumSubInFlow,SumSubInFlowYest Real :: inflowdepth Real :: Runoff,RunoffS,Runoffd Real :: SWYest, PondEvap !Used for Error checking Real :: ROOTDEP Real :: RootDepIncr Real :: SW(MaxLayer) Real :: VirtKSMac(MaxLayer) !Effective Sat. Flow for virtical flow Real :: Flow(MaxLayer),FlowD(O:MaxLayer+l),FlowU(O:MaxLayer+l),FlowR Real :: SWDif(MaxLayer),SWDifD(MaxLayer),SWDifU(MaxLayer),SWDifR(MaxLayer) Real :: RedCOU, RedCOR Real :: ES Real :: EP end type element_type contains !Modified Water balance routine from Joe Ritchie. ! Soil -- all the permenant physical properties associated with a soil ! (LL, DUL, etc.) 1 Element -- all the dynamic properties of a soil in a given ! geographical element (SW, ES, EP, etc.) 1m I ! PROGRAM SOIL WATER BALANCE SUBROUTINE ElemCSWB(DYNAMIC, DOY, YR, SOLAR, TMAX, TMIN, RAIN, & LAT, LONG, ALT, Press, dTT, LAI, RID, SWSTR, Elem, & Soil,CurrSoil,CurrElem,CROPEC,RTLMR,RGRATE, & Elements, Soils) 143 1 DEFINITIONS OF VARIABLES ARE LISTED AT THE END OF THE PROGARM INTEGER, intent(in) :: Dynamic, DOY, Yr, Elem, Soil REAL, intent(in) :: Solar, TMax, TMin, Rain, Lat, Long, Alt, Press, d'IT, LAI Real, intent(in) :: RLD(MaxLayer) REAL, intent(out) :: SWStr type(soil_type), intent(inout) :: CurrSoil,Soils(MaxSoils) type(element_type), intent(inout) :: CurrElem,Elements(MaxElems) Real, intent(in) :: CROPEC,RTLMR,RGRATE Integer :: I,L,LNU Integer :: Month Integer :: DownNum,DownElem,DownSoil Real :: Error !Used for Error checking Real :: Precip,WInf Real :: a(12), b Real :: PBSWabove,PESWaboveFAC(MaxLayer) Real :: MEANDEP, CUMDEP, acoef, ncoef Real :: runoff_slope Real :: Porosity, WatFilPor, halfpesw, swdryfac, swwetfac, coldfac Real :,:TMEAN LATHEAT, GODPG, RadCom, VPD, VPDAY, VP, AEROCOMS, AEROCOMC Real :: EO, EPO, ESO, ESD, PotentialSurfaceEvap Real :: SWN(MaxLayer), SWTI, SWTD, SWT Real :: FLOWEXC Real :: RootPres,Coef_Wat_Ext Real :: RootAct(MaxLayer) Real :: CumDif Real :: Depth,SatDP,DownAmt,KSSat,KSSatDown, KSEff, LateralFlow Real :: DeltaXDeltaZ Real :: DepthFact, RadSlope, ZLayr, TotSubSurfOutFlow I*************************INITIAIJZATION****************************** * ' File reading will be handled outside of this routine for 1 spatial application 1 IF (DYNAMICEQJNIT) THEN !********************************************************************** 1 ' OPEN(UNII‘=2,FILE=’SOIL.DAT’,STATUS=’UNKNOWN') 1 OPEN(UNIT=4,FILE=’CROP.DAT’,STATUS=UNKNOWN’) ' OPEN (UNIT =3 ,FILE=trim(FN AIME)” ’.PRN’, STATUS=UNKNOWN’) 144 1 1 1 1 1 1 OPEN (UNIT=5 ,FILE=trim(FNAME)//’.PRF’, STATUS=UNKNOWN’) OPEN (UNIT: 1 5 ,FIle=tri m(FNAME)// ’.Err’, STATUS=UNKNOWN ’) READ(4,*) READ(4,*) READ(4,*)CROPEC,RGRATE 1C *******************NHMLIZWG AND READING SOILS mE*************** 0‘ 0‘ a- I- a- READ(2,*) READ(2,*) READ(2,*)NLAYR,Slope READ(2,*) Compute PondMax from Slope CurrElem.PondMax = Max(1./(0.2 + 0.35*CurrElem.Slepe),CurrElem.PondMax) Compute the efective virtical saturated flow conductivity First convert Slope from % SIOpe to radians RadSlope = CurrElem.Slope * 100.0/90.0 * 0.00175 ZLayr = 0 Do L = 1,CurrSoil.NLayr At zero depth there is no difference bewteen KSMacro and KSMacVirt at two meters the correction factor will be sin(slope) ZLayr = ZLayr + CurrSoi1.DLayr(L) DepthFact = 1.0 - 0.5 * (ZLayr/ 100.0) CurrElem.VirtKSMac(L) = CurrSoil.KSMacro(L) * (DepthFact + & (1 - DepthFact) * sin(RadSlOpe)) EndDo Set coefficients for Runoff, currently N. Central USA eqn. from Gerakis b=-O.30 Do I=1,12 a(I)=MIN(-0.015,-0.15*sin(2.*3.1417*(I+2.)/12.)-0.05) ENDDO PESWabove = 0. CUMDEP=0. DO L=l,CurrSoil.NLayr CUMDEP=CUMDEP+CurrSoiLDLayr(L) CurrSoil.SWAD(L)=0.44*CurrSoil.SWDUL(L)**2 !Air dry soil water contents WRITE(* ,10)CUMDEP,CurrSOil.DLayr(L),CurrSoil.SWLL(L),CurTSoil.SWDUL(L),Cu rrSoil.SWSAT(L). & 145 CurrSoil.KSMacro(L),CurrSoil.RI-IF(L),CurrSoil.InitSW(L),CurrSoil.BD(L),CurrSoil.S WAD(L) 1 Set factors for root water extraction based on upper soil water PESWabove=PESWabove + max(0.,CurrSoil.InitSW(L)- CurrSoil.SWLL(L))*CurrSoil.DLayr(L) PESWaboveFACa.)=MIN( l .,MAX(0.,exp(2.2-0.62*PESWabove)-0.07)) RootAct(L) = 0 11ndicates that roots have not yet been active,layer L 1 Set root water uptake coefficients for each layer IF(L.LT.5)THEN CurrSoil.RUCO(L)=CurrSoil.RI-IF(L)*0.30 ELSE CurrSoil.RUCO(L)=CurrSoil.RHF(L)*O. 10 ENDIF 1 Compute coefficients for soil-limited evaporation from depths MEANDEP = (CUNHDEP + (CUMDEP - CurrSoil.DLayr(L)))l2 ncoef=-2.2-l.04*CurrSoil.SWDUL(L) acoef=3.20*CurrSoilSWDUL(L)-0. 15 CurrSoil.FlowUCo(L)=acoef*MEANDEP* *ncoef ENDDO CurrSoil.FlowUCo(1)=CurrSoil.FlowUCo(1)*(0.82-4.7*(0.45- CurrSoil.SWDUL(1))**2) 1adjust top layer CurrElem.ROOTDEP = 0.0 CurrElem.RootDepIncr = 0.0 CurrElem.PondY=0.0 CurrElem.EMDAY=0 CurrElem.Pond = 0.0 !AG CurrElem.SumSubInFlow = 0.0 CurrElem.SumSubInFlowYest = 0.0 CurrElem.SatDepth = 9999 CurrElem.SatLayer = CurrSoil.nlayr + 1 do L=CurrSoil.nlayr,l,-1 if (CurrElem.sw(L) .lt. CurrSoil.swdul(L)) then exit else CurrElem.SatDepth = CurrSoil.Z1ayr(L) - & CurrSoil.Dlayr(L) CurrElem.SatLayer = L endif enddo 146 1 1 1 1 1 1 1 1 WRITE(3,’(A32)’)HEAD WRITE(S,’(A32)’)HEAD WRIT E(3, l 8)LAT,LONG,ALT,CROPEC,CurrElem.PondMax WRITE(3,30) WRITE(S,31) WRITE(3,15)(SW(I),I=1,12) I****************************RATE CALCULATIONS“**************************** ELSEIF (DYNAMIC.EQ.RATE) THEN I*********************************************************************** *** 1 CurrElem.PondEvap = 0.0 1BDB CurrElem.SWYest = sum(CurrElem.sw(l:CurrSoil.nlayr)*CurrSoil.DLayr(1:CurrSoil.nlayr)) 1BDB !*********** Runoff Rate calCUIationS************************************ PRECIP=RAIN*0.1 1CONVERTS RAIN FROM min/day to cm/day Added for Spatial Precip = Precip + CurrElem.surfinflow 1BDB CurrElem.runoffS = 0.0 !AG winf = 0.0 !AG if (precip .gt. 0.0) then !AG 1 Slope of curve of rain vs. runoff: !AG call nailuj (doy,yr,month) !AG runoff_slope = MIN(1.0, exp (MAX (-1e1, a(month)* & !AG CurrElem.VirtKSMac(l) + b))) !AG CurrElem.Runoffs = max(runoff_slope * (precip + CurrElem.Pond - CurrElem.PondMax),0.0) !AG !AG winf = min(precip - CurrElem.Runoffs, CurrElem.VirtKSMac(1)) !AG CurrElem.Pond = CurrElem.Pond + precip - CurrElem.Runoffs - winf !AG if (CurrElem.Pond .gt. CurrElem.PondMax) then !AG CurrElem.Runoffs = CurrElem.Runoffs + CurrElem.Pond - CurrElem.PondMax CurrElem.Pond = CurrElem.PondMax !AG endif !AG elseif (CurrElem.Pond .gt. 0.0) then !AG winf = min(CurrElemPond, CurrElem.VirtKSMac(1)) !AG CurrElem.Pond = CurrElem.Pond - winf !AG endif !AG 147 if (CurrElem.Pond .lt. -1e-5) then !AG write (*,*) ’ **** Negative ponding = ’, CurrElem.Pond !AG endif !AG CurrElem.FLOWD(O)= winf !AG CurrElem.Runoffd=0. 1******************RATE OF ROOT GROWTH, DOWNWARD (cm/day)************* IF(LAI.GT.0.)THEN 1******************Initializing Root Depth on Day of Emergence******* CurrElem.EMDay=CurrElem.EMDay+1 IF(CurrElem.E1VIDay.EQ.l)TI-IEN CurrElem.RootDep=7. CurrElem.LR=2 ENDIF Porosity = l. - CurrSoil.BD(CurrElem.LR)/RockBD 1jwj WatFilPor = CurrElem.SW(CurrElem.LR)/Porosity 1jwj halfpesw = CurrSoil.SWLL(CurrElem.LR) + 0.5*(CurrSoil.SWDUL(CurrElem.LR)—CurrSoi1.SWLL(CurTElem.LR)) 1jwj SWDRYFAC=max(O. ,min( 1 .,(CurrElem.SW(CurrElem.LR)- CurrSoil.SWLL(CurrElem.LR)) & /(halfpesw-CurrSoil.SWLL(CurrElem.LR)))) SWWETFAC=MIN(1.,5.*(1.-WatFilPor)) 1jwj coldfac = 1. CurrElem.RootDepIncr = RGRATE * dTT * min(CurrSoil.RHF(CurrElem.LR),swdryfac, & swwetfac,coldfac) ENDIF !**************** POTENTIAL EVAPOTRANSPIRATION CALCULATIONS********** TMEAN=TMAX*0.5+TMIN*0.5 LATHEAT=(25.01-TMEAN/42.3) 1THIS IS FOR CM/DAY UNITS GODPG=(0.598-3.0E-05*PRESS)-0.017*TMEAN-1-l .49E-04*TMEAN* *2 IF(GODPG.GT. l .O)GODPG=1.0 1 THE ABOVE EXPRESSION IS NOT GOOD AT VERY LOW TEMPERATURES RADCOM=( l -GODPG)*SOLAR*0.63 148 VPDAY=0.61 1*EXP(l7.27*T1VIEAN/(T1V[EAN+237.3)) VP=0.61 l*EXP(17.27*TIVIIN/(TIVIIl\l+237.3)) VPD=VPDAY-VP IF(LAI.GT.0)THEN AEROCOMS=GODPG*15.0*VPD**1.5*EXP(-0.75*LAI) ESO=(RADCOM*EXP(-0.40*LAI)+AEROCOMS)/LATHEAT AEROCOMC=GODPG*CROPEC*VPD** l .5 E0 = (RADCOM+AEROCOMC)/LATHEAT EPO=EO*(1 .-EXP(-0.92*LAI)) if (lai .lt. 0.5) then E0 = ESO endif ELSE AEROCOMS=GODPG*15.0*VPD**1.5 ESO=(RADCOM+AEROCOMS)ILATHEAT EPO=0.0 E0 = ESO ENDIF 1 End of LAI if statement for PET calculations CurrElem.FlowR=0. CUMDEP=0. ESD=O. PESWabove=0. SWN = CurrElem.SW 1 SPATIAL: 1 Move water laterally in or out of the saturated zone CurrElem.SumSubInFlowYest = CurrElem.SumSubInFlow CurrElem.SumSubInFlow = CurrElem.SumSubInFlow + & CurrElem.subsurfinflow 1 If water is flowing in start from the top of the water table and 1 move up If (CurrElem.subsurfinflow > 0) then L = CurrElem.SatLayer do while ((L > 0) .and. (CurrElem.subsurfinflow > 0)) SWN(L) = SWN(L) + min(CurrElem.subsurfinflow, & (CurrSoil.SWSAT(L) - SWN(L)» CurrElem.subsurfinflow = CurrElem.subsurfinflow - & min(CurrElem.subsurfinflow,(CurrSoilSWSAT(L) - SWN(L)» L = L -1 enddo 1 If there is no room in the soil to store the water put it into 149 1 the ponded water If (CurrElem.subsurfinflow > 0) then CurrElem.pond = CurrElem.Pond + CurrElem.subsurfinflow CurrElem.subsurfinflow = 0 EndIf 1 If water is to be removed start from the bottom layer and work up ElseIf (CurrElem.subsurfinflow < 0) then L = CurrSoil.NLayr do while ((L > 0) .and. (CurrElem.subsurfinflow < 0)) SWN(L) = SWN (L) - min(-CurrE1em.subsurfinflow, & Max((SWN(L) - CurrSoil.SWDUL(L)),0.0)) CurrElem.subsurfinflow = CurrElem.subsurfinflow + & min(-CurrE1em.subsurfinflow, & Max((SWN(L) - CurrSoil.SWDUL(L)),0.0)) L = L - 1 enddo If (CurrElem.subsurfinflow < 0) then Write(*,*) ’ Outward lateral flow exceeded available water!’ CurrElem.subsurfinflow = 0 EndIf Endif DO L=l,CurrSoil.NLayr 1*******COMPUTE DOWNWARD FLOW (INFILTRATION AND REDISTRIBUTION)****** SWN(L)=SWN(L)+CurrElem.FlowD(L-l)/CurrSoil.DLayr(L) IF(SWN(L).LE.CurrSoil.SWDUL(L))THEN CurrElem.FlowDG.)=0. ELSE SWTI=CurrSoil.SWDUL(L)-I-0.25*(1-EXP(-0.05*CurrE1em.FlowD(L-1))) SWN(L)=MIN(SWTI,CurrSoil.SWSAT(L),SWN(L)) SWTD=CurrElem.SW(L)-0.55*(CurrElem.SW(L)-CurrSoilSWDUL(L)) SWN(L)=MAX(SWN(L),SWTD) CurrElem.FlowD(L)=Cu1rElem.FlowD(L-1)-(SWN(L)- CurrElem.SW(L))*CurrSoi1.DLayr(L) 1 ALTER FLOW DOWN AND WATER CONTENTS ABOVE LAYER FOR RESTRICTED DOWNWARD FLOW IF(CurrElem.FlowD(L).GE.CurrElem.VirtKSMac(L))THEN 150 FLOWEXC=CurrElem.FlowD(L)-CurrElem.VirtKSMac(L)- & (CurrSoil.SWSAT(L)-SWN(L))*CurrSoil.DLayr(L) CurrElem.FlowD(L)=CurrElem.VirtKSMac(L) IF(FLOWEXC.GE.0.)THEN SWN(L)=CurrSoil.SWSAT(L) DO LNU=L-1,l,-l SWT=SWN(LNU)+FLOWEXC/CurrSoil.DLayr(LNU) IF(SWT.GT.CurrSoil.SWSAT(LNU))THEN FLOWEXC=(SWT-CurrSoil.SWSAT(LNU))*CurrSoi1.DLayr(LNU) CurrElem.FlowD(LNU)=CurrElem.FlowD(LNU)-FLOWEXC- & (CurrSoil.SWSAT(LNU)-SWN(LNU))*CurrSoil.DLayr(LNU) SWN(LNU)=CurrSoil.SWSAT(LNU) ELSE CurrElem.FlowD(LNU)=CurrElem.FlowD(LNU)-FLOWEXC FLOWEXC=0. SWN(LNU)=SWT ENDIF CurrElem.SWDifD(LNU)=SWN(LNU)-CurrElem.SW(LNU) ENDDO ELSE . SWN(L)=CurrSoil.SWSAT(L)+FLOWEXC/CurrSoil.DLayr(L) FLOWEXC=0. ENDIF CurrElem.Pond=CurrE1em.Pond+FLOWEXC IF(CurrElem.Pond.GT.CurrElem.PondMax)'I‘I-IEN CurrElem.Runoffd=CurrElem.Pond—CurrElem.PondMax CurrElem.Pond=CurrElem.PondMax ELSE CurrElem.Runoffd=0.0 ENDIF ENDIF ENDIF CurrElem.SWDifD(L)=SWN(L)-CurrElem.SW(L) CurrElem.RUNOFF=CurrElem.Runoffs+CurrElem.Runoffd . THIS IS TOTAL RUNOFF FROM SURFACE AND SUBSURFACE FLOW CAUSES 1*************COMPUTE SOIL LIMITED ROOT WATER [MAKE********************* IF(LAI.GT.0.0)THEN '********************Determine bottom layer with roots*********************** CUMDEP=CUMDEP+CurrSoi1.DLayr(L) IF(CUMDEP-CurrSoil.DLayr(L).LT.CurrElem.RootDep)CurrElem.LR=L 1Layer where roots are growing 151 IC*********************************UP’I‘AKE****************************3|: ****** RootPres = 1.0 IF(L .EQ. CurrElem.LR) then RootPres=MIN( l .0,(CurrElemRootDep-(CUMDEP- CurrSoil.DLayr(CurrElem.LR))) & /CurrSoil.DLayr(CurrElem.LR)) Endif IF(L .GT. CurrElem.LR) RootPres = 0.0 PESWabove=PESWabove+MAX(0.,CurrElem.SW(L)- CurrSoil.SWIL(L))*CurrSoil.DLayr(L) IF(RootAct(L) .EQ. 1) THEN 1 Roots Already have been active in uptake in this layer PESWaboveFAC(L)=MH\l(1 .,MAX(O. , & -exp(—4.8+0.3*PESWabove)+l.03)) Else 1 Roots still have not become active in uptake from this layer yet PESWaboveFAC(L)=MIN( l .,MAX(0., & exp(2.2-0.62*PESWabove)-0.07)) IF(PESWaboveFAC(L)*RootPres .GE. 0.90) RootAct(L) = l Endif Coef_Wat_Ext=CurrSoilRUCO(L)*RootPres*PESWaboveFAC(L) IF(CurrElem.SW(L).LT.CurrSoi1.SWLL(L))THEN CurrElem.SWDifR(L)=0. ELSE CurrElem.SWDifR(L)=Coef_Wat_Ext*MAX(0.,CurrElem.SW(L)- CurrSoil.SWLL(L)) ENDIF ELSEIF(CurrElem.EMDay.GT. 1)THEN 1RESET ROOT GROWTH PARAMETERS WHEN LAI = O CurrElem.EMDay=0. CurrElem.RootDepflD ENDIF CurrElem.FlowR=CurrElem.FlowR+CurrElem.SWDifR(L)*CurrSoil.DLayr(L) !*********CALCULATE SOIL LMTED SOIL EVAPORATION RATE************** CurrElem.SWDifU(L)=CurrSoil.FlowUCo(L)*(CurrElem.SW(L)- CurrSoil.SWAD(L)) IF(CurrElem.SW(L).GT.CurrSoilSWDUL(L))THEN 152 CurrElem.SWDifU(L)=CurrElem.SWDifU(L)*(3 .-2.*(CurrSoil.SWSAT(L)- CurrElem.SW(L))/& (CurrSoil.SWSAT(L)-CurrSoil.SWDUL(L))) ENDIF ESD=ESD+CurrElem.SWDifU(L)*CurrSoil.DLayr(L) ENDDO 1 end of layer do loop IF(CurrElem.RootDep.GT.CUMDEP)THEN CurrElem.RootDep=CUMDEP ENDIF 1*****Calculate actual soil evaporation (cm/day) and water stress factor* if (CurrElem.Pond > 0) then if (CurrElem.Pond > ESO) then PotentialSurfaceEvap = ESO CurrElem.Pond = CurrElem.Pond - ESO CurrElem.PondEvap = ESO else PotentialSurfaceEvap = ESO - CurrElem.Pond CurrElem.PondEvap = CurrElem.Pond CurrElem.Pond = 0 endif else PotentialSurfaceEvap = ESO endif IF(ESD.GT.PotentialSurfaceEvap)THEN CurrElem.RedCOU=PotentialSurfaceEvap/ESD CurrElem.ES=PotentialSurfaceEvap ELSE CurrElem.RedCOU=l . CurrElem.ES=ESD ENDIF SWSTR = 1.0 IF(CurrElem.FlowR.GT.0.)THEN IF(EPO+CurrElem.ES .GT. EO)THEN EPO=EO-CurrElem.ES if (EPO < 0) then EPO = 0 endif ENDIF IF(CurrElem.FlowR.GE.EPO)Tl-IEN CurrElem.RedCOR=EPO/CurrElem.FlowR CurrElem.EP=EPO ELSE 153 CurrElem.RedCOR= l .0 CurrElem.EP=CurrElem.FlowR SWSTR = CurrElem.FlowR/EPO ENDIF ELSE CurrElem.RedCOR=0. CurrElem.EP=0. ENDIF 1*********************INTEGRATION CALCULATIONS******************************** ELSEIF (DYNAMIC.EQJNTEG) THEN !********************************************************************** 1 !*****************Update ROOt Depth******************************* CurrElem.RootDep = CurrElem.RootDep+CurrElem.RootDepIncr !********************************************************************** 1 1 !******CALCULATE FWAL 8011. WATER DEFERENCE****************** CUMDIF=0. CurrElem.FlowU(CurrSoil.NLayr)=0. DO L=CurrSoil.NLayr,l,-l CurrElem.SWDif(L)=CurrElem.SWDifD(L)— & CurrElem.SWDifU(L)*CurrElem.RedCOU— & CurrElem.SWDifR(L)*CurrElem.RedCOR CurrElem.FlowU(L-1)=CurrElem.FlowU(L) + & CurrElem.SWDifU(L)*CurrElem.RedCOU*CurrSoil.DLayr(L) CurrElem.FLOW(L)=CurrElem.FlowD(L)-CurrElem.FlowU(L) CurrElem.SW(L)=CurrElem.SW(L)+CurrElem.SWDif(L) IF(CurrElem.SW(L).LT.CurrSoil.SWAD(L))THEN WRITE(*,19)L,CurrElem.SW(L),CurrElem.SWDifU(L), & CurrElem.SWDifD(L),CurrElem.SWDifR(L),L ENDIF CUMDIF=CUMDIF+CurrElem.SWDif(L)*CurrSoi1.DLayr(L) END DO !************ lateral subsurface flow Added for Spatial ******************* 1 calculate depth of saturated layer if any CurrElem.SatDepth = 9999 CurrElem.SatLayer = CurrSoil.nlayr + 1 do L—-CurrSoil.nlayr,l,-l if (CurrElem.sw(L) .lt. CurrSoil.swdul(L)) then 154 exit else CurrElem.SatDepth = CurrSoil.Z1ayr(L) - & CurrSoil.Dlayr(L) CurrElem.SatLayer = L endif enddo 1 convert depth from cm from the top of the soil to meters above some 1 base level CurrElem.SatDepth = CurrElem.CentreZ - & (CurrElem.SatDepth/100.0) TotSubSurfOutFlow = 0.0 1 write(3,*) ’Sat depth = ’, satdp if (CurrElem.SatLayer .LE. CurrSoil.nlayr) then 1 For all saturated layers above the water table of down hill elements Do L = CurrElem.SatLayer, CurrSoil.nlayr do DownNum = l, CurrElem.ndownslope DownElem = CurrElem.downslope(downnum) if (CurrElem.SatDepth .gt. Elements(DownElem).SatDepth) then deltaZ = CurrElem.SatDepth - Elements(DownElem).satdepth deltaX = sqrt((CurrElem.centrex-Elements(DownElem).centrex)**2+ & (CurrElem.centrey-Elements(DownElem).centrey)**2) Kssat = CurrSoil.KSMacro(L) kssatdown = Soils(Elements(DownElem).soi1).KSMacro(E1ements(DownElem).SatLayer) kseff = 2/(1/Kssat + llkssatdown) LateralFlow = kseff * CurrElem.downfrac(DownNum) * (DeltaZ/DeltaX) * & (CurrElem.Area/Elements(DownElem).Area) 1 If the flow is out of the current element but there is 1 not enough water above DUL to drain If ((LateralFlow . gt. 0) .and. & (LateralFlow .gt. (CurrElem.sw(L)-CurrSoil.SWDUL(L))))then Elements(DownElem).subsurfinflow = & Elements(DownElem).subsurfinflow + & CurrElem.SW(L)-CurrSoilSWDUL(L) TotSubSurfOutFlow = TotSubSurfOutFlow + & CurrElem.SW(L)-CurrSoilSWDUL(L) CurrElem.SW(L) = CurrSoil.SWDUL(L) 155 1 If flow is out of the current element and there is 1 enough water to drain Elself (LateralFlow .gt. 0) then CurrElem.SW(L) = CurrElem.SW(L) - LateralFlow Elements(DownElem).subsurfinflow = & Elements(DownElem).subsurfinflow + LateralFlow TotSubSurfOutFlow = TotSubSurfOutFlow + & LateralFlow 1 If flow is into the current element but there is not enough 1 space to hold the water Elself ((LateralFlow .lt. 0) .and. & (-Latera1Flow .gt. (CurrSoil.SWSat(L) - CurrElem.sw(L))))then _ Elements(DownElem).Subsurfinflow = & Elements(DownElem).subsurfinflow - & (CurrSoil.SWSat(L) - CurrElem.SW(L)) TotSubSurfOutFlow = TotSubSurfOutFlow- & (CurrSoil.SWSat(L) _- CurrElem.SW(L)) CurrElem.SW(L) = CurrSoil.SWSat(L) 1 If flow is into the current element and there is enough 1 space to hold it Elself (LateralFlow .lt. 0) then CurrElem.SW(L) = CurrElem.SW(L) - LateralFlow Elements(DownElem).subsurfinflow = & Elements(DownElem).subsurfinflow + LateralFlow TotSubSurfOutFlow = TotSubSurfOutFlow + & LateralFlow Endlf endif enddo EndDo 1 Old code 1! look up ksmacro of uppermost saturated layer 1 depth 2 0 1 do i=1,CurrSoil.nlayr 1 depth = depth + CurrSoil.dlay'r(i) 1 if (depth > satdp) then 1 CurrElem.outflowlayer = i 156 kssat = CurrSoil.ksmacro(i) exit endif enddo write(*,*) ’Ksat = ’, kssat,’for layer ’, CurrElem.outflowlayer 1 1 1 1 1 1 1 depthlost = 0.0 1 do downnum = l, CurrElem.ndownslope 1 downelem = CurrElem.downslope(downnum) 1 DownSoil = Elements(DownElem).Soil 1 write(*,*) DownSlope connection to ’, downelem 1 11 subsurface flow happens when elevation of saturated layer (2 - satdp) 11 is higher than the elevation of the downslope element’s saturated layer 1 1 deltaZ = (CurrElem.centrez - satdp) - & 1 (Elements(DownElem).centrez - Elements(DownElem).satdepth) 1 write(*,*) ’deltaZ = ’,deltaZ 1 if (deltaZ .gt. 0) then 1 11 set depth of inflow equal from elevation of highest contributing inflow 11 adjusting depth for elevation difference of the two elements if (CurrElem.centrez - satdp .gt. & Elements(DownElem).centrez - & Elements(DownElem).inflowdepth) then Elements(DownElem).inflowdepth = satdp - & (Elements(DownElem).centrez - CurrElem.oentrez) endif 1 1 1 1 1 1 1 1 1 depth = 0 1 do i=1,CurrSoil.nlayr 1 depth = depth + CurrSoil.dlayr(i) 1 if (depth > Elements(DownElem).inflowdepth) then 1 kssatdown = Soils(DownSoil).KSMacro(i) 1 exit 1 endif 1 enddo 1 kseff = 2/(1/kssat + llkssatdown) 1 write(*,*) kssatdown = ’,kssatdown,’ kseff = ’,kseff 1 11 calculate amount of water moving laterally subsurface 1 1 deltaX = sqrt((CurrElem.centrex-Elements(DownElem).centrex)**2-1- & 1 (CurrElem.centrey—Elements(DownElem).centney)**2) 1 SurfSlope = deltaZ / deltaX 157 Slopeangle = atan(SurfSlope) flowdistance = kseff * sin(slopeangle) * & CurrElem.downfrac(downnum) depthlost = depthlost + flowdistance Elements(DownElem).subsurfinflow = & Elements(DownElem).Subsurfinflow + & flowdistance * Currelemarea / elements(downelem).area write(*,*) ’flowdist =’,flowdistance endif enddo CurrElem.sw(CurrElem.outflowlayer) = & CurrElem.SW(CurrElem.outflowlayer) - depthlost .— I- C- 0‘ 0‘ O- O- O- O- I- 0‘ I- 1‘ 0‘ endif ERROR=PRECIP+CurrElem.PondY-CurrElem.Runoffs-CurrElemRunoffd- & CurrElem.Pond-CurrElem.FlowD(CurrSoil.NLayr)+ & CurrElem.SWYeSt- & sum(CurrE1em.SW(l:CurrSoil.NLayr)*CurrSoil.DLayr(1:CurrSoil.NLayr))- & CurrElem.EP-CurrElem.ES-CurrElem.PondEvap-TotSubSurfOutFlow 1BDB IF(ERROR**2.GT.0.0000000001)THEN WRIT'E(15,22)YR,DOY,Elem,ERROR WRITE(*,22)YR,DOY,Elem,ERROR ENDIF CurrElem.PondY=CurrElem.Pond WRIT'E(13,25)YR,DOY,PRECIP,CurrElem.RUNOFF,CurrElem.Pond, & CurrElem.FlowD(CurrSoil.NLayr),CurrElem.ES, & CurrElem.EPLAI,CurrElem.RootDep, SWSTR, & (CurrElem.SW(I),I=l ,12) 1 WRITE(S,26)YR,DOY,PRECIP,CurrElem.Runoffs,CurrElem.Runoffd,CurrElem.Pond,( FLOW(I),I=1,12) 1 1 ENDIF RETURN 10 FORMAT(2F6.0,3F6.3,F6.1,F6.3,3F6.3) 15 FORMAT(8X,’ -------- -cm/day ---------- cm’,8X,12(1x,F7.3)) 18 FORMAT(LAT’,F5. l,’ LONG’,F6. l,’ ALT’,F5.0,’ CROPEC’,F5. 1,’ VPD’ & ’MIN’,F4.2,’ PONDMAX’,F5. l) 19 FORMAT(’ UNACCEPT ABLE RESULT, SW BELOW AIR DRY, SW(L),UPWARD FLOWDOWN FLOW,ROOT UPTAKEDEPTH L’,I3,4F7.3,I3) 158 22 FORMATC YR’,I4, ’DAY’,I4,’ ELEMENT’,I4,’ WATER BALANCE ERROR (cm) =’,F9.6) 25 FORMAT(214,F6.2,6F5.2,F6.0,F6.3,2X,12(lx,F8.7)) 30 FORMAT(’ YR DOY RAIN ROFF POND DRAN ES EP LAI RTDEP SWSTR’& ,’ SW1 SW2 SW3 SW4 SW5 SW6 SW7’& ,’ SW8 SW9 SW10 SW11 SW12) 26 FORMAT(214,4F6.2, 12F7 .4) 31 FORMAT(’ YR DOY RAIN RNOFS RNOFD POND FLOWl FLOWZ FLOW3 F’, & ’LOW4 FLOWS FLOW6 FLOW7 FLOWS FLOW9 FLOWIO FLOW11 FLOW12’) END SUBROUTINE ElemCSWB 1 DEFINITIONS 1 BO Potential evapotranspiration (cm/day) 1 ESD Amount of water available for soil evaporation? (cm/day) 1 E80 Potential soil evaporation (cm/day) 1 ES Actual soil evaporation (cm/day) 1 EP0 Potential plant evaporation (cm/day) 1 EP Actual plant evaporation (cm/day) 1 VPD Estimated vapor pressure deficit for a day (kPa) 1 VPDAY Estimated vapor pressure for a day (kPa) 1 VP Estimated saturation vapor pressure for a day (kPa) 1 CROPEC A crop specific constant--Slope of the aerodynamic component relation- 1 ship in the potential evaporation equation (MJ/kPa) 1 RADCOM The radiation component in the potential evaporation equation (mm/day) 1 AEROCOMS The aerodynamic component in the potential evaporation equation for 1 bare soil surfaces (mm/day) 1 AEROCOMC The aerodynamic component in the potential evaporation equation for 1 crops surfaces (mm/day) 1 GODPG Gamma divided by delta plus gamma in the potential evap. equation 1 PRES The estimated atmospheric pressure based on elevation (kPa) 1 RUNOFFS Daily runoff related to lack of surface infiltration (cm/day) 1 RUNOFFD Daily runoff related to lack of drainage at depth (cm/day) 1 LAI Leaf area index (cm2/cm2) 1 LAI_TAB Table of LAI for given dates (input data) 1 ROOTDEP_TAB Table of root depth for given dates (input data) 1 COEFD Empirical coefficient used to calculate new water content after 1 infiltration (unitless) 1 CUMDEP Cumulative soil depth at bottom of layer (cm) 1 MEANDEP Soil depth at center of layer (cm) 1 SW(L) Soil water content (cm3/cm3) 1 SWT Temporary variable for a new soil water content (cm3/cm3) 1 SWTD Temporary variable for calculating water content during drainage (cm3/cm3) 159 1 SWT'I Temporary variable for calculating water content during infiltration (cm3/cm3) 1 SWDUL(L) Soil water content at the drained upper limit (cm3/cm3) 1 SWSAT(L) Soil water content at saturation (cm3/cm3) 1 RHF(L) Root hospitality factor input from soils file unitless 1 RUCO(L) Root uptake coefficient unitless 1 SWDIFU(L) Soil water content difference resulting from upward flow (vol frac) 1 SWDIFD(L) Soil water content difference resulting from downward flow (vol frac) 1 SWDIFR(L) Soil water content difference resulting from root uptake (vol frac) 1 SWDIF(L) Net soil water content difference from all factors (vol frac) 1 FLOWU(L) Flow rate upward at the layer being evaluated (cm/day) 1 FLOW(L) Net flow rate at the layer being evaluated (cm/day) 1 FLOWEXC Excess water that cannot flow downward (cm/day) 1 FLOWUCO(L) Coefficient for each soil layer used to calculate evaporation rate - unitless 1 RUCO(L) Coefficient for each soil layer used to calculate root uptake - unitless 1 REDCOU Reduction coefficient to reduce upward flow to potential soil evap. 1 REDCOR Reduction coefficient to reduce root uptake to potential plant evap. 1 ROOTDEP Depth of rooting cm 1 INFILT Amount of infiltration during a given period (mm). Temp.var. 1 KSMACRO(L) Saturated hydraulic conductivity with macropores (cm/day). 1 MAXRAIN Maximum rate of rainfall (cm/day) 1 NUMSTEPS Number of steps in a rainstorm 1 PINF Daily precipitation infiltrated (cm) 1 PIP Precipitaion infiltrated from ponding (cm) 1 POND Amount/depth of water held in ponding on the surface (cm) 1 PONDMAX Maximum amount of water that can be stored in ponding (cm) 1 RAIN Daily measured rain from weather input file (mm). 1 is assumed to infiltrate and the ponding routine is bypassed (mm). 1 PPRECIP Amount of precipitation for a given period (cm). 1 PRECIP Sum of RAIN and irrigation (cm) convert from m. 1 RAINDUR Length of rain Storm (hrs) 1 RUNOFFS Daily runoff (cm). 1 TIME Current time into the storm (days) 1 TIMESTEP Amount of time per step in a storm (days) 1 TPAMT For a given period amount of water that could infltrate (cm) 1 TRUNOFF Amount of water runoff during a given period (cm). 1 LR Depth increment in depth counter (L) where root depth is located 1 acoef 1 ncoef 1 (HT Thermal time on current day, oC.d (base of 10 0C) 1 BD Bulk density of soil(g/cm3) 1 RockBD Bulk density of the solid material in soil (g/cm3) 1 Porosity Fraction of soil volume that is pore space 1 WatFilPor Fraction of soil pore space filled with water 1 HalfPESW Soil water content when half of available water remains 1 RootDepIncr Daily growth rate of roots, vertical direction (cm/d) 160 1 coldfac Factor reducing root growth in a layer due to cold temperature (0-1) 1 Slope Average slope of the soil surface (%) 1 RootPres Indicator of root presence in a layer, fraction (0 to l) 1 PESWabove volume of plant avaialble water above a layer (for computing root activity) 1 PESWaboveFAC Factor (0-1) to indicate root activity in water uptake 1 Coef__Wat__Ext Relative rate at which water is extracted from a layer (l/d) 1 RootAct Indicator that roots have been taking up water, each layer 1 end module elemwatbal 161 REFERENCES Abbott, M.B., J.C. Bathrust, J.A. Cunge, P.E. O’Connell, and J. Rasmussen. 1986. An introduction to the European System-System Hydrologique Europeen, “SHE”, 2, Structure of a physically-based, distributed modeling system. J. Hydrol. 87: 61-77. Adams, R.M., Rosenzweig, R.M. Peart, J.T. Ritchie, B.A. McCarl, J .D. Glycer, R.B. Curry, J .W. Jones, K]. Boote, and L.H. Allen, Jr. 1990. Global Climate change and US. Agriculture. Nature 345: 219-224. Addiscott, T.M., and R. J. Wagenet. 1985. Concepts of solute leaching in soils: A review of modeling approaches. J. Soil Sci. 36: 411-424. Barnes, E.M., P.J. Pinter Jr., B.A. Kimball, G.W. Wall, R.L. LaMorte, DJ. Hunsaker, F. Adamsen, S. Leavitt, T. Thompson, and J. Mathius. 1997. Modification of CERES- Wheat to accept leaf area index as an input variable. Written for Presentation at the 1997 ASAE Annual International Meeting, Sponsored by ASAE, Minneapolis, MN, August 10-14, 1997. Paper No. 973016. Batchelor, W.D., and J .O. Paz. 1998. Process-Oriented Crop Growth Models as a Tool to Evaluate Spatial Yield Variability. In: Proceedings of the First International Conference on Geospatial Information in Agriculture and Forestry, Vol. 1. pp. 1198- 205. 1-3 June 1998, Lake Buena Vista, Florida, USA. ERIM International. Bell, J .C., R.L. Cunningham, and M.W. Havens. 1992. Calibration and validation of a soil-landscape model for predicting soil drainage class. Soil Sci. Soc. Am. J. 56: 1860-1866. Bell, J .C., R.L. Cunningham, and M.W. Havens. 1994. Soil drainage class probability mapping using a soil-landscape model. Soil Sci. Soc. Am. J. 58: 464-470. Beven, K.J., and MJ. Kirky. 1979. A physically-based variable contributing area model of basin hydrology. Hydrol. Sci. Bull., 24, 43-69. 162 Blackmer, A.M., and SE. White. 1996. Remote sensing to identify spatial patterns in optimal rates of nitrogen fertilization. In: P.C. Robert, RH. Rust, and WE. Larson (eds.), Proceedings of the Third International Conference on Precision Agriculture, Minneapolis, MN, 23-26 June 1996. ASA Miscellaneous Publication, pp. 33-41. ASA, CSSA, and SSSA, Madison, WI. Booltink, H.W.G., and J. Verhagen. 1997. Using decision support systems to optimize barley management on spatial variable soil. In: M.J. Kropff, P.S. Teng, P.K. Aggarwal, J. Bouma, B.A.M. Bouma, J.W. Jones, and H.H. van Laar (ed.), Application of systems approaches at the field level, pp. 219-233. Kluwer Academic Publishers, Dordrecht, The Netherlands. Boote, K.J., J.W. Jones, G. Hoogenboom, and NB. Pickering. 1998. CROPGRO model for grain legumes. In: G.Y. Tsuji, G. Hoogenboom, and PK. Thornton (eds.), Understanding Options for Agricultural Production, pp. 99-128. Kluwer Academic Publishers, Dordrecht, Netherlands. Braga, R.P., J .W. Jones, and B. Basso. 1999. Weather Variability in Site-Specific Management Profitability: A Case Study. In: P.C. Robert, R.H. Rust, and WE. Larson (ed.), Precision Agriculture, pp. 1853-1863. ASA-CSSA-SSSA, Madison, Wisconsin, USA. Brown, D.G. 1994. Predicting vegetation types at treeline using topography and biophysical disturbance variables. Journal of Vegetation Science, 5: 641-656. Cambardella, C.A., T.S. Colvin, D.L. Karlen, S.D. Logsdon, E.C. Berry, Radke, T.C. Kaspar, T.B. Parkin and DB. Jaynes. 1996. Soil property contributions to yield variation pattern. In P.C. Robert, R.H. Rust, and WE. Larson (eds.). Proceedings of the Third International Conference on Precision Agriculture. ASA, CSSA, SSSA, Inc., Madison, WI. pp 417-424. Carlson, TN, and DA. Ripley. 1997. On the Relation between NDVI, Fractional Vegetation Cover, and Leaf Area Index. Remote Sens. Environ. 62: 241-252. Carter, JR. 1988. Digital representations of topographic surfaces. Photogrammetric Engineering and remote sensing, 54: 1577-1580. 163 Cassel, D.K., L.F. Ratliff, and IT. Ritchie. 1983. Models for estimating in-situ potential extractable water using soil physical and chemical properties. Soil. Sci. Soc. Am. J. 47:764-769. Chong, S.-K. 1983. Calculation of sorptivity from constant rate rainfall infiltration measurement. Soil Sci. Soc. Am. J. 47: 627-630. Cochrane, T.A., and DC. Flanagan. 1999. Assessing water erosion in small watershed using WEPP with GIS and digital elevation models. Journal of Soil and Water Conservation, 54: 678-686. Cora, J.E., P.J. Pierce, B. Basso, and IT. Ritchie. 1999. Simulation of within-field variability of corn yield with Ceres-Maize model. In: P.C. Robert, R.H. Rust, and WE. Larson (ed.), Precision Agriculture, pp. 1309-1319. ASA-CSSA-SSSA, Madison, Wisconsin, USA. Dunne, T., and RD. Black. 1970. Partial-area contributions to storm runoff in a small New England watershed, Water Resour. Res. 6: 1296-1311. Dunne, T. 1983. Relation of field studies and modeling in the prediction of storm runoff. J. Hydrol. 65: 25-48. ESRI. 1993. Surface modeling with TIN, in “FrameViewer (TM) 3.1X”, the on-line help manual of ARC/INFO System that was developed and is supported by Environmental Systems Research Institute, Inc. (ESRI). Feddes, R.A., P.J. Kowalik, and H. Zaradny. 1978. Simulation of field water use and crop yield. Pudoc, Wageningen. Freeze, RA. 1972. Role of subsurface flow in generating surface runoff, 2, Upstream source areas, Water Resour. Res. 8: 1272-1283. Gabrielle, B., S. Menasseri, and S. Houot. 1995. Analysis and field-evaluation of the CERES models water balance component. Soil Sci. Am. J. 59: 1403-1412. 164 Gallant, J .C., and JP Wilson. 1996. TAPES-G: a grid-based terrain analysis program for the environmental sciences. Computers and Geosciences, 22: 713-722. Gallant, J .C., and M.F. Hutchinson. 1997. Scale dependence in terrain analysis. Mathematics and Computers in Simulation 43: 313-321. Gallant, J .C. 1999. TERRAE: A new element network tool for hydrological modelling. 2nd Inter-Regional Conference on Environment-Water: Emerging Technologies for Sustainable Water Management, Lausanne, Switzerland, 1-3 September 1999. Gamma Design Software, GS+: geostatistics for the environmental scienc. Version 3.01. Gamma Design Software, Plainwell, MI, p.44, 1994. Gardner, W. 1960. Dynamic aspects of water availability to plants. Soil Sci.: 89-93. Gerakis, A., and J .T. Ritchie. 1998. Simulation of atrizine leaching in relation to watertable management using the CERES model. Journal of Environmental Management 52: 241-258. Gessler, P.E., I.D. Moore, NJ. McKenzie, and P.J. Ryan. 1995. Soil landscape modeling and spatial prediction of soil attributes. International Journal of Geographical Information Systems, 9: 421-432. Gessler, P.E., NJ. McKenize and M.F. Hutchinson. 1996. Progress in soil-landscape modelling and spatial prediction of soil attributes for environmental models. Proceedings of the Third International Conference/Workshop on Integrating GIS and Environmental Modelling. National Centre for Geographic Information and Analysis, Santa Barbara, CA, USA. Hall, GE, and CG. Olson. 1991. Predicting variability of soils from landscape models. In: M.J. Mausbach and LP. Wilding (eds.), Spatial Variabilities of Soils and Landforms, pp. 9-24. SSSA Special Publication No. 28. Hoogenboom, G.J., J .W. Jones, P.W. Wilkens, W.D.Batchelor, W.T. Bowen, L.A. Hunt, N. Pickering, U. Singh, D.C. Godwin, B. Baer, K.J. Boote, J.T. Ritchie, and J .W. White. 1994. Crop models. In: G.Y. Tsuji, G. Uehara, and S. Balas (eds.). 1994. DSSAT v3. Vol. 2-2. University of Hawaii, Honolulu, Hawaii. 165 Horton, R. E. 1933. The role of infiltration in the hydrologic cycle, Eos Trans. AGU, 14: 446-460. Hutchinson, M.F. 1996. A locally adaptive approach to the interpolation of digital elevation models. Proceedings of the Third lntemational Conference/Workshop on Integrating GIS and Environmental Modelling. National Centre for Geographic Information and Analysis, Santa Barbara, CA, USA. Hutchinson, M.F. 1997. ANUDEM Version 4.6. http://cres.anu.edu.au/software/anudem46.html. IBSNAT. 1984. lntemational Benchmark Sites Network for Agrotechnology Transfer. Department of Agronomy and Soil Science, University of Hawaii, Honolulu, HI, USA. ISBRAM. 1991. lntemational Board for Soil Research and Management. Evaluation for sustainable land management in the developing world; Volume I Towards the development of an lntemational framework. IBSRAM Proceedings no. 12. Jackson, RD. and AR. Huete. 1991. Interpreting vegetation indeces. Preventive Veterinary Medicine 11: 185-200. Jasinski, M.F. 1996. Estimation of subpixel vegetation density of natural regions using satellite multispectral imagery. IEEE Trans. Geosci. Remote Sens. 34: 804-813. Johannsen, C.J., M.F. Baumgardner, P.R. Willis, and PG. Carter. 1998. Advances in Remote Sensing Technologies and Their Potential Impact on Agriculture. In: Proceedings of the First International Conference on Geospatial Information in Agriculture and Forestry, Vol. 1. pp. I6-11. 1-3 June 1998, Lake Buena Vista, Florida, USA. ERIM lntemational. Jones, A.J., L.M. Mielke, C.A. Bartles, and CA. Miller. 1989. Relationship of landscape position and properties to crop production. J. Soil Water Conserv. 44: 328-332. 166 Jones, J .W., K]. Boote, G. Hoogenboom, S.S. Jagtap, and G.G. Wilkerson. 1989. SOYGRO V5.42, Soybean crop growth simulation model user’s guide. Florida Agricultural Experiment Station Journal N ° 8304. University of Florida, Gainesville, Florida, USA. Jones, J .W., and J .T. Ritchie. 1991. Crop growth models. In: G]. Hoffman, T.A. Howell, and K.H. Soloman (eds.), Management of farm irrigation systems. ASAE, St. Joseph, MI. Kirkby, M.J., A.C. Irneson, G. Bergkamp, and L.H. Cammeraat. 1996. Scaling up processes and models. Journal of Soil and Water Conservation, 51: 391-396. Kovacs, G.J., T. Németh, and J .T. Ritchie. 1995. Testing Simulation Models for the Assessment of Crop Production and Nitrate Leaching in Hungary. Agricultural Systems 49: 385-397. Kumler, MP. 1994. An intensive comparison of triangulated irregular networks (TIN s) and digital elevation models (DEMS). Cartographica, 31: 1-99. Laflen, J.M., W.J Elliot,. D.C., Flanagan, CR, and MA. Nearing. 1997. WEPP- Predicting water erosion using a process-based model. Journal of Soil and Water Conservation, 52: 96-103. Lafolie, F.L., Bruckler, A.M. deCockborne, and C. Laboucarie. 1997. Modeling the water transport and nitrogen dynamics in irrigated salad crops. Irrigation Science. 17: 95- 104. Liu, Q., and A. Huete, 1995. A feedback based modification of the NDVI to minimize canopy background and atmosphere noise. IEEE Trans. Geosci. Remote Sens. 33: 457-465. Maune, D. 1996. DEM extraction editing matching and quality control techniques. In: C.W. Greve (ed.), Digital Photogrammetry: An Addendum to the Manual of Photogrammetry. Bethesda, MD: American Society for Photogrammetry and Remote Sensing. 167 Miller, CL, and RA. Laflamme. 1958. The digital terrain model — theory and application. Photogrammetric Engineering 24: 433-42. Moore, I.D., E.M. O’Laughlin, and G]. Burch. 1988. A contour-based topographic model for hydrological and ecological application. Earth Surf .Processes Landforms, 13: 305-320. Moore, I.D., and RB. Grayson. 1989. Hydrologic and Digital Terrain Modelling using vector elevation data, EOS Trans. AGU 70: 1091. Moore, I.D., and RB. Grayson. 1991. Terrain-Based catchment Partitioning and Runoff Prediction Using Vector Elevation Data. Water Resour. Res. 27: 1177-1191. Moore, I.D., and M.F. Hutchinson: 1991. Spatial Extension of Hydrologic process Modelling. Proc. Int. Hydrology and Water Resources Symposium. Institution of Engineers-Australia 91/22, pp. 803-808. Moore, I.D., R.B. Grayson, and AR. Ladson. 1991. Digital Terrain Modelling: A review of hydrological, geomorphological and biological applications. In: K.J. Beven and ID. Moore (eds.), 1991. Terrain Analysis and distributed Modelling in hydrology. J .Wiley & Sons. Moore, I.D., and PE. Gessler, G.A. Nielsen and G.A. Peterson. 1993. Terrain Analysis for Soil Specific Crop Management. In: Soil Specific Crop Management. ASA- CSSA-SSSA. Moore, I.D., J .C. Gallant, and L. Guerra. 1993. Modelling the spatial variability of hydrological processes using GIS in HydroGIS 93: Application of Geographic Information Systems in Hydrology and Water Resources Management. Proceedings of Vienna Conference, April 1993. IAHS publ. No. 211. p. 161-169. Moore, I.D., A. Lewis, and J .C. Gallant. 1993. Terrain attributes: estimation methods and scale effects. In: A.J. Jakeman, M.B. Beck, and M.J. McAleer (eds.), Modeling Change in Environmental Systems, pp. 191-214. New York: John Wiley and Sons. Moore, I.D., P.E. Gessler, G.A. Nielsen, and G.A. Peterson. 1993. Soil attribute prediction using terrain analysis. Soil Sci.Soc. Amen]. 57: 443-452. 168 Moran, M.S., Y. Inoue, EM. Barnes. 1997. Opportunities and Limitations for Irnage- Based Remote Sensing in Precision Crop Management. Remote Sensing of the Environment Vol. 61 : 319-346. Mueller, TC. 1998. Ph.D. Dissertation, Department of Crop and Soil Science, Michigan State University, East Lansing, MI, USA. Nix, HA. 1983. Minimum Data Sets for Agrotechnology Transfer. Proceedings of the lntemational Symposium on Minimum Data Sets for Agrotechnology Transfer. 21-23 March, ICRISAT AP 502324. India. O’Laughlin, EM. 1986. Prediction of surface saturation zones in natural catchments by topographic analysis, Water Resour. Res. 22: 229-246. Onstad, C.A., and D.L. Brakensiek. 1968. Watershed simulation by stream path analogy. Water Resour. Res. 4: 965-971. Onstad, C. A. 1973. Watershed flood routing using distributed parameters, In: E.S. Shultz, V.A. Koelzer, and K. Mahmood (eds.), Floods and Droughts, pp. 418-428. Water Resources Publications, Fort Collins, CO. Passioura, J .B. 1996. Simulation Models: Science, Snake Oil, Education or Engineering? Agronomy Joum. 88: 690-695. Paz, J .O., W.D. Batchelor, T.S. Colvin, S.D. Logsdon, T.C. Kaspar, and D.L. Karlen. 1997. Analysis of Water Stress Effects Causing Spatial Yield Variability in Soybeans. Trans. ASAE 41: 1527-1534. Paz, J .O., W.D. Batchelor, B.A. Babcock, T.S. Colvin, S.D. Logsdon, T.C. Kaspar, and D.L. Karlen. 1999. Model-based technique to determine variable rate nitrogen for corn. Agric. Syst. 61: 69-75. Pierce, P.J., D.D. Wamcke, and M.W. Everett. 1995. Yield and nutrient availability in glacial soils of Michigan. p. 133-151. In P.C. Robert et al. (ed.) Proc. Of the Int. Conf. on Site Specific Management for Agricultural Systems, 2““. Bloomington/Minneapolis, WI. 169 Pierce, P.J., N .W. Anderson, T.S. Colvin, J .K. Schueller, D.S. Humberg and NB. McLaughlin. 1997. Yield Mapping. In: F.J. Pierce and E.J. Sadler (ed.), The State of Site Specific Management for Agriculture, p. 211-243. ASA Misc. Publ., ASA, CSSA, and SSSA, Madison, WI. Pierce, P.J., and P. Nowak. 1999. Aspects of Precision Agriculture. Advances in Agronomy, 67: 1-86. Pierce, F.J., and DD. Warncke. 2000. Soil and Crop Response to Variable-Rate Limiting for Two Michigan Fields. Soil Sci. Soc. Am J. 64: Price, J.C. 1992. Estimating Vegetation Amount from Visible and Near Infrared Reflectances. Remote Sens. Environ. 41: 29-34. Priestley, C.H.B., and R.J. Taylor. 1972. On the assessment of the surface heat flux and evaporation using large-scale parameters. Monthly Weather Review 100: 81-92. Refsgaard, J.C., S.M. Seth, J.C. Bathurst, M. Erlich, B. Storm, G.H. J¢rgensen, and S. Chandra. 1992. Application of the SHE to catchments in India, Part 1. General results. J. Hydrol. 140: 1-23. Richardson, C.W., and J.T. Ritchie. 1973. Soil Water Balance for Small Watershed. Transactions of the ASAE. 16: 72-77. Ritchie, J .T. 1972. Model for predicting evaporation from a row crop with incomplete cover. Water Resour. Res. 8: 1204-1213. Ritchie J .T., D.C. Godwin, and S. Otter-Nacke. 1985. CERES-Wheat: A simulation model of wheat growth and development. Texas A&M University Press, College Station, TX. Ritchie, J .T. 1985. A user-oriented model of the soil water balance in wheat. In: W. Day and R.K. Atkin (ed.), Wheat growth and modeling, pp. 293-305. Series A: Life Sciences Vol. 86. Plenum Press, NY. 170 Ritchie J .T. 1986. Using Simulation Models for Predicting Crop Performance. Symposium on the Role of Soils in Systems Analysis for Agrotechnology Transfer. August 1986. ISSS Congress. Hamburg, Germany. Ritchie, J .T. 1987. Using crop models as a decision support system to reduce nitrate leaching. In: F.M. D’Itri and LG. Wolfson (ed.), Rural groundwater contamination, pp. 179-192. Lewis Publishers, Chelsea, Michigan. Ritchie J .T. and J .Crum. 1989. Converting soil survey characterization data into IBSNAT crop model input. In: J. Bouma and A.K. Berget (eds.), Land Qualities in Space and Time, 155-167.Wageningen, The Netherlands. Ritchie, J .T., and U. Singh, BC. Godwin and L. Hunt. 1989. A user’s guide to CERES Maize-V.2.10- Muscle Shoals, AL. lntemational Fertilizer Development Center. Ritchie J .T. and BS Johnson. 1990. Soil and Plant factors affecting evaporation. In: B.A. Stewart and DR. Nielsen (eds.), Irrigation of Agricultural Crops. Madison, WI. Amer. Soc. of Agronomy. Ritchie J .T., and D.S. Nesrrrith. 1991. Temperature and Crop Development. In: Hanks and Ritchie (eds.), Modelling plant and soil systems. Agronomy 31: 5-29. Ritchie, J .T. 1991. Specifications of the ideal model for predicting crop yield. In: R.C. Muchow and J .A. Bellamy (ed.), Climatic Risk in Crop Production: Model and Management for the serni-arid tropic and subtropic, p. 97-122. Proc. Int]. Symp. S. Lucia, Brisbane, Queensland, Australia. July 2-6, 1990. CAB lntemational Wallingford, U.K. Ritchie, J.T. 1998. Soil water balance and plant water stress. In: G.Y. Tsuji, G. Hoogenboom, and PK. Thornton (eds.), Understanding Options for Agricultural Production, pp. 41-54. Kluwer in cooperation with ICASA, Dordrecht/Boston/London. Ritchie, J .T., A. Gerakis,. and A. Suleiman . 1999. Simple Model To Estimate Field- Measured Soil Water Limits. Tran. ASAE 42: 1609-1614. 171 Rosensweigh, C., and D. Hillel. 1998. Climate Change and the Global Harvest: Potential Impacts of the Greenhouse Effect on Agriculture. (Cary, N .C.: Oxford University Press). Sadler, E.J., and Russell G. 1997. Modeling Crop Yield for Site Specific Management. p. 69-80. In: E]. Pierce and E.J. Sadler. (ed.), The State of Site Specific Management for Agriculture. ASA Misc. Publ., ASA, CSSA, and SSSA, Madison, WI. Shanholtz, V0. and TM. Younos. 1994. A soil water balance model for no-tillage and conventional till systems. Agricultural Water management 26: 155-168. Singh, P., G. Alagarswamy, G. Hoogenboom, P. Pathak, S.P. Wani, and SM. Virrnani. 1999. Soybean-Chickpea rotation on Vertic Inceptisol, 2, Long-term simulation of waterbalance and crop yields. Field Crops Research 63: 225-236. Speight, J .G. 1974. A parametric approach to landforrn regions. In: E.H. Brown and D.L. Linton (eds.), The Progress of Geomorphology, 7: 213-230. The British Geographers. Sudduth, K.A., S.T. Drummond, S.J. Birrell, and NR. Kitchen. 1996. Analysis of spatial factors influencing crop yield. In: P.C. Robert, R.H. Rust, and WE. Larson (ed.), Proceedings of the Third lntemational Conference on Precision Agriculture, pp. 129- 140. ASA, CSSA, SSSA. Madison, WI. Swaney, D.P., J .W. Jones, W.G. lBogess, G.G. Wilkerson, and J .W. Mishoe. 1983. Real- time Irrigation Decision Analysis using Simulation. Transactions of the ASAE 26: 562-568. Tsuji, G.Y., G. Uehara, S. Balas. 1994. DSSAT: A Decision Support System for Agrotechnology Transfer. Vol. 3. University of Hawaii, Honolulu, HI, USA. USDA. 1996. Soil survey Laboratory Methods Manual. Soil Survey Investigation Report No.42. Vers. 3.0. Washington, DC: USDA. Verhagen, A., H.W.G. Booltink, and J. Bouma. 1995. Site-specific management: Balancing production and environmental requirements at farm scale. Agric. Syst. 49: 369-384. 172 Verhagen, J ., and Bouma J. 1997. Modeling Soil Variability. p. 55-68. In: P.J. Pierce and E.J. Sadler. (ed.), The State of Site Specific Management for Agriculture. ASA Misc. Publ., ASA, CSSA, and SSSA, Madison, WI. Vertessy, R.A., J .T. Hatton, P.J. O’Shaughnessy, and M.D.A. Jayasuriya. 1993. Predicting water yield froma amountain ash forest catchment using a terrain analysis based catchment model. J. Hydrol. 150: 665-700. Wang, M. and AT. Hjelmfelt. 1998. DEM Based Overland Flow Routing Model. Journal of Hydrologic Engineering, Vol. 3, No. 1: 1-8. Weibel, R., and M. Heller. 1991. Chapter 19: Digital Terrain Modeling. In Goodchild, M.F. and D. Rhind (eds.), Geographic Information Systems, Principles and Applications. New York: Taylor and Francis. White, I., and P. Broadbridge. 1988. Constant rate rainfall infiltration: A versatile nonlinear model, 2, Applications of solutions, Water Resour. Res. 24: 155-162. White, I., M.J. Sully, and MD. Melville. 1989. Use and hydrological robustness of Time- to-incipient-ponding. Soil Sci. Soc. Am. J. 29: 1343-1346. Wiegand, C.L., A. J. Richardson, D.E. Escobar, and AH. Gerbermann. 1991. Vegetation Indeces in Crop Assessments. Remote Sens. Environ. 35: 105-119. World Resources Institute. 1992. World Resources, 1992-1993. Oxford University Press. Xu, T. 1999. Mapping Soil Properties with Landscape Analysis. Ph.D. Dissertation, Centre for Resource and Environmental Studies, Australian National University, Canberra, ACT, Australia. Yates, D.N. 1996. WatBal: An Integrated Water Balance Model for Climate Impact Assessment of River Basin Runoff. Water Resources Development, 12: No.2. 121- 139. 173 Young, R. A., CA. Onstad, D.D. Bosch, and WP Anderson. 1996. AGNPS: A nonpoint source pollution model for evaluating agricultural watersheds, Journal of Soil and Water Conservation, 44, 168-193. 174 a ..___Dr_ w r. . __;. "7111111111111111111