Thermal Lift Below Sea Level

(1)

Thermal lift below sea level

Walter van Dijk

December 2020

Abstract

Solar radiation heats the Earth’s surface in a diurnal cycle. This influences the temperature gradient in the air above it. If air at the surface is sufficiently warmer than the air at higher altitudes, vertical upwards movements of air can occur, creating thermal updrafts. Some bird species can use the thermal lift obtained from these updrafts to fly more efficiently. Thermal updrafts can occur on various scales and altitudes. However, with increasing proximity to the Earth’s surface, localised variation in land use becomes more important. High-resolution data is therefore essential to simulate thermal updrafts near the surface. This report describes a set of models focused on estimating detailed variation in thermal updrafts close to the surface, where the effect of surface heating is relatively large and the effect of turbulent mixing is relatively low. This is done to aid research on flight patterns and energy expenditure of birds that can make use of localised atmospheric variation. To estimate the occurrence of thermal updrafts, sensible heat flux can be used as a proxy, especially at lower altitudes. The sensible heat flux can be calculated by solving the energy balance at the surface of the earth, which consists of partitioning energy from radiation into warming the soil, warming the air and evaporating water. The amount of energy that goes into warming the air is effectively the sensible heat flux. To solve the energy balance, a soil heat flow submodel is used to calculate the energy flux to the soil, and a crop factor submodel is used to calculate the energy used for evaporation. Solar radiation is obtained from the ERA5 reanalysis data set. The soil heat flow submodel calculates the movement of energy through the soil, accounting for variation in conductivity and moisture content in both parts. The crop factor submodel separates soil evaporation and plant transpiration, also accounting for soil moisture content. The modelling effort described here has resulted in an hourly estimation of sensible heat throughout the growing season of 2018 at a 5 m spatial resolution in North Holland. This data shows a diurnal pattern in surface heating, as well as seasonal variation. It is clear that variation in land use has a significant effect on the sensible heat that is produced. In the research area, the urban landscape produces the most sensible heat, while wet grasslands produce the smallest sensible heat flux. On croplands the sensible heat flux depends significantly on the growth stage and crop type. Areas of open water show different dynamics due to mixing, which often results in a relatively high sensible heat flux at night. The results presented here can, in combination with data from GPS-tracked birds, be used to investigate various aspects of bird flight behaviour.

(2)

Acknowledgements

During this research I had the privilege of work-ing with state of the art data and technology. Additionally, the aid of various people and insti-tutions have been essential for the completion of this thesis.

Firstly, I would like to thank Joachim Hunink for providing the hydrological data. This data was essential for the dynamic representation of the landscape.

I would also like to thank Wageningen Univer-sity & Research for the LGN2018 land use map and the crop factors from SIMGRO. The land use map has been especially important, since it allowed for modelling on a very high spatial res-olution.

Additionally, I am grateful I could use the Lisa system to run my code. Because of the computa-tional requirements, it was not realistically pos-sible to produce output for the whole research area on a normal computer. Therefore, using the Lisa system was essential for producing the large amount of results described here.

Finally, I would like to thank my supervisor prof. Willem Bouten for the guidance and coun-selling throughout the project. Your personal ex-pertise on both hydrology and movement ecology has been invaluable during this project. Your ex-citement in the the combination of these fields of research has also been contagious and greatly in-creased my own motivation for the project.

1 Introduction

Bird flight is a complex type of locomotion. There are various intrinsic and extrinsic factors that influence this behaviour, including physiol-ogy, the environment, the weather and the need for nutrition (Alerstam, 1991). Consequently, there is much variation in the flight patterns of birds, between species and individuals. How-ever, the various types of flight birds perform can be generalised in two main types: soaring and flapping (Dhawan, 1991). Soaring greatly reduces the energetic cost of flight (Baudinette & Schmidt-Nielsen, 1974). Some bird species de-pend on soaring flight, while others use it op-portunistically or don’t use it at all (Pennycuick, 1972).

The birds that use soaring flight are not always able to do so. Some type of favourable atmo-spheric condition is required (Hedenstr¨om, 1993). For a bird to sustain an altitude or climb in al-titude while soaring, air needs to be rising. Ris-ing air in the atmosphere can be caused by vari-ous phenomena, including: orographic lift, ther-mal lift and convergence zones. Orographic lift is caused by obstructions forcing horizontal air movements upwards, thermal lift is caused by a difference between air temperatures at the sur-face and higher in the atmosphere, and in a con-vergence zone the collision of air masses forces air upwards (Stull, 1988). In this research the focus is on thermal lift, and specifically how it varies throughout the landscape.

The lesser black-backed gull (Larus fuscus) is known to use a wide variety of flight modes (Klaassen et al., 2012), and it can adapt its flight behaviour to meteorological conditions (Shamoun-Baranes & van Loon, 2006). Addi-tionally, previous research (Shamoun-Baranes et al., 2016) has indicated that gulls use thermal lift at low altitudes. This makes it an interesting species for studying how environmental variation in thermal lift affects flight patterns.

Because thermal updrafts are created by diur-nal surface heating, variation in the surface influ-ences the distribution of thermal updrafts in the atmosphere. The effect of variation in the surface decreases with altitude, because at higher

(4)

alti-tudes air becomes more mixed (Garratt, 1994). Consequently, the strength of thermal updrafts close to the surface cannot be derived from medium-range weather models because they lack the spatial resolution to include small scale land-scape variation (Shamoun-Baranes et al., 2016).

Thermal updrafts are formed in the atmo-spheric boundary layer, which is the part of the atmosphere directly influenced by surface pro-cesses (Garratt, 1994). There is a clear diurnal influence of surface heating on turbulence and updrafts in the boundary layer, which influence its thickness (Garratt, 1994; Yang et al., 2004; Jacobson, 2005). Thermal updraft are created when the temperature difference between air at the surface and air at a higher level is more than the environmental lapse rate, which is the tem-perature gradient in a stable atmosphere (Jacob-son, 2005). To estimate the velocity of these ther-mal updrafts, convective velocity scale w∗can be

used as a proxy (Stull, 1988). w∗ can be

cal-culated as (Deardorff & Willis, 1974; Holtslag & Boville, 1993)

w∗= (gHzi/T )1/3 (1)

where g is gravitational acceleration, H is the sensible heat flux at the surface, zi is the

bound-ary layer height and T is the mean potential tem-perature.

Boundary layer height and potential tempera-ture only have a limited spatial variation, and can therefore be derived from medium-ranged weather forecasts, such as the ERA5 data set (Dee et al., 2011) produced by ECMWF (Euro-pean Centre for Medium-Range Weather Fore-casts). This leaves sensible heat flux H as vari-able that needs to be estimated in order to cal-culate w∗. Although H is also available in the

ERA5 data set, the spatial resolution is not high enough to aid the investigation of flight behaviour of gulls.

H is one of the fundamental components of the energy balance at the earth’s surface. This balance partitions energy from solar and thermal radiation into warming of the soil (ground flux), warming of the air (sensible heat flux) and evap-oration of water (latent heat). H is the energy from radiation used for heating up air and can be

estimated as the residual of energy after evapora-tion and ground flux have been subtracted from radiation.

Various previous efforts have been made to ob-tain values for H from remote sensing (Friedl, 2002) and modelling (Choudhury & Monteith, 1988; Baldocchi et al., 2000). A large proportion of these studies were done within the context of modelling the water balance. As a result those studies often focused on obtaining estimates of the latent heat flux, rather than H. However, because this also requires solving the energy bal-ance at the earth’s surface, these studies can still be relevant here.

