INTERCEPTION AND THE SENSIBLE HEAT FLUX High-resolution modelling of interception’s influence on the sensible heat flux using limited high-resolution data.

(1)

Image source: http://www.uva-bits.nl/project/flight-energetics-of-urban-lesser-black-backed-gulls-a-case-study-for-bio-inspired-flightuav-technology/

INTERCEPTION

AND THE SENSIBLE

HEAT FLUX

High-resolution modelling of interception’s influence on

the sensible heat flux using limited high-resolution data.

ABSTRACT

This paper presents the interception module for a high-resolution sensible heat flux model for the purpose of microscale analysis of the migratory behavior of soaring birds. Using high-resolution precipitation and land-use data, all other parameters can be interpolated from low-resolution model outputs while still producing realistic model results. This allows data from low-resolution GCM’s to be incorporated in small-scale modelling using high-resolution data only for the most influential parameters for the model. After a sensitivity, model dynamic and spa-tial analysis, this paper concludes the described model to be plausibly correct, although there is not enough data available to quantitatively verify its results.

Maarten Mol

Bachelor thesis Earth sciences 2015/2016 Supervisors:

prof. dr. ir. W. Bouten dr. J.Z. Shamoun-Baranes

(2)

1

1. Introduction

As climate models or global circulation models (GCM’s) have grown in complexity and accuracy, most of the focus of these models has been to explain the earth’s climate as a whole. Employing large grid sizes and time steps in order to reduce the already massive computational power required for such a task. These models capture the earth’s system as a whole but, as a result, often do not allow for predictions on a microscale level (Hewitson & Crane, 1996).

However, there are types of research that do require microscale analyses of atmospheric conditions. One such application is to increase the understanding of bird soar/flap-patterns through the microscale mod-elling of the sensible heat flux. Sensible heat flux is one of the main components in the calculation of convective processes, which are key factors in migratory behaviour of soaring birds (Vansteelant et al., 2015; Mandel et al., 2011). This research project, a joint effort of Casper Borgman, Daniel Kooij, Maarten Mol and Bart Sweerts, is undertaken with the UVA-Bits program from the University of Amsterdam in mind and aims to provide a high-resolution model for the sensible heat flux to enable the analysis of soaring bird behaviour in relation to atmospheric convective processes.

Perhaps there is a way to circumvent this consequence of the current limitations on computing power by using large output data in yet another model, to allow for higher resolution analyses or even predictions. This is done based on a division between slow and fast processes where precipitation, for example is a highly variable, event-driven process and thus needs a high-resolution dataset in order to implement correctly. Whereas temperature, or atmospheric pressure do not differ significantly over distance and can therefore be interpolated based on data from these GCM’s (Haylock et al., 2008).

This research project aims to provide a proof of concept that it is indeed possible to create a relatively accurate model of the sensible heat flux, based on low-resolution data broken down using limited high resolution data as a denominator. This specific study within the research will focus on the modelling of rain interception and interception evaporation as it influences the sensible heat flux.

(4)

3

2. Methods

2.1 Theoretical Framework

2.1.1 Sensible heat flux

The earth’s surface and atmosphere maintain a constant interaction through the exchange of energy. With sensible heat flux as one of the main terms in the earth’s energy balance, it is necessary to solve for every element in order to get an accurate final result. The sensible heat flux can be defined according to the following equation:

(1 − 𝑟)𝑆 ↓ + 𝐿 ↓ = 𝐿 ↑ + 𝑆𝐻𝐹 + 𝜆𝐸 + 𝐺 (W) (eq. 1) The left-hand side of the equation contains the albedo (r) from which the absorbed solar radiation (1-r)S↓ is calculated and the longwave radiation (L↓). The right-hand side contains emitted longwave radiation (L↑), sensible heat flux (SHF), latent heat (λE) and soil heat exchange by conduction (G) (Bonan et al., 2002). The parameter directly influencing thermal lift is SHF, which can be quantified by calculating and subtracting the other energy fluxes from the absorbed solar irradiation.

Latent heat can then again be divided into four parts: transpiration, ground evaporation, canopy evapo-ration and open water evapoevapo-ration (Shuttleworth, W., 1993 as described in Savenije, H., 2002). In our case, open water can be excluded from our model for the purpose of the land-based model intended for this research. While ground evaporation and transpiration are discussed in the evapotranspiration mod-ule by Bart Sweerts, this paper concerns itself with canopy evaporation and the rain interception pro-cesses that are required to accurately model it.

2.1.2 Rain interception

Savenije (2004) provides an overview of hydrological models that do not incorporate interception pro-cesses and shows how this leads to errors in the estimation of the soil moisture stock. Interception great-ly impacts the water balance by limiting the fraction of precipitation actualgreat-ly infiltrating into the ground (Linhoss & Siegert, 2016), which in turn affects both the water and energy available for regular evapo-transpiration. This model makes the assumption that incoming net radiation will be used first to drive canopy evaporation processes before transpiration or ground evaporation processes take place (Tolk et al., 1995).

Interception itself is concerned with the state variable of canopy storage, which is affected by precipita-tion as an inflow and evaporaprecipita-tion as outflow. This research will use a Rutter-like model (Rutter et al., 1972) with the exception that drainage is implicit within interception, rather than an explicit flow varia-ble. The storage for each next step is calculated through as follows:

