PHYSICAL-MATHEMATICAL MODELS IN THE CONTEXT OF INTEGRATING TEACHING, RESEARCH AND EXTENSION IN THE AGRICULTURAL SCIENCES

         !"#$% &''$ (''                                             !"#$%&'()&*& *)++*&, -    .!/0/  10 ,()21340"(#3.56   (%7,  ,  8.&(),( 1  8 3 -  9":!;8 9()(<.!/0/  10  1 BOOK PHYSICAL-MATHEMATICAL MODELS IN THE CONTEXT OF INTEGRATING TEACHING, RESEARCH AND EXTENSION IN THE AGRICULTURAL SCIENCES: PRECIPITATION, RUNOFF AND EVAPOTRANSPIRATION Nelson Mario Victoria Bariani Cassiane Jrayj de Melo 2023 2 SUMMARY 1 INTRODUCTION 5 1.1 Presentation 5 1.2 History of the SWAT Model 8 2. METHODOLOGY 11 3. PRESENTATION OF RESULTS 15 3.1 Precipitation 16 3.1.1 Precipitation as a Physical Variable - Units and Relationships 17 3.1.2 Modelling Precipitation Intensity in the SWAT Model 19 3.1.2.1 Variables and Parameters Involved in the Representation of a Storm 20 3.1.2.2 Interpreting the relationships (equations) between the variables in the physical model 22 3.1.3.1 Programming Rainfall Intensity in Scilab 25 3.1.3.2 Normalised Intensity Distribution (page 931) 29 3.1.3.2.2 Scilab programming of the Normalised Intensity 32 3.1.3.2.3 Calculation of Total Precipitation and its Duration 36 3.1.4 Conclusion 40 3.1.5 Precipitation: Calculating Daily Values 40 3.1.6 CALCULATED PRECIPITATION: Computational definition of the occurrence of a wet or dry day. 42 3.1.7 Results 56 3.1.8 Complete rain generation algorithm 63 3.1.9 Maximum amount of precipitation every half hour 74 3.1.10 Maximum monthly half-hourly rainfall 74 3.1.11 Daily value of the maximum half-hourly rainfall fraction 77 3.1.12 Conclusion on Precipitation Generation Models 82 3.2 Surface Flow Modelling 82 3.2.1 Introduction - Basic Concepts and Relationships 83 3.2.1 Physical process associated with surface runoff 84 3 3.2.3 Factors influencing surface runoff 85 3.2.4 Quantities associated with surface runoff 85 3.2.5 Estimating surface runoff 86 3.2.6 Other methods 89 3.2.7 Water balance 90 3.3 Flow Modelling According to the SWAT Model Production of Daily Flow Values 93 3.3.1 Calculation of the flow of Arroio Olaria and Arroio Cambaí using the Curve Number - SCS method 93 3.3.2 Description of the application of the SCS Method 98 3.3.3 Programming Flow Calculations using the SCN Method in Scilab 100 3.3.4 Conclusion 105 3.4 Maximum Flow Rate 106 3.5 Concentration Time 107 3.5.1 Time of Concentration of Surface Runoff on the Ground 108 3.5.2 Calculating Flows within the Time of Concentration 115 3.5.4 Programming the Previous Equations in the SCILB Environment 119 3.6 Runoff Delay and Transmission Losses 121 3.6.1 Delay 121 3.6.2 Transmission Losses 121 3.6.3 Conclusions 127 3.7 Evapotranspiration 128 3.7.1 Priestley-Taylor, Hargreaves and Penman-Monteith methods 128 3.7.1.1 Introduction 128 3.7.1.2 Concepts and Methods 129 3.7.1.3 HNET Net Radiation Calculation 132 3.7.1.3.1 Extraterrestrial radiation 133 3.7.1.4 Solar Radiation with a Cloudless Sky 139 3.7.1.5 Daily Solar Radiation 141 3.7.1.6 Solar Radiation per Hour 141 3.8 Temperature 147 4 3.9 Water Temperature, Wind Speed and Variable Nomenclature in SWAT 152 3.9 Real Evapotranspiration 169 3.10 Transpiration 171 3.11 Soil Sublimation and Evaporation 172 3.12 Sublimation 173 3.13 Soil Water Evaporation 174 3.14 Daily Temperature and Solar Radiation Calculations 177 3.14.1 Statistical Simulation of Daily Solar Radiation and Maximum and Minimum Temperatures in Southern Brazil 177 3.15 Elevation Ranges, Climate Change and the Incorporation of Weather Forecasting 180 3.15.1 Elevation bands 180 3.16 Climate customisation 184 3.17 Climate Change 184 3.18 Solar radiation and temperature 188 3.18.1 Daily Residuals 188 4 FINAL CONCLUSIONS 207 11 2. METHODOLOGY The content of the SWAT model has been worked on in various academic curricula with increasing intensity since 2017. In the Physics and Hydrology classes, the 634-page theoretical documentation was separated into topics that were assigned to each student. The method was called the parallel approach, where each student is in charge of analysing a specific topic and presenting a report on it. The figure shows an example of a table with the division of topics, which is the starting point for the work. Table 1. Division of SWAT Model Themes for Parallel Study The activity was launched and managed from virtual environments such as Moodle, which was mainly used until 2019, or Google Classroom, which has been more powerful and modern since 2020. In order to prioritise and encourage students to carry out the research work requested, after the initial launch of the activity and the discussion of some examples corresponding to the beginning of the activity by some students, it was organised in the form of a test, because 12 in the academic environment students are conditioned to dedicate their best efforts to test situations. The following table shows an example of the Physics exam text. PHYSICAL TEST JUNE 2022 (300 points) THE TEST WILL BE OPEN FOR SEVERAL WEEKS, BUT EACH STUDENT WILL HAVE TO SHARE THEIR PROGRESS THROUGHOUT THE PROCESS. GRADES WILL BE UPDATED. 1. ORGANISE THE INFORMATION AVAILABLE: a. FROM THE SWAT THEORETICAL DOCUMENTATION IN PORTUGUESE (634 PAGES) SAVE IN PDF (PRINT, PDF) THE PAGES OF THE SUBJECT OF YOUR INTEREST. DO THE SAME FOR THE ENGLISH DOCUMENTATION (CHECK THE PAGES) b. IN THE ENGLISH DOCUMENTATION CUT OUT THE FIGURES AND EQUATIONS OF YOUR WORK USING THE SNIPPING TOOL IN WINDOWS. c. PREPARE A TEXT DOCUMENT CONTAINING THE AVAILABLE THEORETICAL INFORMATION AND FIGURES, WITH ADDITIONAL EXPLANATIONS DERIVED FROM YOUR RESEARCH. d. BUILD A TABLE CONTAINING INFORMATION ABOUT THE VARIABLES THAT ARE WORKED 13 ON IN YOUR RESEARCH SUBJECT. INCLUDE: DEFINITION/CONCEPTUALISATION, UNITS USED IN THE SWAT MODEL AND OTHER TYPICAL UNITS, SOME TYPICAL VALUES FOR THAT VARIABLE, EXAMPLES OF UNIT CONVERSION USING THESE TYPICAL VALUES AND THE CHAIN RULE. 2. DESCRIBE THE RELATIONSHIPS BETWEEN THE VARIABLES IN THE FORM OF EQUATIONS AND TRY TO INTERPRET THEM. WRITE THE EQUATIONS USING SCILAB SYNTAX (SYMBOLS +,-,*, /, ^) 3. PROGRAM THE EQUATION INTO SCILAB AND MAKE GRAPHS EXPLAINING THE BEHAVIOUR OF THE VARIABLE ANALYSED. AT AN EARLY STAGE, ANY VALUES CAN BE USED, JUST TO TEST HOW THE CALCULATIONS WORK, BUT LATER, REASONABLE VALUES FOR THE PROGRAMMED VARIABLES, SUITABLE FOR OUR REGION, MUST BE FOUND IN THE TECHNICAL-SCIENTIFIC LITERATURE (CALIBRATION). 4. CARRY OUT A DIARY AND GENERAL FINAL REPORT EXPLAINING THE RESULTS, WHICH CAN BE IN THE FORM OF AUDIO OR VIDEO. Generally, the subjects and objectives of the work were initially quite complex for most of the students, but the first students' deliveries were used to develop all the steps required in detail in class, with 14 everything recorded on videos available in the virtual environment. In this way, with the support of Information and Communication Technologies, in the style of pedagogy applied in distance learning (DE), subjects were gradually dealt with in class, paving the way for progress in individual research. This process took place over several semesters, allowing us to evolve in our understanding of the different physical processes and their programming. In the Physics subject, the emphasis was on describing and understanding the variables of the system under study, their units, concepts and relationships with other variables, programming the equations in the Scilab environment. In the subject of Hydrology, which is taught in the semester immediately after Physics, most of the students took up the work they had done in Physics in order to make more progress with the programming and emphasise the relationship between the subject under consideration and the hydrological cycle, which is the fundamental aspect behind the SWAT model. In the Interdisciplinary Laboratory Topics (I and II) courses, some more complex programming was carried out, which served as an example for the development of programming in the Scilab environment. Also, some processing of the SWAT model software could be carried out in these courses, as well as in the Remote Sensing courses in the Agronomy course in Itaqui, and Geotechnologies Applied to the Preparation of Agricultural Reports and Expertise in the Aquaculture course in Uruguaiana, which are related to geoprocessing. 15 In Statistics, we worked on the climate model used in SWAT, the Weather Generator (WXGEN), which is analysed in an independent volume. The coursework contributed to building a more general overview of the SWAT model and the results of its processing, as well as helping to organise the developments produced during the courses. Also in this sense, some trainees from the Integrated Interdisciplinary Laboratory (LABii) contributed to the progress of certain tasks and even to participation in classes. 16 3. PRESENTATION OF RESULTS In the rest of this book, and in other volumes, we will be presenting examples of the developments that have been made, which have helped to systematise the understanding of the main physical models, their variables and concepts, including calculations that encourage a deeper understanding of each particular model. In this way, we are getting closer to the major objective of gaining an in-depth understanding of the large SWAT model, taking advantage of the curricular spaces available and suitable for this purpose. The physical-mathematical models studied, the programming algorithms suitable for the Scilab mathematical environment prepared, and the results obtained will be presented for each topic. 3.1 Precipitation Precipitation is the most important variable in the hydrological cycle, as it deposits water from the atmosphere on the earth's surface, triggering all the other phenomena related to the movement of water in the catchment area. It can occur meteorically: rain, hail, snow, or hidden: dew, fog and frost. For the purposes of this study we will be using rain (precipitation) within the interval from the start of precipitation to its interruption, the period that constitutes the so-called storm, distributed over time according to the characteristics of the regional climate. Precipitation is generally the only form of water entering a river basin, and quantifying it is necessary for dimensioning water supply, irrigation, hydraulic works, flood control, soil erosion, among others. 17 There are three primary causes of rain formation and they all have to do with the rise of warm, humid air masses in the atmosphere. This air mass rises to a lower temperature level where the dew point can be reached or exceeded. The different causes of precipitation formation and occurrence are: orographic precipitation, convective precipitation and frontal precipitation. The behaviour of precipitation over time, i.e. the shape of the precipitation intensity graph as a function of minutes or hours, will be related to the type of precipitation, so there are various types of behaviour of these graphs. 3.1.1 Precipitation as a Physical Variable - Units and Relationships The first step in analysing a physical model for a given system is a conceptual understanding of the variables, their units, and the relationships with other variables. The system considered in this book is the Water Response Unit (HRU), which is a region within a river basin that has uniform topography, soil type and land cover, i.e. of a single type, which makes it possible to carry out single calculations, as the equations apply to the entire region. The concept is similar to the agricultural plot used in Brazil, which is an area defined by its uniformity, making it possible to carry out agricultural operations without interruptions caused, for example, by obstacles or very different slopes. Table 1 - Relations and/or equalities of the variables involved Main Variab le Relationship and/or equality of the Units Involved Quantities Involved Units of measurement Derived Quantity Relationships 18 Precipit ation mm / hr Volume/Time (millimetre (mm) / hour (hr)) Length and Time m3 /s / km2 Volume/Area (cubic metre (m3 ) / second / square kilometre (Km2 )) Length, Time, Area 1m =3 1000 litres Volume = Capacity Length x Area = Volume 1dm =3 1 litre 1cm =3 1 millim etre 60 min = 1 hour 3600seconds=60minutes = 1hour Time 60 sec = 1 min 60 seconds = 1 minute Metric Conversions Table - VOLUME k m3 h m3 da m3 m 3 d m3 c m3 m m3 Correspondence to unit of CAPACITY 1 00 0 00 0 000 Equivalent to 1,000L (1m3 =10 0 1 00 0 000 Equals 1L (1dm3 =1L) 1 000 Equals 1ml (1cm3 =1ml) 19 Example of using the CHAIN RULE to convert units: 1mm/h = 1mm/hr x 1m/1000mm x 1hr/3600sec , transformation of mm/hr into m/s So we have the equivalence of rainfall intensity in mm/h to m/s, needing the area of the region affected to transform it into volume; for each square metre, 1mm/hr is equivalent to 1/3600/1000 m3 /s = 1/3600/1000 m3 /s x 1000 litres/m3 /s = 1/3600 litres/s = 2.78x10-4 litres/s x 1000 ml/1litre = 0.278 ml/s. We can see in the example above that a rainfall with an intensity of 1 mm/h is equivalent to receiving 0.278 ml of water per second over an area of 1 m2 , approximately 1 ml every 4 seconds, which is a very weak rainfall. This type of exercise is basically worked on in the Physics subject, where students have their first contact with variables, units and physical models, but already in an applied way, in Agronomy. 3.1.2 Modelling Precipitation Intensity in the SWAT Model This quantity, precipitation, will be represented mathematically in the SWAT model by means of relationships with other variables and parameters that make it possible to characterise the type of precipitation that occurs in the region under study. In other words, the mathematical relationships that represent precipitation in the model can be adapted to each climatic region by adjusting the parameters and eventually the functions used. The process of adaptation is called calibration. 20 3.1.2.1 Variables and Parameters Involved in the Representation of a Storm The following table describes the variables or parameters involved in the physical representation of precipitation during a single rainfall event, which we call a storm, according to the SWAT model. As the ultimate goal is to programme the equations that represent the precipitation phenomenon, there are several steps in describing the variables. On the one hand, we must use the notations common in scientific texts, using sub-indices and supra-indices, as well as Greek letters, to name variables and parameters. On the other hand, as in algorithms programmed in mathematical environments we only have common alphanumeric characters, written inline, we must find a corresponding name for the variables that can be used in the mathematical programming environment, as well as rewriting the equations using the symbols for addition (+), subtraction (-), multiplication (*), division (/) and potentiation (^). In our case, the mathematical environment used is Scilab, but there are several other equivalents available (R, Otave, Matlab). Table 2 - Variables used in the equations according to their common nomenclature 27 was calibrated using the rainfall recorded by the Itaqui automatic station, data every 15 minutes, shown in the previous figure. Table : Programming the equation in SCILAB2 //Exponential equations for rainfall page 921 (representation of a rain storm) T_sub=[0:1:10]/60;//time from start of rain to maximum, at 1min intervals, in hours T_des=[10:1:30]/60; //time from peak rainfall intensity to end, every 1 min, in hours T_peak=0.17;//time from start to peak intensity, in hours T_dur=0.5// total storm time, in hours delta1=0.33;// pp growth coefficient in hours (before peak) delta2=0.42//coefficient of decrease of pp in hours (after peak) i_mx=40 // maximum rainfall mm/h i_T_sub=i_mx*exp((T_sub-T_peak)/delta1); i_T_des=i_mx*exp((T_peak-T_des)/delta2); T=cat (2,T_sub,T_des) //UNIT THE UP AND DOWN VECTOR i_T=cat (2,i_T_sub,i_T_des) xtitle('RAINFALL TEMPLATE','TIME (hr)','PRECIPITATION (mm)') plot2d(T,i_T) //Graph of rainfall intensity as a function of time It is important to emphasise to current undergraduate students that the equation will generally be written in its traditional form in texts, as a 28 mathematical equation, with variables often represented by Greek letters, also using sub-indices and supra-indices for differentiation. However, this notation is not suitable for the mathematical programming and calculation environment, whereby the Greek letters are replaced by their name and the sub-indices or supra-indices are indicated by the underline character (_). The following table illustrates these changes. SYMBOLS AND LEGENDS USED IN THE EQUATION AND GRAPH CONSTRUCTION Correspo nding symbols EQUATION CONCEPTUALISATION IN SCILAB/SWAT UNI TS SW AT SCIL AB i i_T precipitation intensity over time mm/ hr i(T) imx i_mx maximum intensity (peak rainfall during the storm) mm/ hr T T_su b rainfall increase time to maximum time since the storm began hr T_de s time from the peak of rainfall intensity until the end hr Tpea k T_pe ak maximum value, (peak rainfall intensity) hr Tdur T_du r storm duration Hr 29 δ1* delta 1 equation coefficient for rainfall intensity BEFORE reaching peak intensity Hr δ2 * delta 2 equation coefficient for rainfall intensity AFTER reaching peak intensity Hr - i_T_s ub Calculation of the formula itself (increase in rainfall) mm/ hr - i_T_ des Calculation of the formula itself (decrease in rainfall) mm/ hr - T=cat () T=cat(2,i_T_sub, i_T_des) Unites the up and down vector of the graph mm/ hr - plot2 d(T,i _T) Graph generator command - 3.1.3.2 Normalised Intensity Distribution (page 93 )1 The rainfall intensity distribution given in the equation above (1:3.3.1) can be normalised to eliminate the units of both intensity and time. This is convenient for obtaining what we can call the pure form of the storm graph, which can easily be adapted to any duration or maximum intensity that occurs during a storm. 1.3.2.1 Concepts of normalised intensity To obtain the normalised intensity, î , each rainfall intensity value is divided by the average intensity of the storm, iave , although another approach would be to divide by the maximum value imx . 30 To obtain the normalised time,  ෡ , all the time values, T, from the start of the storm are divided by the duration of the storm, Tdur , This quotient will give the time during the storm expressed as a function of the total duration of the storm (0.0-1.0). The normalised distribution of storm intensity is: The relationship between the original equation coefficients and the normalised equation coefficients is: Where: 31 δ1 This is the equation coefficient for rainfall intensity before reaching peak intensity (hr). d1 This is the normalised equation coefficient for rainfall intensity before reaching peak intensity. δ2 This is the equation coefficient for rainfall intensity after the peak intensity has been reached (hr). d2 This is the normalised equation coefficient for rainfall intensity after the peak intensity has been reached. The equations, written in the mathematical language Scilab, are as follows: increase in rainfall: i_T_norm_sub=i_mx*exp((T_norm_subT_peak)/d1); decrease in rainfall: i_T_norm_des=i_mx*exp((T_peakT_norm_des)/d2); Sy m bo ls Variable Name CONCEPTUALISATION OF EQUATIONS 1:3.3.2 , 1:3.3.3 and 1.3.3.4 UN ITS iav e i_ave average rainfall intensity in a storm mm /hr 32 T T_sub and T_des time since the start of the storm. Separated into ascent time and descent time hr Td ur T_dur duration of the storm hr î i_norm normalised value of rainfall intensity over time dim ensi onle ss  ෡ t_norm time during the storm expressed as a function of the total duration of the storm. (0.0 - 1.0) 3.1.3.2.2 Scilab programming of the Normalised Intensity //Normalised exponential rainfall equations (representation of a rain storm) T_dur=0.5// total storm time, in hours T_norm_sub=[0:1:10]/60/T_dur// rainfall increase time to maximum, at 1min intervals in fraction of total storm time T_norm_des=[10:1:30]/60/T_dur; //time from the peak of rainfall intensity to the end at 1 min intervals, in fraction t.t.t. t_peak_norm=0.17/T_dur;//time from start to peak intensity, in hours i_ave=5 //average rainfall intensity in a storm, mm/hr delta_1=0.33// pp growth coefficient in hours (before peak) delta_2=0.42//coefficient of decrease of pp in hours (after peak) d1=t_peak_norm/4.605 // This is the normalised equation coefficient for rainfall intensity before reaching peak intensity. 33 d2=(1-t_peak_norm)/4.605 //Is the normalised equation coefficient for rainfall intensity after peak intensity is reached. i_mx=8.6 // maximum rainfall mm/h i_mx_norm=i_mx/i_ave i_T_norm_sub=i_mx_norm*exp((T_norm_sub-t_peak_norm)/d1); i_T_norm_des=i_mx_norm*exp((t_peak_norm-T_norm_des)/d2); T=cat (2,T_norm_sub,T_norm_des) //UNIT THE UP AND DOWN VECTOR i_T=cat (2,i_T_norm_sub,i_T_norm_des) xtitle('NORMALISED RAINFALL TEMPERATURE','WEATHER (FACTION OF TOTAL DURATION)','PRECIPITATION (FACTION OF MAXIMUM PP)') plot2d(T,i_T) //Graph of rainfall intensity as a function of time 34 3.1.3.2.3 Precipitation: time generated for peak intensity. Total precipitation and duration. The time it takes for the peak intensity to occur can be modelled and can be created as a function of a random component and parameters calibrated according to regional characteristics. The variables involved in this simulation are described below. TABLE OF VARIABLES Name / Symbol Units/value Concept Relationships / Equations t_peak fraction of total storm duration (f.d.t.t) Value: [0 - 1] time normalised to peak intensity t_peak=t_peakM*(t_peak L+(rnd1*(t_peakUt_peakL)*(t_peakMt_peakL))^0.5)/t_peak_m ean t_peakM f.d.t.t Value: [0 - 1] Average time to peak intensity. t_peakL f.d.t.t Value: [0 - 1] minimum time to peak intensity rnd1 no units random 35 [0-1] number between 0.0 and 1.0 t_peakU f.d.t.t Value: [0 - 1] maximum time to peak intensity t_peak_ mean f.d.t.t Value: [0 - 1] average of t_peakM and t_peakU t_peak=t_peakM*(t_peakL+(rnd1*(t_peakU-t_peakL)*(t_peakMt_peakL))^0.5)/t_peak_mean PROGRAMMING IN SCILAB //Precipitation: time generated for peak intensity t_peakM=0.5 //Average time for peak intensity in fraction of total storm duration (f.d.t.t). t_peakL=0.2 //minimum time for peak intensity in f.d.t.t rnd1=rand() //random number between 0.0 and 1.0 t_peakU=0.7 //maximum time for peak intensity in f.d.t.t t_peak_mean=(t_peakM+t_peakU)/2 //average of t_peakM and t_peakU in f.d.t.t //time normalised to peak intensity 36 t_peak=t_peakM*(t_peakL+(rnd1*(t_peakU-t_peakL)*(t_peakMt_peakL))^0.5)/t_peak_mean Example result: rnd1 = 0.7560439 t_peakU = 0.7; t_peak_mean = 0.6 ; t_peak = 0.4472991 in units of a fraction of the total time of the storm; this means that the peak would occur at a time equal to 44.7 per cent of the total duration of the rain. This time can be introduced into the simulation of a storm's rainfall intensity, and its calibration according to regional reality takes place by adjusting the value of the average, maximum and minimum times that are put into the algorithm, and constitute statistically based information that describes regional characteristics. 3.1.3.2.3 Calculation of Total Precipitation and its Duration As rainfall intensity is the volume of rain per unit of time and per unit of area, the total rainfall accumulated over time per unit of area can be calculated using the integral of the function i(t), i.e. intensity as a function of time. The integral is conceptually the area under the graph of i(t) as a function of t, i.e. the product of i(t)*t, which equals the sum of precipitation over time, which is the variable called total precipitation. The value of the integral, which is the total rainfall, can be found numerically (by adding squares or trapezoids) or analytically, when there is a mathematical representation of the rainfall intensity over time. 43 HYDROLOGY PRECIPITATION TEST 1) Describe the concept and units of the hydrological variable PRECIPITATION. 2) Describe the concept of watershed, water response unit (HRU) and agricultural plot, and explain, with an example, how to convert a precipitation value in mm to m3 /s for each of them. Based on your understanding of the text in Annex 1, explain the next questions in your own words: 3) How useful is it to be able to simulate rainfall using a computer method, and thus be able to define which days of a given month will hypothetically be rainy (Wet) or dry (Dry)? 4) Explain what the probabilities Pi (D/W) , Pi (W/W) , Pi (D/D) , Pi (W/D) are and how to calculate them, giving examples of how to calculate some of the monthly probabilities taken from the 2015 Uruguaiana spreadsheet. 5) Explain the fundamentals of the Markov procedure on how to make the computational definition of the occurrence of a wet or dry day using the previous probabilities. 44 6) Use a code (algorithm, programme) in the mathematical programming language Scilab to define the days of the months of 2015 as wet (rainy) or dry. Based on the results, explain how the programme works. See Annex 2. 7) For each day defined as rainy, explain how a rainfall value that is reasonable for our region can be generated. To do this, look at the Scilab algorithm in Appendix 3 and try to explain how it works. 8) Apply the algorithms in items 6 and 7 to generate simulated rainy days for a month in 2015, and compare them with the real values measured by the weather station. 9) Write a summarised text explaining the work carried out and your conclusions. 45 ANNEX 1 - THEORETICAL BASIS For precipitation that is not read from meteorological data, or to fill gaps in measured data, precipitation can be created using the model developed by Nicks (1974). The model has two aspects: i) defining which days will be rainy; ii) defining the volume of rainfall. First-order Markov chains are used to define rainy days. It's easy to see that the occurrence of a wet, rainy day is related to what happened the day before. If the previous day was rainy, here called wet, it is natural to increase the chance that the current day will also be wet, without it being so important what happened in the previous history. This type of behaviour is described statistically by so-called Markov chains. Once you have the probabilities of rainfall events, which can be obtained by analysing meteorological records, you can program the criteria for defining a given day as wet (rainy) or dry (without precipitation). Once we have the Markov probabilities, in the first-order Markov model, the probability of rain on a certain day is conditional on the wet or dry status of the previous day. A wet day is defined when we have 0.1 mm of rain or more, although it is necessary to check that the sensitivity of the sensors of the automatic weather stations used does not introduce noisy rainfall data, often values of up to 0.2 mm, which in this case start to appear with high frequency. From the meteorological data for a given month, the model user needs to calculate the probability Pi (W/W) of having a wet day on day i given a wet day on the previous day i-1; and also the probability Pi (W/D) of a wet day on day i given a dry day on day i-1. This probability is 46 calculated for each month of the year of interest, which can be done by counting the days of rain after rain and rain after dry day directly in the meteorological archives. From this input data, the probability of transition from one type of day to another can be derived: Pi (D/W)=1-Pi (W/W) , probability of a dry day on day i given a wet day on day i-1. Pi (D/D)=1-Pi (W/D) , probability of a dry day on day i given a dry day on day i-1. For the Scilab programming environment we will use the notation: P_i_D_W=1-P_i_W_W // probability of a dry day on day i given a wet day on day i-1. P_i_D_D=1-P_i_W_D // probability of a dry day on day i given a dry day on day i-1 In order to define a certain day in the simulation as dry or wet, the swat model generates a random number between 0.0 and 1.0, which will represent the random probability (of the day changing in relation to the previous day) that any given day has (in other words, that day being considered in the simulation would turn out to be the opposite of the previous day). This random number is compared with the correct (Markov) probabilities according to what happened on the previous day. If the previous day was wet, then it is compared with wet-dry and wet-wet, and one of the following three possibilities can occur: ΁ΐΒΝΖΒΥ؄΁ΐΈΐ͵΁ΐΈΐΈ΅ΙΖΟΖΩΥΕΒΪΨΚΝΝΓΖΨΖΥ ΁ΐΈΐ͵؄΁ΐΒΝΖΒΥ΁ΐΈΐΈ΅ΙΖΟΖΩΥΕΒΪΨΚΝΝΓΖΨΖΥ ΁ΐΈΐ͵؄΁ΐΈΐΈ؄΁ΐΒΝΖΒΥ΅ΙΖΟΖΩΥΕΒΪΨΚΝΝΓΖΕΣΪ 47 If the previous day was dry, then it is compared with wet-dry and wet-wet, and one of the following three possibilities can occur: ΁ΐΒΝΖΒΥ؄΁ΐ͵ΐΈ΁ΐ͵ΐ͵΅ΙΖΟΖΩΥΕΒΪΨΚΝΝΓΖΕΣΪ ΁ΐΈΐ͵؄΁ΐΒΝΖΒΥ΁ΐΈΐΈ΅ΙΖΟΖΩΥΕΒΪΨΚΝΝΓΖΕΣΪ ΁ΐ͵ΐΈ؄΁ΐ͵؄΁ΐΒΝΖΒΥ΅ΙΖΟΖΩΥΕΒΪΨΚΝΝΓΖΨΖΥ To obtain the Markov probabilities for rainfall, you can use a file of measurements from meteorological stations; to see the model used for this work click here. The original file, from the INMET weather station in Uruguaiana, contains only what is in the TIME_DATA tab for various climatological variables, and the other tabs were created by processing this data. The hourly data was converted to daily data in the DAILY_DATA tab. The rainfall was placed on its own in the Pp_dia tab, showing only the rainy days and the total rainfall in mm of rain. This tab is used to count the rainy days after a rainy day and the rainy days after a dry day, which are needed to calculate the Markov probabilities. These probabilities are calculated on the PARAM_SWAT+day tab, in the PR_WD and PR_WW columns, and then copied to the Prob_Markov tab, where the missing probabilities (dry-wet and dry-dry) are completed. Note in the spreadsheets that the calculations correspond to the following ratios: PR_WW = WET DAYS AFTER WET DAY/TOTAL WET DAYS PR_WD = DRY DAYS AFTER WET DAYS/TOTAL DRY DAYS PR_DW = 1-PR_WW 48 PR_DD = 1-PR_WD Using the above formulas, the probability values for each month in the Uruguaiana records for 2015 can be found in the Prob_Markov tab and are reproduced here: Markov probability table (W-wet day; D-dry day) PR_WD PR_WW P_D_W P_D_D January 0.31 0.67 0.33 0.69 February 0.20 0.50 0.50 0.80 March 0.18 0.56 0.44 0.82 April 0.07 0.33 0.67 0.93 May 0.29 0.40 0.60 0.71 June 0.33 0.50 0.50 0.67 July 0.41 0.50 0.50 0.59 August 0.18 0.79 0.21 0.82 September 0.17 0.33 0.67 0.83 October 0.18 0.56 0.44 0.82 November 0.30 0.40 0.60 0.70 December 0.41 0.50 0.50 0.59 This table is the basis for building an algorithm in Scilab that can, using a random number, decide whether a day in the simulation will be considered wet or dry, as shown in the annex. 49 ANNEX 2 - SCILAB ALGORITHM FOR DEFINING RAINY DAYS //Definition of wet days by Markov chains - Uruguaiana 2015 data //Probability of dry day after wet day for each month PR_WD=[0.31 0.20 0.18 0.07 0.29 0.33 0.41 0.18 0.17 0.18 0.30 0.41]; //Probability of wet day after wet day for each month PR_WW=[0.67 0.50 0.56 0.33 0.40 0.50 0.50 0.79 0.33 0.56 0.40 0.50]; //Probability of wet day after dry day for each month PR_DW=[0.69 0.80 0.82 0.93 0.71 0.67 0.59 0.82 0.83 0.82 0.70 0.59]; //Probability of dry day after dry day for each month PR_DD=[0.33 0.50 0.44 0.67 0.60 0.50 0.50 0.21 0.67 0.44 0.60 0.50]; //Random probability for 31-day, 30-day and 28-day months P_Aleat31days=rand(1,31); P_Aleat30days=rand(1,30); P_Aleat28days=rand(1,28); for i=1:12 DecWD31=P_Aleat31days-PR_WD(i); DecWW31=P_Aleat31days-PR_WW(i); DecDW31=P_Aleat31days-PR_DW(i); DecDD31=P_Aleat31days-PR_DD(i); for j=1:31 if DecWD31(j) <0 | DecWW31(j)<0 then //Compares Probab. previous rainy day type_dayW(i,j)=1; //1 represents rainy day else type_dayW(i,j)=0; //0 represents dry day 50 end if DecDW31(j) <0 | DecDD31(j)<0 then //Compare Probab. previous day dry type_dayD(i,j)=0; else type_dayD(i,j)=1; end end end //Monthly results based on Prob previous rainy day //To see the results for a month, delete the final semicolon pp_Uru_jan_2015W=tipo_diaW(1,:); pp_Uru_fev_2015W=tipo_diaW(2,:); pp_Uru_mar_2015W=tipo_diaW(3,:); pp_Uru_abr_2015W=tipo_diaW(4,:); pp_Uru_mai_2015W=tipo_diaW(5,:); pp_Uru_jun_2015W=tipo_diaW(6,:); pp_Uru_jul_2015W=tipo_diaW(7,:); pp_Uru_ago_2015W=tipo_diaW(8,:); pp_Uru_set_2015W=tipo_diaW(9,:); pp_Uru_out_2015W=tipo_diaW(10,:); pp_Uru_nov_2015W=tipo_diaW(11,:); pp_Uru_dez_2015W=tipo_diaW(12,:); //Monthly results based on Prob previous dry day //To see the results for a month, delete the final semicolon pp_Uru_jan_2015D=tipo_diaD(1,:); pp_Uru_fev_2015D=tipo_diaD(2,:); 51 pp_Uru_mar_2015D=tipo_diaD(3,:); pp_Uru_abr_2015D=tipo_diaD(4,:); pp_Uru_mai_2015D=tipo_diaD(5,:); pp_Uru_jun_2015D=tipo_diaD(6,:); pp_Uru_jul_2015D=tipo_diaD(7,:); pp_Uru_ago_2015D=tipo_diaD(8,:); pp_Uru_set_2015D=tipo_diaD(9,:); pp_Uru_out_2015D=tipo_diaD(10,:); pp_Uru_nov_2015D=tipo_diaD(11,:); pp_Uru_dez_2015D=tipo_diaD(12,:); //Select day W or D depending on the previous day //Final result of rainy (1) or dry (0) days for each month //January pp_Uru_jan_2015(1)=pp_Uru_jan_2015D(1); for k=2:31 if pp_Uru_jan_2015(k-1)==1 then pp_Uru_jan_2015(k)=pp_Uru_jan_2015W(k); else pp_Uru_jan_2015(k)=pp_Uru_jan_2015D(k); end end pp_Uru_jan_2015; //February pp_Uru_fev_2015(1)=pp_Uru_fev_2015D(1); for k=2:31 52 if pp_Uru_fev_2015(k-1)==1 then pp_Uru_fev_2015(k)=pp_Uru_fev_2015W(k); else pp_Uru_fev_2015(k)=pp_Uru_fev_2015D(k); end end pp_Uru_feb_2015; //March pp_Uru_mar_2015(1)=pp_Uru_mar_2015D(1); for k=2:31 if pp_Uru_mar_2015(k-1)==1 then pp_Uru_mar_2015(k)=pp_Uru_mar_2015W(k); else pp_Uru_mar_2015(k)=pp_Uru_mar_2015D(k); end end pp_Uru_mar_2015; //April pp_Uru_abr_2015(1)=pp_Uru_abr_2015D(1); for k=2:31 if pp_Uru_abr_2015(k-1)==1 then pp_Uru_abr_2015(k)=pp_Uru_abr_2015W(k); else pp_Uru_abr_2015(k)=pp_Uru_abr_2015D(k); end end pp_Uru_abr_2015; 59 sigma_mon=14.28 //standard deviation of daily rainfall, rainy days only, for the month g_mon= 1.87 //asymmetric coefficient for daily precipitation in the month SND_day=cos(6.283*rand())*sqrt(-2*log(rand())) //normal random standard deviation calculated for the day //Calculation of daily precipitation R_day=mu_mon+2*sigma_mon*((((SND_dayg_mon/6)*(g_mon/6)+1)^3-1)/g_mon) The normal standard deviation for the day is calculated using the following equation: where rnd1 and rnd2 are random numbers between 0.0 and 0.1, which in Scilab notation is: SND_day=cos(6.283*rand())*sqrt(-2*log(rand())) //SCILAB PROGRAMME PRECIPITATION MODEL ASYMMETRIC DISTRIBUTION NICKS pp_mes=210.4 //total precipitation in the month days_pp=15 // number of days with precipitation in the month mu_mon= pp_mes/dias_pp //average daily rainfall (mmH20) sigma_mon=14.28 //standard deviation of daily rainfall, only rainy days, for the month 60 SND_day=cos(6.283*rand(1,31)).*sqrt(-2*log(rand(1,31))) //normal random standard deviation calculated for the day g_mon= 1.87 //asymmetric coefficient for daily precipitation in the month R_day=mu_mon+2*sigma_mon*((((SND_dayg_mon/6)*(g_mon/6)+1)^3-1)/g_mon) days=[1:1:31] xtitle('SIMULATION PRECIPITATION ASYMMETRIC DISTRIBUTION NICKS','TIME (days)','PRECIPITATION (mm)') plot2d(days,R_day) Examples of RESULTS: 61 The exponential distribution is offered as another alternative to the asymmetric distribution, as it requires fewer inputs and is often used in areas that do not have data on precipitation events. daily precipitation in the exponential distribution is calculated using the following equation: In this equation Rday consists of the amount of rainfall in one day (mmH20) , μmon is the average daily rainfall for the month, rnd1 is a random number between 0.0 and 1.0 and rexp is an exponent that must be appropriate between 1.0 and 2.0. As the value of rexp is increased, the number of extreme precipitation events during the year will increase. testing this equation in the USA has shown that a value of 1.3 generates satisfactory results. //SCILAB PROGRAMME PRECIPITATION MODEL EXPONENTIAL DISTRIBUTION mu_mon= 14.03 // average daily rainfall (mmH20) rexp= 1.3 //exponent to be appropriate in the range [1.0-2.0] rnd1=rand(1,31) //random numbers for 31 days //Calculation of daily precipitation R_day=mu_mon*(-log(rnd1))^rexp days=[1:1:31] xtitle('EXPONENTIAL DISTRIBUTION PRECIPITATION SIMULATION','TIME (days)','PRECIPITATION (mm)') plot2d(days,R_day) 62 RESULTS: We can conclude that precipitation can be simulated coherently by both models, and that the algorithms created in Scilab work correctly, thus allowing the availability of precipitation data that responds to regional characteristics. Using these algorithms also prepares the user for working with the SWAT model, as it provides detailed contact with the variables and equations used. To this end, we have listed the variables used in the SWAT model below. Input variables in the swat model that belong to precipitation generation. PCPSIM Precipitation input code 1=measured 2=generated PR_W(1,mon) Pi (W/D): probability of a wet day after a dry day in the month 63 PR_W(2,mon) Pi (W/W): probability of a wet day after a wet day in the month IDIST Rainfall distribution code 0=asymmetric 1=exponential REXPrexp: exponent value (mandatory if IDIST=1) PCPMM(mon) average amount of precipitation falling in the month (mmH20) PCPD(mon) number of days of precipitation in the month (μmon =PCPMM/PCPD) PCPSTD(mon) σmon : Standard deviation for daily precipitation in the month (mmH2O) PCPSKW(mon) gmon : slope coefficient for daily precipitation in the month . 3.1.8 Complete rain generation algorithm At the end of the algorithm to define the rainy days by Markov probability, the following code is appended, which in addition to calculating the magnitude of the rainy days, eliminates the values of the days that do not have rain predicted by the model. 64 //SCILAB PROGRAMME PRECIPITATION MODEL ASYMMETRIC DISTRIBUTION NICKS pp_mes=210.4 //total precipitation in the month days_pp=15 // number of days with precipitation in the month mu_mon= pp_mes/dias_pp //average daily precipitation (mmH20) sigma_mon=14.28 //standard deviation of daily rainfall, only rainy days, for the month SND_day=cos(6.283*rand(1,31)).*sqrt(-2*log(rand(1,31))) //normal random standard deviation calculated for the day g_mon= 1.87 //asymmetric coefficient for daily precipitation in the month R_day_Uru_jan_2015=mu_mon+2*sigma_mon*((((SND_dayg_mon/6)*(g_mon/6)+1).^3-1)/g_mon); for k=1:31 //For each day of the month if pp_Uru_jan_2015(k) ==1 then pp_Uru_jan_2015(k)=R_day_Uru_jan_2015(k); else pp_Uru_jan_2015(k)=0; end end pp_Uru_jan_2015 days=[1:1:31] xtitle('SIMULATION PRECIPITATION ASYMMETRIC DISTRIBUTION NICKS','TIME (days)','PRECIPITATION (mm)') plot2d(days,pp_Uru_jan_2015) //Definition of wet days by Markov chains - Uruguaiana 2015 data //Probability of dry day after wet day for each month 65 PR_WD=[0.31 0.20 0.18 0.07 0.29 0.33 0.41 0.18 0.17 0.18 0.30 0.41]; //Probability of wet day after wet day for each month PR_WW=[0.67 0.50 0.56 0.33 0.40 0.50 0.50 0.79 0.33 0.56 0.40 0.50]; //Probability of wet day after dry day for each month PR_DW=[0.69 0.80 0.82 0.93 0.71 0.67 0.59 0.82 0.83 0.82 0.70 0.59]; //Probability of dry day after dry day for each month PR_DD=[0.33 0.50 0.44 0.67 0.60 0.50 0.50 0.21 0.67 0.44 0.60 0.50]; //Random probability for 31-day, 30-day and 28-day months P_Aleat31days=rand(1,31); P_Aleat30days=rand(1,30); P_Aleat28days=rand(1,28); for i=1:12 DecWD31=P_Aleat31days-PR_WD(i); DecWW31=P_Aleat31days-PR_WW(i); DecDW31=P_Aleat31days-PR_DW(i); DecDD31=P_Aleat31days-PR_DD(i); for j=1:31 if DecWD31(j) <0 | DecWW31(j)<0 then //Compares Probab. previous rainy day type_dayW(i,j)=1; //1 represents rainy day else type_dayW(i,j)=0; //0 represents dry day end if DecDW31(j) <0 | DecDD31(j)<0 then //Compare Probab. previous day dry type_dayD(i,j)=0; 66 else type_dayD(i,j)=1; end end end //Monthly results based on Prob previous rainy day //To see the results for a month, delete the final semicolon pp_Uru_jan_2015W=tipo_diaW(1,:); pp_Uru_fev_2015W=tipo_diaW(2,:); pp_Uru_mar_2015W=tipo_diaW(3,:); pp_Uru_abr_2015W=tipo_diaW(4,:); pp_Uru_mai_2015W=tipo_diaW(5,:); pp_Uru_jun_2015W=tipo_diaW(6,:); pp_Uru_jul_2015W=tipo_diaW(7,:); pp_Uru_ago_2015W=tipo_diaW(8,:); pp_Uru_set_2015W=tipo_diaW(9,:); pp_Uru_out_2015W=tipo_diaW(10,:); pp_Uru_nov_2015W=tipo_diaW(11,:); pp_Uru_dez_2015W=tipo_diaW(12,:); //Monthly results based on Prob previous dry day //To see the results for a month, delete the final semicolon pp_Uru_jan_2015D=tipo_diaD(1,:); pp_Uru_fev_2015D=tipo_diaD(2,:); pp_Uru_mar_2015D=tipo_diaD(3,:); pp_Uru_abr_2015D=tipo_diaD(4,:); pp_Uru_mai_2015D=tipo_diaD(5,:); pp_Uru_jun_2015D=tipo_diaD(6,:); 67 pp_Uru_jul_2015D=tipo_diaD(7,:); pp_Uru_ago_2015D=tipo_diaD(8,:); pp_Uru_set_2015D=tipo_diaD(9,:); pp_Uru_out_2015D=tipo_diaD(10,:); pp_Uru_nov_2015D=tipo_diaD(11,:); pp_Uru_dez_2015D=tipo_diaD(12,:); //Select day W or D depending on the previous day //Final result of rainy (1) or dry (0) days for each month //January pp_Uru_jan_2015(1)=pp_Uru_jan_2015D(1); for k=2:31 if pp_Uru_jan_2015(k-1)==1 then pp_Uru_jan_2015(k)=pp_Uru_jan_2015W(k); else pp_Uru_jan_2015(k)=pp_Uru_jan_2015D(k); end end pp_Uru_jan_2015; //February pp_Uru_fev_2015(1)=pp_Uru_fev_2015D(1); for k=2:31 if pp_Uru_fev_2015(k-1)==1 then pp_Uru_fev_2015(k)=pp_Uru_fev_2015W(k); else pp_Uru_fev_2015(k)=pp_Uru_fev_2015D(k); end end 68 pp_Uru_feb_2015; //March pp_Uru_mar_2015(1)=pp_Uru_mar_2015D(1); for k=2:31 if pp_Uru_mar_2015(k-1)==1 then pp_Uru_mar_2015(k)=pp_Uru_mar_2015W(k); else pp_Uru_mar_2015(k)=pp_Uru_mar_2015D(k); end end pp_Uru_mar_2015; //April pp_Uru_abr_2015(1)=pp_Uru_abr_2015D(1); for k=2:31 if pp_Uru_abr_2015(k-1)==1 then pp_Uru_abr_2015(k)=pp_Uru_abr_2015W(k); else pp_Uru_abr_2015(k)=pp_Uru_abr_2015D(k); end end pp_Uru_abr_2015; //May pp_Uru_mai_2015(1)=pp_Uru_mai_2015D(1); for k=2:31 if pp_Uru_mai_2015(k-1)==1 then pp_Uru_mai_2015(k)=pp_Uru_mai_2015W(k); else 75 here R0.55sm(mon) is the smoothed maximum half-hourly rainfall for a given month (mm H20), R0.5x(mon) is the extreme maximum half-hourly rainfall for the specified month (mm H20), R0.5x(mon-1) and R0.5x(mon+1) are the extreme maximum half-hourly rainfall for the month before and after, respectively. This is necessary to smooth out the great variability that can occur in monthly rainfall values. Once the maximum attenuated half-hourly precipitation is known, the representative fraction of the half-hourly precipitation is calculated using the equation In this fraction α0.5mon is the average fraction of half-hourly precipitation for the month, adj0,5α consists of an adjustment factor, R0.5sm(mon) is the amount of half-hourly precipitation attenuated for the month (mm H20) , μmon is the average daily precipitation for the month, yrs is the number of years of precipitation data used to obtain the extreme values of monthly half-hourly precipitation and dayswet is the total number of wet days in the month. The adjustment factor is included to allow users to modify estimates of fractions of half-hourly precipitation and peak flow for the runoff rate. Adapting the nomenclature to Scilab syntax, we have: adj_0p5_alpha= 1//adjustment factor 76 R_0p5_mon_bef=41.60 //mm extreme maximum half-hourly rainfall for the previous month R_0p5_mon=24.4 //mm extreme maximum half-hourly rainfall for the month in question R_0p5_aft=30.4 //mm maximum extreme half-hourly rainfall for the following month R_0p5sm_mon=(R_0p5_mon_bef+R_0p5_mon+R_0p5_aft)/3 //mm amount of half-hourly precipitation attenuated for the month u_mon= //average daily precipitation for the month yrs=1 //number of years of precipitation data used to obtain the extreme values of half-hourly monthly precipitation days_wet=15 //total number of wet days in the month mu_mon= 14.03 alpha_0p5mon=adj_0p5_alpha*(1exp(R_0p5sm_mon/(mu_mon*log(0.5/(yrs*days_wet))))) RESULTS: adj_0p5_alpha = 1. ; R_0p5_mon_bef = 41.6; R_0p5_mon = 24.4; R_0p5_aft =30.4; R_0p5sm_mon =32.133333; yrs =1; days_wet =15; mu_mon =14.03; alpha_0p5mon =0.4900229 Conclusion: The algorithms prepared in Scilab provided an average halfhourly rainfall fraction of approximately 50% , meaning that the maximum half-hourly rainfall would be half the daily rainfall value. This estimate seems reasonable, and can be adjusted using the adjustment factor, which was taken to be = 1 in this simulation. As 1-hour rainfall 77 data was used, the values may be overestimated, and it is possible to lower the adjustment factor to 0.5-0.7, for example, which may be more realistic. The value of the monthly average fraction of maximum rainfall calculated above can be used to calculate the maximum half-hourly rainfall on every day of the month, or a fraction can be generated for each day, as explained in the next section. 3.1.11 Daily value of the maximum half-hourly rainfall fraction The user of the model analysed here has the option of using the maximum monthly or daily half-hourly precipitation. In the case of the SWAT model, the ised-det variable in the basin input file defines which option the user prefers. In the case of choosing to generate daily fractions, the triangular distribution can be used for this purpose, but it should be noted that the randomness of the triangular distribution used to generate daily values can cause an exaggerated oscillation in the value of the maximum half-hourly rainfall. Especially for small areas or micro-basins, the variability of the triangular distribution can be unrealistic. The distribution used to generate the daily maximum half-hourly rainfall fraction needs four inputs the monthly average half-hourly rainfall fraction, maximum value for the half-hourly rainfall fraction allowed in the month, minimum value for the half-hourly rainfall fraction allowed in the month and a random number between 0.0 and 1.0 78 The maximum fraction of half-hourly precipitation, or upper limit of the triangular distribution, can be calculated from the daily amount of precipitation with the equation where α0,5U is the largest half-hour fraction that can be generated on a given day, and Rday is the precipitation on a given day (mmh20) . In Scilab notation, we have the equivalent equation: alfa_0p5U=1-exp(-125/(R_day+5)) The minimum half-hour fraction or upper limit of the triangular distribution, α0.5L is set at 0.20083. Scilab notation: alpha_0p5L=0.02083 Once the maximum and minimum fraction values have been established, the triangular distribution uses one of two sets of equations to generate the maximum half-hourly rainfall fraction for the day. The equation to be used is defined by comparing a randomly generated value (0-1) and a comparison parameter: if where the right limb can be written: Cond=(alfa_0p5mon-alfa_0p5L)/(alfa_0p5U-alfa_0p5L) 79 then if Cond=(alpha_0p5mon-alpha_0p5L)/(alpha_0p5U-alpha_0p5L) we have the condition: rnd1=rand() if rnd1<=Cond then alfa_0p5=alfa_0p5mon*(alfa_0p5L+(rnd1*(alfa_0p5Ualfa_0p5L)*(alfa_0p5mon-alfa_0p5L))^0.5)/alfa_0p5mean else alfa_0p5=(alfa_0p5mon*(alfa_0p5U-(alfa_0p5Ualfa_0p5mon)*(alfa_0p5U*(1-rnd1)-alfa_0p5L*(1-rnd1))/(alfa_0p5Ualfa_0p5mon))^0.5)/alfa_0p5mean if then where α0.5 is the maximum fraction of half-hourly precipitation per day, α0.5mon is the average of the maximum fraction of half-hourly precipitation for the month, rnd1 is a random number generated each day, α0,5L is the minimum fraction of half-hourly precipitation that can be generated, α0,5U 80 is the largest fraction of half-hourly precipitation that can be generated and α5.0mean is the average of α0.5moL α0.5mon, and α0.5U . alfa_0p5mon=0.5 //average monthly fraction for pp max ½ hour R_day=30 // mm H2O precipitation on a certain day alfa_0p5U=1-exp(-125/(R_day+5)) alpha_0p5L=0.02083 alfa_0p5mean=(alfa_0p5U+alfa_0p5L+alfa_0p5mon)/3 Cond=(alfa_0p5mon-alfa_0p5L)/(alfa_0p5U-alfa_0p5L) rnd1=rand() if rnd1<=Cond then alfa_0p5=alfa_0p5mon*(alfa_0p5L+(rnd1*(alfa_0p5Ualfa_0p5L)*(alfa_0p5mon-alfa_0p5L))^0.5)/alfa_0p5mean else alfa_0p5=(alfa_0p5mon*(alfa_0p5U-(alfa_0p5Ualfa_0p5mon)*(alfa_0p5U*(1-rnd1)-alfa_0p5L*(1-rnd1))/(alfa_0p5Ualfa_0p5mon))^0.5)/alfa_0p5mean end RESULTS: alpha_0p5mon = 0.5 ; R_day = 30. ; alpha_0p5U = 0.9718843 ; alpha_0p5L = 0.02083 ; alpha_0p5mean = 0.4975714 ; Cond =0.5038303 ; rnd1 = 0.2113249 ; alpha_0p5 = 0.3327756 Conclusion: The maximum rainfall fraction can be generated for each day in a way that is consistent with the algorithms prepared, so it is possible to calculate the maximum half-hourly rainfall for each day of rain generated by the model, multiplying it by the total daily rainfall. The 81 fraction was calibrated using the regional meteorological parameters represented by the precipitation variables (average, maximum, monthly) starting with R. For application purposes in the SWAT model, once the programming has been done in Scilab, it is important to note the nomenclature of the variables used in SWAT. Table of input variables in the swat model attributed to the generation of maximum half-hourly precipitation ISED_DET Code governing the calculation of the maximum daily halfhourly rainfall: 0generate daily value 1use maximum monthly halfhourly rainfall value RAINHHMX(mon) R0,5x : extreme half-hourly rainfall for the month (mmh20) ADJ_PKR adj0,5a : peak rate adjustment factor PCPMM(mon) average amount of precipitation falling in the month PCPD(mon) dayswet : average number of rainy days in the month ( hmon = PCPMM/PCPD) RAIN_YRs yrs: number of years of data used to obtain values for RAINHHMX PRECIPITATION Rday amount of rain falling on a given day (mmh20) 82 3.1.12 Conclusion on Precipitation Generation Models As precipitation is the dominant variable in practice for water input into agricultural and environmental systems, its representation is of fundamental importance for modelling aimed at managing these systems. Once this input has been produced, it is possible to use simplified models to obtain the other water balance variables, or you can proceed to apply physically-based mathematical models that represent them in an analogous way. Wind speed may also be required depending on the method used to calculate evapotranspiration, e.g. Penman-Monteith. Once daily precipitation values are available, it is of interest to generate daily values for solar radiation, maximum and minimum temperature, in order to carry out evapotranspiration calculations and any processes to be modelled within the geographical system under study that depend on these variables. These processes will be studied in other sections. 3.2 Surface Flow Modelling 83 In order to analyse the agricultural and environmental management of water, and its interactions with the soil and chemical products, we need to know the flows at the different points in the river basins or systems under study, which are strategic for defining the effective quantity of products in circulation or displacement or accumulation within the system. The most appropriate systems for calculation are those that allow the use of a specific and determined set of parameters for adjusting the equations, which happens in spatially uniform areas. These uniform calculation systems are represented by HRUs, water response units, a concept similar to what in agricultural management practice is called a plot. The modelling of runoff within these areas must take place at two levels: daily and subdaily. At the daily level, we will use the CN number method. At the subdaily level, calculations linked to the so-called rational method will be used. 3.2.1 Introduction - Basic Concepts and Relationships Surface runoff corresponds to the segment of the hydrological cycle relating to the movement of water over the surface of the ground. It is of fundamental importance for the design of engineering works that are dimensioned to withstand the maximum flows resulting from surface runoff. Evaluating the hydrological cycle, it is expected that part of the total volume of rainfall will be intercepted by vegetation, while the rest will reach the soil surface, moistening the soil aggregates and reducing their cohesive forces. As the rain continues to fall, the aggregates disintegrate into smaller particles. The amount of destructured soil increases with the 84 intensity of the rainfall and the speed and size of the drops. As well as releasing particles, which clog the soil pores, the impact of the drops also tends to compact the soil, sealing its surface and consequently reducing the capacity for water infiltration. Water puddling in depressions on the soil surface only begins to occur when the intensity of precipitation exceeds the rate of infiltration, or when the soil's capacity to accumulate water is exceeded. Once the surface retention capacity has been exhausted, the water will begin to run off. Associated with surface runoff is the transport of soil particles, which are only deposited when the speed of surface runoff is reduced. As well as suspended soil particles, surface runoff carries chemical nutrients, organic matter, seeds and pesticides which, as well as causing direct damage to agricultural production, also cause pollution of watercourses. 3.2.1 Physical process associated with surface runoff Estimates of maximum runoff flows are often necessary in both agricultural and urban catchments. The first step in determining the design discharge is to calculate the fraction of precipitation that becomes surface runoff. Applying empirical methods to predict the runoff resulting from precipitation can be considered a first approximation that should be corrected later, based on an assessment of the system in operation. In basins without instrumentation, determining surface runoff is more difficult and less accurate than in instrumented basins. However, the implementation of models that calculate surface runoff can be a great initial stimulus to provoke the instrumentation of measurements that can progressively improve flow forecasts and calculations. 91 exceeds the demand, Q increases. On a local scale, in the case of a crop, the purpose of the water balance is to establish the variation in storage and, consequently, the availability of water in the soil. By knowing the soil's humidity or how much water it stores, it is possible to determine whether the crop is suffering from water deficiency, which is closely linked to the crop's yield levels. Water balance components for natural conditions Considering a control volume of soil, the BH has the following components, described below. ● Inputs - Rain - Dew - Surface Flow - Subsurface Flow - Hair Rising ● Outputs - Evapotranspiration - Surface run-off - Sub-surface flow - Deep drainage Water Balance Equation The water balance can be carried out for a layer of soil, a stretch of river or a river basin. Understanding these components depends on various factors such as: precipitation, potential evapotranspiration, soil conditions and land use, underground geology. The hydrographic basin is the best space for evaluating water behaviour, as it has defined the input space, the basin, and the output location, the 92 river section that defines the hydrographic basin. In a basin, the water balance is determined by: S (t + 1) = S(t) + (P - E - Q). dt Where, - S (t+1) and S(t) - amount of water at time t+1 and t - P - precipitation in the basin area in the interval - E - actual evapotranspiration in the time interval in the basin - Q - outflow in the time interval dt. - dt - differential time, the time interval between calculations, ideally very small (infinitesimal), to produce great accuracy and precision in the calculations. When the evaluation period (dt) is very long, the variations tend to offset each other, and the difference in storage (S) can be considered negligible, especially for large basins, so the equation can be simplified as follows: ΁Ͷ΂ Where, - P - precipitation in the basin area in the interval - E - actual evapotranspiration in the time interval in the basin - Q - surface and underground runoff REFERENCES GARCEZ, Lucas Nogueira; ALVAREZ, Guilhermo Acosta. Hydrology / 2nd edition. São Paulo : Edgard Blucher, 1988. PINTO, Nelson L. de Sousa; HOLTZ, Antonio Carlos Tatit; Martins, José Augusto; GOMIDE, Francisco Luiz Sibut. Basic Hydrology / 1st edition. Rio de Janeiro : Edgard Blucher, 1976. DE CARVALHO, Daniel; DA SILVA, Leonardo (2006). Hydrology [PDF File]. Available at: 93 http://www.ufrrj.br/institutos/it/deng/leonardo/downloads/APOSTILA/H IDRO-Cap7-E S.pdf 3.3 Flow Modelling According to the SWAT Model Production of Daily Flow Values The model used here to model runoff in uniform areas called HRUs, water response units, equivalent to agricultural plots, is the socalled SCS curve number procedure (SCS, 1972), also known as the CN number method, developed by the US Soil Conservation Service (SCS). For larger, non-uniform areas, an average number weighted by the subarea corresponding to each uniform region can be used. As this application was developed in an academic environment, within the scope of Physics, Hydrology and Interdisciplinary Laboratory Topics classes at the Itaqui Campus of the Federal University of Pampa, the development of the theme will be based on activities carried out in this context. 3.3.1 Calculating the flow of Arroio Olaria and Arroio Cambaí using the curve number method - SCS Taking advantage of the students' already conditioned energy to react in test situations, the following exercise was put into the form of a test in the Physics and Hydrology subjects: 94 EXERCISE 1 : a) USING THE DATA IN THE TABLE, CALCULATE THE VOLUME OF WATER PRECIPITATED IN THE CAMBAÍ AND OLARIA ARROY BAYS, IN ITAQUI, RS, Brazil, AFTER A 40 mm RAINFALL THAT DURED 50 MINUTES, FOLLOWED BY ANOTHER 30 mm RAINFALL DURING 15 MINUTES. EXPLAIN THE CONCEPTS INVOLVED AND THE METHODOLOGIES FOR CONVERTING UNITS TO LITRES AND TO THE INTERNATIONAL SYSTEM. b) INTERPRET THE PARAMETERS IN THE TABLE AND THEIR UNITS. Table I - Physical characteristics of the Olaria and Cambaí stream catchments. (BARIANI, 2016) Physical characteristics Olaria watershed Cambaí watershed Drainage area (km )2 14,368 157,519 Urban area (km )2 7,065 2,677 Perimeter (km) 17,961 53,235 Coefficient of compactness, Kc 1,337 1,197 Shape factor, Kf 1,158 0,379 River order 4 5 Drainage density (km/km )2 3,627 1,07067 Average length of surface runoff (km) 0,069 0,220 95 Sinuosity of the watercourse 1,009 1,186 Maximum slope (%) 6 12 Average slope (%) 3 3 Minimum slope (%) 0,8 0,8 Maximum altitude (m) 97 114 Average altitude (m) 65 70 Minimum altitude (m) 40 46 Figure : Olaria Stream Basin (smaller) in yellow and Cambaí Stream Basin in pink 96 EXERCISE 2 : IN ORDER TO CALCULATE THE RUNOFF CAUSED BY THE PREVIOUS RAINFALL, THE SWAT MODEL USES THE "SCS NUMBER RUNOFF EQUATION" (SEE PAGE 99 OF THE THEORETICAL DOCUMENTATION). a) EXPLAIN THE VARIABLES THAT MAKE UP THE EQUATION: ܳ௦௨௥௙ ൌ ሺோ ೏೔ೌିூೌሻమ ሺோ೏೔ೌିூೌାௌሻ b) CONSIDER THAT S REPRESENTS THE MAXIMUM POTENTIAL RETENTION OR STORAGE OF WATER IN THE SOIL AND THAT Ia IS THE INITIAL RETENTION DUE TO WATERLOGGING, VEGETATION OR INFILTRATION. ASSUME THAT Ia =0.2S REPLACE IN THE EQUATION OF ITEM a) AND GET THE FORMULA FOR Qsurf AS A FUNCTION OF Rdia and S. c) NOW ASSUMING THAT THE RETENTION PARAMETER ࡿ ൌ ૛૞Ǥ ૝ሺ૚૙૙૙ ࡯ࡺ െ૚૙ሻSUBSTITUTE AND FIND THE FORMULA FOR Qsurf AS A FUNCTION OF CN d) EXPLAIN THE MEANING OF THE CN PARAMETER AND HOW IT CAN BE OBTAINED FROM THE TABLES IN THE THEORETICAL DOCUMENTATION (TABLES HERE). e) CREATE A TABLE CONTAINING THE VALUES OF THE CN PARAMETER, WITH VALUES BETWEEN 0 AND 100, APPLICABLE TO THE CAMBAÍ AND OLARIA BASINS. IT WILL BE USED TO PREDICT THE RESULTING RETENTION AND RUNOFF FOR TYPICAL LAND USES 97 AND TYPES. TO FIND OUT THE PERCENTAGES FOR EACH LAND USE CLICK HERE. % %$&, $ $XV H Hʙ km2 Rice Indu stry/ Urba n Field Nati ve Vege tatio n / Ripa rian Fore st Hydr ogra phy CN NO. AVE RAG E PON D. Cam baí Pott ery CN % *CNii f) CALCULATE THE POINTED AVERAGE OF THE CN VALUE FOR EACH BASIN AND FROM THEM APPLY THE EQUATION FOR CALCULATING THE SURFACE DRAINAGE Qsurf FOR THE CAMBAÍ AND OLARIA BASINS, FOR THE PRECIPITATION INDICATED IN ITEM a). 98 3.3.2 Description of the application of the SCS Method Models that simulate the water cycle make it possible to predict the environmental conditions and possible impacts of different land uses through the use of physical and mathematical representations of reality. In this study, the theoretical foundations incorporated into the SWAT model (Soil and Water Assessment Tool) were used to analyse surface water runoff in two watersheds near the city of Itaqui, RS, Brazil. To obtain the surface runoff flow (Qsurf) in the basins studied, the following variables were considered: daily precipitation (Rday), initial abstraction (Ia), retention parameter (S) and SCS curve number (CN variable). The relationship between the variables was derived from the water balance in each basin and is represented by the equation: Qsurf=(Rday-Ia)2/(RdayIa +S). This equation was modified so that it only depended on daily precipitation and the CN parameter; this parameter was obtained for each basin as the weighted average of the different existing land uses. To determine the different land uses and the corresponding areas, as well as to visualise the regions, the layout of the basins and the drainage network, some geoprocessing programmes were used, such as INPE's SPRING GIS and Google tools (Maps, Google Earth, GEE). To choose the CN numbers corresponding to each region, the US Department of Agriculture's SCS number tables available in the SWAT Model's theoretical documentation were used, suitable for a 5% slope and later corrected to the average slope of the basins studied. As a result of the analyses, the following considerations were made: i) In the Rice, Native Field, Hydrography, Industry/Urban and Riparian Forest/Dense Vegetation categories, the Olaria stream basin had land use percentages of 13.9%, 27.8%, 1.3%, 99 46.3% and 1.2% respectively; the Cambaí stream basin had 84.6%, 5.8%, 3.8%, 1.5% and 3.0%. ii) The soils found correspond to types C and D of the SCS tables, with slow or very slow infiltration rates. iii) The antecedent humidity condition considered was type III, humid, in terms of water retention capacity, due to the clayey and floodplain nature of the predominant soils in the basins studied. iv) The average value of the CN number calculated for each micro-basin was: CNOlaria= 87 and CNCambaí= 80. v) The flows obtained for the conditions studied varied in the order of 60 to 300 litres/sec for the Olaria stream, and an order of magnitude more for the Cambaí stream. A continuation of this work could be to analyse the temporal variation of the CN number as a function of the time of year, influenced by rice crop cycles and seasonal rainfall. The results obtained could be used to help implement the SWAT model in the basins considered or in future soil and water analysis projects, serving as a basis for study. Cambaí Total (km^2) 157. 519 Pottery Total (km^2) 14.3 68 % %$&,$XVHʙ Rice Fiel d Nati ve Hydro graph y Indu stry /Urb an Ma ta Cili ary Veget ation dense Cambaí (km^2) 134. 9.16 6.00 2.47 2.4 3.45 100 00 4 Pottery (km^2) 2.00 4.00 0.19 8.01 0.0 1 0.15 % %$&,$XVHʙ Rice Fiel d Nati ve Hydro graph y Indu stry /Urb an Ma ta Cili ary Veget ation dense Cambaí (%) 85.0 7 5.82 3.81 1.57 1.5 5 2.19 82.49 Pottery (%) 13.9 2 27.8 4 1.35 55.78 0.0 7 1.05 87.41 CN Table 82 83 98 91 74 73 Figure: Spreadsheet for calculating the number of LU weighted by the area of each type of land use in the Olaria and Cambaí stream basins considered. 3.3.3 Programming Flow Calculations using the SCN Method in Scilab //Surface Flow Calculation - Pgs. 122 to 132 SWAT //Paste in https://cloud.scilab.in/ on the left and run Area_ola=14.368E6 //m^2 Area_cam=157.519E6 //m^2 //Daily precipitation in January 2019 107 qpeak Maximum flow rate, cubic metres per second, m3 /s C Flow coefficient, dimensionless i Rain intensity, millimetres of rain per hour, mm/h Area Sub-basin area, square kilometres, km2 3.6 Unit conversion factor (1h/3600s*1000mm/1m*1000000m /1km )22 3.5 Concentration time The flow calculation equation of the rational method contains, in a hidden way, the time during which the rainfall occurs, as this is what will define the intensity of the rainfall. The same volume of rainfall, whether it falls in 1 hour or 24 hours, will have completely different effects on the runoff values in the area under consideration. Processes such as erosion, sedimentation, nutrient transport and many others will depend critically on these flow values. For this reason, it is very important to define the time to be used to calculate the flow. The time of concentration is the amount of time from the start of a rainfall event until the entire area of the sub-basin is contributing to the flow at the outlet. In other words, the time of concentration is the 108 time for a drop of water to flow from the most remote point in the subbasin to the sub-basin outlet. The time of concentration is calculated by adding the time of surface flow (the time it takes for the flow from the most remote point in the sub-basin to reach the channel) and the time of channel flow (the time it takes for the flow in the upstream channels to reach the outlet): t_conc=t_ov+t_ch where tconc is the time of concentration for a sub-basin (h), tov is the time of concentration for surface runoff on the ground (h) and tch is the time of concentration for channel runoff (h). tconc Time of concentration for a sub-basin, hours tov Time of concentration for surface runoff on the ground, hours tch Time of concentration for flow in the channel, hours 3.5.1 Time of Concentration of Surface Runoff on the Ground The runoff concentration time, tov , can be calculated using the equation: t_ov=L_slp/(3600*v_ov) where Lslp is the length of the sub-basin slope (m), vov is the runoff 109 velocity (m s-1 ) and 3600 is a unit conversion factor. tov Runoff concentration time, hours Lslp Length of sub-basin slope, metres vov Runoff velocity, metres per second 3600 Unit conversion factor The runoff velocity can be calculated using the Manning equation, considering a 1 metre wide strip below the sloping surface: v_ov=q_ov^0.4*slp^0.3/n^0.6 where qov is the average runoff rate (m3 s-1), slp is the average slope of the sub-basin (m m-1 ), and n is the Manning's roughness coefficient for the sub-basin. Assuming an average flow of 6.35 mm/h and conversion units of qov Average runoff rate, cubic metres per second slp Average slope in metres of elevation per horizontal metre n Manning's roughness coefficient for the sub-basin 110 Substituting equation 2:1.3.5 into equation 2:1.3.3 we get CONCENTRATION TIME FOR THE FLOW IN THE CHANNEL The flow concentration time of the channel tch can be calculated using the equation: where Lc is the average length of the sub-basin drainage channel 111 (km), vc is the average channel velocity (m s-1 ), and 3.6 is a unit conversion factor. Lc Average length of the sub-basin drainage channel, km Vc Average channel velocity, m/s 3,6 Unit conversion factor The average flow duration of the channel can be calculated using the following equation: where L is the length of the channel from the furthest point from the sub-basin outlet (km) and LCEN is the distance along the channel to the sub-basin centroid (km). Assuming that Lcen = 0.5 L, the average flow length of the channel is L Channel length from the furthest point from the sub-basin outlet, km LCEN Distance along the channel to the centroid of the sub-basin, km The average velocity can be calculated from the Manning equation assuming a trapezoidal channel with 2:1 side slopes and a 10:1 widthdepth bottom. 112 where vc is the average velocity of the channel (m s-1 ), qch is the average flow rate of the channel (m s3-1 ), slpch is the slope of the channel (m m-1 ), and n is the Manning's roughness coefficient for the channel. vc Average channel velocity, m/s qch Average channel flow rate, cubic metres per second slp ch Channel slope, vertical metres per horizontal metre n Manning's roughness coefficient for the channel To express the average flow of the channel in units of mm/h, the following expression is used: where qch is the average flow of the channel (mm h-1 ), Area is the area of the sub-basin (km2 ), and 3.6 is a unit conversion factor. The average flow of the channel is related to the unit flow of the source area (unit source area = 1 ha), q*ch . qch Average channel flow in mm of water per hour, mm/h 113 q*c h Unit flow rate of the source area (unit source area = 1 ha), mm/h Area Sub-basin area, in km2 3,6 Unit conversion factor q*0 Flow rate of the unit area (of 1 ha) measured in millimetres per hour, mm/h Area Sub-basin area, in square kilometres, km2 100 Unit conversion factor Assuming that the unit flow rate is 6.35 mm/h and substituting equations 2:1.3.11 and 2:1.3.12 into 2:1.3.10 we have substituting equations 2:1.3.9 and 2:1.3.13 into 2:1.3.7 we get 114 where tch is the time of concentration for the channel flow (h), L is the length of the channel from the furthest point of the sub-basin outlet (km), n is the Manning's roughness coefficient for the channel, Area is the area of the sub-basin (km2 ), and slpch is the slope of the channel (m m ). -1 tch Concentration time for channel flow in hours L Channel length from the furthest point of the sub-basin outlet, km n Manning's roughness coefficient for the channel, dimensionless Area Sub-basin area, km2 slpch Channel slope, vertical metres per horizontal metre Table 2.1-4: Values of Manning's roughness coefficient, n, for channel flow (Chow, 1959). 1 115 Although some of the assumptions used in the development of equations 2:1.3.6 and 2:1.3.14 may seem undemanding, the values of the times of concentration obtained generally give satisfactory results for homogeneous sub-basins. Since equations 2:1.3.6 and 2:1.3.14 are based on hydraulic considerations, they are more reliable than purely empirical equations. Once we know the time of concentration, during which precipitation occurs, we can proceed to find the maximum flows within the basin, which is presented in the next section. 3.5.2 Calculating Flows within the Time of Concentration Flow coefficient (C) - The flow coefficient is ratio of the flow or runoff rate, Q_peak, produced at peak time, and the volume of incoming precipitated water (i. Area), which The runoff coefficient theoretically has values 116 is the daily precipitation, Rday . The coefficient varies depending on the storm and is calculated using the equation. C= Q_surf /R_day between 0 and 1. In practice the values vary between 0.05 and 0.5 for most basins. Rainfall intensity: Rainfall intensity is the average rate of rainfall during the time of concentration. Based on this definition, it can be calculated using the equation I= R /ttcconc or I= R_tc/T_conc - concentration time - tc proportional to the 24-hour period - short storm I=R_tc/T_conc (concentration time) - I is the rainfall intensity (mm/h). - R_tc is the amount of rain that drops during the concentration time (mm H2O). - T_conc is the time of concentration for the sub-basin (h). Rtc = αtcԫ΃day (tc proportional to 24-hour period). - Rtc is the amount of rain that falls during the concentration time (mm H20). - αtc is the fraction of daily αtc,min = Rtc /RdayΚԫΥconc Κԫ = t /24conc - The maximum values of αtc occur for storms of short 123 Q'   Amount of runoff generated in the sub-basin on a given day. mm 2 (millimetre water column) Bar Q   Surface runoff stored or delayed from the previous day. mm 2 (millimetre water column) Bar surla g Runoff delay coefficient. -------------------- ------------ -------------- ------------ t   Time of concentration of the sub-basin. h (hour) min (minute) vol   Flow volume after transmission losses. 3 (cubic metre) 3 (cubic kilometre) a Regression intercept for a channel of length L and width W. 3 (cubic metre) 3(cubic kilometre) b Regression slope for a channel of length L and width W. 3 (cubic metre) 3(cubic kilometre) vol   Flow volume before transmission losses. 3 (cubic metre) 3(cubic kilometre) vol ℎ Threshold volume for a channel of length L and width W. 3 (cubic metre) 3(cubic kilometre) 124 q   Peak rate after transmission losses. 3Ȁ (cubic metre per second) 3Ȁ ℎ(cubic kilometres/ hour) dur  Duration of the flow. h (hour) s (seconds) q   Peak rate before accounting for transmission losses. 3Ȁ (cubic metre per second) 3Ȁℎ Area Area of a sub-basin. 2(square kilometres) 2(squar e centimetre) q   Maximum flow rate. 2Ȁ (square metre per second) 2Ȁℎ (quad km/hour) k Decay factor. െ1 െ1(m etre kilometre) ି1(centi metre) a Unit channel regression intersection. 3 (cubic metre) 3(cubic kilometre) b Regression slope of a unit channel. 3 (cubic metre) 3(cubic kilometre) kℎ Effective hydraulic mm/h m/s (metre 125 conductivity of the channel alluvium. (millimetre per hour) per second) L Length. km (kilometre) cm (centimetre ) W Width. m (metre) mm (millimetre ) The algorithm, written in the mathematical computer language Scilab, for taking into account the delay and transmission losses in surface runoff, is shown below. //Surface runoff delay equations surlag=[1:2:12]; //runoff delay coefficient t_conc = 12; // concentration time for the sub-basin (h) Q_stor = 50; //Stored or delayed runoff from the previous day (mm H2O) Q_surfgen = 50; // amount of runoff generated on a given day Q_surf = (Q_surfgen + Q_stor)*(1-exp((-surlag)/t_conc)); //amount of runoff discharged into the main channel on a given day //Transmission loss equations //Calculation of unit channel parameters K_ch = 25; // effective hydraulic conductivity of the channel alluvium (mm/h) 126 Area = 20; //area of the sub-basin (km^2) q_peak = 28.32; //maximum flow rate (m^3/s) vol_Qsurfi = 41938; //initial flow volume without losses (m^3) Q_surf = 20.4; // Surface runoff (mm H20) dur_flw = Q_surf*Area/(3.6*q_peak); //flow flow duration mprintf('Flow duration = %f',dur_flw); k_r = -2.22*log(1-(2.6466*K_ch*dur_flw/vol_Qsurfi)); //decay factor (m^-1*km^-1) a_r = -0.2258*K_ch*dur_flw; //intercept of regression (m^3) of the unit channel; b_r = exp(-0.4905*k_r); //regression slope of the unit channel mprintf('Channel parameters unit:\n a_r = %f \n b_r =%f\n k_r =%f',a_r, b_r, k_r); L = 8.05; //length of the channel from the furthest point from the sub-basin outlet (km) W = 64; //average flow width (m) b_x = exp(-k_r*L*W); // slope of the regression for a channel of length L and width W a_x = a_r*(1-b_x)/(1-b_r); //intercept regression for a channel of length L and width W (m^3) mprintf('Channel parameters: \n a_x = %f \n b_x =%f\',a_x, b_x); vol_thr = -a_x/b_x; //threshold volume for a channel of length L and width W (m^3) mprintf('Channel threshold volume =%f',vol_thr); q_peaki = q_peak; //peak rate before transmission losses (m^3/s) q_peakf = 1/(3600*dur_flw)*(a_x-(1-b_x)*vol_Qsurfi)+b_x*q_peaki; //peak rate after transmission losses (m^3/s) 127 mprintf('Peak transmission rate =%f',q_peakf) if vol_Qsurfi > vol_thr then vol_Qsurff = a_x+(b_x*vol_Qsurfi); //flow after transmission losses (m^3) else vol_Qsurff = 0; end mprintf('Flow volume after losses =%f',vol_Qsurff); 3.6.3 Conclusions The calculation of surface runoff was presented here using algorithms that work with the CN number method, as an example for two basins in the city of Itaqui, RS, Brazil. These algorithms are used to calculate daily flows from precipitation data, either obtained daily from meteorological stations or through simulations created from monthly average data. The calculation of the time of concentration was also developed, which is very important because it is within this time that maximum runoff events usually occur, influencing various processes such as erosion. The calculation of the delay in the flow reaching the main channel, which occurs in large basins, has also been carried out. These are initial algorithms, of great didactic value for understanding flow models, but they will continue to be developed to improve their connection with others. In the same way, losses that occur during the transmission of the flow began to be dealt with in this edition. 128 Once we have this initial basis for calculating precipitation and runoff, we now move on to the phenomenon of evapotranspiration 3.7 Evapotranspiration 3.7.1 Priestley-Taylor, Hargreaves and Penman-Monteith methods 3.7.1.1 Introduction Considering an agricultural plot with uniform soil type and cover, topography and management, we have a favourable basic unit for defining physico-mathematical models that represent the hydrological cycle, thus promoting more accurate water and soil management. In particular, the hydrological cycle is of fundamental importance in the management of agriculture and the environment. Evapotranspiration is a variable whose concept includes all the processes that remove water from the earth's surface to put it into the atmosphere. In the region under study, these include evaporation from the plant canopy, transpiration and evaporation from the soil. While precipitation is the main cause of water entering the system, evapotranspiration is the primary cause of water leaving it, especially when we consider large tracts of land. For this reason, its knowledge and estimation is of great importance on a daily scale in order to optimise water management in agricultural activities. Physicalmathematical models can be used to monitor or predict evapotranspiration in a given region. However, these models involve variables and relationships that represent the situation of the atmosphere and the earth's surface, which need to be well understood and their magnitude assessed. 129 3.7.1.2 Concepts and Methods Potential evapotranspiration (PET) was a concept originally presented by Thornthwaite (1948) as part of a climate classification project. He defined PET as the rate at which evapotranspiration would occur from a large area evenly covered with ground vegetation and with access to an unlimited supply of soil water and which was not exposed to advection or heat storage effects. Because the rate of evapotranspiration is strongly influenced by a number of characteristics of the vegetative surface, Penman (1956) redefined PET as "the amount of water transpired... by a short green crop, shading the ground completely, of uniform height and never short of water". Penman used grass as his reference crop, but later researchers (Jensen, et al., 1990) suggested that alfalfa at a height of 30 to 50 cm may be a more suitable choice. Numerous methods have been developed to estimate PET. Three of these methods have been incorporated into the SWAT model, and are also used in this work: the Penman-Monteith method (Monteith, 1965; Allen, 1986; Allen et al., 1989), the Priestley-Taylor method (Priestley and Taylor, 1972) and the Hargreaves method (Hargreaves et al., 1985). As the aim of these models is to generate daily PET values for later use in other calculations, it is easy for the user who prefers to apply a different method, or to take advantage of data calculated by meteorological stations, to use this potential evapotranspiration data as input for subsequent calculations related to the hydrological cycle. The three PET methods considered vary in the number of inputs required. The Penman-Monteith method requires solar radiation, air temperature, relative humidity and wind speed. The Priestley-Taylor 130 method requires solar radiation, air temperature and relative humidity. The Hargreaves method only requires air temperature. We can see that PET calculations have developed by adapting to different situations of input data availability. We will now describe the variables, units and relationships that form part of the calculations made, as an initial part of the preparation of the algorithms that make it possible to carry out the calculation, which requires in-depth knowledge of them. Each of these variables will later be important for calculating potential evapotranspiration in the different methods. The approach used is aimed not only at the structure of the calculation, but also at carrying out a learning process through successive repetition of the calculations, at different levels of depth and accuracy, starting, when appropriate, with the most general calculations, and then delving into the complexity of the parts of the calculation, i.e. the individual calculations of the variables used in the general calculation. We'll start with an example calculation using the most general and complete equation, which is the Penmann-Monteith equation. Penmann-Monteith equation As a first contact with the P-M equation, we see in its expression in mathematical language a set of variables with which we must familiarise ourselves, linked by mathematical operators. As each variable, in turn, maintains a relationship with others, then the 131 equation can appear in other forms, which in many situations can make calculations easier. In this example, the variable cp is replaced by others, arriving at a new form of the equation. Where λ is the latent heat of vaporisation (MJ/kg), E is the depth evaporation rate (mm/d), Et is the maximum transpiration rate (mm/d), λE is the latent heat flux density (MJ/d/m2 ), Δ is the slope of the saturation vapour pressure-temperature curve, de/dt (kPa/°C), Hnet is net solar radiation (MJ/d/m2 ), G is the density of heat flow into the soil (MJ/d/m2 ), ρair is the air density (kg/m3 ), cp is the specific heat at constant pressure (MJ/kg/°C), e0z is the air saturation vapour pressure at height z (kPa), ez is the air water vapour pressure at height z, γ is the psychrometric constant, rc is the plant canopy resistance and ra is the air layer diffusion resistance aerodynamic resistance For an initial understanding of the equation, and initial contact with the order of magnitude of the values of the variables and their units, it is important to programme an algorithm in Scilab that performs this calculation, which is presented below. However, each of the variables is also obtained through relationships with other variables, so we will proceed to analyse each of them independently, and then combine all the algorithms into one. // PENMAN-MONTEITH PROGRAMMING IN SCILAB lambda_e= 1 //M (mega) J (joules) m^-2 d^-1 latent heat flux density ( L=q/m) 132 E=1 //mm d^-1 depth evaporation rate T_av=25; //degrees C - average daily air temperature e_0=exp((16.78*T_av-116.9)./(T_av+237.3)); delta_PT=4098*e_0/(T_av+237.3)^2 delta=delta_PT //de/dT (kPa °C^-1) slope of the saturation vapour pressure-temperature curve H_net=11 // Net Radiation (liquid in question of what's left) J.m^-2.d^-1 G=0 //MJ m^-2 d^-1(megaJoules divided square metre divided day) density of heat flow into the ground p_air =1 //kg m^ 3 (kg per cubic metre) Air density C_p=1.005 //MJ kg^-1 °C^-1 specific heat at constant pressure e_0_z=94 //(kPa) air saturation vapour pressure at height z e_z= 0.6108 //(kPa) air water vapour pressure at height z. Y=1 //(kPa .°C^-1) psychrometric constant r_c=1 //s m^-1 - plant canopy resistance r_a=1 //s m^-1 diffusion resistance of the air layer K_l = 8.64 * 10^4 // coefficient needed to ensure that the two terms in the numerator have the same units P=1 //kPa atmospheric pressure lambda_E=(delta*(H_net-G)+p_air*C_p*(e_0_ze_z)/r_a)/(delta+Y*(1+r_c/r_a)) 3.7.1.3 Calculating Net Radiation H NET 139 where Isc is the solar constant (4.921 MJ m-2 h-1 ), EO is the eccentricity correction factor of the Earth's orbit, and ω is the angular velocity of the Earth's rotation (0.2618 rad h-1 ), the time of rising, TSR , is defined by the equation δ is the solar declination in radians, and ϕ is the geographical latitude in radians. If multiply all the constants together and you get 3.7.1.4 Solar Radiation under Cloudless Skies The insolation received at a given location on the Earth's surface can vary, and one of the factors responsible for this variation is cloudiness (Silva, 2011). Clouds cause a variation in the intensity of solar radiation incident on the surface, and this variation in radiation explains the variability of incident radiation in a region (Bastos et al., 2002). In addition, only part of the solar radiation reaches the earth's surface due to atmospheric transmissivity (Pereira and Oliveira, 2011). Therefore, before reaching the ground, the characteristics of solar radiation (intensity, spectral and 140 angular distribution) are affected. These changes depend on the thickness of the atmospheric layer, also identified by a coefficient called Air Mass (AM), the sun's zenith angle, the Earth-Sun distance and atmospheric and meteorological conditions, influencing the performance of solar panels (Silva, 2011). Gonçalves (2013) carried out energy production simulations of solar panel installations at Belém International Airport over the course of a year, stipulating the energy generated according to the solar radiation received in certain periods and showing that there is a direct relationship between the solar radiation received by the solar panels and the electrical energy produced by them. When solar radiation enters the earth's atmosphere, some of the energy is removed by scattering and absorption. The amount of energy lost is a function of the transmittance of the atmosphere, the composition and concentration of the air constituents at the site, the length of the path that the radiation travels through the air column, and the wavelength of the radiation. Due to the complexity of the process and the precision of the information needed to accurately predict the amount of radiant energy lost during passage through the atmosphere, the SWAT model makes a broad assumption that approximately 20 per cent of extraterrestrial radiation is lost as it passes through the atmosphere under cloudless skies. Using this assumption, the maximum possible solar radiation, HMX , at a given location over the earth's surface is calculated as: where the maximum possible solar radiation, HMX , is the amount of 141 radiation that reaches the Earth's surface under a clear sky (MJ m-2 d ).-1 In this way, we have the maximum solar radiation that reaches the Earth's surface, which will be greater than the effective daily radiation that reaches the ground, Hday , due to the presence of clouds and atmospheric effects. The intensity of solar radiation that reaches a horizontal surface varies due to the attenuation it suffers as it passes through the atmosphere, due to the presence of clouds, dust, pollution and others. Naturally, on a cloudy day, the intensity of solar radiation will be lower and consequently the module's performance will suffer. The opposite is true on a clear day or with a cloudless sky (Marques; Pereira; Assis, 2000). 3.7.1.5 Daily Solar Radiation For the use of models that depend on solar radiation, daily solar radiation data can be obtained directly from a data file measured at a meteorological station, or it can be created by statistical models based on descriptive parameters of the regional climate, also obtained from previous meteorological records. 3.7.1.6 Solar Radiation per Hour I0 is the extraterrestrial radiation for one hour, centred around the hour angle wt 142 //Solar Radiation per Hour - Scilab Programming I_sc= 4.921//MJ/ m^2/ h Cte.Solar E0=0.0167 //Excentricity orbit T delta=23.27*%pi/180 //solar declination fi=29*%pi/180 //Latitude in radians w=0.2618 //rad.h-1 ang. vel. Earth t=[0:1:11] //solar time I0=I_sc*E0*(sin(delta)*sin(fi)+cos(delta)*cos(fi)*cos(w*t)) plot2d(t,I0) Results: // Solar radiation per hour I_frac= 14.90 // (MJ/m2/d) H_day=16 //(MJ/m2/d)Solar radiation reaching the ground on a given day I_hr= I_frac*H_day Result: I_hr = 238.4 MJ/m2/h // Fraction of daily solar radiation, I_frac 143 delta=23.27*%pi/180 //solar declination fi=29*%pi/180 //Latitude in radians w=0.2618 //rad.h-1 ang. vel. Earth t=[0:1:11] //solar time // Daily Net Radiation alpha= //shortwave reflectance or albedo H_day= //Shortwave solar radiation reaching the earth MJ/m /d2 H_b= //Longwave net radiation MJ/m /d2 H_net=(1-alpha)*H_day+H_b //Net radiation MJ/m /d2 // Long wave radiation //H_R is the radiant energy MJ/m /d2 epsilon= // is the emissivity sigma=4.903E-9 //MJ/m2 /d/K^-4 Stefan-Boltzmann constant T_K= //Average air temperature in Kelvin (273.15+Celcius) H_R=epsilon*sigma*T_K^4 144 // Long wave net radiation H_b f_cld= //cloud cover adjustment factor epsilon_a= //atmospheric emittance epsilon_vs= //vegetative or soil emittance H_b=f_cld*(epsilon_a-epsilon_vs)*sigma*T_K^4 //Cloud cover adjustment factor f_cld MJ/m^2/day a= //constant b= //constant H_day= //MJ/m^2/day , solar radiation reaching the earth on a given day H_MX= //MJ/m^2/day , maximum possible solar radiation reaching the surface on a given day f_cld= a*H_day/H_MX-b //MJ/m^2/day //Cloud cover adjustment factor Emittances can be combined and expressed as a function of the vapour pressure of the day, in KPa epsilon_net=epsilon_a-epsilon_vs=-(a_1+b_1*sqrt(e)) Combining equations results in a general equation for net longwave radiation 145 H_b=-(a*H_day/H_MX-b)*(a_1+b_1*sqrt(e))*sigma*T_K^4 Using the coefficients from Doorenbos and Pruitt (1977), the equation is as follows: H_b=-(0.9*H_day/H_MX-0.1)*(0.34+0.139*sqrt(e))*sigma*T_K^4 The values of the coefficients proposed in the SWAT model are shown in the following table. In this way, solar radiation calculations are resolved, based on the daily values measured or generated, and the maximum possible radiation calculated with the physical models of the Earth's orbit. Table of Variables Used Varia ble Value Definition 146 Isc rate of total incident solar energy of all wavelengths in a unit area UA distanc e from the sun 1UA = 1.496 x 108 km Astronomical unit I0n Daily extraterrestrial irradiance incident on a normal surface (MJ m-2 h-1) E0 Earth eccentricity correction factor I0 Daily extraterrestrial irradiance incident on a horizontal surface (MJ m-2 h-1) H0 Daily extraterrestrial radiation Assuming that EO remains constant during the time interval of one day 147 H0 or and converting dt from time to the hour angle, the equation can be written TSR e Tss Sunrise time on a solar day (h) e Sunset time on a solar day (h) H0 If we multiply all the constants together we get HMX Maximum possible solar radiation (MJ m2 d-1) 3.8 Temperature The influence of temperature on numerous physical, chemical and biological processes is undeniable. Temperature can be inserted into the models by means of measurement data files, or it can be generated with 148 monthly average values extracted from regional records. Soil and water temperatures are derived from air temperature. Below we present the Scilab algorithms for the equations, which contain explanations of the variables. Hourly air temperature T_hr T_av= //degrees Celsius - average temperature of the day T_mx= //degrees Celsius - maximum temperature of the day T_mn= //degrees Celsius - minimum temperature of the day T_hr=T_av+(T_mx-T_mn)/2*cos(0.2618*(hr-15)) Soil temperature (Carslaw and Jaeger, 1959) - T_soil_z_dn //T_soil_z_dn is the soil temperature in degrees Celsius at depth z (mm) and on the day of the year dn T_AA= //degrees Celsius , average annual temperature A_surf= //degrees Celsius , amplitude of the temperature fluctuation on the surface dd= //mm , damping depth w_tmp= //angle frequency In the model adopted, the calculation of soil temperature is based on parameters from the same day and the previous day, as explained in the following algorithm. z= //mm , depth in the ground