1
Rice Cultivation on the edge of the Sahara;
Modeling the hydrological effects of irrigated rice at Office du Niger
2 Cover photo: Rice cultivation at Office du Niger
Source: http://www.hornn.biz/dadodi/images/stories/373.jpg Martijn Vermeer Studentnumber: 10243879 Cornelis Krusemanstraat 73-3 1075 NJ Amsterdam 06-25371919 martijnvermeer92@gmail.com Bachelor Thesis Earth Sciences Faculty of Science University of Amsterdam
Supervisor: dhr. prof. dr. ir. W. Bouten Amsterdam 30-06-2014
3
Summary
The Inner Niger Delta is a large wetland area in the semi-arid Sahel region in central Mali. The water level and corresponding inundated area in delta depends on the peak discharge of the Niger and Bani rivers (Liersch et al., 2013). The seasonal flooding of the delta is vital for its 1.5 million inhabitants as it enables the cultivation of non-irrigated rice. On top of that, the delta is of high ecological importance, with 3 to 4 million resident or migratory birds and a high biodiversity (Zwarts et al., 2009). The presence of irrigated vegetation, mainly rice, upstream of the delta affects the water balance in the Delta. Most of this irrigated vegetation is localized in an 83.000 ha irrigation scheme operated by Office du Niger. Food demand is growing in Mali resulting in plans to expand the irrigated areas upstream to 200.000-300.000 ha by the year 2020. This study investigates the influence of this irrigated rice cultivation on the discharge of Niger. For this purpose the hydrological HBV-model is extended with an irrigated rice module. It was found that due to irrigation at Office du Niger the peak discharge is currently reduced with 41m3/s, increasing to 102-158m3/s in 2020 depending on the extent of the expansions. The average reduction of the discharge throughout the year is even higher, reaching a maximum in the end of the dry season. In this period the discharge of the Niger is at lowest and without the construction of dams no water would be available for irrigation. Therefore Office du Niger indirectly causes an even higher peak discharge reduction as dams are forced to withhold water during the wet season that can be used for irrigation in the dry season. This reduction of the peak discharge could have serious consequences for the inhabitants and ecology of the Inner Niger Delta.
4
Table of Contents
Summary ...3 1. Introduction ...6 2. Method ...8 2.1. Model description ...8 2.2. Data preparation ...9 2.3. Calibration...92.4. Irrigated rice implementation ... 10
2.4.1. Create irrigated catchments ... 10
2.4.2. Evapotranspiration ... 10 2.4.3. Irrigation ... 11 2.5. Scenarios ... 11 3. Results ... 12 3.1. Water balance ... 12 3.1.1. Inputs ... 12 3.1.2. Storages ... 13 3.1.3. Outputs ... 13 3.1.4. Balance ... 15 3.2. Scenarios ... 16 3.2.1. Irrigation ... 16 3.2.2. Discharge... 16 4. Discussion ... 18 4.1. Water balance ... 18 4.2. Scenarios ... 18 4.3. Input data ... 19 4.4. Model ... 19 5. Conclusion ... 20 6. Acknowledgements ... 21 7. Literature ... 22 8. Appendix... 24 8.1. Data preparation ... 24 8.1.1. ArcGIS ... 24
8.1.2. Conversion ArcGIS to Matlab ... 25
5
8.1.4. Lumping ... 30
8.2. Calibration... 32
8.3. Irrigated rice implementation ... 33
8.3.1. Create Irrigated Catchments ... 33
8.3.2. Interpolation evapotranspiration ... 35
8.3.3. Conversion and interpolation irrigation ... 35
8.3.4. Irrigation Module ... 36
6
1. Introduction
The Inner Niger Delta is a large wetland in the semi-arid Sahel region in central Mali covering an area of approximately 36,000 km2, consisting mostly of floodplains. (Zwarts et al., 2009). The water level and coupled inundation of the floodplains is determined by the discharge of the Niger and Bani Rivers that enter the delta from the south-west. The catchment of these rivers, the upper Niger catchment, has a large seasonal variability in precipitation caused by the west-African monsoon. As a result the discharge of the Niger and the Bani strongly fluctuates between the wet and dry season (Liersch et al., 2013). The peak flood caused by the high precipitation in wet season (July-October) reaches the delta in October and the northern plains even two months later when water levels in the south are already declining.
The Inner Niger Delta is inhabited by 1.5 million people and this population is expected to grow to 3.5 million in 2050. These people provide for their livelihoods by fishing, farming and stock farming, strongly depending on the delta as seasonal inundation enables cultivation of rice and bourgou (Liersch et al., 2013). On top of that, the delta is of high ecological importance, with 3 to 4 million resident or migratory birds and a high biodiversity (Zwarts et al., 2009).
Just upstream of the delta an irrigation scheme of 83,000 ha is operated by Office du Niger, a semi-autonomous government agency. The irrigation scheme itself is also called Office du Niger and is the largest irrigated area in Mali, producing 810,000 tons of paddy rice corresponding to 40% of the national production (Vandersypen et al., 2006). During the wet season almost the entire area is used for rice cultivation, in the dry season only 26,000 ha is irrigated of which 11.000 ha rice (van Dijk, 2008). Together with the population, food demand is growing in Mali, resulting in plans to extend Office du Niger to 200,000-300,000 ha by the year 2020 (Zwarts, 2012). The irrigation at Office du Niger influences the discharge of the Niger and therefore the water balance in the Inner Niger delta. Several studies (Kuper et al., 2003; Zwarts et al., 2005; Zwarts, 2010) have been carried out that investigated the influence of irrigation at Office du Niger on the discharge of the Niger. On average
7
8.3% of the discharge of the Niger is diverted to the irrigation scheme (Zwarts et al., 2005). The amount of discharge that is diverted varies seasonally between 60m3/s in January and 130m3/s in October (table 1). Zwarts (2010) examined that an extension of the irrigation scheme to 280.000 ha leads to an increase in maximal water extraction from 130m3/s at present to nearly 500m3/s. As a result the loss of the inundated area in the Inner Niger Delta, which is currently 13%, could increase to 20% or even 27%. These studies however simply subtract the amount of water that is diverted to the irrigation scheme from the discharge. This is not realistic as the water does not disappear, only part of the water is lost due to evapotranspiration and the remainder stays in the system being stored in the ground or becoming runoff. This research will build on these previous studies by incorporating the irrigated rice in a hydrological model. An important aspect of simulating irrigated rice is determining the evapotranspiration. For rice evapotranspiration can be estimated by multiplying the reference evapotranspiration with the crop factor. The best way to calculate reference evapotranspiration in the Sahel region is with the use of Penman-montieth equation (Reas, Sy & Feyen, 1995). The crop factor of paddy rice varies between 1.0 and 1.3 depending on the climate and time in the season; Vandersypen et al. (2006) applied a crop factor of 1.2 for Office du Niger. The change from natural vegetation to irrigated rice changes the hydrology and eventually affects the water balance of Inner Niger Delta.
The overall aim of this research is to assess the impact of changing land-use, from natural vegetation to irrigated rice, at Office du Niger on the amount of discharge entering the Inner Niger Delta. The findings of the research can contribute to map the changing water balance of the Inner Niger Delta, which affects its inhabitants and ecology.
The report starts with a description of the research method (section 2), second the result are displayed (section 3). Subsequently, these results are discussed (section 4) and finally conclusions are presented (section 5).
Table 1: The amount of water taken per month, on average (m3/s), for irrigation by Office du Niger near Markala between 1989 and 2009 (“mean”) and during five years with the lowest river flow (1989-1993; "dry”) and five years with the highest flow (1994, 1995, 2001, 2003, 2008; “wet”) (Zwarts, 2010).
8
2. Method
The method consist of five parts, first a general description of the HBV-model, second the preparation of the input data for the HBV-model (2.2), third the calibration of the model (2.3), fourth the implementation of irrigated rice (2.4) and finally the different scenarios that are modeled are discussed (2.5). Section 2.1, 2.2 and 2.3 were written and conducted in close cooperation with Berend-Christiaan Wijers, since both the research of Wijers (2014) and this research use the HBV-model in the same research area.
2.1. Model description
The HBV-model is a semi-distributed precipitation runoff model originally developed by the Swedish Meteorological and Hydrological Institute (SMHI) in 1970 (Buiteveld et al., 2005). The HBV-model has been adapted to the Rhine by Buiteveld et al. (2005) and has been rewritten to Matlab code by Spaaks (2011). The HBV-model is used in many countries to model river discharges. The river Basin that is simulated in the HBV-model is subdivided into different catchments (figure 2). Per catchment the water flowing out of the catchment is calculated using the Snow & Rain module, the Soil & vegetation module and the triangular filter. Finally the river discharge is calculated from the different catchment flows.
The Snow & Rain module determines the fraction of precipitation which becomes rain and snow. As our research area is in Mali where temperatures are high all precipitation becomes rain.
The rain enters the Soil & Vegetation module and is stored in the soil as soil moisture. Soil moisture can either evaporate or become part of the upper zone storage, which is the water above the groundwater table. From the upper zone water can percolate to the lower zone, the groundwater. The amount of water in the upper zone determines the quick flow, the overland flow. Similarly the amount of water in the lower zone determines the base flow, the groundwater flow. Together these flows are the total flow. The triangular filter determines during which time step water will flow. Each catchment has a
Figure 3: General working of the HBV model
Figure 2: Upper Niger basin subdivided into catchments and the Niger and Bani rivers.
9
unique shape and size and therefore also has a different retention time of water. Larger catchments will spread the total flow from a single time step over multiple time steps whereas smaller catchments will take less time steps to spread the total flow over. At last the total flow of each catchment is summarized with the total flows of all upstream catchments to get the river discharge for each catchment. Figure 3 shows the general working of the HBV-model.
2.2. Data preparation
Input data for the model can be subdivided into four main categories forcings, initial states, parameters and other catchment data. Forcing data consists of precipitation, temperature and potential evapotranspiration data. This data is already available in gridded form for entire Africa on a resolution of roughly 5 by 5 km (Kamps, 2013). The HBV-model, however, requires data per catchment; therefore the gridded data was lumped using a catchment grid. Initial state data consist of storage of water in interception, soil, upper zone, lower zone, liquid snow and solid snow. The storages for snow do not apply for Mali as temperatures never drop below zero. Like the forcings, the other initial states are available in gridded form and were lumped. For all catchments the initial parameters values, before the calibration process, were set to the average parameter values of the HBV-Rhine model. An exception is the parameter maxbas (pMAXBAS), which determines the flow delay. pMAXBAS is linearly related to catchment size, as a larger catchment corresponds with higher flow delay. Finally catchment data about: forestation, area, connectivity and midpoints is required as input. In ArcGis a catchment grid, connectivity, area size and forestation were obtained using a DEM of the research area (U.S. Geological Survey, 2013). First catchments were created and the connectivity between them. At this point the exact research area could be defined. The confluence of the Niger and Bani rivers was chosen as the most downstream catchment. All catchments that do not flow into that catchment were deleted. Subsequently, the area size of the catchments was calculated and forestation was determined using a Bing satellite basemap. A more detailed description of the steps taken in ArcGIS is available in the appendix (8.1.1). In Matlab the catchment grid, connectivity, area size and forestation obtained in ArcMap were converted to .mat files and midpoints were calculated. The script for this can be found in the appendix (8.1.2). Furthermore, Matlab was used to clip the forcings and initial states data of entire Africa to our research area and to lump this data in order to retrieve only one averaged value per catchment. The clip and lump scripts can be found in the appendix (8.1.3 & 8.1.4.).
2.3. Calibration
As calibrating is a time consuming job and the time for this research was very limited it was decided to calibrate only a few parameters. The parameters that were used in the calibration process are pMAXBAS, pK4, pPERC, pKHQ, pBETA and pFC. pMAXBAS was chosen to calibrate because it determines the delay of the flows, this is important for the timing of the peak discharge in the wet season. pPERC determines the percolation to the lower zone storage and pK4 the amount of water in the lower zone that becomes base flow. These parameters are important to calibrate the discharge during the dry season, when precipitation is insignificant and the discharge is determined by the base flow. pKHQ is for controlling the quick flow from the upper zone storage whereas pFC and PBETA determine respectively the capacity of the soil moisture and the flow from soil moisture to the upper zone. For calibration the mean monthly discharges for the Niger at Markala were used (Zwarts, 2010).The model calculates these mean monthly discharges at Markala as well; subsequently it summates the deviation of the modeled discharges from the discharges in the literature to retrieve a
10
simple model error. This model error was minimized during several calibration cycles. The final values can be found in the appendix (8.2).
2.4. Irrigated rice implementation
The implementation of irrigated rice consists of three parts. First separate catchments are created for the irrigated areas, second evapotranspiration is calculated and third the irrigation logic itself is implemented.
2.4.1. Create irrigated catchments
In the HBV-model each catchment has its own input data and calculations are separated for the catchments. As irrigated areas require different input and calculations than non-irrigated areas, the irrigated area needs to be separated from the rest of the catchment in which it is located. In figure 3 the division of a catchment into two separate catchments, one non-irrigated and one irrigated is visualized. The existing catchment becomes the non-irrigated catchment and a new catchment is created to become the irrigated catchment. All the input data for the irrigated catchment needs to be added and for the non-irrigated catchment some data needs to be updated. The size of the area of the irrigated catchment can be specified in the beginning the model. This area is subtracted from the non-irrigated catchment. For the connectivity it was
decided to let water of the non-irrigated catchment flow into the irrigated catchment which in turn flows into the catchment into which the water of the original catchment was flowing (figure 4). Forestation in the irrigated catchment is always set to zero as the entire catchment consists of irrigated rice. All other data for forcings, initial states, parameters and midpoints are copied from the original catchment. In the appendix (8.3.1) the code with which this part was done can be found.
2.4.2. Evapotranspiration
The potential evapotranspiration of irrigated rice can be calculated by multiplying the reference evapotranspiration by the crop factor (pKC). As the quality of the original reference evapotranspiration data is uncertain (discussion 4.3.), reference evapotranspiration was recalculated for the Office du Niger catchment. For this purpose the ET0 calculator software of the FAO was used (Raes, 2009). As input the program requires temperature, air humidity, wind speed and radiation data and the coordinates and altitude of the measurement station. Mean monthly data for Segou, a city located next to Office du Niger, was used as input (Hong Kong Observatory, 2003; World Climate, 2010). Subsequently, the program uses this data to calculate the mean monthly reference evapotranspiration. This was interpolated to daily data in matlab (appendix 8.3.2). In the HBV-model
11
the reference evapotranspiration data for catchment 4, in which Office du Niger is located, is substituted by the new reference evapotranspiration data. Crop factors are also added to the model, for the irrigated rice catchments a crop factor of 1.2 is applied and for the other catchments a crop factor of 1. In the soil vegetation module the potential evapotranspiration is calculated by multiplying the reference evapotranspiration by the crop factor. The creation crop factors is done when the catchments are separated, the code can be found in the appendix (8.3.1).
2.4.3. Irrigation
Data about the monthly average amount of water in m3/s that is diverted to the Office du Niger irrigation scheme in the period 1989-2009 (table 1) is used. As the size of Office du Niger in this period is approximately 83.000 ha this data can be converted to the amount of irrigation in mm/day. In the appendix (8.3.3) this conversion script is attached. In the HBV-model directly after the soil vegetation module an irrigation module is implemented. The module sets irrigation for all irrigated catchments equal to the irrigation for Office du Niger and for all other catchments irrigation is zero. The irrigation is subtracted from the discharge and added to the soil moisture. If the storage capacity of the soil is exceeded the soil moisture storage is filled to the maximum and the rest of the water is irrigated to the upper zone storage. If negative discharges are undesirable the module can check whether the amount of water
that was used for irrigation is not higher than the discharge of the river. If this is the case irrigation is reduced to the exact amount of the discharge. A new overview of the general working of the HBV-model is shown in figure 5. The code of the irrigation module can be found in the appendix (8.3.4).
2.5. Scenarios
Four scenarios were studied in this research to investigate the impact of irrigated rice at Office du Niger on the Discharge of the Niger entering the Inner Niger Delta. In all four scenarios no dams are present.
(S0) The baseline scenario is the natural situation in which the Office du Niger irrigation scheme is not present.
(S1) Scenario 1 is the current situation in which the Office du Niger covers an area of approximately 83.000 ha.
(S2) Scenario 2 is the Office du Niger with minimal future expansions to a 200.000 ha irrigation scheme in 2020.
(S3) Scenario 3 is the Office du Niger with maximal future expansions to a 300.000 ha irrigation scheme in 2020.
12
3. Results
The results consist of two sections; first the water balance of Office du Niger independent of its size will be examined and compared to the natural situation in which no irrigation is present (3.1). Second the effects of the different scenarios on the discharge of the Niger will be presented (3.2).
3.1. Water balance
The water balance of Office du Niger can be expressed in the following equation: Irr + Prec = Δstorages + evap + flow
This equation can be divided in three parts, the inputs (Irrigation and Precipitation), the storages (soil, upper zone and lower zone) and the outputs (evapotranspiration and outflow). The amounts of water in this subsection will be expressed in height (mm); therefore the amounts will be independent of the area that Office du Niger is covering. As a result values are the same for S1, S2 and S3, these values will be compared with the reference situation S0 in which no irrigation present.
3.1.1. Inputs
In Figure 6 the inputs are plotted over time. The combined water input for Office du Niger is highest from July till November; this corresponds with the main growing season for rice. In the end of the growing season and during the harvest (November and December) the water input rapidly declines as irrigation drops and the dry season starts. In the dry season irrigation remains relatively low as less rice is cultivated. The yearly total water input for Office du Niger is 3591 mm; only 300 mm of this amount is comprised of precipitation.
13
3.1.2. Storages
Irrigation water is stored in the soil moisture; if the soil moisture full capacity is reached the water is stored in the upper zone. In figure 8 the storage of water in soil, upper zone and lower zone is plotted for Office du Niger and the natural/reference situation in which there is no irrigation. For Office du Niger the soil moisture storage becomes completely filled around July as the irrigation is increasing. During the wet season from July to October the rainfall and high irrigation cause the upper zone volume to fill. In October the rainfall stops and the irrigation starts to decrease. This causes the upper zone storage to become empty in the end of November. At this point soil moisture is still at its maximum, in December however it starts to decline as irrigation reaches a minimum. In the dry season less water is irrigated, causing the soil moisture to reach a minimum in April. The lower zone storage is filled during the wet season as water percolates from the upper zone, during the dry season there is no water in the upper zone so there is no
percolation possible and the lower zone slowly empties again. In the reference situation the soil, upper and lower zone are completely empty during the dry season. In the wet season rainfall causes the soil moisture to fill up slightly. The upper zone stays empty, as the little amount of water that reaches it directly percolates to the lower zone.
3.1.3. Outputs
The two outputs are evapotranspiration and outflow. In figure 8 the actual and potential evapotranspiration for Office du Niger and the reference situation are plotted. The evapotranspiration depends on soil moisture storage, if the soil moisture storage is at 75% of its maximum capacity or more evapotranspiration is equal to the potential evapotranspiration. Therefore for Office du Niger potential evapotranspiration is reached the entire year except for a short period, from March till June, when soil moisture is lower than 75% of its maximum capacity. This period corresponding to the end of the dry season has the highest potential evapotranspiration as air humidity is low and temperatures are high. In the reference situation the potential evapotranspiration is lower; this is because the potential evapotranspiration is calculated by multiplying the reference evapotranspiration by the crop factor. As stated in the method the crop factor for Office du Niger is 1.2 and for the reference situation 1. The actual evapotranspiration in the
Figure 7: Soil, upper zone and lower zone storages for Office du Niger and the reference situation.
14 reference situation is near zero
in the dry period as soil moisture is exhausted. However, evapotranspiration increases in the wet season as more water is stored in the soil and becomes available for evapotranspiration. Total yearly evapotranspiration for Office du Niger is 2664 mm compared to 287 mm in the natural reference situation. The Outflow consists of two parts as well, namely the base flow and the quick flow. Base flow is fed by the lower zone storage and the quick flow is fed by the upper zone storage. In figure 9 the outflow for Office du Niger and the reference situation is visualized. For Office du Niger the quick flow is only present in the wet season when there is water in the upper zone, during this period the base flow steadily increases as water percolates to the lower zone. After the wet season the upper zone quickly empties causing the quick discharge to drop to zero. Subsequently the lower zone isn’t fed by the upper zone anymore and starts to empty as well. As a result the base flow slowly drops during the dry season. The total outflow shows a distinct peak during the wet season reaching approximately 7mm/day and plummets to almost zero in the end of the dry season. For the reference situation only the total flow is plotted as Quick flow is not
present. The reason for the absence of the quick flow is that the upper zone is empty throughout the entire year. As the Total flow is the sum of the quick and base flow the base flow is exactly the same
Figure 9: Evapotranspiration for Office du Niger and the reference situation.
15
as the total flow and therefore disregarded as well. The total flow in the natural (reference) situation can almost be neglected; some base flow occurs in the wet season when a small amount of water reaches the lower zone. The total yearly outflow for Office du Niger is 926 mm, without irrigation (reference) this is 12 mm.
3.1.4. Balance
The water balance resulting from all the discussed in- and output is illustrated in figure 10. The storages were left out, as the model was run for several years with the same data. The storages reach an equilibrium in which they have the same value in the start and the end of the year. The water balance on the right is the reference balance subtracted from the Office du Niger balance. This balance therefore shows the change in the water balance when land-use changes from natural vegetation to irrigated rice. For Office du Niger the extra input in the form of irrigation
evaporates for 72% and the other 28% flows back into the Niger.
Figure 10: Water Balance for Office du Niger, reference situation and the change between these them
16
3.2. Scenarios
This section consists of two parts, first the amount of water that is diverted for irrigation is evaluated and second the effects on the discharge of the Niger are shown.
3.2.1. Irrigation
In figure 11 the amount of water that is diverted by Office du Niger is shown for the different scenarios. As the water diversion is linearly related to the irrigated area the shape of the three irrigation curves are the same. A clear peak is visible in September and October, corresponding with the middle of the main growing season for rice. In the current situation (S1) this peak irrigation is around the 130m3/s. With minimal expansion (S2) this is increases to 310 m/s and with maximal expansion to 470m3/s. After harvesting the rice in the beginning of the dry season the amount of water diverted for irrigation drops to respectively 60m3/s (S1), 130m3/s (S2) and 200m3/s (S3).
3.2.2. Discharge
The Discharge of the Niger entering the Inner Niger Delta in al scenarios is plotted in figure 12. In natural situation (S0), which was taken as baseline scenario, no dams are present and the discharge is extremely low in the dry season. The minimum is in May when the discharge becomes 25m3/s. Because the S0 discharge is so low in the dry season the discharge drops
Figure 11: The amount of water diverted for irrigation at Office du Niger
17
below zero in all three irrigation scenarios. Indicating that the amount of water used for irrigation is higher than the discharge of the river. The average discharge of the Niger in the four scenarios are respectively 738m3/s (S0), 678m3/s (S1), 595m3/s (S2) and 523m3/s (S3).
Figure 13 shows the net discharge reduction due to irrigation at Office the Niger. This is the discharge entering the Inner Niger Delta in S1, S2 and S3 subtracted from the discharge entering the delta in S0. Most discharge is lost in the end of the dry season, at this point in time irrigation is already relatively high and outflow is at its lowest. In the wet season the large amount of precipitation resulting in a high quick flow and increased base flow partly compensate for the water that is
consumed for irrigation. In August the net discharge loss increases again as more water is consumed for irrigation, the peak halfway October corresponds with the irrigation peak. Subsequently irrigation drops together with the discharge loss and another minimum is reached in December. From that point till the start of the next wet season the discharge loss increases again as irrigation increases en base flow decreases. The average discharge reduction is 59m3/s for S1, 143m3/s for S2 and 215m3/s for S3.
18
4. Discussion
This section consists of four parts. First the water balance will be evaluated (4.1); second the effects of the scenarios on the discharge of the Niger will be discussed and interpreted (4.2). At last the input data (4.3) and the model itself (4.4) will be discussed.
4.1. Water balance
The inputs and outputs of the water balance (section 3.1) for Office du Niger and the reference situation are approximately equal and the change in storages is near zero. This indicates that the water balance is correct, meaning that no water is created or destroyed in the model. Of the water that is diverted by Office du Niger 72% evaporates, the other 28% flows back into the Niger. Total evapotranspiration during the main growing season is roughly 1100 mm. Zwart & Leclert (2010) stated that evapotranspiration in this period is between 700-900 mm. The potential evapotranspiration rates on the other hand are consistent with the literature. Since rice at Office du Niger grows under a layer of water it is likely that actual evapotranspiration is equal to potential evapotranspiration as is the case in the simulation. The difference in evapotranspiration therefore indicates disagreement in literature rather than an error in the model. The loss of more than two third of the irrigation water due to evapotranspiration causes a reduction of the discharge of the Niger, the magnitude of this reduction depends on the size of Office du Niger. In the next paragraph the discharge reduction of the 4 modeled scenarios will be discussed.
4.2. Scenarios
The baseline scenario (S0) simulates the natural flow of the Niger in which Office du Niger is completely absent. The average discharge of the Niger entering the Inner Niger Delta in S0 is 738m3/s and the peak discharge is 2.981m3/s. As calibration was done on the basis of this scenario these values correspond with the literature. In the current situation (S1) in which Office du Niger covers an area of 83.000 ha on average 87m3/s is used for irrigation, this reduces the average discharge with 59m3/s. However, the reduction of the peak discharge is only 41m3/s. During this period the combination of a relatively low potential evapotranspiration and a high water input of precipitation and irrigation, causes large amounts of water to flow back into the Niger via groundwater and overland flow. In S0 almost no water flows back into the Niger as practically all precipitation evaporates. As relatively so much water flows back into the Niger in S1 and not in S0 the reduction of the peak discharge is lower than the average discharge reduction. This also applies for S2 and S3, in which Office du Niger covers an area of respectively 200.000 and 300.000 ha. In S2 on average 208m3/s is diverted for irrigation, leading to an average discharge reduction of 143m3/s and peak
discharge reduction of 102m3/s. In S3 on average 313m3/s is diverted for irrigation, leading to an average discharge reduction of 215m3/s and peak discharge reduction of 158m3/s. These scenarios indicate that the peak discharge will be reduced by 102-158m3/s by the year 2020 due to water consumption by Office du Niger. For S2 and S3 the amount of water used for irrigation was linearly related to the size of the irrigated area, by doing this it is assumed that the water efficiency stays approximately the same between the current situation and 2020. The resulting amounts of water intake by Office du Niger in both S2 and S3 are similar to values Zwarts (2012) expects. In S1, S2 and S3 the discharge of the Niger becomes negative during the dry season, caused by the fact that no dams were incorporated in the model. Currently there are two dams present in the upper Niger basin, the Selingue dam upstream and the Markala dam next to Office du Niger. The construction of another dam (Fomi dam) is already planned; this dam is implemented in the HBV-model by Wijers
19
(2014). Further research could be done by incorporating the Selingue and Markala dams in the HBV-model and combining this with the Office du Niger and Fomi extensions. In this way the entire human impact on the discharge of the Niger can be simulated for the current situation and future scenarios.
4.3. Input data
Some critical points must be mentioned concerning the input data of the HBV-model. For this research only data of the year 2010 was used and duplicated to simulate multiple years. For example the total yearly input of precipitation in the catchment in which Office du Niger is located is 300 mm. On average yearly total precipitation is approximately 600 mm, so for the Segou region in which Office du Niger is located 2010 was a relatively dry year. This does not have major impact as 2010 was a quite normal year concerning precipitation over the entire research area. For this research the parameters were kept constant over the different catchments, in reality there is large variability in parameters over space. However, this is compensated by the process of calibration in which the actual discharge is approximated. Finally the reference evapotranspiration data was of very low quality and therefore it was decided to use the recalculated reference evapotranspiration for Office du Niger for the entire Upper Niger Basin. Further improvement to this data could be the use of different meteorological station and interpolation over the entire basin. The lack of precision in reference evapotranspiration data is not a very disastrous as potential evapotranspiration rates are rarely reached, except for the Office du Niger catchment for which the data is very precise.
4.4. Model
Hydrological models, such as the HBV-model, are very useful tools for the prediction of river discharges, however they remain a simplification of the reality and therefore results should be carefully interpreted. The HBV-model has been successfully applied in basins all over the world (Liden & Harlin, 2000 & Buiteveld et al., 2005 & Jia & Sun, 2012). Some of the limitations of the HBV-model that are relevant concerning the upper Niger basin will be discussed. First the model uses saturation excess overland flow instead of Hortonian overland flow, or infiltration excess overland flow. Saturation excess overland flow occurs when the soil becomes saturated, Hortonian flow occurs when precipitation exceeds infiltration. Hortonian flow is the main contributor to overland flow in arid regions, especially when intense showers occur as is the case in the upper Niger basin. Not using Hortonian flow, therefore leads to an underestimation of quick flow. Second the model does not incorporate lakes, floodplains or marshes. The northern part of our research area is located in the Inner Niger Delta and therefore entirely consists of lakes, floodplains and marshes. Large quantities of water evaporate from these areas, which is not included in the model. As the Upper Niger Delta is located downstream of Office du Niger this is not a problem for discharges at ON. Including the delta in the model could also be an area of further research.
20
5. Conclusion
Currently Office du Niger covers an area of approximately 83.000 ha and the amount of water diverted for irrigation purposes is on average 87m3/s. 72% of this irrigation water that is diverted by ON is lost due to evapotranspiration, 28% finds its way back in the Niger again via groundwater or overland flow. As a result, the average discharge of the Niger is reduced by 59m3/s in the current situation due to irrigation at Office du Niger. It is intended to expand the irrigated area to roughly 200.000 ha or even nearly 300.000 ha by the year 2020, to feed the growing Malian population or possibly to export rice to China. In these scenarios the amount of diverted water will increase to respectively 208m3/s and 313m3/s on average, resulting in an average discharge reduction of 143-215m3/s with respect to the natural situation in which no irrigation is present. Currently, the natural peak discharge is reduced with only 41m3/s, increasing to 102-158m3/s in 2020 depending on the extent of the expansions. The absolute discharge reduction is the largest in the end of the dry season, when discharge of the Niger is already at its lowest. Without the construction of dams, no water would be available to feed the irrigation scheme during this period. The total reduction of the peak discharge is therefore much higher as large quantities of water are stored in reservoirs during the wet season and used for irrigation in the dry season. The peak discharge determines the extent of the inundated area in the Inner Niger Delta. These inundated areas are vital for the inhabitants of the delta as it enables cultivation of non-irrigated rice and bourgou and stimulates fishery. On top of that, the flooding of the delta draws millions of migrating birds each year and therefore is of high ecological importance. The reduction of peak discharge could have serious consequences for the inhabitants and ecology of the delta.
21
6. Acknowledgements
Special thanks to my supervisor Willem Bouten for the guidance and numerous feedback moments throughout the bachelor project. Furthermore I would like to thank Martijn Kamps for his input data and help with the data preparation. At last I would like to thank Sofie te Wierik and
22
7. Literature
Buitelveld, H., Wilke, K. & Krahe, P., 2005. Hydrological Modelling in the River Rhine Basin Part III - Daily HBV Model for the Rhine Basin. bundesanstalt für gewässerkunde.
Climate-Charts.com, 2010. Segou, Mali: Climate, Global Warming, and Daylight Charts and Data. [Online] World Climate Available at: http://www.climate-charts.com/Locations/m/M161272.php
[Accessed 17 June 2014].
Hong Kong Observatory, 2003. Climatological Information for Segou, Mali. [Online] Available at:
http://www.hko.gov.hk/wxinfo/climat/world/eng/africa/w_afr/segou_e.htm [Accessed 17 June 2014].
Jia, Q.Y. & Sun, F.H., 2012. Modeling and forecasting process using the HBV model in Liao river delta. Procedia environmental sciences, 13, pp.122-28.
Kamps, M., 2013. The influence of climate change and population growth on the discharge of the Nile and water distribution in the Nile Basin. Amsterdam.
Kuper, M., Mullon, C., Poncet, Y. & Benga, E., 2003. Integrated modelling of the ecosystem of the Niger river inland delta in Mali. Ecological Modelling, 164(1), pp.83-102.
Liden, R. & Harlin, J., 2000. Analysis of conceptual rainfall–runoff modeling performance in different climates. elsevier, 238(3-4), pp.231-47.
Liersch, S. et al., 2013. Vulnerability of rice production in the Inner Niger Delta to water resources management under climate variability and change. Environmetnal Science & Policy, 34, pp.18-33. Raes, D., 2009. ETo Calculator. [Online] FAO Available at:
http://www.fao.org/nr/water/docs/referencemanualeto.pdf [Accessed 15 June 2014].
Raes, D., Sy, B. & Feyen, J., 1995. Water use in rice schemes in the Senegal river Delta and Valley. Irrigation and Drainage Systems, 9(2), pp.117-28.
U.S. Geological Survey, 2013. Global Data Explorer. [Online] Available at:
http://gdex.cr.usgs.gov/gdex/ [Accessed 20 May 2014].
van Dijk, G.J., 2008. Potential water conflicts in Mali: identification of competing claims for surface water and the impact of climate change in the West Niger basin.
Vandersypen, K. et al., 2006. Irrigation performance at tertiary level in the rice schemes of the Office du Niger (Mali): Adequate water delivery through over-supply. Agricultural water management, 83(1), pp.144-52.
Wijers, B., 2014. Controlling the river; Simulating hydrological effects of the future Fomi dam in the Upper Niger river basin. University of Amsterdam.
Zwart, S.J. & Leclert, L.M., 2010. A remote sensing-based irrigation performance assessment: a case study of the Office du Niger in Mali. Irrigation science, 28(5), pp.371-85.
23
Zwarts, L., 2010. Will the Inner Niger Delta shrivel up due to climate change and water use upstream. Altenburg & Wymenga ecologisch onderzoek.
Zwarts, L., 2012. Water crisis in the Inner Niger Delta (Mali) Causes, consequences, solutions. Altenburg & Wymenga ecologisch onderzoek.
Zwarts, L., Beukering, P.V., Kone, B. & Wymenga, E., 2005. The Niger, a lifeline: effective water management in the Upper Niger Basin. Lelystad: RIZA-Rijkswaterstaat.
Zwarts, L., Bijlsma, R.G., van der Kamp, J. & Wymenga, E., 2009. Living on the edge. Wetlands and birds in a changing Sahel. Utrecht: KNNV uigerverij.
24
8. Appendix
8.1. Data preparation
This section contains detailed information and scripts that were used to prepare the input data for the HBV-model.
8.1.1. ArcGIS
This section contains a detailed description of the steps and tools (table 2) needed to create a catchment grid and corresponding catchment information about forestation, area size and connectivity.
A 3 arc second digital elevation model (DEM) of the upper Niger catchment was retrieved from … In ArcGIS the DEM was used to create a catchment grid. First the sinks of the DEM are filled, to make sure that all water eventually reaches the edge of the DEM. Subsequently a flow direction grid is made using the filled DEM; each cell in this grid contains the direction of the neighbouring cell with the lowest altitude. Next a flow accumulation grid is created from the flow direction grid, containing the number of cells that flow into each cell. With the flow accumulation grid streams can be defined, if a value is higher than a certain threshold the cell becomes a stream, otherwise not. Out of the resulting steam grid a stream link grid is made, in this grid each of the stream segments, linking two stream junctions is given a unique identifier. The stream link grid together with the flow direction can be used to make the catchment grid. Cells in this grid belong to the same catchment if they flow into the same stream segment. The threshold chosen by the stream definition therefore eventually determines the number of catchments. From the catchment grid and the stream link grid catchment polygons and streamlines are created. At this point the exact research area could be defined. The confluence of the Niger and Bani rivers was chosen as the most downstream catchment. All
streamlines and corresponding catchments that do not flow into this catchment are deleted. To link the data the catchment polygons and the streamlines the attribute tables of the two shape files are joined on their Grid-ID’s. An extra field is added to the attribute table indicating the next
downstream streamline and thus catchment as the attribute tables were joined. Another field is added to the attribute table in which the area of the catchment polygons is calculated in square meters instead of square units of the coordinate system. Finally one more field is added by hand containing a Boolean indicated whether a catchment is forested or not. Forestation was determined by laying the catchments over the Bing satellite base map. The attribute table containing
connectivity, area and forestation for all catchments was exported as a text file and further processed in Matlab.
The created catchment grid is needed to lump the meteorological forcing and initial state data. However, before lumping is possible the resolution and the size of the catchment grid and the data must be exactly the same. First the catchment grid is rescaled to the same resolution. Second, the forcing and initial state data of entire Africa must be clipped to the same size as the catchment grid. One grid for entire Africa map is geo referenced in ArcGIS and the catchment grid is laid over it. By doing this it can be checked whether the resolutions are the same. On top of that the row and column numbers of the Africa grid which compromise the catchment grid and thus research area can be retrieved. These row and column numbers are used for clipping in Matlab. Appendix … contains a table with all the used ArcGIS tools
25
8.1.2. Conversion ArcGIS to Matlab
The scripts below were used to convert the data retrieved in ArcGIS to .mat files, to create midpoints, to change the ArcGIS ID’s and to select the catchments within the research area. First the main script is shown and then the three used function files.
%% Title: creation HBV input
% Function: convert connectivity, area & forestation from ArcGIS to % matlab, calculate midpoints, create new ID's ascending % from 1, throw away catchments outside research area % Author: Vermeer, M. % Date: 22-5-2014 % Matlab version: R2012b %% Initialisation % Cleaning environment clc clear all Tool Description
Fill Sinks (Arc Hydro) Fill sinks for an entire DEM (grid).
Flow Direction (Arc Hydro) Create flow direction grid from a DEM grid.
Flow Accumulation (Arc Hydro) Create flow accumulation grid from a flow direction grid.
Stream Definition (Arc Hydro)
Create a new grid (stream grid) with cells from a flow accumulation grid that exceed used-defined
threshold. Stream Segmentation (Arc Hydro)
Create a stream link grid from the stream grid (every link between two stream junction gets a unique identifier).
Catchment Grid Delineation (Arc Hydro)
Create a catchment grid for segments in the stream link grid. It identifies areas draining into each stream link.
Catchment Polygon Processing (Arc Hydro) Create catchment polygons out of the catchment grid.
Drainage Line Processing (Arc Hydro) Create streamlines out of the stream link grid.
Find Next Downstream Line (Arc Hydro)
Find the HydroID of the next downstream linear feature class and store it in the NextDownID field of the feature. The directionality is based on the digitized direction. Connectivity is established by the physical connection of the linear features (does not require hydro network).
Calculate Areas (Spatial Statistics) Calculates area values for each feature in a polygon feature class.
Resample (Data Management) Alters the raster dataset by changing the cell size and resampling method.
26 close all
% loading
[subcatchments refinf bbox] = geotiffread('catchments61.tif'); [ID,gridID,hydroID,nextDownHydroID,area,isForest] =
textread('catchmentdata.txt','%d %*d %d %*d %*d %*d %*d %d %*d %d %*f %*f %f %d','delimiter',';','headerlines',1);
% variables & initialization
[row col] = size(subcatchments);
subcatchments = double(subcatchments); numberOfCatchments = 57;
%% Calculations
% Subcatchments
% Clean subcatchments (select the used subCatchments and change from GridID to ID)
subcatchments = cleanSubcatchments(subcatchments, gridID);
% As subcatchment 1 is lost due to upscaling it will be added by hand
subcatchments(1,145) = 1; % Connectivity flowsInto = connectivity(hydroID,nextDownHydroID); % Area subArea = round(area/1000000)'; % Forest isForest = isForest'; % Mid points
[midPosX, midPosY] = midpoints(subcatchments);
%% Output
% saving
save('K:\Bachelor scriptie\Matlab\midpoints.mat', 'midPosX', 'midPosY') save('K:\Bachelor scriptie\Matlab\subcatchments.mat', 'subcatchments') save('K:\Bachelor scriptie\Matlab\connectivity.mat', 'flowsInto') save('K:\Bachelor scriptie\Matlab\area.mat', 'subArea')
save('K:\Bachelor scriptie\Matlab\forest.mat', 'isForest')
function [ subcatchments ] = cleanSubcatchments(subcatchments, GridID)
% cleanSubcatchments: Select the used subCatchments and change from ArcGIS % GridID to ID
% Author: Vermeer, M.
[row col] = size(subcatchments); ID = 1:length(GridID)';
for i = 1:row for j = 1:col
val = subcatchments(i,j); index = find (GridID==val);
27 if isempty(index) subcatchments(i,j) = NaN; else subcatchments(i,j) = ID(index); end end end end
function [ flowsInto ] = connectivity( hydroID,nextDownHydroID )
% connectivity: the in GIS retrieved hydroID and nextDownHydroID are used to
% calculate the standard flowsInto ID's % Author: Vermeer, M. numberOfCatchments = length(hydroID); NextDownID = zeros(numberOfCatchments,1); array = [hydroID,nextDownHydroID,NextDownID]; for i = 1:numberOfCatchments index = find(array(:,1)==array(i,2)); if index array(i,3) = index; end end flowsInto = array(:,3)'; end
function [ midPosX, midPosY ] = midpoints( subcatchments )
% midpoints: Calculate midpoints of the subcatchments % Author: Wijers, B. editted by Vermeer, M.
[row col] = size(subcatchments); y = 1:row;
x = 1:col;
[X Y] = meshgrid(x,y);
% calculate x y coordinates for centroid
for i = 1:numberOfCatchments
% calculate mean x / y coordinates
cY = mean(Y(subcatchments==i)); cX = mean(X(subcatchments==i));
%round to whole number
midPosY(i) = round(cY); midPosX(i) = round(cX); end
midPosY = midPosY'; midPosX = midPosX';
28 end
8.1.3. Clipping
The forcings and initial states data for entire Africa is clipped in two steps. First the data is clipped to a rectangular area that exactly fits the dimensions of the research area. Second the data is clipped to shape of the research area, meaning that all cells outside the research area are converted to NaNs. The two following scripts were used for clipping the data.
%% Title: Clipping 1
% Function: Clipping data to a rectangular area in which the % research area is located
% Author: Wijers, B. editted by Vermeer, M. % Date: 22-5-2014 % Matlab version: R2012b %% Initialisation % Cleaning environment close all clear all clc
% assuming the to be clipped files to be 1540x1411 % row and column numbers determined in ArcGIS
rows = 497; %row start
rowe = 614; %row end
cols = 162; %col start
cole = 315; %col end
% load data
% Gather information on variables to be calculated
% Creating struct with all variables inside loaded mat file
vars = load('Temperature_interp.mat');
%% calculations
%# get the fieldnames stored in the structure
flds = fieldnames(vars);
%Loop over all variables and clip only the research area (see %initialisation)
for i=1:length(flds) %Variable Name
VarName = [flds{i} '_clip']; %Clipped area of the variable
clip = vars.(flds{i})(rows:rowe,cols:cole); %Adding new data to struct evap with a new field
temp.(VarName) = clip; end
%% output
29
%% Title: Clipping 2
% Function: Clipping data to research area % Author: Wijers, B. editted by Vermeer, M. % Date: 22-5-2014 % Matlab version: R2012b %% Initialisation % Cleaning environment close all clear all clc %load data load('evap_clip.mat'); load('initial_africa_clip.mat'); load('rain_africa_clip.mat'); load('temp_africa_clip.mat');
%load research area
load('subcatchments.mat');
%create template
template = logical(subcatchments >=0); template = double(template);
template(template==0) = NaN;
%Gather name list
fldsevap = fieldnames(evap); fldstemp = fieldnames(temp); fldsrain = fieldnames(rain); fldsinits = fieldnames(inits); %% Calculations % Clip evaporation for i = 1:length(fldsevap)
VarName = [fldsevap{i} '_res'];
clip = evap.(fldsevap{i}) .* template; evap_res.(VarName) = clip;
end
% Clip temperature
for i = 1:length(fldstemp)
VarName = [fldstemp{i} '_res'];
clip = temp.(fldstemp{i}).* template; temp_res.(VarName) = clip;
end
% Clip precipitation
for i = 1:length(fldsrain)
VarName = [fldsrain{i} '_res'];
clip = rain.(fldsrain{i}).* template; rain_res.(VarName) = clip;
end
% Clip initial sates
30 VarName = [fldsinits{i} '_res'];
clip = inits.(fldsinits{i}).* template; inits_res.(VarName) = clip;
end
%% Output
save('J:\Matlab\data\temp_clip_res.mat','temp_res'); save('J:\Matlab\data\rain_clip_res.mat','rain_res'); save('J:\Matlab\data\evap_clip_res.mat','evap_res'); save('J:\Matlab\data\inits_clip_res.mat','inits_res');
8.1.4. Lumping
The following script was used to lump the forcings and initial states data.
%% Title: Lumping
% Function: Lumping forcings & initial states data % Author: Wijers, B. & Vermeer, M.
% Date: 22-5-2014 % Matlab version: R2012b %% Initialisation % Cleaning environment clear all close all clc % Load data load('evap_res_fix.mat'); load('rain_res_fix.mat'); load('temp_clip_res.mat'); load('inits_clip_res.mat'); load('subcatchments.mat'); % Gather info fldsevap = fieldnames(evap_res_fix); fldstemp = fieldnames(temp_res); fldsrain = fieldnames(rain_res_fix); fldsinits = fieldnames(inits_res);
% Remove nans for unique
subcatchment2 = subcatchments; subcatchment2(isnan(subcatchment2)) = []; catchments = unique(subcatchment2); %% Calculations % Lump precipitation for j = 1:length(fldsrain); currgrid = rain_res_fix.(fldsrain{j}); for i = 1:length(catchments)
currcatch = double(logical(subcatchments == catchments(i))); currcatch(currcatch == 0) = NaN;
31 select = currcatch .* currgrid;
lumpedval = (nanmean(select(:))) * 1000; obsPrec(j,i) = lumpedval; end end % Lump temparature for j = 1:length(fldstemp); currgrid = temp_res.(fldstemp{j}); for i = 1:length(catchments)
currcatch = double(logical(subcatchments == catchments(i))); currcatch(currcatch == 0) = NaN;
select = currcatch .* currgrid; lumpedval = nanmean(select(:)); obsTemp(j,i) = lumpedval; end end % Lump evaporation for j = 1:length(fldsevap); currgrid = evap_res_fix.(fldsevap{j}); for i = 1:length(catchments)
currcatch = double(logical(subcatchments == catchments(i))); currcatch(currcatch == 0) = NaN;
select = currcatch .* currgrid;
lumpedval = (nanmean(select(:))) * -1000; obsEref(j,i) = lumpedval; end end % Lump interception currgrid = inits_res.(fldsinits{1}); for i = 1:length(catchments)
currcatch = double(logical(subcatchments == catchments(i))); currcatch(currcatch == 0) = NaN;
select = currcatch .* currgrid; lumpedval = nanmean(select(:)); sInterceptInit(i) = lumpedval; end
% Lump soil moisture
currgrid = inits_res.(fldsinits{3}); for i = 1:length(catchments)
currcatch = double(logical(subcatchments == catchments(i))); currcatch(currcatch == 0) = NaN;
select = currcatch .* currgrid; lumpedval = nanmean(select(:)); sSoilInit(i) = lumpedval;
end
% Lump upper zone
32 for i = 1:length(catchments)
currcatch = double(logical(subcatchments == catchments(i))); currcatch(currcatch == 0) = NaN;
select = currcatch .* currgrid; lumpedval = nanmean(select(:)); sUpperZoneInit(i) = lumpedval; end
% Lump lower zone
currgrid = inits_res.(fldsinits{5}); for i = 1:length(catchments)
currcatch = double(logical(subcatchments == catchments(i))); currcatch(currcatch == 0) = NaN;
select = currcatch .* currgrid; lumpedval = nanmean(select(:)); sLowerZoneInit(i) = lumpedval; end
% Create empty initial snow
numberOfCatchments = 57; % Create timeseries obsErefTime = [0:364]'; obsPrecTime = [0:364]'; obsTempTime = [0:364]'; sSnowLiquidInit = zeros(1,numberOfCatchments); sSnowSolidInit = zeros(1,numberOfCatchments); %% Output
save('K:\Bachelor scriptie\Matlab\initial-state.mat', 'sInterceptInit', 'sLowerZoneInit', 'sSnowLiquidInit', ...
'sSnowSolidInit', 'sUpperZoneInit', 'sSoilInit');
save('K:\Bachelor scriptie\Matlab\precipitation.mat', 'obsPrec', 'obsPrecTime');
save('K:\Bachelor scriptie\Matlab\reference-evaporation.mat', 'obsEref', 'obsErefTime');
save('K:\Bachelor scriptie\Matlab\temperature.mat', 'obsTemp', 'obsTempTime');
8.2. Calibration
In the table 3 the final values of the parameters are stored. The final values are a multiplication of the original value, the averaged value from the HBV-Rhine model. The pMAXBAS value is the value for the largest catchment.
parameter value pMAXBAS 11 pPERC pPERC * 3.5 pK4 pK4 * 0.84 pKHQ pKHQ * 0.5 pBETA pBETA * 1 pFC pFC * 1.4
33
8.3. Irrigated rice implementation
This section contains the scripts used to implement the irrigated rice.
8.3.1. Create Irrigated Catchments
The following script was used for creating separate catchments for the irrigated areas and for the addition of the crop factors.
%% Title: Initialize irrigation
% Function: creating sperate catchments for irrigated areas % Author: Vermeer, M.
% Date: 17-6-2014 % Matlab version: R2012b
irrigatedArea = zeros(1,nSubs);
%irrigatedArea(3) = 300; % specialize which subcatchment is irrigated and how large that area is
irrigatedArea(4) = 830; % specialize which subcatchment is irrigated and how large that area is
irrSubs = find(irrigatedArea ~= 0); nIrrSubs = length(irrSubs); % forcings obsPrec = [obsPrec,obsPrec(:,irrSubs)]; obsTemp = [obsTemp,obsTemp(:,irrSubs)]; obsEref = [obsEref,obsEref(:,irrSubs)]; % initial state sSnowSolidInit = [sSnowSolidInit,sSnowSolidInit(irrSubs)]; sSnowLiquidInit = [sSnowLiquidInit,sSnowLiquidInit(irrSubs)]; sInterceptInit = [sInterceptInit,sInterceptInit(irrSubs)]; sSoilInit = [sSoilInit,sSoilInit(irrSubs)]; sUpperZoneInit = [sUpperZoneInit,sUpperZoneInit(irrSubs)]; sLowerZoneInit = [sLowerZoneInit,sLowerZoneInit(irrSubs)]; % parameters pALPHA = [pALPHA,pALPHA(irrSubs)]; pBETA = [pBETA,pBETA(irrSubs)]; pCEVPFO = [pCEVPFO,pCEVPFO(irrSubs)]; pCFLUX = [pCFLUX,pCFLUX(irrSubs)]; pCFMAX = [pCFMAX,pCFMAX(irrSubs)]; pCFR = [pCFR,pCFR(irrSubs)]; pEPF = [pEPF,pEPF(irrSubs)]; pFC = [pFC,pFC(irrSubs)]; pFOCFMAX = [pFOCFMAX,pFOCFMAX(irrSubs)]; pHQ = [pHQ,pHQ(irrSubs)]; pICFI = [pICFI,pICFI(irrSubs)]; pICFO = [pICFO,pICFO(irrSubs)]; pK4 = [pK4,pK4(irrSubs)]; pKHQ = [pKHQ,pKHQ(irrSubs)]; pLP = [pLP,pLP(irrSubs)]; pMAXBAS = [pMAXBAS,pMAXBAS(irrSubs)]; pPERC = [pPERC,pPERC(irrSubs)]; pRFCF = [pRFCF,pRFCF(irrSubs)]; pSFCF = [pSFCF,pSFCF(irrSubs)]; pTT = [pTT,pTT(irrSubs)]; pTTI = [pTTI,pTTI(irrSubs)]; pWHC = [pWHC,pWHC(irrSubs)];
34
% no interception in irrigated catchments
for i = nSubs+1:nSubs+length(irrSubs) pICFI(i) = 0;
end
% crop factor or crop coefficient
pKC = ones(1,nSubs+length(irrSubs)); for i = nSubs+1:nSubs+length(irrSubs) pKC(i) = 1.2; end % forestation isForest = [isForest,zeros(1,nIrrSubs)]; % connectivity flowsInto = [flowsInto,flowsInto(irrSubs)]; for i = 1:length(irrSubs) flowsInto(irrSubs(i)) = nSubs + i; end
% activate if you want the irrigated subcatchment before the non-irrigated % catchment in the connectivity
% for i = 1:length(irrSubs)
% indexUpSubs = find(flowsInto == irrSubs(i)); % flowsInto(indexUpSubs) = nSubs + i;
% end
% flowsInto = [flowsInto,irrSubs];
% sub area
subArea = subArea - irrigatedArea;
subArea = [subArea,irrigatedArea(irrSubs)]; % midpoints midPosX = [midPosX;midPosX(irrSubs)]; midPosY = [midPosY;midPosY(irrSubs)]; % subcatchments [row,col] = size(subcatchments); for i = 1:row for j = 1:col
index = find(irrSubs == subcatchments(i,j)); if ~isempty(index)
if round(rand(1)) == true %if midPosX(nSubs + index) < j
subcatchments(i,j) = nSubs + index; end
end end end
isIrrigated = [zeros(1,nSubs),ones(1,nIrrSubs)]; nSubs = nSubs + nIrrSubs;
35
8.3.2. Interpolation evapotranspiration
The following script was used to interpolate the mean monthly reference evapotranspiration to daily data.
%% Title: Interpolation reference evapotranspiration
% Function: Interpolating mean monthly reference evapotranspiration % data to daily data
% Author: Vermeer, M. % Date: 10-6-2014 % Matlab version: R2012b %% Initialisation % cleaning environment clear all close all clc %% Calculations
% monthly reference evapotranspiration Office du Niger [mm/day]
monthlyRefEvapON = [5.5 6.4 7.1 7.3 7.0 6.2 5.2 4.7 5.0 5.7 5.7 5.2]'; % linear interpolation xMonthly = [-14 16 45 75 106 136 167 197 228 259 289 320 350 381]; yMonthly = [monthlyRefEvapON(12);monthlyRefEvapON;monthlyRefEvapON(1)]'; xDaily = 1:365; yDaily = interp1(xMonthly,yMonthly,xDaily); refEvapON = yDaily'; %% Output
save('E:\Bachelor scriptie\Matlab\refEvapON.mat', 'refEvapON')
8.3.3. Conversion and interpolation irrigation
The following script was used to convert the irrigation data from m3/s to mm/day and, subsequently, to interpolate the mean monthly reference evapotranspiration to daily data.
%% Title: Interpolation irrigation
% Function: Interpolating mean monthly irrigation data to daily % data % Author: Vermeer, M. % Date: 10-6-2014 % Matlab version: R2012b %% Initialisation % cleaning environment clear all close all clc %% Calculations
36
% monthly irrigation rate in m3/s for area of 83.000 ha
monthlyIrrigation =
[58.3,62.5,71.1,74.5,87.0,95.8,95.7,97.1,121.4,131.4,88.7,54.6]'; irrigatedArea = 83000;
% m3/s to mm/day
monthlyIrrigation = monthlyIrrigation * 60 * 60 * 24 * 1000 / irrigatedArea / 10000; % linear interpolation xMonthly = [-14 16 45 75 106 136 167 197 228 259 289 320 350 381]; yMonthly = [monthlyIrrigation(12);monthlyIrrigation;monthlyIrrigation(1)]'; xDaily = 1:365; yDaily = interp1(xMonthly,yMonthly,xDaily) irrigation = yDaily'; %% Output
save('E:\Bachelor scriptie\Matlab\irrigation.mat', 'irrigation')
8.3.4. Irrigation Module
The following script is the irrigation module. Note: for this research negative discharge was allowed.
%% Title: Irrigation Module
% Function: Irrigating water from discharge to soil % Author: Vermeer, M. % Date: 15-6-2014 % Matlab version: R2012b function [TotalFlow1,SoilMoisture8,UpperZoneVolume9,Irrigation1]... = hbvIrrigationModule(TotalFlow,SoilMoisture6,UpperZoneVolume8,... Irrigation,isIrrigated,fDischargeTotalLagged,lagWeights,subArea,... flowsInto,pFC,timeStep);
Irrigation1 = Irrigation * isIrrigated;
TotalFlow1 = TotalFlow - Irrigation1 * timeStep;
SoilMoisture7 = SoilMoisture6 + Irrigation1 * timeStep;
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % deactivite to allow negative discharge
% calculate total discharge
[fDischargeTotalDelayed,fDischargeTotalLagged] =...
hbvTriangularFilter(fDischargeTotalLagged,TotalFlow1,lagWeights); [fDischarge3TotalDelayed,fDischarge3TotalDelayedAcc] =...
hbvAccFlow(fDischargeTotalDelayed*1e-3/86400,subArea*1e6,flowsInto);
% discharge back to mm/day
fDischarge3TotalDelayedAcc =...
37
CannotBeIrrigated = abs(min(0,fDischarge3TotalDelayedAcc));
TotalFlow1 = TotalFlow1 + CannotBeIrrigated;
SoilMoisture7 = SoilMoisture7 - CannotBeIrrigated; Irrigation1 = Irrigation1 - CannotBeIrrigated;
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
ToUpperZone = (SoilMoisture7>pFC).*(SoilMoisture7-pFC) / timeStep; SoilMoisture8 = SoilMoisture7 - ToUpperZone * timeStep;
UpperZoneVolume9 = UpperZoneVolume8 + ToUpperZone * timeStep;
8.4. Year duplication
The following script was used to let the HBV-model run for multiple years with the same input data.
%% Title: Year duplication
% Function: copy the data of one year to simultate multiple years % with the same data
% Author: Wijers, B. editted by Vermeer, M. % Date: 12-6-2014 % Matlab version: R2012b function [obsErefTime,obsPrecTime,obsTempTime,obsEref,obsPrec,obsTemp,irrigation,pot EvapON,plotEnd,timeEnd]... = yeardupe(nryear,timeEnd,obsEref,obsPrec,obsTemp,irrigation,potEvapON,plotEn d)
% Creating temporarily variables.
obsEreftemp = obsEref; obsPrectemp = obsPrec; obsTemptemp = obsTemp;
tempIrrigation = irrigation; tempPotEvapON = potEvapON;
% Dependant on nryear duplication of data
for i = 1:nryear-1;
obsEreftemp = [obsEreftemp; obsEref]; obsPrectemp = [obsPrectemp; obsPrec]; obsTemptemp = [obsTemptemp; obsTemp];
tempIrrigation = [tempIrrigation; irrigation]; tempPotEvapON = [tempPotEvapON; potEvapON]; end
% Final update of the Forcing
obsEref = obsEreftemp; obsPrec = obsPrectemp; obsTemp = obsTemptemp;
irrigation = tempIrrigation; potEvapON = tempPotEvapON;
38
% Forcing time adjustment
timeEnd = (timeEnd *nryear)+(nryear-1);
% Plot year adjustment
plotEnd = plotEnd*nryear;
% Time adjustment
obsTempTime = [0:1:timeEnd]'; obsPrecTime = [0:1:timeEnd]'; obsErefTime = [0:1:timeEnd]';