The goal of this research is to assess how en-ergy partitioning at the earth’s surface can be modelled on a high enough spatial and temporal resolution, to simulate the low-altitude sensible heat flux. All within the context of supporting the analysis of flight behaviour of lesser black-backed gulls. This includes the representation of a diurnal heating pattern in time and relevant landscape features in space. Model predictions can be compared to GPS tracks of the gulls to improve our understanding of how this species uses thermal lift.

In this research a set of models are used to sim-ulate the processes underlying thermal updrafts. This simulation aims to investigate the effects of local environmental variation on formation of thermal updrafts that facilitate low-altitude soar-ing flight.

2 Methods

Because gulls can make use of local environmen-tal variation in the landscape (Shamoun-Baranes & van Loon, 2006), estimating the spatial and temporal dynamics of H on a high resolution is the focus of this research. The modelling effort described here is fundamentally based on solving the energy balance at the surface of the earth. This balance describes the partitioning of radi-ation into warming of the soil, warming of the air and latent heat (Budyko, 1961). This can be formulated as

(5)

where H is sensible heat flux, Rnis net radiation,

G is soil heat flux and λE is latent heat flux. Rn

is the sum of solar (Rs) and thermal radiation

(Rt). Rn can be retrieved from meteorological

data, such as the hourly ERA5 data set (Dee et al., 2011), and the model calculates G and λE every hour, which allows for an hourly estima-tion of H. It is important to consider the effect of surface albedo here, since this determines the effective amount of available energy. In case of the ERA5 data set, radiation data corrected for albedo is available. This corrected data is used in this research.

Except for the soil temperature, all model pa-rameters are derived from input data every time step. Therefore, only the soil temperature needs to be initialised. This is done by using soil tem-perature measurements of KNMI. The deepest available measurements over a period of 39 years are averaged to obtain a temperature of 11.21◦C for the lowest soil layer in the model, at a depth of 3 m. The other layers are then assigned a temperature by linearly interpolating to the air temperature.

The temperature at 3 m depth below the sur-face is assumed to remain constant for the du-ration of the simulation. This is the boundary condition at the lower edge of the system. At the upper edge, radiation and air temperature are the boundary conditions. These variables are derived from the ERA5 data set. Between the boundary conditions, the dynamic parameters such as land use, crop factors and soil water content influence the energy balance. This results in a represen-tation of the spatial and temporal dynamics of H.

This research aims to provide a functional set of models to accurately predict the variation in atmospheric conditions gulls encounter during flight. A comparison to H derived with SEBAL (Bastiaanssen et al., 1998; Jacob et al., 2002) is made to test whether the estimations are suffi-ciently similar.

In this section, first the theory behind calcu-lating G and λE is explained, after which it is showed how H is derived from the results. Sec-ondly, the research area is described. Addition-ally, further information is given on the different

data sets that are used and how they are com-bined. And lastly, the comparison with SEBAL is described.

2.1 Model Description

2.1.1 Soil Heat Flux

Calculating the energy flux through the soil (G) consists of finding the amount of radiation ab-sorbed by the soil, and calculating the conduc-tive energy transfer between the vegetation, soil surface and deeper soil. The amount of energy absorbed by the soil can be derived from the amount of solar radiation reaching the soil sur-face, which is limited by vegetative cover. The model uses leaf area index LAI and solar zenith θ to calculate the land cover fraction τ by us-ing a version of Beer’s law, as described in Monsi & Saeki (1953). This is further explained in ap-pendix A.1.1. The absorbed energy A is then calculated as:

A = α ∗ τ ∗ Rs (3)

where α is the albedo of the soil surface, τ is the land cover fraction and Rs is the incoming solar

radiation.

The energy from solar radiation absorbed by the soil is summed with the net long-wave radia-tion and added to the heat content of the soil sur-face layer. From the heat content and the ther-mal capacity of the soil, the temperature of the soil can be derived. The difference in tempera-ture between the air and the soil surface and the difference in temperature between the deeper soil and the soil surface determines how much energy is moved by conductive transfer. Heat flow from the air to the soil surface is influenced by aero-dynamic resistance of the landscape, while heat flow through the soil is calculated with the energy continuity equation (eq. 4) (Campbell, 1985).

CδT δt = λ

δ2_T

δz2 (4)

The formulation of Johansen (1977) is used to ap-proximate thermal conductivity of the soil based on soil type and soil water content, further de-scribed in appendix A.1.1. Soil heat capacity is

(6)

calculated as the sum of the heat capacities of dry soil and the amount of water that is present in the soil.

The energy flux is calculated for 10 layers: a vegetation layer, soil surface layer and eight soil layers, down to a depth of 3 m. Flux between layers is calculated from the difference in heat content, conductivity and layer thickness as seen in equation 4. Conductivity for the soil layers is determined by soil type and soil water content, while the conductivity of the vegetation layer de-pends on the aerodynamic resistance of the land-scape. The energy flux between the soil surface and the vegetation layer is G in the energy bal-ance (eq. 2).

Water bodies are seen as fully saturated soils that consist only of ’pore space’. The surface layer is extended to 0.5 meter to account for the deeper penetration of solar radiation in water compared to soils. To account for the mixing of water, the conductivity variable calculated in equation 11 (appendix A.1.1) is substituted for an effective thermal diffusivity parameter that includes molecular heat transfer, turbulent mix-ing and any other process influencmix-ing heat trans-fer in the water column, as explained in Dut-ton & Bryson (1962). This value is set to 1000 W m−1K−1, which is in the order of values found by Lozzi-Koˇzar & Koˇzar (2017). This is a sig-nificant simplification of reality because the mix-ing is not homogeneous throughout the water col-umn. However, because water bodies are not a focus of this model, this simplification is deemed sufficiently realistic.

2.1.2 Latent Heat Flux

The energy used for evaporating water is calcu-lated with the crop factor approach (eq. 5) (Allen et al., 1998). In such methods a reference evap-oration is calculated for a standard type of land use, often well-watered grassland with a specific height, which is multiplied by a crop factor to account for the difference in evaporation com-pared to the reference. This allows calculation of the evaporation of land use types with a sin-gle parameter. A crop factor varies throughout the year to account for vegetative development

(Yoder et al., 2005).

Eact= Kc∗ Eref (5)

When using a single crop factor, the contri-butions of plant transpiration and soil evapora-tion are combined. Instead of a single crop fac-tor per land use type, the model uses the dual crop factor method as described by Allen et al. (1998) is used in the model (equation 6). Both crop factors, Kcb for vegetation and Ke for soil

evaporation, include a evaporation reducing func-tion, determined by water availability, soil type, wind speed, and relative humidity. This allows for the implementation of the effect of limited water availability on evapotranspiration originat-ing from plants and soil separately. The primary mechanisms of these evaporation reducing func-tions are explained in appendix A.1.2, the calcu-lation of Kcband Keadjusted for water stress is

further explained in Allen et al. (1998).

Eact = (Kcb+ Ke) ∗ Eref (6)

Various methods to calculate the reference evapotranspiration exist, with various levels of complexity (Yoder et al., 2005). It is shown by De Bruin & Lablans (1998) that in the Nether-lands the Makkink approach is an adequate substitute to the internationally used Penman-Monteith method. The main advantage of the Makkink method is that it only requires in-coming solar radiation Rs and air temperature