𝑆𝑡+1 = (𝑆𝑡+ 𝐼)/(1 + 𝐸0

𝑐 ∗ 𝑑𝑡) (mm) (eq. 2)

Which is a backwards integration of:

𝑑𝑆/𝑑𝑡 = 𝐼 − 𝐸 (mm) (eq. 3)

Where the left-hand side is the change in canopy storage, and the right-hand side is made up of inter-ception (I) and Evaporation (E). Interinter-ception is defined as

(5)

4

𝐼 = 𝑎 ∗ 𝑃 for S+I ≤ c (mm/dt) (eq. 4)

𝐼 = 𝑐 − 𝑆 for S+I > c (mm) (eq. 5)

In these equations, a stands for the interception fraction and P for the precipitation. When Storage (S) would exceed the canopy storage capacity (c), the interception is the difference, or the remaining amount of water that can be contained within the canopy.

Evaporation is defined as a process which is dependent on the potential evaporation, multiplied by the fraction of stored water over storage capacity:

𝐸 = 𝐸0∗ 𝑆

𝑐∗ 𝑑𝑡 (mm) (eq. 6)

2.1.3 Canopy evaporation

For the purpose of this module, canopy evaporation will be treated as open water evaporation. This is a simplification, as evaporation from interception can be much higher as a result of advection from border-ing areas. However, this effect has only been proposed as the cause of higher evaporation rates of inter-cepted water (Wallace & McJannet, 2006), and not conclusively proven or analysed so as to allow for inclusion into a model. Therefore, potential evaporation will be calculated using a version of Penman-Monteith potential evaporation rate parametrised to open water, as described in Finch et al. (2001).

𝐸 =1 𝜆

𝑅_𝑛+𝜌𝐶_𝜌 (𝑒_𝑠− 𝑒_𝑎)/𝑟𝑎

∆+ 𝛾 (mm/s) (eq. 7)

Where λ is the latent heat of water, Rn the incoming radiation. Then, air density (ρ) is multiplied by air specific heat (Cρ) and the difference between saturated vapour pressure (es) and actual vapour pressure (ea). This is then divided by the aerodynamic resistance (ra) and then divided by the slope of saturation vapor pressure over temperature (∆) and the psychrometric constant (γ). For the calculation of es, ea, ∆ and γ, the equations described in Allen et al. (1998) were used. However, ra presents a problem, as most literature estimates it based on observations (Deguchi et al., 2006; Liu et al., 2007). Since this is not pos-sible for this study, it will use a simplified method of estimating ra developed by the agriculture research council (CRA):

𝑟𝑎= 𝐾/𝑢2 (s/m) (eq. 8)

(6)

5

2.2 Input data

The input data consists of low-resolution wind speed, temperature and dew temperature from the ERA-Interim reanalysis datasets and high-resolution precipitation data from the KNMI

RAD_NL25_RAC_MFBS_5min dataset. Aside from this, the land-use parameters of canopy openness and canopy storage capacity, which are discussed in the paper of Daniel Kooij, and net radiation, which is covered by Bart Sweerts, were also used as input. Table 1 provides an overview of the datasets and their resolution.

Initially, all datasets were converted to either ETRS89 (EPSG: 4258) or WGS 1984 (EPSG: 4326), as, for the purposes of modelling for the Netherlands, these projections only have a difference of 0.5 meters, which is hardly significant when modelling at a scale of 100x100m grid sizes. However, the program used to project these datasets, Arcmap, exports these two projections differently, leading to an inconsistency in the overlap of the different datasets. Therefore, the datasets which were processed using Arcmap have all been converted to ETRS89.

Source Data provided Spatial resolution Temporal resolution

Data format Units

The Satellite Appli-cation Facility on Climate Monitoring (CM SAF)

-Solar radiation 0.05 x 0.05 (de-grees) 1-hourly .nc (netcdf) [MJ] ERA-interim (rea-nalysis) by ECMWF - Temperature - Dewpoint tem-perature - Wind speed 0.0125 by 0.0125 (degrees) 6-hourly .nc (netcdf) [K] [K] [m s-1_] Copernicus land monitoring services - Land-cover 100 by 100 (m) - .csv [-]

KNMI Rad_NL25 - Precipitation 1 by 1 (km) 5 minutes (aggregated to hourly)

(7)

6 The raw precipitation data is set in a custom projection created through the use of a program called proj.4 (OSGeo/proj.4., n.d.). In order to get the coordinates in WGS84, the coordinates of each point in the custom projection has been saved to a txt file in an [X, Y] format. After this, the C-program “CS2CS.c” was compiled and applied to the dataset using the parameters included in the precipitation metadata. Furthermore, the terms: “+to +proj=latlong +datum=WGS84”, the paths to the file containing original [X Y] data and the destination file were included. After saving this to a .txt-file and extracting the relevant coordinates (see appendix B), this resulted in a file with the WGS 1984 coordinates for each point of the precipitation dataset. This data has then been set as X and Y fields in Excel, following the precipitation for each hour and loaded into ArcMap using the add xy data tool. This point data has been converted to a raster using the arcpy point to raster tool with a cell size of 0.015° and the assignment type “Mean”. This assignment type was chosen so as to minimize alteration of the water balance. Since precipitation is expressed in mm, or kg m-2_{, it cannot simply be summed when merging grid cells as this would strongly}

