HIGH RESOLUTION MODELLING OF EVAPOTRANSPIRATION AND
SENSIBLE HEAT FLUX IN THE NETHERLANDS
Explorative research for a global thermal convection model for predicting bird
behaviour
ABSTRACT -‐ Many bird species have evolved traits to increase flight efficiency to help cover the long migration routes taken in search of resources required for breeding or survival. For birds migrating over land, soaring flight is the principle adaptation seen in a number of species. Soaring flight is supported by rising air, which is largely caused by thermal lift. For this reason, a model that accurately describes thermal lift through calculated proxies (SHF, w*) at a high spatial and temporal resolution and across large spatial domains contributes significantly to the study of migratory birds. As thermal lift is part of the energy balance of the surface, it is dependent on atmospheric and surface conditions. Latent heat plays a major role in the energy balance of the Earth’s surface and consequently in the occurrence of atmospheric conditions suitable for soaring flight. As the latent heat flux is largely dependent on soil moisture content, this research sets out to incorporate existing modelled –and measured hydrological data in the Netherlands into a thermal lift model to produce output of high spatial and temporal resolution. Furthermore, performance of the model under different hydrological input is evaluated. This is done in order to determine the spatial and temporal scale at which hydrological data is required to accurately model thermal lift. The results from this research fit in a larger effort to successfully model migrating bird behaviour.
BSc. Thesis by: Bart Sweerts
Supervisors: Prof. dr. ir. Willem Bouten & dr. Judy Shamoun-‐Baranes Reviewer: dhr. Wouter Vansteelant MSc.
Date: July 7th 2016
Table of contents
1. Introduction ... 2
2. Theoretical Framework ... 4
Energy balance and thermal convection ... 4
Joint model ... 5
3. Methods ... 5
Data collection ... 5
Data processing ... 6
Evapotranspiration model description ... 7
Parameterization and computation of soil moisture ... 9
Sensitivity analysis ... 11
Analysis area and period ... 12
4. Results ... 12
Model output with data from NHI ... 12
Model output with data from PCR-‐GLOBWB ... 14
Variable sensitivity ... 16
5. Discussion ... 17
Model output with data from NHI ... 17
Model output with data from PCR-‐GLOBWB ... 18
Variable sensitivity ... 18
6. Conclusion ... 19
5. References ... 20
Appendix ... 23
Appendix 1 – Variable conversion ... 23
Appendix 2 – MATLAB evapotranspiration script ... 24
Appendix 3 – MATLAB net solar radiation script ... 28
1. Introduction
Most birds display migratory behaviour on varying temporal and spatial scales to exploit resources for breeding and survival (Newton, 2008). To deal with these relatively long migration routes, birds have evolved a range of traits to increase flight efficiency (Hedenström, 2002, Norberg, 2012). For birds migrating over land, soaring flight is the principle adaptation seen in a number of species (Hedenström, 1993). As some of these birds are not able to sustain flapping flight for an extended period of time, they are dependent on thermal lift that supports soaring flight. Consequently, much of migrating birds behaviour can be linked to atmospheric conditions (Vansteelant et al., 2014). For this reason, studies on avian migratory behaviour are intrinsically linked to research on atmospheric dynamics (Kunz et al., 2008).
Extended work has been done on the influence of atmospheric dynamics on animal migration. This includes a considerable amount of models that simulate migration of a variety of avian animals based on meteorological data (Shamoun-‐Baranes et al., 2010). Shamoun-‐Baranes et al. state that these models are either data-‐driven or concept-‐driven, but always include atmospheric condition data as input. Furthermore, they note that data quantity and quality, the
model framework and communication and collaboration with other researchers or institutes are important factors to consider when integrating atmospheric conditions into migration models.
Tracking data for migrating birds is collected on high spatial and temporal resolution and across large spatial domains (Bouten et al. 2013). Furthermore, atmospheric conditions and the resulting bird behaviour show large gradients over time and space. For this reason, developing a high-‐resolution thermal convection model that covers large areas contributes significantly to bird behaviour research by linking bird-‐tracking data to atmospheric conditions.
To model thermal convection with the use of existing data, the energy balance of the Earth’s surface has to be considered. When solar radiation meets the Earth’s surface, the energy is roughly divided into four components:
1 − 𝑟 𝑅 ↓ = 𝑅𝑛𝑙 ↑ + 𝑆𝐻𝐹 + 𝜆𝐸 + 𝐺 (1)
With incoming solar radiation (R↓), albedo (r), net longwave radiation (Rnl↑), sensible heat flux (SHF), latent heat (λE) and soil heat exchange by conduction (G) (Bonan, G., 2002). SHF has in combination with boundary layer height (BLH), a proxy for thermal depth, been used to calculate the convective scale velocity (w*), which acts as a proxy for thermal lift velocity (Bohrer et al, 2011, Shamoun-‐Baranes et al, 2016).
In moist conditions like the Netherlands, latent heat is the most important component of the energy balance of the surface, as solar radiation is primarily used to heat and evaporate moisture (Bonan, 2002). Since latent heat flux is largely dependent on soil water content, hydrological models provide important input for the modelling of thermal lift. The soil moisture content parameterization currently used in meteorological models often displays high degrees of uncertainty. For this reason, there is a need for a method that uses existing hydrological data from varying sources to accurately model thermal lift parameters used to predict soaring bird behaviour.
The aim of the research is twofold. First of all, it aims to effectively calculate evapotranspiration rates to determine SHF on a high spatial and temporal resolution and across a large domain. For this, daily hydrological data will be used from the complex and high spatial resolution (1km) Nationaal Hydrologisch Instrumentarium (NHI). Second, it aims to determine the possibilities for a global sensible heat flux model using the coarser 5’ (geographical minutes) PCRaster Global Water Balance model (PCR-‐GLOBWB). In both cases, emphasis is placed on the impact of various input variables, parameters and model mechanics to determine what is required for an accurate high-‐resolution sensible heat flux model.
2. Theoretical Framework
Energy balance and thermal convection
There exists a constant exchange of energy between the Earth’s surface and the atmosphere (Figure 1). As net radiation reaches the Earth’s surface, it is divided in the latent-‐heat flux (lE), soil heat flow (G) and the sensible heat flux (H or SHF).
Fig 1. Surface energy budget. Obtained from USGS (2016)
The latent heat flux is the energy that is used to evaporate water and can be further divided in the soil evaporation, transpiration and interception evaporation flux. The energy flow directly influencing thermal lift is SHF, which can be quantified by calculating and subtracting the other energy fluxes from the absorbed solar irradiation.
Sensible heat is the energy that is used to raise the temperature of the air (Stull, 1988). As warm air is lighter than cold air, the airs buoyancy rises resulting in thermal lift. However, there are more factors that influence soaring flight behaviour. An important factor is the boundary layer height (BLH). The BLH is the border of the lowest layer of the atmosphere that is heavily influenced by surface characteristics. Energy fluxes concerning the Earth are in large confined to this layer. Above this layer is the free atmosphere, where air moves with little to no influence from the surface. Shamoun-‐Baranes et al. (2016) find that BLH provides a good proxy for where soaring birds fly and BLH is a useful predictor of ground speed in soaring birds engaged in migration (Vansteelant et al. 2015). However, other studies propose to use a proxy for thermal lift velocity provided by a scaling coefficient, known as convective scale velocity (w*) (Stull, 1988; Bohrer et al, 2011). An adapted calculation by Shamoun-‐Baranes et al. (2016) includes several atmospheric parameters (eq. 2).
𝑤∗ = !"#$ ! . !"# !!! ! ! (2)
With gravitational constant (g), boundary layer height (BLH), temperature (T), sensible heat flux (SHF), density of air (ρ) and specific heat of air (Cp). In this research, the main output is SHF. However, in combination with atmospheric variables different proxies like w* can be calculated with relative ease.
Joint model
To determine values for SHF, the other energy fluxes in equation 1 have to be quantified. In order to do this, a model will be written that contains all parts of the energy balance equation. The computed energy fluxes are latent heat of evapotranspiration, latent heat of interception evaporation and soil heat flow. Additionally, the parameterization of the modules is computed separately on a 100m spatial resolution. Figure 2 provides a general overview of the total model. The model will consist of a main script that contains the initialization and calls the other script to calculate the energy fluxes. The initialization consists of data loading and processing. From there, the different energy fluxes are calculated in separate modules. These are returned to the main script to produce SHF, which can then be used to calculate other thermal lift proxies.
Fig 2. General overview of the total model. The module discussed in this research in bold
3. Methods Data collection
For the computation of evapotranspiration, data from a variety of sources was obtained (Table 1). First off, the ERA-‐interim reanalysis provides the atmospheric conditions. Secondly, two hydrological models, the Nationaal Hydrologisch Instrumentarium (NHI) by Deltares (De Lange et al., 2014) and the PCRaster Global Water Balance model (PCR-‐GLOBWB) by the physical
geography department of the University of Utrecht (Van Beek et al., 2008), provide the soil moisture content. The hydrological models are discussed in more detail in the methods. Thirdly, land-‐use maps were acquired from the CORINE land cover programme. Finally, radiation data was acquired from The Satellite Application Facility on Climate Monitoring (CM-‐SAF).
Table 1.
List of used data with source and resolution
Source Data used
Spatial resolution Temporal resolution
ERA-‐interim (reanalysis) by ECMWF
-‐ Temperature at 2m -‐ Dewpoint temperature -‐ Windspeed at 10m -‐ Total cloud cover
0.0125° (degrees) 6-‐hourly National Hydrologisch Instrumentarium (NHI) by Deltares
-‐ Soil moisture content -‐ Soil moisture deficit -‐ Rootzone depth -‐ Soil properties
250m Daily
PCR-‐GLOBWB (model) by University of Utrecht
-‐ Soil moisture content 5’ (geographical
minutes)
Daily CORINE land cover
programme -‐ Land-‐use 100m -‐
The Satellite Application Facility on Climate Monitoring (CM SAF)
-‐ Surface incoming shortwave radiation 0.05° (degrees) 1-‐hourly
Data processing
Some of the data required conversion or processing to be used in the model calculations.
Formulas related to evapotranspiration rates use net solar radiation (Rn). To determine net radiation at the surface (Rn), equation 3-‐5 are used.
𝑅!= (1 − 𝑎)𝑅!− 𝑅!" (3) 𝑅!" = 𝜎 𝑇! 0.34 − 0.14 𝑒! 1.35 𝑅! 𝑅!"− 0.35 (4) 𝑅!"= 0.75 𝑅! (5)
With net radiation at the surface (Rn), albedo (a), surface incoming shortwave radiation (Rs), net longwave radiation (Rnl), hourly Stefan-‐Boltzman constant (𝜎 = 2.043 10-‐10 MJ m-‐2 hour-‐1), average hourly temperature (T), actual saturation vapour pressure (ea) (equation 9), clear sky solar radiation (equation 5) and extra-‐terrestrial radiation (Ra). Hourly values for extraterrestrial radiation were determined according to the method provided by Allen et al. (2006) (Annex 1, eq. 1 -‐ 4).
Formulas for the computation of potential evapotranspiration (ET0) require windspeed at 2m. Therefore, ERA-‐interim windspeed values at 10m are converted using the FAO conversion curve (Allen et al., 1998) (Annex 1, figure 1).
A third issue arose in the computation of evaporation rates. The hydrological models used in this research do not separately compute the soil moisture content in the top 100mm of the soil. Therefore, this soil moisture content is estimated by determining if the moisture deficit exceeds 100mm. Any soil moisture present in the top 100mm is then used to compute top-‐layer soil moisture content.
To be able to use all variables in calculations, interpolation in space and time was required. For this, the simplest technique of ‘nearest neighbour’ interpolation was used. This technique was chosen because it makes the least assumptions regarding the newly formed query points and their assigned value (Hartkamp et al., 1999).
Evapotranspiration model description
The latent heat flux module calculates latent heat flux (W m-‐2) by multiplying evapotranspiration with the latent specific heat of water. To determine evapotranspiration rates the module uses a crop coefficient – reference evapotranspiration (Kc -‐ ET0) procedure. A reference evapotranspiration (ET0) is computed for a reference crop like grass or alfalfa and is then multiplied by an empirical crop-‐specific coefficient (Kc) ranging from [0 to 1.4] to produce crop potential evapotranspiration (ETc). To calculate evapotranspiration rates, the American Society of Civil Engineers (ASCE) standardized reference evapotranspiration equation is used, which is an adapted version of the FAO-‐56 Penman-‐Monteith equation (Allen et al., 2005; Walter et al., 2001): 𝐸𝑇!= 0.408𝛥 𝑅!− 𝐺 + 𝛾𝑇 + 273 𝑈𝐶𝑛 ! 𝑒!− 𝑒! 𝛥 + 𝛾 1 + 𝐶𝑑 𝑈! 6
With slope vapour pressure (𝛥) (equation 7), net solar radiation on the surface (Rn), soil heat flux density (G), psychometric constant (𝛾), hourly mean temperature at 2m (T), hourly mean wind speed at 2m (U2), saturation vapour pressure (es) (equation 8) and actual vapour pressure (ea) (equation 9) (Allen et al., 1998).
𝛥 = 4098 (0.6108 𝑒 !".!" ! ! ! !"!.!) (𝑇 + 237.3)! (7) 𝑒! = 0.6108 𝑒 !".!" ! ! ! !"#.! (8)
𝑒! = 0.6108 𝑒 !!".!" !!"# ! !"#.!!"# (9)
With temperature at 2m (T) and dewpoint temperature at 2m (Tdew).
Furthermore, parameters Cn and Cd can be adjusted to account for different vegetation heights, nighttime calculations and to allow for 1-‐hour time-‐step calculations (Table 2) (Walter et al., 2001).
Table 2.
Values for Cn and Cd in equation 1. Adapted from Walter et al. (2001).
Calculation Time Step Short Reference ET0s Tall Reference ET0t
Cn Cd Cn Cd
Daily 900 0.34 1600 0.38
Hourly during daytime 37 0.24 66 0.25
Hourly during nighttime 37 0.96 66 1.7
However, Walter et al. (2001) state that short vegetation has an approximate height of 0.12m and tall vegetation an approximate height of 0.50m. As this method was specifically developed for agricultural needs, it does not take into account natural vegetation that can grow much higher than the tall vegetation class. Nonetheless, for natural vegetation tall reference parameters are used to compute potential evapotranspiration. Secondly, nighttime parameters are used whenever net radiation equals or is lower than 0.
To determine actual evapotranspiration (ETact), ET0 is multiplied by a crop coefficient Kc. However, determining values for Kc requires several steps. First off, values for Kc vary over time under the influence of the seasons and growing conditions. In the FAO dual crop coefficient method Kc is divided in a transpiration component (Kcb) and a soil evaporation component (Ke). Values for Kcb follow a seasonal curve with peak values occurring during the growing season (figure 2). Secondly, values of Kcb are reduced by multiplication with a water stress coefficient Ks (eq. 10).
𝐾!= 𝑇𝐴𝑊 − 𝐷!
1 − 𝑝 𝑇𝐴𝑊 (10)
Where TAW is total available water in the rootzone at field capacity, Dr is water shortage relative to field capacity and p is the fraction of TAW a crop can extract without suffering water stress. Values for p were obtained from the FAO (Allen et al., 1998).
Evaporation from the soil is predicted by estimating the amount of energy that is able to penetrate the canopy and reach the ground. This is done by subtracting the current vegetation
coefficient (Kcb) from the maximum value the vegetation coefficient can have at optimal growing conditions (Kc max) (eq. 11).
𝐾! = 𝐾! !"#− 𝐾!" (11)
When the soil is wet, evaporation is only limited by the empirically derived maximum value for total crop coefficient Kc max. When the soil moisture contents drop beneath ideal levels, a reducing coefficient Kr is introduced (eq. 12). For this coefficient, soil moisture content in the top 100mm of the soil is considered, from which evaporation occurs.
𝐾! = 𝑇𝐸𝑊 − 𝐷!,!!!
𝑇𝐸𝑊 − 𝑅𝐸𝑊 (12)
Where total evaporable water (TEW) is the amount of water in the top 100mm of soil at field capacity, readily evaporated water (REW) is the amount of water that can be extracted through evaporation before water stress and De,j -‐1 is the cumulative moisture depletion from the soil surface layer at the previous day. Combining the crop-‐coefficients and potential evapotranspiration rates eventually produces actual evapotranspiration ETc act in millimetres (eq. 13).
𝐸𝑇! !"# = 𝐾!𝐾!"+ 𝐾!𝐾! 𝐸𝑇!= 𝐾! !"#𝐸𝑇! (13)
To convert to an energy flux, the evapotranspiration [mm] is multiplied by the latent evaporation heat of water (l = 2257000 J kg-‐1) and divided by the amount of seconds in 1 hour (3600) to produce the latent heat flux (λE) [W m-‐2]. Additionally, the latent heat flux can be subtracted from the net solar radiation to estimate sensible heat flux (SHF).
Parameterization and computation of soil moisture
The NHI is used to calculate the impact of precipitation and other water flows on water levels in the Netherlands. In this fashion it is used as a decision making tool for water management strategies. Therefore, for the scope of this research the output from the NHI is considered adequate. However, since the NHI only spans the Netherlands, a different hydrological input is required for the model to be used on larger scales. For this, the global PCR-‐GLOBWB model is used. To evaluate whether hydrological input from the PCR-‐GLOBWB may be used for the purpose of high resolution modelling of evapotranspiration, the parameterization, resolution and accuracy of the soil moisture computation has to be determined.
This section starts with a description of the parameters used in the evapotranspiration model. After this, the computation of the soil moisture content in the NHI and PCR-‐GLOBWB is briefly discussed.
Parameters
Not only the soil moisture content, but also the parameterization of the NHI model is not available on a global scale. Therefore, it is important to determine the differences between the parameters used in the NHI and the parameters used in the global PCR-‐GLOBWB model. Table 3 contains the model-‐specific parameters used in the evapotranspiration model. Furthermore, it contains the parameters that were empirically derived on a 100 by 100m scale for this study. A more detailed description of their computation can be found in Kooij (2016).
Table 3.
Model-‐specific and 100 by 100m parameters. Model-‐
specific Parameter
NHI
Resolution NHI computation PCR-‐GLOBWB resolution PCR-‐GLOBWB computation Rootzone
depth 250m Empirically derived from 250x250(m) land-‐use map of the Netherlands.
5’ (geographical
minutes) Obtained from Global Crop Water Model (Siebert and Döll, 2010) Total
available water (TAW)
250m Derived from rootzone depth and 250x250(m) soil type map of the Netherlands.
5’ (geographical
minutes) Derived from rootzone depth and aggregated (5’) soil type map.
Readily available water (RAW)
250m Fraction of TAW, depends on
soil type 5’ (geographical minutes) Fraction of TAW, depends on soil type Total
evaporable water (TEW)
250m Derived from soil type and evaporation layer depth (100mm)
5’ (geographical
minutes) Derived from soil type and evaporation layer depth (100mm)
Readily evaporable water (REW)
250m Fraction of TEW, depends on soil type
5’ (geographical minutes)
Fraction of TEW, depends on soil type
Parameter (100 by 100m) Computation Description
Crop depletion factor (p) Assigned to CORINE land-‐use according to Allen et al. (1998)
Fraction of water a crop can extract from the rootzone without suffering from water-‐stress Maximum vegetation
coefficient (Kc max)
Derived from climate conditions and CORINE land-‐ use map
The coefficient for a vegetation type under optimal growing conditions
Current vegetation coefficient (Kcb)
Derived from leaf area index
(LAI) and Kc max The current vegetation coefficient
Soil moisture – NHI
The NHI contains five hydrological models that combine to solve the water balance in the Netherlands. Concerning soil moisture content, the Soil Vegetation-‐Atmosphere water Transfers (SVAT) model computes the vertical water transfer between the groundwater and the
atmosphere. On a daily basis and a 250m resolution, evapotranspiration and unsaturated transfer of water (infiltration, percolation and capillary rise) are computed to form the essential processes that determine water flow and content in the unsaturated soil layer (Van Walsum and Supit, 2012; Van der Bolt et al., 2015). Moisture content levels in the unsaturated zone are then combined with rootzone depth and TAW to compute soil moisture content and soil moisture deficit in the rootzone.
Soil moisture – PCR-‐GLOBWB
PCR-‐GLOBWB computes moisture content levels for two stacked soil layers and an underlying groundwater layer, where the upper soil layer is defined as the rootzone layer. Exchange between the layers consists unsaturated transfer. Additionally, water is extracted from the upper soil layer through evapotranspiration. In general, the mechanics of PCR-‐GLOBWB are similar to those of the NHI. However, the resolution is considerably lower and the parameterization is more simplified.
Sensitivity analysis
Due to the complex nature of the model and the large amount of parameters and variables involved, the model has many degrees of freedom. Consequently, one could argue that such a model can be made to produce virtually any desired behaviour (Hornberger and Spear, 1981). For this reason, a sensitivity analysis (SA) is conducted to determine the variation in model output that can be attributed to the different input variables (Saltelli, 2002; Saltelli et al., 2000). Through the SA, it is possible to determine which variables have the highest influence on the output and thus require the most attention concerning quality of input data and parameterization in the model. Due to the limited scope of this research, only the input variables were analysed.
The first step of a SA is to analyse the input variables on scale and shape. This includes the range of the input variation and the form of its probability density function (Saltelli et al., 2008). Secondly, the input variables are visualized by plotting them against the model output. Through inspection, these plots can yield general response trends and to some extent their strength (Saltelli et al., 2008). The third step is to statistically analyse and quantify the influence input variables exert on the output. For this, regression curves are fitted to the data. Resulting from the regression, the explained variance (R2) and residual mean square error (RMSE) are analysed. High values for R2 and low values for RMSE point to strong response of the output to the variable in question (Judd et al., 2011). Unlike done in the standard procedures of a SA, this research does not use variable values randomly drawn from the variables’ respective probability
density functions but conducts a similar analysis on actual model input data. For the analysis, soil moisture data from the NHI is used.
Analysis area and period
The research aims to provide insight in the question whether a high-‐resolution model can be extrapolated on European or Global scale. For this reason, any analysis of the model output in the Netherlands should cover a large variety of lower and upper boundary conditions. However, due to constraints in computing power and time, a selection is made concerning the analysed input and output data. Consequently, the analysis focuses on an area with relatively wet polders (Flevoland) and dry sandy soils (Veluwe) (Fig. 3A)). Furthermore, it spans the last week of June 2012, a period that contained warm sunlit days as well as cooler cloudy days (Fig. 3B)). Additionally, in this period the growing season has sufficiently started to show the effects of different vegetation types.
Fig 3. A) Soil moisture deficit (NHI), B) Analysis period hourly net radiation and temperature
4. Results
The model results are discussed in three parts. In the first part, the model output using soil moisture data from the NHI is discussed in terms of probable outcome, temporal and spatial gradients and systematic errors. In the second part, a similar procedure is followed with output resulting from PCR-‐GLOBWB soil moisture input. In the third part, the sensitivity analysis provides the influence of the different input variables on the model output.
Model output with data from NHI
To determine whether the model generates probable output in general, the energy balance was computed for the last week of June 2012 (Fig. 4).
Fig 4. Average hourly energy balance of the land surface of the Netherlands.
Figure 4 shows net radiation peaks between 150-‐450, evapotranspiration 75-‐200, transpiration 50-‐150 and evaporation 25-‐100 Wm-‐2. Additionally, subtracting evapotranspiration from net radiation produces the sensible heat flux peaking between 80-‐220 Wm-‐2.
However, sensible heat flux does not just vary through time but also through space. To gain insights in this, Figure 5 shows the average daytime sensible heat flux for June 21st. Furthermore, it contains the energy balance in a grid cell with high and low soil moisture content for June 21-‐23rd.
Fig 5. A) Average daytime sensible heat flux June 21st, high sensible heat flux (red) and low sensible heat
flux (blue). Energy balance June 21-‐23 in B) dry and C) wet cell. With soil moisture data from NHI.
There are clear differences between areas with high sensible heat flux like the Veluwe, Flevoland and urban areas and areas with low sensible heat flux. In terms of spatial resolution, differences in surface characteristics are easily recognizable in the sensible heat flux. Furthermore, Figure 3B-‐C illustrate that large differences exist in energy balance dynamics between relatively wet and dry grid cells. Wet grid cells show highly dynamic behaviour, with significant fluctuation in
the percentage of solar radiation that is expended as sensible heat flux. Dry grid cells on the other hand can be considered relatively inactive, only producing a small amount of transpiration.
Model output with data from PCR-‐GLOBWB
Figure 6 shows the energy balance in the Netherlands in the last week of June 2012. The slightly higher values of net solar radiation are caused by differences in clipping extent. Furthermore, evapotranspiration values seem to account for a larger portion of the energy balance resulting in lower values for SHF. Nonetheless, the magnitude of the energy fluxes is similar to those previously found in the computation with soil moisture data from the NHI.
Fig 6. Average hourly energy balance of the land surface of the Netherlands with soil moisture content
data from PCR-‐GLOBWB.
Fig 7. Average daytime sensible heat flux June 21st, computed using soil moisture data and
However, the spatial differences in computed SHF yielded less promising results. Figure 7 shows that computed SHF differs greatly from SHF computed using data from the NHI. First off, the low spatial resolution of PCR-‐GLOBWB is clearly visible in the results, with lower resolution spatial gradients in SHF as a result. For some areas the results can still be generally interpreted concerning lower and higher values for SHF, but for determining high-‐resolution border values these results are not suitable. Secondly, Figure 7 shows that there are large areas with negative average values for SHF occurring throughout The Netherlands. Analysing soil moisture deficit data as computed by PCR-‐GLOBWB (Figure 8) shows that those areas have a negative soil moisture deficit.
Fig 8. Soil moisture deficit (PCR-‐GLOBWB). Negative values
Variable sensitivity
Figure 9 shows input variables plotted against model output with linear and quadratic curves fitted through the data. Table 4 contains values for explained variance (R2) and residual mean square error (RMSE).
Fig 9. Scatterplots of input variables A) net solar
radiation B) temperature C) dewpoint temperature D) soil moisture deficit E) windspeed against model output. With linear (yellow) and quadratic (purple) regression lines.
Table 4.
Correlation between input variables and computed evapotranspiration over the Netherlands (N = 18000). R2 and RMSE denote the coefficient of determination and residual mean square error respectively.
Variable Polynomial R2 RMSE
Net radiation (Rn) Linear Quadratic 0.7349 0.7558 0.0426 0.0409 Temperature (T) Linear Quadratic 0.5814 0.5982 0.0534 0.0523 Dewpoint
temperature (D) Linear Quadratic 0.1033 0.1652 0.0961 0.0927
Soil moisture deficit
(Sd) Linear Quadratic 0.7274 0.7891 0.0031 0.0027
Windspeed (U2) Linear
Quadratic 0.2057 0.2057 0.0694 0.0694
From figure 7 it is obvious that radiation (A) and temperature (B) have a strong positive effect on the actual evapotranspiration. Furthermore, the soil moisture deficit (D) has a clear negative effect on the actual evapotranspiration. Windspeed (E) has a mild positive effect on evapotranspiration. The effect of dewpoint temperature (C) is the least pronounced, especially when computed through a quadratic regression.
The regression results of net radiation and soil moisture are similar in terms of explained variance, R2. However, the smaller variation in output response to the soil moisture deficit causes a significantly lower RMSE. Therefore, soil moisture deficit has the strongest relation to actual evapotranspiration. Temperature also has fairly high R2 and low RMSE values and observation of the scatterplots indicates a strong correlation with net radiation. Windspeed yields a fairly low R2 and average RMSE. R2 and RMSE are left unchanged when computing a quadratic relation, meaning the relation between windspeed and evapotranspiration is strongly linear. Dewpoint temperature has a significantly lower R2 and higher RMSE compared to the other variables. The negative relation is expected as higher dewpoint temperatures coincide with higher relative humidity, which negatively impacts evapotranspiration (Wallace & Hobbs, 2006).
5. Discussion
Model output with data from NHI
The energy-‐balance values seen in Figure 4 are similar to values produced by studies measuring and modelling the components of the surface energy balance (Meijninger et al., 2002; Wilson et al., 2002; Jia et al., 2003; Santos et al., 2009).
Looking at the spatial differences and wet versus dry grid cell in Figure 5, the results are largely realistic, as evaporation from the top layer quickly drops during prolonged dry periods.
However, the topsoil layer is also the first to receive water from the atmosphere, a characteristic that is not included in the hydrological models. Another systematic error of the model becomes evident when analysing the energy balance in Figure 5B. In the absence of a significant evapotranspiration flux, all energy is immediately expended as sensible heat flux. In reality, a portion of the radiation is used to heat the soil as soil heat exchange conduction (G). However, this variable was set to 0 in the model. The effect of this energy flux would be most pronounced during the start and end of the day, when it respectively extracts and adds energy to the system.
Model output with data from PCR-‐GLOBWB
Large areas showed negative average values for SHF. In terms of model mechanics, this is caused by the fact that according to moisture data from PCR-‐GLOBWB, these areas have a negative water deficit meaning they contain more water than the total amount present at field capacity. The consequence of this unlikely phenomenon is that evaporation rises dramatically, drowning out most spatial differences in SHF related to vegetation and soil type. In reality this may occur after periods of extreme rainfall in poorly drained areas. This issue can either be caused by faulty computation of soil moisture content or by wrong parameterization. Either way, it makes the soil moisture data from PCR-‐GLOBWB unreliable in the computation of SHF.
Variable sensitivity
The results found in the sensibility analysis were largely expected. However, the extent of the effect of soil moisture was less expected. In the Netherlands fairly wet conditions are generally dominant. Furthermore, the analysis period was not particularly dry or moist. Therefore, it was expected that the effect of soil moisture content on actual evapotranspiration would be less pronounced than observed. However, this effect may be caused by a variety of factors.
First off, there exists a strong correlation between soil type, vegetation and soil moisture deficit. For this reason, soil moisture deficit may be strongly correlated to other parameters related to evapotranspiration rates rather than being a cause on its own. To fully understand the differences in output and their cause, the parameters used in the model have to be analysed. However, the computation and therefore the analysis of these parameters are beyond the scope of this research.
Secondly, it is important to note that under wet conditions evaporation plays a large role in the total evapotranspiration (Fig. 5C). The current computation for the soil moisture content related to evaporation (top 100mm) is basic and lacks vital dynamics. It is computed on a daily basis from data obtained from the hydrological models and remains constant over a 24-‐hour computation period. However, in reality evaporation extracts and precipitation adds water to
the top layer throughout the day thereby influencing evaporation rates. Also, rootzone extraction, water percolation and dew formation further influence the top layer moisture content on a small timescale. To successfully capture the influence of moisture on evapotranspiration on a 1-‐hour resolution, a top layer water balance that is computed per hour should be added to the model.
6. Conclusion
With the available data and formulas it is possible to model sensible heat flux on a 100m spatial and 1 hour temporal resolution in the Netherlands. The results are probable concerning the order of magnitude and gradients over different surface characteristics. However, the output was heavily influenced by soil moisture, a characteristic of the model confirmed in the sensitivity analysis. Since soil moisture data is only computed on a daily basis, this characteristic flattens the energy balance dynamics throughout the day.
Furthermore, running the model with input data from PCR-‐GLOBWB proved to be problematic. The large influence of soil moisture caused the low resolution of PCR-‐GLOBWB to resonate strongly in the sensible heat flux output. Furthermore, issues related to the parameterization of PCR-‐GLOBWB affected the output. However, it is not entirely clear whether this results from errors in or insufficient study of the parameterization in question.
A possible solution to these problems consists of two efforts. First off, adding a 1-‐hour dynamic topmost soil moisture layer to the model would help to capture fluctuations in evaporation throughout the day. A similar layer can be constructed for the rootzone, with basic formulas for root uptake, percolation and capillary rise to catch the hourly fluctuations in moisture content. The second effort consists of creating high-‐resolution parameterization for areas outside of the Netherlands. With the availability of land-‐use and soil characteristics on a high resolution, parameterization can be computed in a similar fashion to the NHI to force the low-‐resolution soil moisture input.
With the input of realistic soil moisture content data, the evapotranspiration module produces an acceptable base for more complex computation of SHF. With the addition of soil heat exchange by conduction and interception latent heat of evaporation, the occurrence and size of the SHF can be modelled in more detail. However, more research is required to apply this method to larger spatial domains, as the current parameterization is incapable of coping with the soil moisture data provided by PCR-‐GLOBWB. Consequently, the base evapotranspiration most likely drowns out any further complexity in SHF dynamics produced by the soil heat exchange and interception evaporation module.
5. References
Allen, R.G., Pereira, L.S., Raes, D., Smith, M. (1998). Crop evapotranspiration – Guidelines for computing crop water requirements FAO Irrigation and drainage papers, 56. FAO – Food and Agriculture Organization of the United Nations, Rome.
Allen, R.G., Pereira, L.S., Smith, M., Raes, D., Wright, J.L. (2005). FAO-‐56 Dual crop coefficient method for estimating evaporation from soil and application extensions. Journal of Irrigation and Drainage Engineering, 1(2), pp. 2-‐13.
Van Beek, L.P.H. and M.F.P. Bierkens (2008), The Global Hydrological Model PCR-‐GLOBWB: Conceptualization, Parameterization and Verification, Report Department of Physical Geography, Utrecht University, Utrecht, The Netherlands, http://vanbeek.geo.uu.nl/suppinfo/vanbeekbierkens2009.pdf
Van der Bolt, F., Van Walsum, P., Veldhuizen, A. (2015). Toelichting op berekenen van verdamping in MetaSWAP in LHM.
Bohrer, G. et al. (2011). Estimating updraft velocity components over large spatial scales: contrasting migration strategies of golden eagles and turkey vultures. Ecology Letters, 15(2), pp. 96-‐103.
Bonan, G. (2002). Ecological Climatology. New York: Cambridge University Press
Bouten, W., Baaij, E.W., Shamoun-‐Baranes, J., Camphuysen, C.J. (2013). A flexible GPS tracking
system for studying bird behaviour at multiple scales. J. Ornithology, 154, pp. 571-‐580.
De Lange, W.J., Prinsen, G.F., Hoogewoud, J.C., Veldhuizen, A.A., Verkaik, J., Oude Essink, G.H.P., van Walsum, P.E.V., Delsman, J.R., Hunink, J.C., Massop, H.Th.L., Kroon, T. (2014). An operational multi-‐scale, multi-‐model system for consensus-‐based, integrated water management and policy analysis: The Netherlands Hydrological Instrument. Environmental Modelling & Software, 59, pp. 98-‐108.
Hartkamp, A.D., De Beurs, K., Stein, A., White, J.W. (1999). Interpolation techniques for climate variables. Mexico, DF (Mexico): CIMMYT. Series: CIMMYT NRG-‐GIS Series.
Hedenström, A. (1993). Migration by Soaring or Flapping Flight in Birds: The Relative importance of Energy Cost and Speed. Philosophical Transactions of the Royal Society B, 342(1302).
Hedenström, A. (2002). Aerodynamics, evolution and ecology of avian flight. Trends in Ecology & Evolution, 17(9), pp. 415-‐422.
Hornberger, G.M., Spear, R.C. (1981). An approach to preliminary analysis of environmental systems. Journal of Environmental Management, 12, pp. 7-‐18.
Jacovides, C.P., Kontoyiannis, H. (1995). Statistical procedures for the evaluation of evapotranspiration computing models. Agricultural Water Management, 27, pp. 365-‐371.
Jia, L., Su, Z., van den Hurk, B., Menenti, M., Moene, A., De Bruin, H.A.R., Yrisarry, J.J.B., Ibanez, M., Cuesta, A. (2003). Estimation of sensible heat flux using the Surface Energy Balance System
(SEBS) and ATSR measurements. Physics and Chemistry of the Earth, Parts A/B/C, 28(1-‐3), pp. 75-‐88.
Judd, C.M., McClelland, G.H., Ryan, G.S. (2011). Data analysis: A model comparison approach. New York: Routledge
Kunz, T. et al. (2008). Aeroecology: probing and modelling the aerosphere. Integrative & Comparative Biology, 48(1), pp. 1-‐11.
Meijninger, W.M.L., Hartogensis, O.K., Kohsiek, W., Hoedjes, J.C.B., Zuurbier, R.M., De Bruin, H.A.R. (2002). Determination of area-‐averaged sensible heat fluxes with large aperture scintillometer over heterogeneous surface – Flevoland field experiment. Boundary-‐Layer Meteorology, 105(1), pp. 37-‐62.
Newton, I. (2008). The migration ecology of birds. Oxford: Academic press.
Norberg, U.M. (1990). Vertebrate flight: mechanics, physiology, morphology, ecology and evolution. Berlin, Springer-‐Verlag.
Rutter, A.J., Kershaw, K.A., Robins, P.C., Morton, A.J. (1971). A predictive model of rainfall interception in forests, 1. Derivations of the model from observations in a plantation of Corsican pine. Agricultural Meteorology, 9, pp. 367-‐384.
Saltelli, A. (2002). Sensitivity Analysis for Importance Assessment. Risk Analysis, 22(3), pp. 579-‐ 590.
Saltelli, A., Chan, K., Scott, M. (2000). Sensitivty analysis, New York: John Wiley & Sons.
Saltelli, A., Ratto, M., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Taratola, S. (2008). Global Sensitivity Analysis: The Primer, Chichester, West Sussex: John Wiley & Sons.
Santos, C.A.C., Silva, B.B., Rao, T.V.R., Neale, C.M.U. (2009). Energy balance measurements over a banana orchard in the Semiarid region in the Northeast of Brazil. Pesquisa Agropecuaria Brasileira, 44(11).
Shamoun-‐Baranes, J., Bouten, W., van Loon, E., Meijer, C., Camphuysen, C. (2016). Flap or soar? How a flight generalist responds to its aerial environment. Philosophical Transactions of the Royal Society B.
Shamoun-‐Baranes, J., Bouten, W., van Loon, E. (2010). Integrating meteorology into research on migration. Integrative & Comparative Biology, pp. 1-‐14.
Siebert, S., Döll, P. (2012). Quantifying blue and green virtual water contents in global crop production as well as potential production losses without irrigation, Journal of Hydrology, 384, pp. 198-‐217.
Stull, R.B. (1988). An Introduction to Boundary Layer Meteorology. Dordrecht: Kluwer Academic Publishers.
Vansteelant, W., Bouten, W., Klaassen, R., Koks, B., Schlaich, A., van Diermen, J., van Loon, E., Shamoun-‐Baranes, J. (2015). Regional and seasonal flight speeds of soaring migrants and the role of weather conditions at hourly and daily scales. Journal of Avian Biology, 46(1), pp. 25-‐39.