Ta. Because of this simplification, availability

of Makkink-based crop factors from SIMGRO (Van Walsum et al., 2004) and the similarity of the Makkink method to the Penman-Monteith method in the Netherlands (De Bruin & Lablans, 1998), the Makkink method is used in the model to obtain values for reference evaporation:

Eref =

C λ

s

s + yRs (7) where Eref is the reference evaporation, λ is the

latent heat of vaporization of water, s is the slope of the saturation water vapour-temperature curve at Ta and C is an empirical constant found

(7)

multiplied with a crop factor the respective ac-tual evaporation value is derived (equation 5 and 6). Through this method it is possible to ad-dress the difference in evaporation between vari-ous land use types in a relatively simplistic man-ner, without requiring many parameters.

2.1.3 Sensible Heat Flux

After the calculation of G and λE, H can be derived from equation 2 as the result of sub-tracting λE and G from Rn. The model is

one-dimensional, thus the resulting H is not influ-enced by the surrounding cells.

The H can subsequently be used in equation 1 to calculate the convective velocity scale w∗.

This comprehensible measure of the influence of energy partitioning on the atmospheric condi-tions can be used to analyse the relevance for bird flight.

2.2 Research Area

The research area concerns North Holland, a 4,092 km2 _{coastal province in the Netherlands.}

It is a relatively flat area, with coastal dunes on the western side. A significant proportion of the area has a surface below sea-level and is arti-ficially kept dry. The area has a wide variety of land uses including grassland, crop cultivation and towns.

North Holland is chosen because there are GPS-tracked gulls in this area and there is good data availability as input for the model. Addi-tionally, wet and dry areas are present, which makes choosing reference locations required for validation with SEBAL possible (Jacob et al., 2002). Figure 1 and 2 show the land use in the research area in 2018.

The area is modelled for the growing season of 2018, from March to October. These months contain the highest variation in crop stages, and is the period with the most lesser black backed gull activity in the area (Baker, 1980). In the search area, 2018 was a very dry year, which re-duces the influence interception evaporation may have on the model results. The progressive dry-ing of the soil throughout the year also increases

the visibility of the effect a low soil water content has on H.

2.3 Data Sets

Various data sets were used to run the model for North Holland during 2018 (table 1). The land use data from LGN2018 (Rip & Hazeu, 2019) is used to classify the area in land use types at 5m resolution (figure 1 and 2). This detailed land use map forms the foundation of the model, because it determines the spatial scale. All other data sets are spatially interpolated to this map.

The temporal resolution of the model is related to the temporal resolution of the meteorological variables, obtained from the ERA5 global reanal-ysis data set (Dee et al., 2011). This data set has an hourly temporal resolution, which translate to the model’s time steps of one hour. The vari-ables that are derived from this data set include short- and long wave radiation, temperature at 2 m, dew point at 2 m and wind velocity at 10 m. Using the meteorological data and crop fac-tors derived from SIMGRO (Van Walsum et al., 2004), reference evaporation was calculated with equation 5.

The soil type is derived from ESDAC soil tex-ture data (Panagos, 2006). This provides infor-mation on the field capacity and wilting point of the soil, used to calculate the evaporation re-duction. Soil types are also used for the conduc-tivity and heat capacity in the heat flux model (Johansen, 1977). Because this data has a spatial resolution of 1 km, some of the finer details in the LGN2018 land use map are not represented here, however the variation in soil type is represented in the data.

Soil water data is derived from Nederlands Hy-drologisch Instrumentarium (NHI) (De Lange et al., 2014). This data is used to calculate the soil conductivity and evaporation.

Because LAI values are different every year due to variation in growing conditions, it is inac-curate to assign standard values of LAI to land cover classes. Therefore, as input for the model, processed LAI data from the Copernicus data set (Baret et al., 2016) is used. This data has a spatial resolution of 300 m, which implies that

(8)

Figure 1: The land use map used in this research (Rip & Hazeu, 2019). It has a spatial resolution of 5 m, and is made to represent the situation of 2018.

(9)

Figure 2: Legend for the land use map (figure 1). These are the original Dutch names for the various land use classes (Rip & Hazeu, 2019).

(10)

Name Information Spatial grid Temporal scale Source

LGN2018 Land use 5m year Rip & Hazeu (2019) NHI Soil saturation 100m day De Lange et al. (2014) ERA5 Meteorology 0.10/0.25 deg hour Dee et al. (2011) European Soil Database Soil texture 1km - Panagos (2006) LAI300 Leaf area index 300m dekad Baret et al. (2016)

Table 1: Data sets used as input for the model.

simply using this data as it is would result in major edge-effects due to the 5 m spatial reso-lution of the model. To overcome this problem, the LAI for every land use class is statistically determined.

First a statistical summary is made of the LAI for every class of the land use map. The inter-polation method used to derive values is nearest neighbour, because it is expected that a large ho-mogeneous field will have the same LAI through-out the whole field, thus being accurately repre-sented in a 300 m grid cell. The mode of the LAI is subsequently assigned to the land use class to represent every cell of that type. This works well for land use classes that consist of large homoge-neous areas, such as croplands, because most of the LAI values are derived from 300m cells that are completely within this land use type. How-ever, classes that mainly consist of narrow spatial features have an unrealistic value. Consequently, the LAI of roads is inaccurate because it is al-ways the average of the roads and the 300m next to it. Therefore roads and water are assumed to have no vegetative cover and the LAI is set to 0 for these classes.

The data sets used in the model are varied in their spatial and temporal resolution. How-ever, they all capture their own relevant varia-tion. Having more detailed input data is there-fore not expected to increase the model accuracy.

2.4 Combining Data

To perform calculations at the LGN2018 resolu-tion, all data was converted to the same spatial scale. This conversion consists of interpolation to the spatial scale of LGN2018, or assigning parameter values according to land use classes.

For soil texture, interpolation was done with the nearest neighbour method because the soil tex-ture is a discrete parameter. Linear interpola-tion was used for all continuous input parame-ters. To make the data sets compatible, some pre-processing was necessary. This is explained in appendix A.3.1.

Complementary to the interpolation, some data is adjusted with the high-resolution land use map. From the land use map, impermeable surfaces and areas with open water can be accu-rately located. The specific physical attributes of these areas are incorporated in the soil and hydrological data. Water bodies are classified as having only pore space and a constant saturation of 100%. Impermeable surfaces are considered to have no saturation. Additionally, the field capac-ity and wilting point of impermeable surfaces are set to zero.

The temporal aspect of model inhibits a sim-ilar process, every time step information on all parameters is required. To reduce unnecessary duplication of data, some of the model variables are updated on a daily basis instead of hourly. As some of the variables only show variation on a day to day scale, updating daily captures the dynamics adequately. This is further explained in appendix A.1.3.

There are computational limits to the resolu-tion a model can have. Because the spatial reso-lution of LGN2018 is used, there are 176 million grid cells in North Holland. In combination with 5880 time steps of 1 hour, the amount of calcu-lations required to produce the results presented here are approaching the current limitations to what is possible within reasonable time. Calcu-lations were structured as efficiently as possible, to limit the computational resources required to

(11)

generate results. This included the use of parallel computing. The methods that were used to do this are further explained in appendix A.3.2.

2.5 Testing with SEBAL