increase the total amount of water entering the system.

After this, it was projected to ETRS89 using the arcpy project raster tool. Finally, it was exported to an ascii format in order to allow it to be read in Matlab. The python-scripts used for the ArcGIS conversion are included in appendix C. Also, the script for the final loading of the precipitation data in Matlab can be found in appendix D.

2.3 Modelling approach

The interception model has been written in Matlab because of its capability of handling large matrices, which is convenient in the handling of large geographical datasets. For more details of the complete module, please see appendix A. In order to allow for a complete integration of the interception module, the evapotranspiration module and the ground heat flux module, the interception module has been writ-ten as a function, which returns the “consumed” lawrit-tent heat for the calculation of the SHF, and the cano-py storage, which is then used as an input value for the next time step.

The SHF-calculation in this module is rather simplistic. As interception is only a part of the SHF, it is only a preliminary value, and set as the difference between radiation and used latent heat. In the full model, the SHF will be calculated more precisely, as the full latent heat and ground flux will also be added. For the purpose of this research, the model will run for May and June 2012. This was chosen so as to get both a wet period and a dry period to run the model for and to be able to see possible differences in circumstances.

2.4 Analysis

The final output of the sensible heat flux model contains spatial (x/y), value-based (z) and temporal (t) values. This presents a problem for the adequate analysis and visualisation of the data and therefore it is important to define the statistics and processes we are actually interested in. This subsection will cover the methods used for the analysis of relevant model characteristics. For the statistical testing, an alpha of 0.05 will be used to determine significance.

(8)

7

2.4.1 Model validation

Firstly, it is important to know whether the model behaviour is correct, in the sense that it functions in accordance with the used equations and its results are representative of the actual processes involved. As there is little to no data available to compare the model outputs to, this needs to be done qualitative-ly, describing the dynamics within the model or within a single grid cell and verifying whether the behav-iour seems logical and reasonable.

2.4.2 Sensitivity analysis

Secondly, it is relevant to know how big the influences of input values are, so that in future research, one can decide which factors should be measured with more accuracy in order to provide reliable results. For this purpose, a sensitivity analysis is needed. This was done by taking the empirical minimum, mean and maximum values within all the input data. These three values were then used as input values for the interception module and modelled for 10

hours or time steps. Please see table 2 for the exact values.

The final analysis was done by grouping the three unique values for each parameter while keeping all other values at their respective means. After this, a t-test was performed on each time series compared to the overall mean of the parameter-set output, in order to de-termine how representative the output is for the three combined values, with the alterna-tive hypothesis that they are significantly dif-ferent. The script for the analysis can be found

in appendix E. Table 2: Parameter values for the sensitivity analysis

2.4.3 Interception’s influence on SHF

Secondly, for the purpose of the overall study of the SHF, it is important to determine whether intercep-tion actually has a significant impact on the SHF. If there is hardly a significant impact on the SHF, it would merely be wasted computational power to include it into the full model.

This test will be done using a pairwise t-test on the net radiation and the modelled sensible heat flux, treating all values as part of a population and thereby, seeing whether the impact of interception and interception evaporation is significant. Although the latent heat is always subtracted from the net radia-tion, one might consider a one-sided test. However, the latent heat can become negative at night, being effectively “added” to net radiation. But perhaps more importantly, since we cannot make the assump-tion of normality within both samples, the t-test would not be robust against a one-sided t-test of a non-normal sample.

Parameter: Min: Mean: Max: Canopy Cover Fraction 0.34 0.85 0.99 Canopy Storage Capacity 1.74 2.56 5.85 Precipitation 0 0.11 133.84 Net Radiation -3.13 *105_{4.15 * 10}5_{2.13 * 10}6 Temperature 1.88 14.24 28.86 Dewpoint Tem-perature -1.05 9.65 20.92 Wind speed 0 3.57 11.80 Relative canopy saturation 0 0.5 1.0

(9)

8

3. Results

3.1 Model output and validation

3.1.1 Spatial analysis

Running the interception module with all data for the Netherlands, the main dynamics seem to be in order. As can be seen in figure 1, the canopy storage gets higher in locations where forests are, such as the Veluwe, or the forests in Noord-Brabant.

Figure 1: Canopy storage over the Netherlands on day 3 hour 8 clearly shows high canopy storage in forested areas.

Also, as can be seen in figure 2, when canopy storage is around the same value over the Netherlands and net radiation starts coming in, the area around the Randstad and Amsterdam produces a lot more latent heat, as the canopy is relatively more saturated.

Figure 2: Canopy storage and latent heat from interception evaporation on day 4, hour 6 shows the areas with a lower canopy storage capacity.

Figure 2 also shows an unfortunate side-effect from using different datasets and using them as input for the same model. Around the Dutch borders, there appears to be a higher latent heat flux than for the rest of the country. This is not because of an actual physical process, but because the precipitation and land-use datasets do not align perfectly, thereby rendering the model less reliable.

(10)

9

3.1.2 Model Dynamics

