A hydrologic model in semi-arid urban areas with urban irrigation and water management

(1)

MSc Computational Science

Master Thesis

A hydrologic model in semi-arid urban

areas with urban irrigation and water

management

by

Maarten W.J. van den Ende

10412939

September 20, 2019

42EC

Daily Supervisor:

Dr. em. prof. Hilary McMillan

Examinator:

Dr. Valeria

Krzhizhanovskaya

Assessor:

Dr. Mike Lees

(2)

Abstract

In this research a hydrological model has been created for semi-arid urban areas with significant urban irrigation and water management. A framework around the extensive WRF-Hydro model was created to simulate urban irrigation. Land-cover, impervious surface satellite data, soil characteristics and temporal model output are used to daily calculate urban irrigation. Calibrating irrigation on spatial evapotranspiration data significantly increased the models performance during arid time periods. Calibration on streamflow is performed as well as a sensitivity analysis on two main model parameters and on urban irrigation. While urban irrigation does affect streamflow, the conclusion that urban irrigation significantly impacts flooding could not be made.

(3)

1 Introduction

In the last century urbanization has increased immensely all over the world. In 2008 more than half of the total population lived in cities, and this is expected to increase to more than 80 percent by 2030 [1]. These urban areas need to provide their inhabitants with adequate amounts of water and need to be able to dispose of all water used by them.

Unfortunately, more than half of all cities with more than one hundred thousand in-habitants are located in watersheds stressed for water [2]. In the United States, cities in California like Los Angeles and San Diego, as well as cities in Arizona like Phoenix require much more water than the local environment is able to provide. This requires them to im-port water from different areas, changing the local hydrological water balance. The planted vegetation in these cities is often not adapted to the local climate and needs significant watering. In addition to this these lawns are often over irrigated by up to 40 percent during summer [3]. This urban irrigation can greatly exceed the natural precipitation in these (semi-)arid areas and therefore significantly impact the hydrological water balance.

The predictions of the California Pacific Institute shows that the urban water use in California will increase by 1.5 million acre feet (MAF) to an incredible 5.8 MAF [4] and following current trends total water use could reach 12 MAF in 2030, as can be seen in Figure 1. To understand the effects of this increase of stress on the urban water balance, a better understanding of the effects of urban irrigation on the hydrological water balance is increasingly important.

Figure 1: California total water use. Courtesy of Pacific Water Institute [4].

(5)

In order to better understand the effects of urban irrigation in semi-arid urban areas, this research will answer the question:

”To what extent can the performance of WRF-Hydro be improved by the incorporation of urban irrigation in semi-arid urban areas and what are its effects on the streamflow of the main watershed river?”

To answer this, a framework is created that adds a model for irrigation onto the WRF-Hydro model. It uses multiple spatial data-sources and high resolution satellite imagery as well as feedback from the model itself to calculate the amount and spatial distribution of urban irrigation as accurately as possible.

(6)

2 Literature review and related work

Related work can be split into two subjects, the general study of urban influence on the hydrological cycle and land surface interactions, and the specific research on the effects of irrigation. The latter has significantly less research done in the case of urban irrigation.

2.1 Urban water balance and irrigation

Quite some studies on the urban water balance have been conducted, although almost all research ends with the conclusion that more and more in-depth research is required for better understanding of this complex system. Taking all irrigation into account, both urban and agricultural, Lawston shows that irrigation has an impact on soil temperature, decreasing it during day and increasing during night [5]. It also shows that irrigation repartitions surface sensible and latent heat fluxes. Focusing on the urban recharge of groundwater [6], this review shows that different approaches on irrigation, holistic or component wise, give different results and challenges. They show how complicated these systems are and emphasize that urban recharge is an under-researched topic.

Other papers researching the water and heat balance of urban areas focus on the surface urban energy [7] and urban heat island effect [8] [9], showing that urban irrigation counter-acts urban heat island effects. The effects of irrigation on cooling the soil temperature could even be masking global warming, as recently global irrigation has increased significantly and its effects on decreasing soil temperature are often not taken into account [10] [11].

Research on agricultural irrigation shows that irrigation-induced soil moisture memory lasts from a couple of weeks to a month [12]. It shows that parametrization of irrigation is important for climate studies and to improve understanding of the role of human activities on land-air interactions.

2.2 Urban irrigation

Less well researched is urban irrigation, both in its effects and in accuracy of its magni-tude. Almost all studies that do take urban irrigation into account use an overly simplified irrigation scheme, and most of them use an old or very simplified land surface model.

An example of such very simplified simulation of irrigation is the private model Aqua-cycle. While they do include urban wastewater and storm water, they simulate irrigation by forcing the top soil layers in the model to always be above a minimum value in urban area [13]. This value is set such that the soil moisture is always high enough to sustain plant growth. The earlier mentioned paper on surface urban energy [7] includes a basic irrigation scheme splitting into automatic irrigation and manual irrigation. The amount irrigated depends in their model on four site specific parameters as well as time since last irrigation event.

(7)

to further investigate hydrological processes in urban areas [14]. They focus on discovering best management practices in the Los Angeles area, while using the same CIMIS data as used in this research for validation to see how much irrigation is being done. They focused on storm runoff and its protection. Their irrigation methods is quite extensive and can be adjusted depending on the irrigation method. They calculate the amount of water needed by the urban vegetation as P W R = ET0KcropD with Kcropthe so called ’crop coefficient’, a

factor defined by the amount of water a plant need. This factor is used to go from reference evapotranspiration ET0 to the plant specific evapotranspiration ET . These terms will be

explained more in the Theory section. D is the time in between irrigation events. Then they either choose to run the irrigation event for the needed time to get the required amount of water in the soil, or they calculate the soil water deficit. This is done by looking at the soil moisture in the top soil layer and calculating the deficiency. The net needed water is defined by

P W Rnet= (P W R − SM C ∗ ds)Ks (1)

with SM C the soil moisture and ds the depth of the roots in the soil, which depends

on the vegetation in their model. Ks is called the stress factor but they do not go deeper

into its meaning or values. Their main result was that the soils physical properties was very important in the size of the storm runoff, and that according to their model, the amount of water used for urban irrigation could be decreased by 53% based on the plants water demands.

Research on the effects extreme cases of irrigation in the south-east of the United-States, where soil moisture of the top one meter soil layer is set to 100%, shows that while evap-oration is seen to be increased, decrease in soil temperature prevents a positive feedback loop [10]. This shows that there is a limit to the increase in evapotranspiration caused by irrigation, and that the atmosphere is able to stabilize fluxes on a large-scale.

A very recent paper on the effects of irrigation on urban environments uses the Common Land Model as well as the same USGS Land Cover satellite data, while splitting the urban area in three classes. Their focus is on land surface temperature, which does indeed decrease with added irrigation. They note that in their simulations, the timing of irrigation did not impact the hydrological fluxes much, but that added irrigation greatly increases the modelling capabilities of explaining the land surface temperature differences between high and low water use areas with similar land cover classification [15].