A method to test the model is required to assess its accuracy compared to other models. There-fore, the variation and average values of H in the model are compared to sensible heat flux maps constructed from remote sensing. The remote sensing based method to derive H is Surface En-ergy Balance Algorithm for Land (SEBAL) (Bas-tiaanssen et al., 1998; Jacob et al., 2002). This provides momentary estimations of H, which can identify discrepancies with the model at that point in time.

SEBAL is a surface energy balance algorithm that estimates sensible heat with infrared remote sensing imagery (Bastiaanssen et al., 1998). It is meant to work on cloudless days and requires availability of wet and dry areas within the re-search area. All of these requirements can be fulfilled with available MODIS satellite images, though only a limited amount of cloudless mea-surements are available.

Inputs for the SEBAL algorithm include: wind speed at blending height, air temperature, leaf area index and surface temperature. Wind speed and air temperature are provided by the ERA5 data set, vegetation height is estimated from the land use type in the LGN2018 land-use map, leaf area index is derived from the LAI300 data and the LGN2018 land-use map, and the surface tem-perature is calculated from MODIS infrared im-agery.

In SEBAL, H is computed using the following equation

H = ρ ∗ cp∗ δT rah

(8) where ρ is air density, cp is air specific heat

con-tent, δT is the air temperature difference between two altitudes, and rah is the aerodynamic

resis-tance to heat transport. This equation inhibits an inverse problem because the values for rah

change according to atmospheric stability. Espe-cially in dry conditions, H can influence stability

of the atmosphere. SEBAL applies the Monin-Obukhov theory in an iterative process to solve this problem. This is explained in great detail in Bastiaanssen et al. (1998). In appendix A.2 additional information on the implementation of SEBAL in the research area is given.

3 Results

A simulation of sensible heat over North Holland during the growing season of 2018 has been per-formed. The results include a data set indicating the sensible heat flux (H) in the research area for every hour between 1 March and 30 October 2018 with a spatial resolution of 5m.

Figure 3 shows the average H of the simulation period. The actual obtained data set consists of 5880 maps, one for every hour between 1 March and 30 October 2018. Because of the high spa-tial and temporal resolution, this data set (146 GB) cannot be included in a digital appendix. Instead, a video of the simulation is included to give an accurate impression of the results, albeit in a lower resolution.

Although figure 3 gives a brief indication of the spatial variation in the model results, a lot of the dynamics of the behaviour of H are lost when only analysing the average values. Because the model has an hourly temporal resolution, there is also a temporal dimension that can be explored. Figure 4 and 5 show two weeks of data. The mean of 100 randomly chosen locations with the same type of land use is represented in every line. The locations used for figure 4 are the same as the locations used for figure 5.

It is important to note that figures 4 and 5 only show a limited subset of the available data. The figures give a fair impression of the distribution of H between common land use types, but there is a lot more variation than displayed. There are over 40 land use classes used in the model. Al-though some classes are more similar than oth-ers, all have a distinct behaviour regarding pro-duction of sensible heat. Additionally, there can be variation within a given land use class due to differences in soil type, water availability and weather.

(12)

(13)

Figure 4: A week of data showing the daily patterns in sensible heat over six land use classes from the 28 April to 5 May. The land use classes are each represented by 100 randomly chosen locations. Tick marks on the x-axis indicate the start of a day at midnight.

Figure 5: A week of data showing the daily patterns in sensible heat over six land use classes from 28 July to 4 August. The land use classes are each represented by 100 randomly chosen locations. Tick marks on the x-axis indicate the start of a day at midnight.

(14)

Figure 6: Sensible heat, cover fraction and soil saturation near Anna Paulowna. The upper images show the situation at midday on 1 May, while the lower images show the situation at midday on 1 August. At the bottom the land use of the area is shown, along with a legend.

(15)

Figure 7: Sensible heat, cover fraction and soil saturation near Grootschermer. The upper images show the situation at midday on 1 May, while the lower images show the situation at midday on 1 August. At the bottom the land use of the area is shown, along with a legend.

(16)

Figure 6 and 7 give an indication of the relative contribution of land coverage and soil saturation to sensible heat at two points in time. These im-ages are not meant to completely dissect the pat-terns in sensible heat to all underlying variables, however they can be useful to gain insight in the spatio-temporal dynamics of the model. Two similar days have been chosen, but no two days are exactly the same. Soil saturation and cover fraction, as shown in the figures, can explain a large proportion of the sensible heat. However other variables such as crop development can also influence production of sensible heat.

The location of figure 6 is Anna Paulowna, a village surrounded by agricultural fields that pro-duce various crops. This figure shows the effect of different agricultural practices. The western fields contain predominantly flower bulbs, while various food crops are cultivated in the south-western area.

Figure 7 shows a more natural grassland area within agricultural fields. The grassland is grazed extensively and has many ditches between the fields. Here the effect of drought can be seen more clearly. The sensible heat over the central grassland is higher on 1 August, while the soil water content and the cover fraction is lower.

It is important to acknowledge the role of crop factors in both figure 6 and 7. The crop factors effectively link vegetative development and wa-ter availability to evaporation. Since the dual crop factor method is used (Allen et al., 1998), it not only indicates how the vegetation influences transpiration, but also the effect of vegetation on evaporation. This makes it a complicated process not easily distilled to a few images.

4 Discussion

4.1 Average Sensible Heat Flux

In figure 3 the average H during the simulated period is depicted. When comparing this with figure 1, a clear relation between sensible heat and land use becomes apparent. Areas with the lowest values of H consist mainly of wet grass-lands. These are areas that almost always have water available for evaporation, and cool down

severely during the night. Water bodies show a higher H because they are constantly subject to mixing, which results in a higher surface temper-ature and consequently a higher H at night. An intermediate H can be found at areas where crops are grown and the dunes at the western coast. Both of these areas are dryer than the grasslands and soils are bare more often. The highest values of H are found at urban areas and roads. This logically follows from their properties: the ’soil’ consists mainly of concrete and asphalt and there is no evaporation1_.

Water bodies also have some interesting varia-tion in the average H. The values are moderately high, which might be counter-intuitive because evaporation is not limited. Most of the energy from radiation is still used for evaporation, but a substantial part of the energy is heating the wa-ter. Because the water is constantly mixing, the surface does not warm up or cool down quickly. This has as a consequence that during the day the surface of water bodies is generally colder than the surface of land, but at night the surface is generally warmer. Moreover, the surface of water is often warmer than the air at night, resulting in a positive H during the night. Therefore, when averaging day and night, water can have a high H.

4.2 Sensible Heat Through Time

The difference in H between land use classes becomes even clearer in figure 4. Firstly, it is shown that the values for roads are large com-pared to the other classes. Another interesting aspect is the behaviour of water, being fairly con-stant compared to the other classes, which have a clear daily pattern. Maize has a high H, com-pared to other vegetated areas. This is because in May maize crops have not undergone much vegetative development yet, while for the other classes the soil is already covered to a larger ex-tent. This is different later in the year, as shown in figure 5, where the flower bulbs show a larger heat flux because these fields are bare.

1_{Urban green areas have a separate land use class,}

(17)

Because the angle of the sun is similar in figure 4 and 5, the amount of solar radiation received at the surface is also approximately the same. Consequently the peaks in H are comparable be-tween these two figures. If a time period closer to the longest day is chosen, there is a larger in-flux of available energy to produce sensible heat, resulting in higher peaks.