Figure 3 displays the model dynamics for a randomly chosen grid cell in the overall model running from 16.00 on the second day until 23.00 the day after. The period displayed was selected since this clearly shows how the interception module functions and a visualisation over the full runtime of the model would not be as clear. At first, Interception takes place, increasing canopy storage. After this, at around 5.00 in the morning, radiation starts coming in, increasing the evaporation and subsequently reducing canopy storage. Finally, as evaporation reaches its peak, the canopy storage becomes so low that it be-comes the limiting factor in evaporation, slowly reducing until the canopy is empty.

Figure 3: Modelled interception dynamics for a single grid cell

3.2 Sensitivity analysis

In the sensitivity analysis, the results of which are displayed in figures 4 and 5. It appears that dewpoint temperature and wind speed are not significant in any way when varied with all the other values at their mean. Precipitation and radiation, however, are significantly different. Given the fact that the flow varia-bles of interception and evaporation are directly driven by precipitation and radiation, respectively, whereas all other parameters can be seen as modifiers of these rates, but do not explicitly drive these processes.

Interestingly enough, although canopy storage capacity has a strong impact on the canopy storage itself, its influence on the latent heat, although present, is not significant. However, this is likely this is due to the way initial saturation was implemented, as it has the value of half the canopy storage capacity. If this is set at an absolute amount, assuming it does not strongly exceed the storage capacity, its significance quickly fades.

It can be seen that temperature and canopy cover fraction do impact the canopy storage and latent heat processes, but the canopy cover fraction does not have a big enough influence to be called significant. Temperature, on the other hand, has a strong impact at its highest value on the canopy storage. Indicat-ing the importance of usIndicat-ing accurate temperature data, particularly in the extremes.

(11)

10

Figure 4: P-values for the effect of different parameter values on canopy storage.

(12)

11

3.3 Interception’s influence on the SHF

The t-test for the impact of the interception module on SHF indicates p-values below 0.02, which would mean its influence is significant (see figure 6). However, there is one specific outlier at time step 457 (day 20, hour 1). But looking at the spatial results for this time step (figure 7), there do not seem to be any major differences between it and all other time steps. This analysis will need to be investigated when the model has been completely integrated, but for now it suffices to say interception, and interception evaporation, has a significant impact on the sensible heat flux.

Figure 6: P-values of the pairwise t-test of radiation and calculated sensible heat flux

(13)

12

4. Discussion

4.1 Results

Unfortunately, the matrices for the land-use parameters and precipitation are not fully aligned, despite both being in ETRS 89 and having been exported to ascii using the same methods. As a result of this, the model can show extreme values for latent heat. This may have impacted the analysis of interception on the sensible heat flux, because these extreme differences might cause the entire sample to skew away from the incoming radiation. It was tried to assign these cells a no-data-value, but this did not appear to have any effect. Figure 8 shows the discrepancy in data availability. A part of this discrepancy can be explained by open terrain without any canopy cover fraction, but the edges are likely the result of a pro-jection error. Because of this error, the model is unlikely to be very accurate on the small scale intended for this research, as values are likely shifting from their original point of measurement.

However, both the single grid cell and spatial dynamics show expected behaviour. Therefore, it can be concluded that, with a better dataset as input, it will function well enough to predict the influence of interception on the sensible heat flux at a high resolution.

Figure 8a: Locations where there is precipitation data, figure 8b: locations where there is land-use data, but no land-use data but no precipitation data.

When applying the interception module from this research to other areas or topics, it is important to know which values take most priority to get from accurate measurements. From the results, it is shown that radiation and precipitation are the most important parameters for the calculation of canopy storage and latent heat. Temperature can also have a strong influence, whereas wind speed and dewpoint tem-perature appear to be negligible in their impact. Wind speed, or the aerodynamic resistance which it influences, may in fact be more important than this analysis suggests, as the equation used for aerody-namic resistance is not quite certain.

Also, the sensitivity analysis as conducted here may not be fully representative, as canopy cover fraction, for example, may have a stronger impact in situations or areas with more, or more intense, precipitation. For this, a more elaborate analysis of the model sensitivity may be required.

(14)

13

4.2 Suggestions for improvement

Apart from the uncertainty around the aerodynamic resistance and the role it plays in interception, sev-eral other factors may be important to include in future research. In this case, interception evaporation was treated as open water evaporation using the Penman-Monteith equation. However, potential evap-oration can be much higher than open water evapevap-oration as a result of advection (van Dijk et al., 2015). In fact, during summer, advection can form over half of the available energy (Ringgaard et al., 2014), making it a significant factor in potential and actual evaporation.

In addition, most of the values in the current model have been the result of interpolation. Although val-ues are not likely to differ highly within different valval-ues, it would naturally be more reliable when based on more accurate measurements than partial degrees.

There is hardly any data available specifically concerning interception evaporation over the whole of the Netherlands. As a result of this, it is not possible to verify the model results using actual data. Such a quantitative analysis would be very interesting, as looking at the difference between modelled results and actual data could indicate for which areas or circumstances the model is less accurate and where the model can be improved. It has not yet been possible to integrate the complete SHF-model over the course of this research, but when this has been done, it will also provide the possibility to analyse its output against actual data.