Their result that the timing of irrigation does not impact the surface fluxes much differs to the research of [16], [17], which results in that a cycle of three days between each irrigation event matches the diurnal fluctuations land surface temperature as well as on evapotran-spiration best. They use a single soil layer version of the Noah land surface model, which is the predecessor of the Noah-MP LSM used in this research. Their irrigation scheme is also the basis of the scheme used in this research. Their focus is on the Los Angeles area. They show that the impact of irrigation is important in modeling studies, especially in arid

(8)

regions. Accurate representation of irrigation impacts mainly long term simulations where the steady state cycles of soil moisture is significantly altered by irrigation.

In contrast to this, [18] came to the conclusion that while commercial irrigation impacts the anthropogenic moisture in southern California, outdoor domestic water use had no or no significant impact on the surface temperatures of the WRF-UCM model during most of the day. The forecasting skills of WRF-UCM did not improve significantly with this added moisture, making their conclusion that more research is needed on the distribution of anthropogenic moisture in the forecasting period.

The Methods section contains more theory and background information on related work, including a technical description of WRF-Hydro, based on [19]. Theory on important hy-drological terms used and research on behaviour in urban irrigation will be discussed there as well.

(9)

3 Case Study Area

As a case study the semi-arid county of San Diego has been chosen. The area includes the urban city of San Diego, as well as the coastline and deeper in land a mountainous region. The watershed of this region provides the main river in this area, the San Diego river. The characteristics of this San Diego river will be used as a reference for the accuracy of the predictions of the model. Figure 2 shows the outline of the region, as well as the model output streamflow locations and USGS streamflow gauges [20] at Mast Road and Fashion Valley.

Figure 2: The San Diego case study area, shown with black outline. The location of artificial streamflow gauges and their channel ID are shown in red dots and the two USGS streamflow gauges shown as black diamonds.

3.1 Meteorological environment

San Diego is known as the city with the perfect weather, with monthly average high tem-peratures ranging between 24.7 and 18.2 degrees Celsius. The average temperature is 17.55 degrees Celsius. During day the sun shines for approximately 70% of the time [21]. Rain-fall during summer is almost nonexistent, with some summers not having a single day of precipitation. Only during late fall and winter months an average of around 45mm of rain is recorded. This monthly precipitation generally falls within a period less than 7 days, indicating that rain often comes in huge storms. In the last century almost every year San Diego has dealt with one or more major floods [22]. Even though an extensive system of

(10)

storm drains is in place, the San Diego river is not always able to handle the discharge of the rain during these events. The San Diego river is infamous for its frequent flooding of Fashion Valley, accounting for significant economic damage. Therefore more knowledge on these flash floods and how to reduce them could help create policies combating them.

3.2 Vegetation

The San Diego climate makes its native vegetation exist mostly of shrubs, chaparral shrubs, some riparian woodlands and torrey pines woodland [23]. The San Diego county has more species of plants and animals groups than any other county in the United States [24], however this vegetation is not well suited for gardens so most of the urban vegetation exists of non-native vegetation. This vegetation needs additional irrigation to survive, especially during the dry summers. Although some residents choose palm trees, cacti and succulents that are well adjusted for this climate, many still have grassy lawns and deciduous trees which need significant irrigation to survive altering the latent heat flux. Even the municipal trees on the side of the roads need to be irrigated by homeowners close by. According to Jarvi et al. [7], residential land types relatively even demand a larger amount of water compared to farmland. In addition to this, even though it is sparse, water is still relatively cheap and urban residents are not knowledgeable on irrigation, many homeowners irrigate between 40 [3] and 50 [14] percent too much.

3.3 Water management

Increased urbanization and the ever increasing water usage, along with global warming and the fact that California depends mainly on depleting snow in the mountains to resupply surface water streams, appropriate water management policies are increasingly crucial for the future of San Diego. Even though its authors express hope and belief that sustainable water management in California is possible, its future is uncertain [4]. Shortage and political conflicts might be expected if California does not meet its environmental, economic and social water needs with new technology, strong management and the right incentives.

In order to model the streamflow of the San Diego river properly, data which included the reservoirs of San Vicente, El Capitan, Loveland reservoir and Sweetwater reservoir was needed and used in the later stages of the research. These reservoirs, located at the base of the mountain ranges, hold and keep most of the water that originates from rainfall in the mountains. Their main source however, is the Colorado river. These reservoirs are part of the ’Emergency Storage Project’ [25], which is able to provide the County of San Diego with drinking (and irrigation) water for a maximum of six months in emergencies by a pipeline system connecting multiple reservoirs and municipalities. Information on the policy of retaining water in these reservoirs is unfortunately unavailable. This means that modeling the hydrological cycle with these reservoirs comes with added difficulties.

(11)

4 Methodology and technical description

A significant amount of work in this research has gone into understanding and working with this WRF-Hydro model: editing input data to be compatible with the WRF-Hydro model, running the model itself and creating a framework to implement the irrigation model. This section starts with theory on WRF-Hydro and gives a short technical description on relevant subjects. It describes the WRF-Model, the physics and methods it uses to calculate the hydrological cycle. It also describes important parameter files it uses and the alterations made to them in this research. Additional precipitation satellite data used here is also described as well as the processing to this data done to make it work with the model.

This section then explains the methodology and theory used for creating the new frame-work, handling of the data, the irrigation model itself as well as the validation, calibration and sensitivity analysis methods.

4.1 WRF-Hydro model description

This research uses the extensive WRF-Hydro modeling system [26]. WRF-Hydro is an open-source community driven model and can be used for many hydrometeorological and hydrological studies. In this study it is used as a stand alone model, however it was originally designed to be coupled with the Weather Research and Forecasting model (WRF) [27].

WRF-Hydro uses the Noah-MP Land Surface Model [28], a multiparameterized land surface model. It incorporates many realistic processes, and parameters like a separated vegetation canopy top and bottom, leaves with prescribed dimensions and even the orien-tation of the vegeorien-tation. It uses two-stream radiation, multiple different soil depths and incorporates shading as well as many other physical processes with the end goal to model all land surface interactions and fluxes as accurately as possible.

WRF-Hydro has a few core features: it is multi-scale and can handle different spatial grids for different input data, it is modularized so its complexity can be adjusted and specific parts can be used, the code is written for the model to run as paralleled as possible and it is capable of handling many data formats as input thanks to pre-processing workflows. On top of this is the advantage that it is open-source, which means that it is free to use. Therefore the model is used for many studies so it is well known in the academic environment.

The model can be used on different scales, from global to regional. Here we use the San Diego county area, so the model will be used using grids and settings fitting for this scale.

The WRF-Hydro model needs to be compiled into an executable file using a FORTRAN compiler and a configuration file to determine some main characteristics of the model. We use the most extensive version, using the spatially variable soil characteristics capabilities of the Noah-MP module. Large file support is turned on as well as diagnostics in order to have all possible output available. Nudging,a way of arbitrarily forcing the streamflows of the model to match desired or measured values is turned off as we have clear channel routing input and want the model to correctly model the streamflows.