The day to day variation shown in figure 4 and figure 5 is the result of changing conditions of the atmosphere, soil and vegetation. The weather of-ten has the largest short-term effect, because it can change rapidly. It however depends on the situation what the relative contributions of all variables are to the patterns in H. As a conse-quence, the H expected from one variable might not reflect the actual H.

Figure 6 illustrates an example of such a situa-tion. In the southwestern area soil water satura-tion has decreased between 1 May and 1 August, though there is also a lower sensible heat flux on 1 August. Normally, a decrease in water results in a higher sensible heat flux since less energy is used for evaporation. In this situation however, land cover plays an important role. The figure shows that the change in sensible heat is similar to the change in the cover fraction. The decrease in sensible heat is caused by the growth of crops, which increase the cover fraction and transpira-tion.

4.3 Comparison with SEBAL

A comprehensive set of models as used in this study is fundamentally difficult to validate as there are no measurements available at an ap-propriate extent and resolution. Satellite images are useful, but they can only provide snapshots of the H energy landscape.

Additionally, algorithms that calculate H from remote sensing data have their own inherent un-certainty. Nevertheless, one validation effort is made using SEBAL with MODIS data, which re-sulted in two snapshots of H to compare with the model. For this comparison the model out-put was aggregated from 5m to 1000m to match the MODIS resolution. The results of this com-parison are shown in figure 8 and 9.

The average H and the distribution through-out the area was relatively similar, as shown in figure 8 and 9. The main differences consist of a slightly higher H in the model and higher H the dune area. The higher H in the model can be caused by three inaccuracies: the model of this report, SEBAL itself (Bastiaanssen et al., 1998), or the estimation of H for the ’hot’ pixel in SEBAL (Voogt & Grimmond, 2000). Which of these contributes the most remains to be seen, though it is likely that the model overestimates H over urban areas, since it assumes that there is no cover and no moisture.

The discrepancy in values for H near the west-ern coast is also difficult to explain. First it is important to note that this is the only part of the research area that borders the sea. The east-ern water body is a lake, and the northeast-ern water is the Wadden Sea, which is an intertidal zone. Although the model uses more categories, two main types of land-use can be identified near the western coast: dunes and urban area. Mainly the dunes in the model have a lower H in figure 8 and a higher H in figure 9, compared to the results from SEBAL.

A major downside of the use of MODIS data to acquire surface temperature for validation is the relatively low spatial resolution. The model output consequently needs to be aggregated from a 5 m grid to a 1000 m grid. Landsat data can be used to derive surface temperature maps as well, and can potentially reduce this discrepancy since its spatial resolution is 100 m. This makes Landsat data an interesting candidate for further validation efforts. It is not used in this report because deriving surface temperature data from Landsat requires more intermediate steps than deriving this information from MODIS data (So-brino et al., 2004).

Remote sensing data, such as MODIS and Landsat, allows for comparison of the spatial context, but omits the temporal aspect of the model. To assess whether the model accurately explains the temporal dynamics of H, measuments with a similar temporal resolution are re-quired. As a topic of future research, the model can be used to simulate an area where H mea-surements are available. One potential location

(18)

Figure 8: Sensible heat flux in W m−2 _{at 06-Jun-2018 10:36 calculated with the model (A) and with with SEBAL}

using MODIS data (B). Results from the model are aggregated to unify the resolution. Below the maps, boxplots of all the pixel values (C) are displayed for both methods and a scatterplot (D) shows the similarity per pixel.

(19)

Figure 9: Sensible heat flux in W m−2 _{at 13-Oct-2018 11:18 calculated with the model (A) and with with SEBAL}

using MODIS data (B). Results from the model are aggregated to unify the resolution. Below the maps, boxplots of all the pixel values (C) are displayed for both methods and a scatterplot (D) shows the similarity per pixel.

(20)

is Cabauw, located approximately 40 km south of the research area. Here, a 200-m meteorologi-cal tower observes atmospheric and land surface processes (Beljaars & Bosveld, 1997; Bosveld et al., 2020).

4.4 Application of the Model

The energy balance at the surface of the Earth is often investigated to model evaporation or mete-orological processes. However, since sensible heat provides insight in the vertical movements of air, the energy balance can be interesting for ecolog-ical research as well.

In this report the data is summarised in figures and text, however exploring the data becomes most interesting if it is put in some context. Al-though the output data consists of only sensible heat flux through time and space, it can be in-vestigated in many ways. Comparing the effect of different types of land-use on the sensible heat production can be interesting, but since sensible heat is influenced by many variables, other rela-tions can be investigated as well.

Essentially, the model was created to investi-gate bird flight behaviour in the context of ther-mal soaring. This can be done in various ways using the results presented here. Some options include: analysing whether flight paths are

in-fluenced by the probability of encountering high sensible heat flux, measuring if a longer distance is flown when much sensible heat is encountered, or analysing the effect of sensible heat on the en-ergetic cost of flight.

Another research topic possible with the re-sults presented here is the use of bird flight data to analyse atmospheric processes and conditions (Treep et al., 2016). It is known that rising air in the atmosphere can exhibit turbulent move-ments. Near the surface diffusive processes are dominant, however with enough movement even-tually turbulent motion takes over. This con-cept is known as eddy diffusion. A key param-eter that governs this process is the turbulent Prandtl number (Li, 2019). This parameter is the ratio between momentum and heat transfer eddy diffusivity. Birds that can use the thermal lift from both thermal columns and diffusion pro-cesses, such as Larus fuscus, can provide a valu-able insight in the separation of these processes. If it is possible to determine until which altitude a bird uses thermal lift from diffusive processes, and from which altitude a bird starts to use ther-mal columns, it might be possible to relate that to the sensible heat footprint and improve our understanding of the turbulent Prandtl number throughout the atmospheric boundary layer.

(21)

References

Alerstam, T. (1991). Bird flight and optimal migration. Trends in Ecology & Evolution, 6 (7), 210–215.

Allen, R. G., Pereira, L. S., Raes, D., & Smith, M. (1998). Fao irrigation and drainage paper no. 56. Rome: Food and Agriculture Organization of the United Nations, 56 (97), e156.

Avdan, U., & Jovanovska, G. (2016). Algorithm for automated mapping of land surface temperature using landsat 8 satellite data. Journal of Sensors, 2016 .

Baker, R. R. (1980). The significance of the lesser black-backed gull to models of bird migration. Bird Study, 27 (1), 41–50.

Baldocchi, D. D., Law, B. E., & Anthoni, P. M. (2000). On measuring and modeling energy fluxes above the floor of a homogeneous and heterogeneous conifer forest. Agricultural and Forest Meteorology, 102 (2-3), 187–206.

Baret, F., Weiss, M., Verger, A., & Smets, B. (2016). Atbd for lai, fapar and fcover from proba-v products at 300m resolution (geov3). INRA: Paris, France.

Bastiaanssen, W. G., Menenti, M., Feddes, R., & Holtslag, A. (1998). A remote sensing surface energy balance algorithm for land (sebal). 1. formulation. Journal of hydrology, 212 , 198–212. Baudinette, R., & Schmidt-Nielsen, K. (1974). Energy cost of gliding flight in herring gulls. Nature,

248 (5443), 83–84.

Beljaars, A. C., & Bosveld, F. C. (1997). Cabauw data for the validation of land surface parame-terization schemes. Journal of climate, 10 (6), 1172–1193.