Finally, for more detailed model results, more computational power will need to be made available. Just this interception model took over 16 hours to model 8 weeks. If the evapotranspiration and ground heat flux models were added, this can easily take several days on an above-average computer. For more spe-cific purposes, such as the previously described UVA-Bits program, it may be advisable to use this model only for a pre-selected area and timespan, as reducing the scope of the model can save a lot of computa-tion time.

5. Conclusion

In conclusion, although the land-use and precipitation dataset do not align, the model shows realistic dynamics which are significant for the effect of interception processes on the sensible heat flux. Thereby, it serves as a proof of concept for high-resolution modelling based on limited high-resolution data. This research has found that radiation, precipitation and temperature are the most important factors to get accurate measurements from, in order to effectively model interception processes. This means that these temporal resolutions should be improved in order to get even more detailed results than hourly time steps. Meanwhile land-use at 100x100m is used as a common denominator for the spatial dimen-sion, which serves to show a significant contrast in the modelled spatial dynamics. As land-use is more variable with distance than most atmospheric processes, this would be the parameter by which the model’s spatial resolution can be increased even further.

(15)

14

6. Literature

- Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998). Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56. FAO, Rome, 300(9), D05109.

- Bohrer, G., Brandes, D., Mandel, J. T., Bildstein, K. L., Miller, T. A., Lanzone, M., ... & Tremblay, J. A. (2012). Estimating updraft velocity components over large spatial scales: contrasting migration strategies of gold-en eagles and turkey vultures. Ecology Letters, 15(2), 96-103.

- Bonan, G. B., Levis, S., Kergoat, L., & Oleson, K. W. (2002). Landscapes as patches of plant functional types: An integrating concept for climate and ecosystem models. Global Biogeochemical Cycles, 16(2).

- Deguchi, A., Hattori, S., & Park, H. T. (2006). The influence of seasonal changes in canopy structure on in-terception loss: application of the revised Gash model. Journal of Hydrology, 318(1), 80-102.

- Donatelli, M., Bellocchi, G., Confalonieri, R., & Habyarimana, E. (n.d.). Aerodynamic resistance. Retrieved June 27, 2016, from http://agsys.cra-cin.it/tools/evapotranspiration/help/Aerodynamic_resistance.html

- Finch, J. W., Hall, R. L., & Environment Agency, Bristol (United Kingdom);. (2001). Estimation of Open Wa-ter Evaporation: A Review of Methods. Environment Agency.

- Haylock, M. R., Hofstra, N., Klein Tank, A. M. G., Klok, E. J., Jones, P. D., & New, M. (2008). A European daily high-resolution gridded data set of surface temperature and precipitation for 1950–2006. Journal of Geo-physical Research: Atmospheres, 113(D20).

- Hewitson, B. C., & Crane, R. G. (1996). Climate downscaling: techniques and application. Climate Re-search, 7(2), 85-95.

- Linhoss, A. C., & Siegert, C. M. (2016). A comparison of five forest interception models using global sensi-tivity and uncertainty analysis. Journal of Hydrology, 538, 109-116.

- Liu, S., Lu, L., Mao, D., & Jia, L. (2007). Evaluating parameterizations of aerodynamic resistance to heat transfer using field measurements. Hydrology and Earth System Sciences Discussions, 11(2), 769-783. - Mandel, J. T., Bohrer, G., Winkler, D. W., Barber, D. R., Houston, C. S., & Bildstein, K. L. (2011). Migration

path annotation: cross-continental study of migration-flight response to environmental condi-tions. Ecological Applications, 21(6), 2258-2268.

- McNaughton, K. G., & Jarvis, P. G. (1983). Predicting effects of vegetation changes on transpiration and evaporation. Water deficits and plant growth, 7, 1-47.

- OSGeo/proj.4. (n.d.). Retrieved July 02, 2016, from https://github.com/OSGeo/proj.4/wiki

- Ringgaard, R., Herbst, M., & Friborg, T. (2014). Partitioning forest evapotranspiration: Interception evapo-ration and the impact of canopy structure, local and regional advection. Journal of Hydrology, 517, 677-690.

- Rutter, A. J., Kershaw, K. A., Robins, P. C., & Morton, A. J. (1972). A predictive model of rainfall interception in forests, 1. Derivation of the model from observations in a plantation of Corsican pine. Agricultural Me-teorology,9, 367-384.

- Savenije, H. H. (2004). The importance of interception and why we should delete the term evapotranspira-tion from our vocabulary. Hydrological Processes, 18(8), 1507-1511.

- Tolk, J. A., Howell, T. A., Steiner, J. L., Krieg, D. R., & Schneider, A. D. (1995). Role of transpiration suppres-sion by evaporation of intercepted water in improving irrigation efficiency. Irrigation Science, 16(2), 89-95. - van Dijk, A. I., Gash, J. H., van Gorsel, E., Blanken, P. D., Cescatti, A., Emmel, C., ... & Montagnani, L. (2015). Rainfall interception and the coupled surface water and energy balance. Agricultural and Forest Meteorol-ogy, 214, 402-415.

- Wallace, J., & McJannet, D. (2006). On interception modelling of a lowland coastal rainforest in northern Queensland, Australia. Journal of Hydrology,329(3), 477-488.

(16)

15

Appendix A: Interception Module:

function [CanopyStore, LatentHeat] = SHF_Interception(a, c, ...

dt, Precipitation, CanopyStore, RNet, Temperature, DewTemperature, ... Windspeed) % a, c

%% Interception Module for SHF-Model % Created by Maarten Mol

% Bachelorthesis Earth Sciences 2016 % Inputs:

% a: Canopy cover fraction % [-] % c: Canopy storage capacity % [mm] % dt: Timestep % [Hour] % Precipitation: Incoming rainfall rate % [mm/hour] % RNet: Net incoming radiation % [MJ/m2] % Temperature: Actual air temperature % [oC] % DewTemperature: Dewpoint Temperature % [oC] % CanopyStore: Stored water in canopy % [mm] %

% Outputs:

% CanopyStore: Amount of water stored in Canopy % [mm] % Throughfall: Amount of water passed on to other modules % [mm] % LatentHeat: Energy used in the evaporation process % [J] %% constants: AirDensity = 1.225; % [kg m-3] AirSpecHeat = 1.005e3; % [J kg-1 K-1] WaterLatHeat = 2450e3; % [J kg-1] MWratio = 0.622; % []

% Possibly changing constants

AtmosphericPressure = 101.3; % [kPa] % [mm] %% Calculating other values

PsychrometricConstant = (AirSpecHeat * AtmosphericPressure)... / (WaterLatHeat * MWratio); % [kPa K-1]

% Vapor Pressure and air moisture content

ActualVaporPressure = 0.6108 * exp(17.27 * DewTemperature ./ ... % [kPa]

(237.3 + DewTemperature));

SaturatedVaporPressure = 0.6108 * exp(17.25 * Temperature ./ ... % [kPa]

(237.3 + Temperature));

% Slope of saturation vapor pressure / temperature

VaporSlope = (4098 * (0.6108 * exp( (17.27 * Temperature) ./... % [kPa K-1]

(Temperature + 237.3) ))) ./ (Temperature + 237.3).^2;

(17)

16 Windspeed(Windspeed < 0.5) = 0.5; K = 208; % Open air AerodynamicResistance = K ./ Windspeed; % [s m-1] %% Interception Modelling % Evaporation

PotEvap = 1 ./ WaterLatHeat .* (RNet + (AirDensity .* AirSpecHeat .* ... % [mm h-1] [kg m-2 h-1]

(SaturatedVaporPressure - ActualVaporPressure) ./ AerodynamicResistance))... ./ (VaporSlope + PsychrometricConstant);

for i = 1:size(CanopyStore, 1) for j = 1:size(CanopyStore, 2) % Interception

Interception(i,j) = min(c(i,j) - CanopyStore(i,j), ... % [mm]

a(i,j) .* Precipitation(i,j) * dt); end

end

%% Integration and Output

CanopyStore = (CanopyStore + Interception) ./ (1 + PotEvap ./ c * dt); % [mm]

LatentHeat = ((CanopyStore ./ c .* PotEvap) .* WaterLatHeat) / 3600; % [W m-2]

(18)

17

Appendix B: First coordinate conversion for reading in arcmap:

%% conversion coordinates:

% Script to read coordinates from output proj4 file,

% reading precipitation data from the knmi RAD_NL25 dataset % writing to xlsfiles small enough to read in ARCGIS

% By Maarten Mol - 27-06-2016 clear all

fileID = fopen('wgsclipcoords.txt');

[E] = textscan(fileID, '%s %s', 'delimiter', '\t'); fclose(fileID);

East = E{1}; North = E{2}; clear E fileID

%% Reading and calculating coordinates for i = 1:length(East) %% Name adjustment for j = 1:length(East{i}) if strcmp(East{i}(j),'d') East{i}(j)= '°'; EDeg(i,:) = str2double(East{i}(1:j-1)); for k = j+1:length(East{i}) if strcmp(East{i}(k),East{2}(5)) EMin(i,:) = str2double(East{i}(j+1:k-1)); for l = k+1:length(East{i}) if strcmp(East{i}(l),'"') ESec(i,:) = str2double(East{i}(k+1:l-1)); end end end end end end for j = 6:length(North{i}) if strcmp(North{i}(j-4:j),'0.000') North{i}(j-4:j)= []; end end for j = 1:length(North{i}) if strcmp(North{i}(j),'d') North{i}(j)= '°'; NDeg(i,:) = str2double(North{i}(1:j-1)); for k = j+1:length(North{i}) if strcmp(North{i}(k),East{2}(5)) NMin(i,:) = str2double(North{i}(j+1:k-1)); for l = k+1:length(North{i}) if strcmp(North{i}(l),'"') NSec(i,:) = str2double(North{i}(k+1:l-1)); end end

(19)

18 end end end end end

%% Final Coordinate Values

Easting = EDeg + EMin./60 + ESec./3600; Northing = NDeg + NMin./60 + NSec./3600; Coordinates = [Easting,Northing];

clear EDeg EMin ESec NDeg NMin NSec %% Precipitation reading

cd('C:\Users\Maarten\Documents\Bachelorproject\Precipitation_final_coords')

MyFolderInfo = struct2table(dir('Precipitation')); Names = table2array(MyFolderInfo(3:end,1));

clear MyFolderInfo

cd('Precipitation')

Time = 1;

for i = 1:length(Names)