(12)

Figure 3: WRF-Hydro modules and their interactions. Courtesy of [19]. 4.1.1 Model Physics description

In this section the physics and workflow of WRF-Hydro is described, with explanations needed input and a focus put on important factors in urban modelling. This information is described in detail in the technical user guide of WRF-Hydro [19]. The coupling between the different modules and their outputs withing WRF-Hydro is shown in Figure 3.

First the 1-dimensional column land surface model is called. It calculates the vertical energy fluxes as well as the moisture and thermal states. It uses a 1-dimensional land surface model (LSM) grid file which can range in resolution from 30 meters to 4 kilometers. In this research a one kilometer resolution is used. Then the subsurface lateral flow is calculated. This happens before the overland flow to account for exfiltration water from fully saturated soil which needs to be added to the infiltration excess variable. After this the water table depth is defined using the top of the saturated soil layer. In this research we use the standard four soil layers. Then the overland flow module is called. This uses diffusive wave equations as well as Manning’s equation to calculate resistance of overland flow momentum. WRF-Hydro is capable of using different sized grids for the land surface model and the routing grid and uses disaggregation-aggregation methods to give the model the important ability to process overland and subsurface flow processes on a scale much smaller than the rest of the land surface model grid. In this research this factor is four, the LSM works on a 1-km grid, while the routing is done on a 250 meter scale grid.

The subsurface module is two dimensional, simulating flow to and from different soil layers in each cell, as well as lateral flow from one cell to another. The subsurface flow is

(13)