Blanken, P. (2015). Land-atmosphere interactions — canopy processes. In G. R. North, J. Pyle, & F. Zhang (Eds.), Encyclopedia of atmospheric sciences (second edition) (Second Edition ed., p. 244 - 255). Oxford: Academic Press. Retrieved from http://www.sciencedirect.com/science/article/pii/B9780123822253001985 doi: https://doi.org/10.1016/B978-0-12-382225-3.00198-5

Bosveld, F. C., Baas, P., Beljaars, A. C., Holtslag, A. A., de Arellano, J. V.-G., & Van De Wiel, B. J. (2020). Fifty years of atmospheric boundary-layer research at cabauw serving weather, air quality and climate. Boundary-Layer Meteorology, 177 (2), 583–612.

Budyko, M. I. (1961). The heat balance of the earth’s surface. Soviet Geography , 2 (4), 3–13. Campbell, G. (1985). Soil temperature and heat flow. Developments in soil science’.(Eds B Singh,

M Gr¨afe) pp, 26–39.

Choudhury, B. J., & Monteith, J. (1988). A four-layer model for the heat budget of homogeneous land surfaces. Quarterly Journal of the Royal Meteorological Society, 114 (480), 373–398.

Deardorff, J., & Willis, G. (1974). Computer and laboratory modeling of the vertical diffusion of nonbuoyant particles in the mixed layer. Advances in Geophysics, 18 , 187–200.

De Bruin, H. (1987). From penman to makkink. In Evaporation and weather: Technical meeting 44, ede, the netherlands 25 march 1987. the hague, netherlands. 1987. p 5-31. 1 fig, 4 tab, 34 ref.

(22)

De Bruin, H., & Lablans, W. (1998). Reference crop evapotranspiration determined with a modified makkink equation. Hydrological Processes, 12 (7), 1053–1062.

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., . . . Vitart, F. (2011). The era-interim reanalysis: configuration and performance of the data assimilation system. Quarterly Journal of the Royal Meteorological Society, 137 (656), 553-597. Retrieved from https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.828 doi: 10.1002/qj.828 De Lange, W. J., Prinsen, G. F., Hoogewoud, J. C., Veldhuizen, A. A., Verkaik, J., Essink, G. H. O.,