Data = ncread(Names{i}, 'image1_image_data');

Precipitation5min(:,:,mod(i,12)+1) = Data(200:550,250:600); if mod(i,12)==0 for j = 1:351 for k = 1:351 PrecipitationHour(j,k,Time) = sum(Precipitation5min(j,k,1:end)); end end Time = Time + 1; end end loc = 1; for i = 1:1464 for j = 1:351 for k = 1:351 ReadPrecip(loc,1,i) = PrecipitationHour(k,j,i); loc = loc+1; end end loc = 1; end for i = 1:1464 PrecipRead(:,i) = ReadPrecip(:,:,i); end write1 = [Coordinates,Precip(:,1:1464/6)]; write2 = [Coordinates,Precip(:,1464/6+1:1464*2/6)]; write3 = [Coordinates,Precip(:,1464*2/6+1:1464*3/6)]; write4 = [Coordinates,Precip(:,1464*3/6+1:1464*4/6)]; write5 = [Coordinates,Precip(:,1464*4/6+1:1464*5/6)]; write6 = [Coordinates,Precip(:,1464*5/6+1:1464)]; xlswrite('Precipitation1.xlsx',write1)

xlswrite('Precipitation2.xlsx',write2) xlswrite('Precipitation3.xlsx',write3) xlswrite('Precipitation4.xlsx',write4) xlswrite('Precipitation5.xlsx',write5) xlswrite('Precipitation6.xlsx',write6)

(20)

19

Appendix C: Python analysis scripts

XY Data – Point to raster

import arcpy import array import string import os import math arcpy.env.workspace = r"C:\Users\Maarten\Documents\ArcGIS\Bachelorproject.gdb" arcpy.env.scratchWorkspace = r"C:\Users\Maarten\Documents\BachelorFiles.gdb"

inFeatures = ["\PData1", "\PData2", "\PData3", "\PData4","\PData5", "\PData6"]

for i in range(0,1464): loc = int(ceil(i / (1464/6)))

In = r"C:\Users\Maarten\Documents\ArcGIS\Bachelorproject.gdb"+inFeatures[loc] Value = "Precip"+str(i+1)

Out = "C:\Users\Maarten\Documents\ArcGIS\BachelorFiles.gdb\Precip"+str(i+1) arcpy.PointToRaster_conversion(In, Value, Out, "MEAN", "", 0.015)

Raster to ascii:

import arcpy arcpy.env = r"C:\users\Maarten\Documents\ArcGIS\BachelorFiles.gdb" for i in range(1,1465): In = r"C:\users\Maarten\Documents\ArcGIS\BachelorFiles.gdb\Precip"+str(i) Out = r"C:\users\Maarten\Documents\ArcGIS\asciis\Precip"+str(i)+".asc" arcpy.RasterToASCII_conversion(In, Out)

(21)

20

Project raster:

import arcpy arcpy.env = r"C:\users\Maarten\Documents\ArcGIS\BachelorFiles.gdb" for i in range(1,1465): In = r"C:\users\Maarten\Documents\ArcGIS\BachelorFiles.gdb\Precip"+str(i) Out = r"C:\users\Maarten\Documents\ArcGIS\BachelorFiles2.gdb\Precip"+str(i) Proj = r"C:\users\Maarten\Documents\ArcGIS\Default.gdb\WGS84LandUse_Clip_PolygonToR" arcpy.ProjectRaster_management(In, Out, Proj,\

(22)

21

Appendix D: Matlab script to load the precipitation from asci files

%% Conversion from ascii % Maarten Mol - 28-06-2016

cd('C:/Users/Maarten/Documents/ArcGIS/Asciis')

for i = 1:1464

Precip(:,:,i) = textread(['precip',num2str(i),'.asc'],

'%357s','headerlines',6); end for i = 1:1464 for j = 1:78183 if strcmp(Precip{j,:,i},'-9999') Precip{j,:,i} = NaN; continue end for k = 1:length(Precip{j,:,i}) if strcmp(Precip{j,:,i}(k),',') Precip{j,:,i}(k) = '.'; end end end end for i = 1:1464 loc = 1; for j = 1:219 for k = 1:357 Precipitation(j,k,i) = str2double(Precip{loc,:,i}); loc = loc + 1; end end end

(23)

22

Appendix E: Matlab analysis script for sensitivity analysis

%% Analysis script for the interception module

% Maarten Mol - Bachelorthesis Earth Science 2015/2016 % University of Amsterdam %% Initialisation close all clear all clc %% loading data % Landuse parameters load('Smax.mat'); load('Kvg.mat'); a = 1 - FracOpenK; a(a == 0) = NaN; c = Smax; c(c == 1) = NaN; % Atmospheric data load('Precipitation.mat'); RNet = ncread('Rn.nc','Rn')*10^6; Temperature = ncread('T2m_1h.nc', 'T2m_1h') - 273.15; DewTemperature = ncread('D2m_1h.nc', 'D2m_1h')-273.15; Windspeed = ncread('U2_1h.nc','U2_1h') * 0.748;

Temperature(:,:,1) = []; % remove first (erroneous) value. DewTemperature(:,:,1) = []; % remove first (erroneous) value.

%% defining parameter values and clearing unnecessary variables dt = 1; P1 = [nanmin(a(:)),nanmean(a(:)),nanmax(a(:))]; % a P2 = [nanmin(c(:)),nanmean(c(:)),nanmax(c(:))]; % c P3 = [nanmin(Precipitation(:)),nanmean(Precipitation(:)), nan-max(Precipitation(:))]; % Precipitation

P4 = [nanmin(RNet(:)), nanmean(RNet(:)), nanmax(RNet(:))]; % Radiation

P5 = [nanmin(Temperature(:)), nanmean(Temperature(:)), nan-max(Temperature(:))]; % Temperature

P6 =

[nan-min(DewTemperature(:)),nanmean(DewTemperature(:)),nanmax(DewTemperature(:))]; % Dewtemperature

P7 = [nanmin(Windspeed(:)), nanmean(Windspeed(:)), nanmax(Windspeed(:))]; % Wind speed

P8 = [0, 0.5, 1]; % Ini-tial saturation fraction of canopy

% preallocating space for speed CanopyStore(1:6561,1:11) = 0; LatentHeat(1:6561,1:11) = 0;

(24)

23 % clearing unnecessary variables

clear Precipitation Smax FracOpenK a c RNet Temperature DewTemperature Wind-speed

%% Iterating for each value of each parameter loc = 1; for I1 = 1:3 % a for I2 = 1:3 % c for I3 = 1:3 % Precip for I4 = 1:3 % Radiation for I5 = 1:3 % Temperature for I6 = 1:3 % Dewtemperature for I7 = 1:3 % windspeed

for I8 = 1:3 % CanopyStore saturation for Time = 2:10 % Hours

% define first timestep as a fraction of Canopy Storage Capacity CanopyStore(loc,1) = P2(I2) * P8(I8);

[CanopyStore(loc, Time), LatentHeat(loc, Time)] ...

= SHF_Interception(P1(I1), P2(I2), dt, P3(I3), CanopyStore(loc, Time-1),...

P4(I4), P5(I5), P6(I6), P7(I7));

end Values(loc,:) = [I1,I2,I3,I4,I5,I6,I7,I8]; loc = loc+1; end end end end end end end end

%% Taking unique values while keeping the others constant CanopyStoreP1 = CanopyStore([1094, 3281, 5468],:); LatentHeatP1 = LatentHeat([1094, 3281, 5468],:); CanopyStoreP2 = CanopyStore([2552, 3281, 4010],:); LatentHeatP2 = LatentHeat([2552, 3281, 4010],:); CanopyStoreP3 = CanopyStore([3038, 3281, 3524],:); LatentHeatP3 = LatentHeat([3038, 3281, 3524],:); CanopyStoreP4 = CanopyStore([3200, 3281, 3362],:); LatentHeatP4 = LatentHeat([3200, 3281, 3362],:); CanopyStoreP5 = CanopyStore([3254, 3281, 3308], :); LatentHeatP5 = LatentHeat([3254, 3281, 3308], :);

(25)

24 CanopyStoreP6 = CanopyStore([3272, 3281, 3290],:); LatentHeatP6 = LatentHeat([3272, 3281, 3290],:); CanopyStoreP7 = CanopyStore([3278, 3281, 3284], :); LatentHeatP7 = LatentHeat([3278, 3281, 3284],:); CanopyStoreP8 = CanopyStore(3280:3282, :); LatentHeatP8 = LatentHeat(3280:3282, :); %% Analysis [h1(1,:), p1(1,:)] = ttest(CanopyStoreP1', mean(CanopyStoreP1(:))); [h1(2,:), p1(2,:)] = ttest(CanopyStoreP2', mean(CanopyStoreP2(:))); [h1(3,:), p1(3,:)] = ttest(CanopyStoreP3', mean(CanopyStoreP3(:))); [h1(4,:), p1(4,:)] = ttest(CanopyStoreP4', mean(CanopyStoreP4(:))); [h1(5,:), p1(5,:)] = ttest(CanopyStoreP5', mean(CanopyStoreP5(:))); [h1(6,:), p1(6,:)] = ttest(CanopyStoreP6', mean(CanopyStoreP6(:))); [h1(7,:), p1(7,:)] = ttest(CanopyStoreP7', mean(CanopyStoreP7(:))); [h1(8,:), p1(8,:)] = ttest(CanopyStoreP8', mean(CanopyStoreP8(:))); [h2(1,:), p2(1,:)] = ttest(LatentHeatP1', mean(LatentHeatP1(:))); [h2(2,:), p2(2,:)] = ttest(LatentHeatP2', mean(LatentHeatP2(:))); [h2(3,:), p2(3,:)] = ttest(LatentHeatP3', mean(LatentHeatP3(:))); [h2(4,:), p2(4,:)] = ttest(LatentHeatP4', mean(LatentHeatP4(:))); [h2(5,:), p2(5,:)] = ttest(LatentHeatP5', mean(LatentHeatP5(:))); [h2(6,:), p2(6,:)] = ttest(LatentHeatP6', mean(LatentHeatP6(:))); [h2(7,:), p2(7,:)] = ttest(LatentHeatP7', mean(LatentHeatP7(:))); [h2(8,:), p2(8,:)] = ttest(LatentHeatP8', mean(LatentHeatP8(:)));