defined by Qnet(i,j)= hi,j X x γx(i,j)+ hi,j X y γy(i,j (2) γx(i,j)= −

wi,jKsati,jDi,j

ni,j βxi,j (3) hi,j= 1 − zi,j Di,j (4) with the flow from one cell to another in the x direction and from saturated soil layer to other layers in the y direction. Ksatis the hydraulic conductivity, n the local power law

exponent, which in the version of WRF used in this study is set to 1. D is the soil thickness, z the depth to the water table w the cell width. βi,j is the difference of the water table

depth between two adjacent cells.

The important parameters for the subsurface flow are the ones in

hydro2dtbl.TBL: saturated hydraulic conductivity, porosity and field capacity.

The overland flow in WRF-Hydro is performed using an explicit fully-unsteady finite difference wave formulation. It takes into account backwater effects and can handle flow on adverse slopes. It uses a simplification of the St.Venant equations of continuity and momentum in a shallow water wave. In our simulations we use the following two-dimensional version of the continuity equation of a wave flowing over land:

∂h ∂t = ∂qx ∂x + ∂qy ∂y = ie (5) (6) With h the depth of the flow and q the unit discharge in x or y direction. ieis the infiltration

excess. q is calculated using the Mannings equation.

qx= αxhβ (7) αx= S_{f x}1/2 nov (8) S_{f x}1/2= Sox− ∂h ∂x (9)

With nova tuneable parameter and β a coefficient set on 3₅. Sf xis the friction slope, which

is calculated using the terrain slope Soxand the change in depth of the water surface above

land.

To accurately calculate overland flow it needs to be calculated at a relatively small scale of 30-300m, as overland waves are typically smaller then the standard grid size of one kilometer. Bigger scale effectively smooths the topography adding more inaccuracy. Smaller scale however adds computation time. The time step for the wave calculation is directly dependent on this grid resolution. If a timestep too big is chosen, the wave might

(14)

have ’skipped’ a grid cell during this time, which leads to erroneous numerical diffusion. To prevent numerical diffusion the Courant Number Cn= c(∆x/∆t) needs to be around one.

To keep stability of the routing this value needs to be less than one. For the grid size of 250m used here the simulations use the recommended timestep of 15 seconds.

WRF-Hydro 5.0.0 is able to use a fine LSM grid, inducing the channels, or a vectorized network of channel reaches. In this research the vectorized network is used matching the National Water Model. The file route_link.nc explicitly describes these linked channels. Using explicitly defined channels is called linked routing. If WRF-Hydro uses linked rout-ing, the Muskinham-Cunge method is used to calculate the hydrologic routing. This is a one-dimensional explicit method, which does result in it not being able to allow local or backwater effects. The method uses the depth of each channel, as well as the trapezoid shape of the channel and the Manning’s Roughness coefficient. All these values are possible to change using the route_link.nc file for each single channel.

Halfway through the research data was requested which included the reservoirs in the San Diego area, which made the Lake routing scheme of WRF-Hydro relevant. It calculates the inflow and uses the size of the lake, the height of the water in the lake and the maximum height the lake can hold water at to calculate the outflow of the lake. This outward flow, which in lakes happens only if it is filled above capacity, is called the Weir flow. If the water level is higher than hmaxthe Weir flow is given by:

Qw(t) = CwLh2/3 (10)

With Cw the Weir coefficient, L the length of the weir and h the height of the water. If the

water body is a reservoir, an orifice flow can be available as well. This represents an opening in the dam, with a specified orifice area and coefficient. These values of the obtained data are arbitrary and will be calibrated to make the streamflow match the observed data better. The orifice flow is given by:

Qo(t) = CoAo

p

2gh (11)

With Co the orifice coefficient and A the area of the orifice, g the gravity acceleration and

h the height of the water.

4.1.2 Preprocessing and input files

This research uses a specific dataset of the San Diego area provided by NCAR [29]. It includes all files created by the pre-processing GIS tools of the WRF model. Here the important files are described as well as changes made to them to improve the models output. In later sections these files could be referenced as more specific alterations or calibrations are made.

The geospatial data is defined in the geo_em.d01.nc file. It is used by the WRF-Hydro preprocessing tools to create the other files in the same grid using ArcGIS.

(15)

Fulldom_hires.nc: The main information that this file contains are: the channel grid, the information of where the channels are located and the flow direction grid, which shows in what direction the surface water flows in the whole research area. Finally it also defines the location of the so called ’gauge points’, where the streamflow of the channels are measured. The location of these points were chosen by NCAR [29] and while many of these locations do not match the location of USGS streamflow gauges, the main San Diego river has two matching gauge locations: Fashion Valley, which lies far downstream and has frequently occurring flash floods, and Mast Road, more upstream and therefore less influenced by urban irrigation. The streamflow of these gauges will be used as validation data.

Route_link.nc: This file, only created when the routing of the channels are not created by WRF-Hydro but given according to the real locations of the channels. This vector file defines the exact locations of the channels and sets which and where channels are connected. It also gives information on the channel depth, slope, roughness etc.

LAKEPARM.nc: Gives the lake parameters for each individual lake. It defines initial water level fraction as well as Weir elevation level, Weir length and coefficient, orifice area coefficient and elevation level. This file will be adjusted by calibrating the model to the streamflows.

The advantage of Noah-MP LSM over Noah LSM and other land surface models is the possibility to define spatially different parameters. This uses three parameter table files: GENPARM.TBL, which gives general global parameters for the model, SOILPARM.TBL, which gives the soil parameters for different soil classes, and the above mentioned MPTABLE.TBL, which gives parameter values for each land use class defined in geo_em.d01.nc used by NOAH-MP.

MPTABLE.tbl: This file gives the parameters for the land surface model. It gives over twenty different parameter for each spatial land cover data given in the geo_em.d01.nc file. Examples are leaf area index, height of canopy, leaf dimension and tree density. Even wood to non-wood ratio and microbial respiration parameter are given. Normally these values are not changed or calibrated for standard model runs. However, normally urban area is defined using the parameter isurban, which blocks the model from running many physical modules and forces for example evapotranspiration to be zero. This is not wanted here as we want to model urban area with some included vegetation patches instead of a fully concrete jungle. Where MPTABLE.tbl gives properties of the land cover, the properties of the soil them-selves are given in the file soil_properties.nc. It defines about fifteen parameters on the high 250 resolution scale, going from soil quartz content to parameters in calculating the runoff as well as maximum and minimum soil moisture contents.

Finally the file wrfinput_d01.nc gives the initial values of the starting state of the model. This goes from initial soil moisture to initial channel water level height and even the initial sea ice content.

(16)

4.1.3 Precipitation data

The precipitation data used comes from the NLDAS-2 dataset [30]. NLDAS stands for National Land Data Assimilation System. This gridded data is extracted from the Regional Climate Data Assimilation System, or R-CDAS, while older data comes from NARR, the North American Regional Reanalysis. Its analysis fields are on a 32km grid with a three hour temporal frequency. This information is then interpolated to the NLDAS 1/8th degree spatial grid as well as temporally disaggregated to an hourly frequency. This interpolation is done using the method of Cosgrove [31]. Next to the precipitation, the NLDAS data also includes many other values, like wind components, air temperature and humidity. In short it includes all information needed by most land surface models. Even though it is spatially interpolated it is currently the standard in precipitation data and has the the highest resolution.

The NLDAS-2 Forcing data is downloaded for the whole United States and comes in a GRIB file. In order to re-grid this data from this full sized GRIB file NCAR has developed a script in the NCL programming language. The script and the programming language are downloaded from the NCAR website [29]. Using the preprocessing geospatial geo_emd01.nc file as input it regrids in two steps. The first one creates a regridding weight file specifying interpolation weights between the source (NLDAS) coordinate data grid and the destination (geo_emd01.nc) coordinate data grid. The second step performs this regridding for all forcing input files, resulting in hourly NCAR .nc files in the same coordinate system as the other WRF-Hydro input on a 1km grid. The steps above were performed multiple times, for the period of January to April and September to December 2016. The downloaded files have a size of 2.5MB each, resulting in over four gigabytes for the second period, while the regridded files have a size of 200kb, giving a total size of 250MB for the second period. 4.1.4 Parameter choices and initialization

The model has two main initialization files, also called namelist files. They are namelist.hrldas and hydro.namelist. namelist.hrldas defines the NOAH LSM options like time to run, starting date, physics options to use and output frequency. hydro.namelist defines the settings for the hydrological part: the routing components, which modules to be included, what kind of output to give and where the input files are.

It is very important to have the grid spacing DXRT in hydro.namelist matching the high resolution input files, which is 250m in our fulldom.hires file. As the land model grid defined in geo_em.d01.nc has a resolution of one square kilometer the aggregation factor needs to be 1000/250 = 4. Furthermore the terrain and channel routing steps are very important to have correct. They are directly dependent on the size of the grids used. As channel flow flood waves can have great flow depth and wave amplitude, they can have a strong accelerations of the propagating wave. These waves cannot pass multiple grids per calculation iteration or numerical diffusion can occur. To prevent this the Courant Number

(17)

Cn= c(_∆x∆t) needs to be close ot one. As our grid spacing is 250 meter this means our time

step should be around 25 seconds. However, the Cn also affects the stability of the routing

routine in general and therefore should be less than 1. As recommended for a 250 meter scale in the Technical Description our timestep is chosen to be 15 seconds.

In this research all hydrological modules are turned on in the hydro.namelist file to have the most powerful and accurate modelling capabilities. Overland and subsurface flows are turned on as well as channel routing using the Muskingham Cunge reach method. Lakes are turned off for the initial exploratory phase and are turned on for the final results. The baseflow bucket model is activated with the exponential bucket method as most extensive groundwater model. More about these modules can be found in the technical description section or in the WRF-Hydro Description [19].

4.1.5 Docker

Recently NCAR has released a Docker module to help with the installation of WRF-Hydro [32]. In this study this Docker container is used as a virtual environment. It comes with Python with many of the required packages and makes it possible to run WRF-Hydro, as well as the model created in this study on any computer. The Docker container with newly created jupyter notebooks including the framework and irrigation model can be provided if requested.

4.2 Framework implementing irrigation model

Implementing custom scripts or models in WRF-Hydro is possible in two ways. The first one is making your own module that can be added when compiling WRF-Hydro. It was chosen to not do this however, as the code would have to be written in Fortran and the model needs to be compiled and thus could give incompatibilities with future versions. In addition to this the documentation of writing modules was completely absent. This is why it was chosen to create a framework around WRF-Hydro, making use of its restarting capabilities to change its input data depending on its previous output. This makes the framework stand-alone from WRF-Hydro and more compatible with future WRF versions as long as input and output data used does not change. However, the disadvantage of this is that every time the irrigation model is called multiple files need to be loaded, increasing running time.

The framework is created using Python and the Jupyter bash magic command to run the bash script calling WRF-Hydro_MP.exe.

Figure 4 shows the UML for this framework. The WRF-Hydro is set to run for a certain amount of time, for example a single day. At the end of its run WRF-Hydro saves its state in a RESTART and HYDRO.RESTART file, for the Noah-MP and the hydrological states respectively. These files contain the initial conditions used in the subsequent run.

The irrigation model obtains the spatial soil moisture content of the top soil layer at the time of the irrigation event from the LSMOUT output file. The model also uses the spatial

(18)

Figure 4: UML of the irrigation model framework, with spatial and human behavioral factors impacting irrigation amount and timing.

(19)

data on the maximum soil moisture content as defined in the soil_properties.nc file. The irrigation model uses this information along with additional higher resolution spatial information from other sources to create the spatial irrigation values, scaled to the one kilometer grid. It is converted to the dimensions used in the precipitation files: [mm]_s . If the irrigation event is set to take multiple hours this total irrigation rate is divided by that value.

Then this irrigation is added to the precipitation rate of the HRLDAS forcing input, the precipitation file. It uses the file(s) of the following hour or hours after the WRF-Hydro interrupted its run. These precipitation files are then used as input for the next WRF-Hydro run. This way the irrigation is simulated as additional precipitation. The advantage of this is that no new WRF-Hydro input file or parameter needs to be created. However, an implication of this is that even though it is possible to calculate irrigation rates on a high spatial resolution, it is reduced to a one kilometer grid when being added to the precipitation input. So even though irrigation in reality occurs only on the garden area, which is for example forty percent, this amount of irrigation is spread out to the full cell size in WRF-Hydro, reducing local garden soil moisture, especially immediately after irrigation when lateral spreading has not occurred yet.

In order to automate the script to call WRF-Hydro using the right settings and initial conditions for each run the initialization files namelist.hrldas and hydro.namelist are adjusted each cycle by the framework. Altering these text files is done using python and regex.

namelist.hrldas includes the starting date and time of the WRF-Hydro run, as well as the amount of time in hours or days the model should run. These are updated initially, while every cycle the date and time are updated to the new ones as well. In addition to this the name of the RESTART file to use is updated every cycle, using the file created by the previous cycle. In addition to this the hydro.namelist reference to the HYDRO.RESTART file is updated, facilitating each run to match the hydrological state of before the irrigation event.

In reality everyone irrigates at a time convenient to them, which results in temporal spatial spread of influx of irrigation water. However, in order to keep computational time minimal it is not feasible to have WRF-Hydro stop and restart at every hour. Therefore it has been chosen to irrigate once every day for the whole area, over a time period of three hours. The time of irrigation is done at midnight, to avoid a three hour peak in evapotranspiration and soil temperature during the day.

4.3 Irrigation Model

This section goes in depth on the newly created irrigation model, starting with some theory and related work and explaining the input data. Then the scheme is explained as well as the scenarios chosen.

(20)

Figure 5: Field capacity and wilting point for different soil types. Courtesy of [33]. 4.3.1 Soil moisture and healthy vegetation

To understand the irrigation model some terms need to be explained about the relationship between soil moisture and vegetation. Plants need a certain soil moisture in order to not dry out and wilt away. This soil moisture is called the wilting point. Field capacity is the value of soil moisture at which the soil that has drained all its excess water ends up at. Figure 5 shows the relationship between field capacity and wilting point for different soils. Added water can not be taken up by the soil and will be drained away. This depends highly on the type of soil, as clay is able to hold a lot more water than for example sand. When irrigating correctly, soil moisture is aimed to be somewhere between the wilting point and field capacity, as it needs to be higher than wilting point for the plants to survive and grow, but adding more water than field capacity will have this water just drained away.

4.3.2 Irrigation Scheme

The basis of the irrigation scheme comes from [16] and [17], as described in the Related Works section. It assumes homeowners irrigate their garden up to a certain point of soil moisture of the top soil layer. This means that if a storm has just occurred, no irrigation will be done, while during the dry summer this irrigation is the sole source of soil moisture. The calculation of amount of irrigation is done as following, for each grid cell:

(21)

SMirr= α · SMmax

DEF = max[(SMirr− SM (t)), 0]

IRR = ρw

∆tDEF ∗ h ∗ Ag

SMirr is the target percentage of Soil Moisture, just after an irrigation event and

SMmax the maximal possible value of moisture the soil type can handle, obtained from

soil_properties.nc. α is the target soil moisture as a fraction of the maximum soil mois-ture saturation. DEF is the resulting deficiency of soil moismois-ture and target soil moismois-ture, and is calculated to the absolute irrigation amount IRR over time with the h height of soil layer used in the Land Surface Model of WRF-hydro, ρw the density of water and ∆t the

time of irrigation. Agis the percentage of area per cell that actually receives irrigation, the

garden area. The method used to calculate this area will be explained in the next section. α is where the decisions and needs of households and vegetation types come into play. It is a unitless value between 0 and 1, with 1 meaning that homeowners irrigate their garden to the point of full saturation. It can be seen as αphys∗ αhuman where αphys depicts physical

needs of the garden in order to keep the plants alive and growing. This value is best be-tween the ’stress point’, also called the allowable depletion, and the field capacity, as can be seen in Figure 5. αhuman depicts the effects of human decision making and can be used to

define different scenarios where residents different irrigation decisions. h, the depth of the soil to be irrigated, is set on ten centimeters, which is the height of the top soil layer used in WRF-Hydro. The soil moisture of this layer, which is ten centimeters high, is used for calculating the needed irrigation.

This way it the assumption is made that when people irrigate they take the top ten centimeters into account for estimating how much they want to irrigate. However, since both h and α are multiplied with each other, if homeowners use a different height of the soil they irrigate, this could be incorporated in a different α, making

Using soil moisture content of the top layer to calculate the amount of irrigation comes with the assumption that people irrigate depending on the soil moisture. Manual irrigation directly takes this into account, as people irrigate up to an arbitrary point of wetness of the soil ??. Sprinklers that are on a schedule without the use of soil moisture meters or weather data irrigate a set flat amount, not taking soil moisture into account. However, after a big storm most homeowners change their irrigation schedule and often stop both sprinkler and manual irrigation as long as soil moisture stays above a certain value.

4.3.3 Satellite data

The most important thing needed to accurately calculate the amount of irrigation water that enters the soil is an accurate estimation of the area of irrigated gardens. To do this we obtain high resolution land cover data to see what area is classified as urban, as only urban

(22)

Figure 6: Skewed impervious urban area .

areas have gardens that are irrigated in the San Diego area. Then the percentage of this urban area that is not impervious is calculated. Impervious areas could be rooftops, roads and terraces, and thus we assume that pervious urban area can be classified as gardens and thus as the area of irrigated soil.

A combination of the 2011 NLCD land cover [34] and NLCD impervious data [35] is used. They have a resolution of 30 meters, the highest resolution available. The land cover data classifies the cells in 20 different classes [36], ranging from shrub, forest types, water as well as three different classes of urban area. The impervious area describes the percentage of impervious area for each cell.

The NLCD data uses the satellite data from the USGS Landsat mission [37], [30]. These are 7 satellites obtaining light in many different spectra. The NLCD uses mainly the data of the Landsat 7 satellite, which has the ”Enhanced Thematic Mapper Plus” sensor which is used for the most accurate classification as it also obtains thermal and panchromatic data [30].

A regression tree algorithm was developed to accurately classify the land cover, with an inaccuracy of lower than four percent [38]. The impervious data was calculated by Yang et al. [35], having an accuracy of around 90 percent.

To use this data in our model it needed to be adjusted significantly. Doing this data manipulation consisted of the following steps, using both ArcGIS Pro and python. First the data was cropped and matched with San Diego modelled area using ArcGIS, taking into account the cells overlapping the border. This resulted in Figure 6, a skewed and flipped grid, as the coordinate systems were different. Using python the data was flipped and realigned to match the square area used as model input. This resulted in a grid of 1219 by 2621 cells of 30 by 30 meters.

(23)

Figure 7: Percentage of area that is classified as garden in a 30 by 30 meter grid. .

Then to calculate the area of gardens all cells that were not classified as urban were set to 100 percent imperviousness. These cells are also set to 100 percent imperviousness to not overshoot the calculated garden area. The cells that were classified as urban were set to their value of imperviousness. This results in 1 − A2_imp = A2_garden, taking into account that non urban pervious area cannot be classsified as garden. This garden area can be seen in Figure 7. Interesting to see is that in downtown San Diego, bottom left, the impervious area is low, while more north, reaching the richer La Jolla area the percentage of garden is higher. It can also be seen that the important distinction between urban garden pervious area and urban (non irrigated) canyon pervious area is made.

Then, as the model uses a 1 kilometer grid for the land surface model, the resolution is reduced to the standard grid of 35 by 77 as shown in Figure 8. Finally the values of these cells are used

4.4 Timescale and Spinup

The model is ran for two different time periods. For the initial exploratory results we have chosen to run the model during winter time: from January to the end March or February 2016. In this period the most storms of the year happen, so most natural precipitation occurs, which is helpful in investigating initial model output. Most of the other months are completely arid and thus do not give a good indication of the quality of the model runs. Later, when investigating the impact of irrigation is of importance we focus on the arid summer, on the first few rainfall events of the year, from September to the first storms in November and the rainy December. We expect these first storms to be the ones most impacted by irrigation, as the garden soils normally completely dry after the arid summer

(24)

Figure 8: Percentage of area that is classified as garden in a 1 square kilometer meter grid. .

now have a significant soil moisture.

The model is run for different timescales depending on the required results, for hourly time traces the scale is a few days up to a month, while daily results can be shown for multiple months. The main factor here is the number of output files to work with. If output is taken hourly for a full month it results in over seven hundred files for each output file. With over five different possible output files and these files being up to 200kb, processing this output is a limiting factor on model running time.

The model reaches a steady state in an arid setting after about a month of warm-up time. This warmed-up state is stored in a RESTART file which is used as a starting point for all runs. However, as irrigation occurs the model needs to reach a new steady state. Figure 9 shows the time it takes to reach this state using α = 0.65.

4.5 Validation Methods

4.5.1 Evapotranspiration

The main validation method used in this research is evapotranspiration.

Evapotranspiration (ET) is the total amount of water lost from the ground to the at-mosphere by the combined processes of plant transpiration and evaporation as well as soil evaporation [39]. A schematic view can be seen in Figure 10. Evapotranspiration is one of the most uncertain terms in the water surface balance, as it is influenced by many soil and hydrological processes. It depends on wind, solar radiation, humidity, temperature, soil characteristics and plant characteristics like root depth and stomatal resistance. Therefore it is a good indicator for assessing land surface models performance [3]. It is also concluded

(25)

Figure 9: Soil moisture during the first few irrigation events. A steady state cycle is reached after about five iterations. Also interesting to see is that after each irrigation event where the soil moisture is increased significantly, it then drops off by a little during the night, caused by infiltration. Then when the day starts it can be seen that the soil moisture drops significantly caused by the increase in evapotranspiration.

that evapotranspiration is highly correlated with irrigation [40] in Canberra, Australia, a semi-arid urban area. This is understandable as in these areas irrigated plants and soil have significantly more water to evaporate than in a non irrigated case, especially during dry summers.

While it is possible to measure ET directly using lysimeters, often they are estimated using other more easily measurable characteristics and analytical and empirical equations. Using a reference crop like grass or alfalfa whose biophysical characteristics are well known, as well as observed wind, solar radiation, humidity and more, a reference ETr or ETo is calculated. ETo can then be converted to the evapotranspiration of other crops using a crop factor: a coefficient Kc [39]. While this does add inaccuracy, it also gives the possibility to

have spatially separated ET based on spatial vegetation data. ET calculation using ETo is shown to be sufficiently accurate especially for vegetation less than one meter tall, which is predominantly the case for the bushy and dry San Diego natural vegetation [42].

The evapotranspiration data used in this research is obtained from the ”California Ir-rigation Management System”, or CIMIS [39]. CIMIS is mainly used as a tool for farmers to manage their water usage more efficiently as it helps in estimating how much irrigation is needed. CIMIS has multiple sensor stations, with each station having a solar radiation meter, soil thermometer, air temperature meter, wind direction meter, wind speed meter and precipitation meter. Using these observations the ET for the case of a closely clipped

(26)

Figure 10: Evapotranspiration includes evaporation from soil and transpiration from vege-tation. It is dependent on characteristics of this vegetation, of the air like humidity, tem-perature and windspeed as well as of the soil. Courtesy of [41]

grass that is well watered and actively growing at these stations is calculated, resulting in ETo.

This station data is then converted to spatial data using a combination of interpolation and additional satellite data. It uses a statewide two kilometer spatial estimated evapo-transpiration of the American Society of Engineers using a Penman-Monteith equation [43] and GOES satellite data. While CIMIS does not give a value of the accuracy, it does state that it is significantly more accurate than taking the closest station data and it is accurate enough for farmers to be used as indication.

On the website of CIMIS one can download historical station data or daily spatial data of last year. After contact with CIMIS I have obtained daily spatial data for all years since 2010.

This spatial data is stored in a .tiff format of 2km grid of the whole of California. Each day is classified in a different band. Again a difficult task was manipulating this data to a netcdf file that is easy to use.

(27)

Figure 11: Land Cover and edited land cover used for ET. The only significant land cover other than urban and shrubs is the orange area, which represents evergreen needle forest.

Using ArcGIS Pro the data was clipped to solely the case study area, taking into account the boundaries. Files with the .tiff extension can have multiple bit depths. Working with a different bit depth than expected results in corrupt displays and possibly fatal errors, so it first had to be converted to a 16 bit depth for ArcGIS to be able to work with it. Then, since ArcGIS does not work well with multiple bands of a .tiff file a GIS python script was created to perform this clipping for each band and store it as a netcdf file. Lastly these netcdf files were put into one single file with the additional dimension of day in the year.

This resulted in a daily ETo, which still needs to be converted to ET using crop coeffi-cients. The left of Figure 11 shows the Land Use Index of the file geo_em.d01x. This is a grid with a spatial resolution of 1 square kilometer, given as initial data by NCAR shows the main land use of each cell. For the sake of simplicity in calculating the ET this was changed to only three main land covers, water, urban and shrubs, as can be seen in the right image of 11. The main difference is that the orange colour disappears. This colour represents evergreen forest. However even though this is in the original data, it does not exist anymore due to recent wildfires.

In the urban area we assume that only the gardens contribute to the evapotranspiration. To calculate the right value the cells classified as urban are multiplied with the percentage of garden in that cell as well as with the Kcurban, which is taken to be 0.8, which is in

accordance with the prevalent turfgrass lawns as well as with home fruit crops and flower beddings [44]. The shrubland is multiplied with 0.4, also in accordance with [44].

Figure 12 shows the mean of the daily ETo as well as the adjusted ET. It is clear that during the summer more is evaporated and evapotranspiration is higher.

(28)

Figure 12: Average daily evapotranspiration of the base and adjusted value of CIMIS in 2016. Note that ET increases during spring and summer.

.

4.5.2 Streamflows

Streamflow data in the San Diego river is obtained from the USGS streamflow gauges [45]. As can be seen in Figure2: USGS has two gauges in the San Diego river. One at Fashion valley which is far downstream and at a location where the river often overflows, and one at Mast road, more upstream and therefore expected to be less impacted by urban irrigation. This data can be downloaded for a specific time period in 15 minute table format. After conversion from [f t]/[h] to hourly data in [m]_[s] it is used in this research as the main validation data of the streamflow.

4.6 Parameter Calibration and Sensitivity Analysis

In this research calibration on multiple parameters and validation output will be per-formed.The paper by [46] performs slave-master parallel parameter calibration and vali-dation of almost all parameters in all tables, as well as for all soil types in MPTABLE. Their results were helpful in deciding the important parameters for calibration in this research, especially in parameters important for reproducing accurate streamflows, stating for exam-ple clearly that they put the RETDEPRTFAC to an extremely low value to be able to obtain streamflows of the same scale as the validation data.

A local calibration and sensitivity analysis on α will be performed, using the average evapotranspiration of the area as validation. This is done in order to obtain a value for α such that modelled and observed urban evapotranspiration match best. The α found here will then be used to most accurately estimate the amount of irrigation being done for the subsequent experiments.

(29)

Furthermore sensitivity analysis of α and on two other main streamflow impacting pa-rameters following [46] is performed, to investigate the importance of urban irrigation on the streamflow of the San Diego river and on possible flash flooding. These parameters are the RETDEPFRAC multiplier and OVROUGHRTFAC multiplier in the file Fulldom_hires.nc. RETDEPFRAC stands for retaining depth fraction, which defines the depth of retaining local water puddles. Decreasing this increases the total amount of water reaching the channels. OVROUGHRTFAC impacts the speed at which water flows overland before reaching a channel. Increasing this decreases the time needed for water to reach the channels after a rain event.

(30)

5 Results

The most important results of this research for policy making would be to understand the effects irrigation has on the overflowing of the San Diego river, while the academic value of this research lies more in the performance of the model in general, which is measured well by accuracy on the evapotranspiration. Therefore, the result section is split up in three main parts. The first section explores the results of WRF-Hydro and gives an indication of the process that this research has gone through in order to work with WRF-Hydro. It shows multiple hurdles that have been overcome and shows the changes made to create a workable model.

Then the second section goes into implementing irrigation. It shows the effects on multiple outputs of an arbitrary estimation of irrigation. It again goes into the process and problems occurring with different irrigation amounts. It shows the effects of this arbitrary irrigation amount on multiple validation data. A calibration of the amount of irrigation is done, using evapotranspiration as the validation data.

The last section explores the results of irrigation on the streamflow of the San Diego river, both upstream of most urban area at Mast Road and downstream close to the coast at Fashion Valley. Their locations are shown in the section on the case study area, in Figure2. A calibration on streamflow without irrigation is performed after which irrigation is added. This, in combination with a sensitivity analysis on irrigation and two main parameters gives an impression of the impact of urban irrigation on streamflow.

5.1 Computational time and efficiency

Before model output results are discussed, the efficiency of the framework is investigated. The model runs take a significant amount of time, especially if ran on a personal laptop. It takes approximately ten minutes to run the model for a single month, and many runs with different parameters needed to be done to get a feel of the model and its important parameters, initialization values and to make sure input data is correct. In estimate a total over one terabyte of output has been created.

The model runs are performed on two machines: a MacBook 13’ 2015, with a 2.7Ghz dual core i5 processor and 8gb of ram and a desktop with an i5 2700k quad-core clocking at 4.2Ghz and 8gb ram. Both machines use an SSD for storage. Note that the docker container is limited to using only two cores.

Table 1 shows the different computing times. Initialisation time, which loads many input files and creates the garden area is more hard drive dependent and less CPU dependent, as both devices take about 40 seconds for it. It is also clear that WRF is capable of fully using the CPU’s power, as during runs it uses two cores at 100%. This shows in the table, making the difference between the laptop and desktop on WRF-Hydro cycle time significant.

Therefore a desktop was used for longer runs as well as for the especially time consuming calibration and sensitivity analysis runs.

(31)

Action Laptop [s] PC [s] Framework Initialisation time 43 41 Framework cycle time 3.3 2.25 WRF processing cycle time 23 13.3

Table 1: Computation times of different model actions for diffferent devices

Using these devices results in a framework cycle taking between 15% and 17% of total cycle time. Framework cycle time decreases by more than 30%, while WRF processing time decreases with more than 40% using desktop over laptop. This means that while the framework does scale with processing power, it does not scale as well as the WRF-Hydro model does. This implicates that the use of a supercomputer, which is often done with hydrological models, is not viable for this framework, as the framework has limited scaling with processing power. This limited scaling can be explained by the fact that the framework uses python, and while all computing processes use numpy, xarray and pandas which are optimized on efficiency, it still is not a compiled script. Furthermore, most of the time spent by the framework is spent on reading and writing files. For example calculation of irrigation takes only 50ms while each read-write action takes about 250ms.

5.2 Initial WRF-Hydro results

First of all we look at the results of the model with no adjustments and with all parameters at their standard value. We first compare the streamflow output at Mast Road and Fashion Valley, the two locations where streamflow observations are taken. This initial result can be seen in Figure 13. This Figure shows two main problems: the modeled streamflow is too ’late’ compared to the observed value by around five hours. The other problem is that especially at the end of the storm the decay is slower than the observed decay resulting in a significant higher total volume. The fact that the streamflow is too late can be caused by two main factors. The speed at which the precipitated water reaches the channels could be too slow, so either or both overland flow and groundwater flow takes too long to reach the streams. The other reason could be incorrect precipitation input data. This data is taken from NLDAS and is, as explained in the previous sections, a product of the 32km NCAR data, spatially interpolated to a one kilometer grid and temporally interpolated from three hours to one hour. All these adjustments can bring errors with them.

In order to find out the reason for this offset in streamflow we compared measured precipitation data with the forcing input data. Initially, one weather station with hourly precipitation data available in San Diego was used. It is the [45] land-based weather station at the San Diego International Airport. The 15-minutes data was downloaded and adjusted to compare with the given precipitation at the gridpoint corresponding with the location of the airport: 32.740291, −117.187002. With the units matched to mm/h this resulted in Figure 14. Here it is clear that the behaviour is similar, but with what looks like an offset of

(32)

Figure 13: Modeled streamflow of San Diego river gauges at Fashion Valley and Mast Road compared with the observed values.

9 hours and different magnitudes. Seven of these nine hour difference can be explained by taking into account time-zones: NLDAS data is in GMT while the airport precipitation is in Pacific Time. The other two hours could be caused by the three-hourly time interpolation as well as local differences with the spatial interpolation of NLDAS. In Figure 15 the NLDAS data is rolled back by nine hours, and as can be seen, the timing of the peaks match perfectly. As the NLDAS data would fit reality more if it were changed by these two hours, we have chosen to implement this for the following parts in the research.

The magnitude of precipitation differs significantly as well. although it is in the same order of magnitude. In addition to this comes the fact that the modeled streamflow mag-nitude matches the observed values quite well, at least in the same order of magmag-nitude, indicating that the total precipitation input in the area matches reality. 16.

However, the discrepancy of precipitation along with streamflow results not matching in time gave the need for more validation of the precipitation data. CIMIS station data was used for its precipitation measurements. CIMIS has three observation stations in the San Diego case study area. Their coordinates were matched with their corresponding position in the 1 kilometer grid. These observations were from September through the end of December 2016. This time frame encapsulates the end of the dry summer, where precipitation is almost zero, to the the first winter storms. In Figure 17 a comparison is shown of the model input precipitation and the measured values. Again the timezone difference and two hour error have been taken into account. The input precipitation magnitude however is significantly closer than the one from USGS, confirming that no change is needed. The main difference

(33)

Figure 14: Comparison of observed precipitation at San Diego Int. Airport with the NLDAS Forcing data. It can be seen that the rainrate peaks behavior is the same but with an offset. of observation and input data can be seen at the end of December, where it is clear that precipitation input is created using interpolation. Here the precipitation input is similar in all three locations, while for example the observed values at San Diego did not measure any precipitation.

5.3 WRF with lakes and streamflow calibration

After seeing that the initial WRF results were far off the observed values, it was hypothesized that the reservoirs in the San Diego area changed the hydrological balance significantly. While peak streamflow height was similar, the decrease of streamflow over time was slower in the model. This indicates that the model modelled significantly more water from upstream, which takes longer to reach the San Diego river gauge points. Therefore it was hypothesised that water management in the form of reservoirs in the San Diego area needed to be included in the mode.

Especially the San Vicente and the El Capitan reservoir which are directly connected to the San Diego river are expected to have a significant influence on the streamflow. These reservoirs hold much of the water that arrives from the mountain ranges upstream. This could be the reason that the streamflows in the validation data drops off faster than the modeled one, the water that flows all the way from the mountain ranges and should arrive at the measurement locations some time after the storms, simply is stopped by these reservoirs and never reaches the San Diego river. To see if this is correct a new dataset with lakes included was obtained. It uses the four main reservoirs in the area, including the El Capitan and San Vicente reservoirs, which are directly connected to the San Diego river.

The initial results of the streamflows with included lakes were farther off the observed values than the version without lakes. A comparison with the no-lakes and the validation data can be seen in Figure 18.

(34)

Figure 15: Comparison of observed precipitation at San Diego Int. Airport with the NLDAS Forcing data rolled back 9 hours. It can be seen that the timing of the peaks match perfectly, but the amplitude is off.

Figure 16: Modeled streamflow of San Diego river gauges at Fashion Valley and Mast Road compared with the observed values, adjusted for the time difference. It is clear that the timing fits the observed data better.

It is clear that there is a consistent six cubic meter per second in the WRF with lakes coming from the reservoirs. The impact of the storms are also very minimal, indicating that the reservoirs take up most of the rainfall from the mountains. However, it is possible that other changed parameters and files in the system have altered this as well.

(35)

Figure 17: Cimis station observations of hourly rainrate compared with their precipitation input file values at their respective locations.

explained in the technical description, the output of the lakes mainly depends on the water level, and the orifice and weir outlets. As these reservoirs are set to never overflow, no output should come from the weir. The main parameters then are the orifice area OrificeA, coefficient OrificeC and height of orifice OrificeH. The initial fraction water depth ifd has an impact similar to the height of the orifice, deciding how much water needs to be added to the lake for the orifice to spill. Other important parameters are in the spatial fulldom_hires.nc file. These are the RETDEPRTFAC, which is a multiplying factor of the maximum retention depth. This is the depth at which the ground holds water before it is routed as overland flow. Decreasing this makes the model have less small ponds and standing water, and have more water routed to the channels. This is expected to help the overall magnitude of the streamflow. The other is OVROUGHRTFAC, which is a multiplier on the Mannings’ Roughness coefficient. This impacts the speed at which the water reaches the channels, and thus impacts the timing of the peak.

A short Markov-Chain Monte-Carlo calibration was performed altering these parameters, using a single parameter for all reservoirs. The result is significant and can be seen in Figure 19. The FAST sensitivity index can be seen in Figure 20 and Table 5.3. The most important parameter is the RETDEPRTFAC, which alters the total volume of water reaching the streams.

(36)

Figure 18: Streamflow comparison of WRF, WRF with lakes and the validation data during the first storms of the year in November and December.

Also important is the OVROUGHRTFAC. Less important appears to be the initial depth of the lakes, which is surprising, as the base streamflow without precipitation is still significantly too high. This might be caused by too few iterations and the Markov Chain getting stuck in a local minimum.

Parameter name initial value calibrated value FAST index

ifd 0.9 0.177210 0.243132

orificeC 1 0.383204 0.874352 RETDEPRTFAC 1 0.481287 0.929871 OVROUGHRTFAC 1 0.221927 0.477100

Before more extensive calibration was performed, more effort was put into phenomeno-logically improving the output. After more research on the reservoirs and their output it was decided that the reservoirs either cannot open their dam or they are in practice never opened to the San Diego river. They do have some ’orifices’, but they are not connected to

(37)

Figure 19: Calibrated results of the streamflow. Maximum objective function converges from the initial −99 to −3.022

the hydrological cycle, as they are used to provide the San Diego population with drinking water. Much of this water reaches the hydrological cycle again in the form of irrigation. In addition to this it was found that while these reservoirs catch all water from the upstream mountain ranges, their main source is their artificial connection with the Colorado River. Therefore the volume of these reservoirs is such that even for big storms the water from the mountain ranges does not impact the reservoir levels enough to cause regular overflowing of the dams. Because of this the reservoirs were chosen to have an output of zero, which was done by setting the OrificeC and WeirC to zero. This resulted in the streamflow shown in Figure 21.

It is clear that, while the base streamflow is now back to zero, the increase caused by precipitation during is significantly too low. A factor that could influence this is errors made in the spatial interpolation of the precipitation model input files as can be seen in Figure 17. It could be the case that in the forcing input most rain falls in the mountainous area, which then is caught by the reservoirs, while in reality more rainfall occurs downrange, west of the mountains and does reach the river. Another source could be that the soil moisture before

(38)

Figure 20: Sensitivity analysis on important streamflow parameters on the streamflow of the San Diego river Fashion Valley gauge.

Figure 21: Streamflow of San Diego River with zero output from the reservoirs. the storm is very low which makes the soil able to uptake much of the precipitated water only to slowly release it back into the hydrological cycle, instead of it resulting in overland flow

A hydrologic model in semi-arid urban areas with urban irrigation and water management

MSc Computational Science

Master Thesis