. . . others (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 , 98–108.

Dhawan, S. (1991). Bird flight. Sadhana, 16 (4), 275–352.

Dutton, J. A., & Bryson, R. A. (1962). Heat flux in lake mendota 1. Limnology and Oceanography, 7 (1), 80–97.

Friedl, M. A. (2002). Forward and inverse modeling of land surface energy balance using surface temperature measurements. Remote sensing of environment , 79 (2-3), 344–354.

Garratt, J. R. (1994). The atmospheric boundary layer. Earth-Science Reviews, 37 (1-2), 89–134. Hedenstr¨om, A. (1993). Migration by soaring or flapping flight in birds: the relative importance

of energy cost and speed. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, 342 (1302), 353–361.

Holtslag, A., & Boville, B. (1993). Local versus nonlocal boundary-layer diffusion in a global climate model. Journal of Climate, 6 (10), 1825–1842.

Idso, S., Jackson, R., Reginato, R., Kimball, B., & Nakayama, F. (1975). The dependence of bare soil albedo on soil water content. Journal of applied meteorology, 14 (1), 109–113.

Jacob, F., Olioso, A., Gu, X. F., Su, Z., & Seguin, B. (2002). Mapping surface fluxes using airborne visible, near infrared, thermal infrared remote sensing data and a spatialized surface energy balance model. Agronomie, 22 (6), 669–680.

Jacobson, M. Z. (2005). Fundamentals of atmospheric modeling. Cambridge university press. Jiménez-Muñoz, J. C., Sobrino, J. A., Skoković, D., Mattar, C., & Cristóbal, J. (2014). Land surface

temperature retrieval methods from landsat-8 thermal infrared sensor data. IEEE Geoscience and remote sensing letters, 11 (10), 1840–1843.

Johansen, O. (1977). Thermal conductivity of soils (Tech. Rep.). Cold Regions Research and Engineering Lab Hanover NH.

Klaassen, R. H., Ens, B. J., Shamoun-Baranes, J., Exo, K.-M., & Bairlein, F. (2012). Migration strategy of a flight generalist, the lesser black-backed gull larus fuscus. Behavioral Ecology, 23 (1), 58–68.

Li, D. (2019). Turbulent prandtl number in the atmospheric boundary layer-where are we now? Atmospheric Research, 216 , 86–105.

(23)

Lozzi-Koˇzar, D., & Koˇzar, I. (2017). Estimation of the eddy thermal conductivity for lake botonega. Engineering Review , 37 (3), 322–334.

Monsi, M., & Saeki, T. (1953). Uber den lichtfaktor in den pflanzen-gesellschaften und seine bedeutung fur die stoffproduktion. Jap. Journ. Bot., 14 , 22–52.

Panagos, P. (2006). The european soil database. GEO: connexion, 5 (7), 32–33.

Pennycuick, C. (1972). Soaring behaviour and performance of some east african birds, observed from a motor-glider. Ibis, 114 (2), 178–218.

Rip, F., & Hazeu, G. (2019). Lgn2018. Geo-Info, 16 (5), 38–39.

Shamoun-Baranes, J., Bouten, W., van Loon, E. 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: Biological Sciences, 371 (1704), 20150395.

Shamoun-Baranes, J., & van Loon, E. (2006). Energetic influence on gull flight strategy selection. Journal of experimental biology, 209 (18), 3489–3498.

Sobrino, J. A., Jim´enez-Mu˜noz, J. C., & Paolini, L. (2004). Land surface temperature retrieval from landsat tm 5. Remote Sensing of environment , 90 (4), 434–440.

Stull, R. B. (1988). An introduction to boundary layer meteorology (Vol. 13). Springer Science & Business Media.

Treep, J., Bohrer, G., Shamoun-Baranes, J., Duriez, O., Prata de Moraes Frasson, R., & Bouten, W. (2016). Using high-resolution gps tracking data of bird flight for meteorological observations. Bulletin of the American Meteorological Society , 97 (6), 951–961.

Van Walsum, P., Veldhuizen, A., Van Bakel, P., Van der Bolt, F., Dik, P. v., Groenendijk, P., . . . Smit, M. (2004). Simgro 5.0. 1; theory and model implementation (Tech. Rep.). Alterra.

Voogt, J. A., & Grimmond, C. (2000). Modeling surface sensible heat flux using surface radiative temperatures in a simple urban area. Journal of Applied Meteorology, 39 (10), 1679–1699. Yang, K., Koike, T., Fujii, H., Tamura, T., Xu, X., Bian, L., & Zhou, M. (2004). The daytime

evolution of the atmospheric boundary layer and convection over the tibetan plateau: observations and simulations. Journal of the Meteorological Society of Japan. Ser. II , 82 (6), 1777–1792. Yoder, R., Odhiambo, L. O., & Wright, W. C. (2005). Evaluation of methods for estimating

daily reference crop evapotranspiration at a site in the humid southeast united states. Applied engineering in agriculture, 21 (2), 197–202.

(24)

A

Appendices

A.1 Description of the Model

In this section more detailed information about the calculation of the different aspects to the en-ergy balance equation are given. First the calcu-lation of soil heat flux and the latent heat flux are discussed. Additionally, the implementation of time in the model is explained.

A.1.1 Soil Heat Flux

To calculate the amount of energy absorbed by the soil, first the amount of radiation reaching the soil surface needs to be discovered. To start, the fraction of soil covered is calculated using Beer’s equation of light extinction, similarly applied as Monsi & Saeki (1953)

τ = exp(−βL) (9) where τ is the fraction covered, L is LAI and β is an extinction coefficient. This coefficient is dependant on leaf architecture and solar zenith, and thus varies with day and season (Blanken, 2015). Choudhury & Monteith (1988) found β could be approximated to be a constant 0.5. In this research however, a constant value does not capture the variation enough because the diur-nal pattern has a large influence of the model re-sults. Therefore, the change in vegetation layer length with a change in angle is incorporated in the function to arrive at

τ = exp(−β L

cos(θ)) (10) where θ is the solar zenith angle. This equa-tion essentially introduces an adjustment of LAI based on the increase in distance the solar radi-ation travels through the vegetradi-ation layer when the sun is closer to the horizon.

Land use class Albedo Soil dry 0.30 Soil wet 0.15 Roads 0.10 Urban area 0.20 Water 0.06

Table 2: Albedo classes

Of the radiation reaching the soil surface, a proportion is reflected. This is determined by the albedo of the soil. Idso et al. (1975) showed a re-lation between soil water and albedo, the albedo values for soil are chosen to be comparable to their findings. A linear relation between soil wa-ter and albedo is assumed. Depending on the soil type, albedo can be somewhat higher or lower than for a loamy soil. However, because the in-fluence of vegetative cover is much larger in most cases, all soils are assumed to have the same re-lationship between soil water and albedo as loam has.

Multiplying albedo with the proportion of land covered and incoming radiation (equation 3) re-sults in the energy absorbed by the soil surface. To calculate the flow of this energy through the deeper layers of the soil, first the conductivity of the soil needs to be known. This conductivity is influenced by the amount of water that fills the pore space in the soil. Due to the attraction be-tween water and soil particles, this is however a more intricate relation than the sum of the con-ductive properties of the soil particles and the available water. Therefore, the formulation of Johansen (1977) is used to approximate thermal conductivity of the soil, as seen in equation 11. This calculation is repeated daily to include the effect of a changing soil water content. It is writ-ten as

k = ke(ksat− kdry) + kdry (11)

where k is conductivity, ke is the Kersten

num-ber, which is a logaritmic function of the degree of saturation, and ksatand kdryare saturated and

dry conductivity, respectively.

Heat capacity of the soil is also derived from soil type and water content, though this can sim-ply be done by adding the heat capacity of the soil particles with the heat capacity of the water that is present. Knowing the heat capacity al-lows for the conversion of energy to temperature, and vice versa. Consequently, the flow of energy through the soil can be calculated with equation 4. This calculation is made every hourly time step. The amount of energy stored in the soil

(25)

layers is the only variable that needs to be saved between time steps, because it is not derived from data, but computed mechanistically.

A.1.2 Latent Heat Flux

Beside the flux of energy in the soil, the energy that is used for evaporation needs to be calcu-lated as well, to leave sensible heat as a resid-ual. The model uses a crop factor method to estimate the energy used for evapotranspiration at the Earth’s surface. This includes calculating reference evapotranspiration, and multiplying it with a crop factor to obtain the actual evapo-transpiration. The reference evapotranspiration is calculated with the Makkink method (eq. 7).

The model uses the dual crop factor method, as described by Allen et al. (1998). Here calcu-lation of the crop factor is split in two parts: the soil evaporation Ke and the plant transpiration

Kcb. Equation 6 shows that these are summed

to obtain the total evapotranspiration. Both Ke

and Kcbhave their own evaporation reducing

pa-rameters: Krand Ksrespectively.

Plant transpiration is calculated by using the values of Kcbwithout water stress from SIMGRO

(Van Walsum et al., 2004) and multiplying them with the related evaporation reducing parameter Ks. The calculation of Ks consists of analysing

the water content of the soil. If the water content is at field capacity, Ks is 1. When the water

content of the soil is reduced, Ksremains 1 until

the readily available water is removed from the system. After that Ks decreases linearly to the

point where the water content is at wilting point and Ksis 0. Readily available water is calculated

with equation 12.

RAW = p ∗ T AW (12) RAW is the readily available water, T AW is the total available water and p is the fraction of water that can be depleted from the root zone without water stress. T AW is the difference between the water content at wilting point and the water con-tent at field capacity, which can be derived from the soil texture data. p normally varies from 0.3 for shallow rooted plants with a high reference evaporation to 0.7 for deeply rooted plants with a

low reference evaporation. A value of 0.5 is com-monly used for many crops (Allen et al., 1998). In this model a value of 0.5 is used everywhere.

Calculating soil evaporation starts with the de-termination of Ke. Since the sum of Keand Kcb

determines the actual evapotranspiration, there is a physical limit to this value. That limit is calculated, as described by Allen et al. (1998), as a function of wind speed, relative humidity and plant height. Hereafter, Ke is calculated as the

difference between this limit and Kcb.

Kr, the evaporation reducing parameter of Ke,

is also calculated according to different stages. These stages are related to total evaporable water (TEW) and the readily evaporable water (REW). TEW is derived from the field capacity and wilt-ing point of the soil with equation 13.

T EW = (θf c − 0.5 ∗ θwp) ∗ Ze (13)

Here Ze indicates the depth up to which the soil

is subject to drying by evaporation. A commonly found depth of 100 mm is used for this parameter in the whole research area.

REW is the depth up to which water can evap-orate without restriction, this is set to a typical value of 10 mm (Allen et al., 1998) throughout the research area. When the soil water content is above REW, Keis 1 since there is no restriction

to evaporation. However, Ke decreases linearly

when REW is depleted, until Kereaches 0 when

the soil water content is at half the wilting point. A.1.3 Temporal Dynamics

Considering the dimension of time, three distinct types of data are used in the model: hourly up-dated, daily updated and non-updated. This se-lection is made based on the amount of variation within a given variable. Soil type does not change much over one year and water content does not change much in one day, but temperature varies to much to assume it to remain constant during the day.

The non-updated variables include land use and soil type. These variables can change on a yearly basis, but since the model only simulated one year it was not necessary to make these vari-ables dynamic.

(26)

Soil water content is updated daily, because this is assumed to capture the relevant variation. As a result, conductivity is also calculated ev-ery day to include the effect of changing soil wa-ter on conductivity of the soil. Beside soil mois-ture, variables related to plant development are also updated daily. These include crop factors with their evaporation reducing factor and the leaf area index.

The remaining variables have a temporal reso-lution of one hour. Most of these are calculated from the ERA5 weather data, which has a reso-lution of one hour. However, there are also vari-ables that are mathematically acquired, such as the solar zenith. For these variables the resolu-tion can be as detailed as needed, but since the output will have an hourly resolution, updating these every hour is sufficient. The surface energy balance is solved for every hour of the simulation, using all variables. This combines the different temporal resolutions in a consistent output for every hour.

A.2 Comparison with SEBAL

To test the accuracy of the model described here, a comparison with SEBAL has been made. SE-BAL is an algorithm that solves the surface en-ergy balance (eq. 2), with as input weather data, land surface characteristics and surface temper-ature. In this report the surface temperature is derived from MODIS.

SEBAL requires the estimation of sensible heat at two reference pixels: one ’cold’ and one ’hot’. The cold pixel is chosen above a lake where it is assumed that the sensible heat equals 0, because all the energy would go into either warming or evaporating the water (Bastiaanssen et al., 1998). Normally the hot reference pixel for SEBAL is taken at a location where there is no water avail-able, but due to the resolution of MODIS (1000 m) this is not possible. Instead, the hot pixel was chosen at the city centre of Amsterdam. It is as-sumed that half the amount of radiation is con-verted to sensible heat at this urban area, which is a crude estimate based on literature (Voogt & Grimmond, 2000).

The variation of sensible heat in the model is in

reasonable agreement with the results from SE-BAL, as can be seen from figure 8 and 9. Al-though the sensible heat is generally higher in the model, areas with relative high sensible heat in the model also show high values in the SE-BAL analysis. The higher average sensible heat in the model can be the result of three inaccura-cies: the model of this report, the SEBAL model (Bastiaanssen et al., 1998), or the estimation of H for the hot and cold pixel in SEBAL (Voogt & Grimmond, 2000).

Two processes that result in additional discrep-ancies can be identified. Both are a result of the absence of horizontal influence in the model. Firstly there is the urban heat island effect. Since the model does not account for horizontal influ-ence, accumulation of heat by large cities is not included. In surface temperature measurements derived from remote sensing the urban heat is-land effect is included because the surface tem-perature is measured directly. Consequently it is likely that the model underestimates the amount of sensible heat produced over cities. However, since the model assumes no evaporation in urban areas, the H in urban areas can also be higher.

The second difference appears in the cells on the Western coast. This in the only place in the research area where there is a border between land and ocean. In figure 8 the sensible heat in the model appears to be higher than in SEBAL at this location, though in figure 9 it is lower. This difference might be due to variation in the difference between ocean temperature and land temperature, and how this influences neighbour-ing cells.

One major shortcoming of the validation effort presented here is the difference in spatial reso-lution. Where the model has a 5 m resolution, MODIS data has a 1000 m resolution. A sur-face temperature map derived from Landsat im-agery might provide a more accurate compari-son, because it has a 100 m spatial resolution. This also makes it possible to use for example a large parking lot for the hot reference pixel in SEBAL. Converting Landsat imagery to a sur-face temperature map is possible (Avdan & Jo-vanovska, 2016; Jim´enez-Mu˜noz et al., 2014), but there are no analysis-ready surface temperature

(27)

maps available, as is the case with MODIS data. Therefore Landsat data is not used in this report, but in further research a more accurate compari-son can be made by performing a SEBAL analysis with Landsat data.

Conversely, it can be argued that a higher res-olution will not necessarily improve the validity of the comparison. Firstly, the temporal aspect is still not tested. Since only a few relevant cloud-less Landsat images are available, it is not possi-ble to derive the diurnal pattern. Secondly, the inaccuracies of estimating sensible heat at the hot and cold pixel remain, and SEBAL has it’s own inherent uncertainty as well. Ultimately, running the model for a time and place where measure-ments are available is a more robust method.

A.3 Practicalities

A.3.1 Reprojection of Data

Because data from multiple sources are used, a way to homogenise the spatial references is re-quired. This model will ultimately be used in analysis of GPS data point in lat-long format, therefore the output data is in lat-long format as well. LGN2018 has a spatial reference that is standard for the Netherlands (RD new), but not all data has the same reference. The soil tex-ture map, for example, is stored in a projection valid for Europe. Because LGN2018 is the map with the highest detail and the correct spatial reference for usage in the Netherlands, all data will be adjusted to fit this map and correspond-ing spatial reference. This is done by calculat-ing the latitude and longitude of every cell in all data sets, which allows interpolating them to the LGN2018 spatial grid.

Calculation of the coordinates of cells is done in ArcGIS. First the map is imported with its projection as a raster. Then the raster is con-verted to a point data set with points at the cen-tre of each raster cell. Finally the point data is exported as a table containing the latitude and longitude for every point. Doing this for all data sets provides a solution for variation in projec-tions because all data will have explicit values for latitude and longitude.

A.3.2 Computational Resources

Applying the model to a limited amount of cells is not computationally intensive. Modern com-puters are more than able to handle the calcu-lations required. However, when the amount of cells are increased it will eventually become too much to handle. In this research an area consist-ing of 176 million cells is analysed for 5880 time steps. In this section the efforts made to reduce the amount of computational resources needed is explained.

Memory requirements become an important aspect when the input data sets with a high tem-poral resolution are interpolated to the LGN2018 spatial resolution of 5m. If all input data are interpolated to this spatial resolution, approxi-mately 100 Terabytes of space would be required to store the input data. Moreover, this data has to be loaded into the memory of the computer to run the model. This is very costly, highly un-practical and not necessary.

To reduce the requirements of the model, non-interpolated data are used as input data, polation will be done when required and inter-polated data will be deleted when it is no longer needed. This is possible because the model can produce hourly output of H using only input data for that hour and some information of the previ-ous hour. All input data before that is no longer required. This reduces the memory requirements by approximately a factor of 5000.

Because the non-interpolated data has a rela-tively small size, this method has the additional benefit that input data for a larger area than is essential can be used. In this research, input data for all of the Netherlands is used, there-fore the model can easily be adjusted to other areas within the Netherlands by using a land use map of the desired area. If a different time pe-riod or an area outside of the Netherlands is re-quired, the respective input data does need to be acquired again. This is often easily achievable, because all data except the land use map is open access.

Reducing the memory requirements with a fac-tor of 5000 still leaves the need for approximately 20 Gigabytes of memory, just to interpolate the

(28)

input data every time step. During the calcula-tions, this amount is increased significantly due to the creation of all variables derived from the input data. However, because the model is one dimensional, the calculations can be parallelised. Dividing the simulation area in subsections re-duces the memory requirements to an acceptable level. This also allows the calculation load to be divided amongst processing units to reduce the

time to run the model.

To acquire the results presented in this report, the computational load was shared over 32 nodes. Each node calculated the sensible heat flux for a part of the research area. After the simulation was completed for all sections the output data has been combined to produce the hourly sensible heat maps for the whole area.

(29)

A.4 Parameter List

Parameter Symbol Unit

Sensible Heat Flux H W m−2

Net Radiation Rn W m−2

Ground Flux G W m−2

Latent Heat λE W m−2

Solar Radiation Rs W m−2

Thermal Radiation Rt W m−2

Leaf Area Index LAI

-Solar Zenith θ degree

Proportion of land covered τ

-Albedo of the soil α

-Solar Radiation Rs W m−2

Soil Heat Capacity C J m−3K−1

Layer Thickness δz m

Actual Evapotranspiration ETact W m−2

Reference Evapotranspiration ETref W m−2

Crop Factor for Vegetation Kcb

-Crop Factor for Soil Ke

-Vegetation Evaporation Reduction Kr

-Soil Evaporation Reduction Ks

-Total Evaporable Water T EW mm Readily Evaporable Water REW mm Total Available Water T AW mm Readily Available Water RAW mm

Field Capacity θf c -Wilting Point θwp -Air Temperature Ta K Albedo α -Conductivity k W m−1K−1 Kersten Number ke

-Conductivity Saturated Soil ksat W m−1K−1

Conductivity Dry Soil kdry W m−1K−1

Wind Speed, Latitudinal v ms−1 Wind Speed, Longitudinal u ms−1

Dew Point Tdp K

Thermal Lift Below Sea Level