Heat Transfer Through Grass: A Diffusive Approach

Heat transport through short and closed vegetation such as grass is modelled by a simple diffusion process. The grass is treated as a homogeneous ‘sponge layer’ with uniform thermal diffusivity and conductivity, placed on top of the soil. The temperature and heat-flux dynamics in both vegetation and soil are described using harmonic analysis. All thermal properties have been determined by optimization against observations from the Haarweg climatological station in The Netherlands. Our results indicate that both phase and amplitude of soil temperatures can be accurately reproduced from the vegetation surface temperature. The diffusion approach requires no specific tuning to, for example, the daily cycle, but instead responds to all frequencies present in the input data, including quick changes in cloud cover and day–night transitions. The newly determined heat flux at the atmosphere–vegetation interface is compared with the other components of the surface energy balance at this interface. The budget is well-closed, particularly in the most challenging cases with varying cloud cover and during transition periods. We conclude that the diffusion approach (either implemented analytically or numerically) is a physically consistent alternative to more ad hoc methods, like ‘skin resistance’ approaches for vegetation and bulk correction methods for upper soil heat storage. However, more work is needed to evaluate parameter variability and robustness under different climatological conditions. From a numerical perspective, the present representation of vegetation allows for both slow and rapid feedbacks between the atmosphere and the surface. As such, it would be interesting to couple the present surface parametrization to turbulence-resolving models, such as large-eddy simulations.


Introduction
In this work, we present a novel view on heat transfer through short and closed vegetation, such as grass. In analogy to heat transfer through the soil, we model this transfer through grass as a diffusive process, with a harmonic analysis method. Thus, the vegetation is treated as a separate homogeneous layer with its own physical properties like the thermal diffusivity and thermal conductivity. After determination of these parameters, the model is used to calculate the temperature profile in the vegetation and heat flux at the atmosphere-vegetation interface.
Despite its simple, homogeneous appearance, grass is inherently a complex threedimensional composite consisting of biomass, air, and water. In fact, we only tend to realize this by close inspection. Given this complexity, it is not surprising that heat transfer through short vegetation layers is often described by simple empirical formulations, for example, by a skin layer in numerical weather and climate models (see, e.g., Viterbo and Beljaars 1995). Although such formulations (resistance law) may indeed account for the insulating effect of vegetation, they do not account for the frequency-dependent damping effect that the vegetation has on the temperature dynamics of the soil below it. For example, short-term variation in temperature at the atmosphere-vegetation interface due to cloudiness will only penetrate through the upper part of the vegetation. Long-term variation, such as the diurnal cycle, will penetrate through the entire vegetation and reach the upper soil. This is illustrated in Fig. 1, which shows the temperature at the top of the vegetation T veg,0 , the top of the soil T soil,δ , and at 0.05-and at 0.10-m depths (in the soil), T soil,5+δ , and T soil,10+δ , during two representative periods at the Haarweg station, the Netherlands (see Sect. 2). Whereas T veg,0 clearly exhibits high-frequency components with a time scale of 10-60 min, such components are already mainly damped out at the top of the soil. Additionally, the vegetation layer induces frequency-dependent phase shifts of the temperature: for example, the delay in the diurnal component is 2.5-3 h for this grass field between T veg,0 and T soil,δ .
Such frequency-dependent behaviour of heat transfer is known to occur in the soil medium, which is often studied apart from the vegetation above it. In that case, heat transfer is simply described by solution of the heat equation using Fourier series, if one assumes homogeneous soil properties with depth (see, e.g., Carslaw and Jaeger 1959;van Wijk and de Vries 1963). This approach (also called harmonic analysis) has indeed been proven successful in previous studies, for example, using observations from the Cooperative Atmosphere-Surface Exchange Study in 1999  in Kansas (van de Wiel et al. 2003), the Negev desert, Israel (Heusinkveld et al. 2004), and the Haarweg site (Jacobs et al. 2008). Analytical expressions for all frequency components to model the soil temperature dynamics were given in these studies after determination of the soil diffusivity and conductivity requiring two observed temperature time series and one soil-heat-flux time series.
In the current study, we investigate if the heat transfer through the vegetation can be modelled similar to that in the soil. We assume that the vegetation layer can be regarded as a quasi-homogeneous medium (sponge layer) with constant diffusivity and conductivity over depth, albeit different than those in the soil. An explicit analytic expression of the heat equation is derived for this two-layered system: a finite vegetation layer on top of the semi-infinite soil layer. The diffusivity and conductivity of both layers are determined by optimization, viz., they follow from fitting observations of temperature and heat flux onto each other using the analytic expression (nonlinear least-squares approach). The resulting expression then describes the frequency-dependent response of both temperature and heat flux in both the vegetation and soil layers as a result of the temperature forcing at the top of the vegetation. Although a similar mathematical expression has been obtained for a two-layered The temperature at the top of the vegetation is estimated from the radiative components and the temperatures in the soil are observed using shielded Pt-100 elements (see Sect. 2). The temperature at the top of soil T soil,δ is modelled using the harmonic model (see Sect. 3.2). All temperatures are given at 10-min resolution, corresponding to the sampling time of the observations soil system (Lachenbruch 1959), we are the first, to our knowledge, to apply such a model to a short vegetation layer.
We will show the potential of the model in two ways. First, the soil temperature at 0.05-m depth (in the soil) as predicted from the temperature at the top of the vegetation is compared to the observed soil temperature. Second, the newly derived surface heat flux at the atmospherevegetation interface is compared to the other main components of the surface energy balance, viz., the total radiation, and the turbulent sensible and latent heat fluxes. It is shown that the closure of the surface energy budget (SEB) is relatively good when using the new model. From a physical point-of-view, the heat storage within the vegetation is now automatically accounted for and needs no separate parameterization. This is because the surface heat flux is directly evaluated at the top of the vegetation. Heat storage in the underlying vegetation layer then directly follows from the vertical flux divergence over that layer. The method captures rapid transitions between day and night, and rapid variation due to changes in cloudiness. This is due to the fact that the novel harmonic method includes phase shifts between the flux of energy and the response of the vegetation, which also occurs in nature: to produce warming, a flux is required. Hence, rapid changes in incoming solar radiation directly alter the surface heat flux, whereas the temperature at the top of the vegetation responds with a short time lag. This is generally not the case with skin-layer formulations or bulk (caloric) methods, which either have instantaneous energy transfer or are 'tuned' to only a few frequencies (Viterbo and Beljaars 1995;Oncley et al. 2007;EBEX). As such, the benefit of the new method will increase with higher temporal-resolution observations. Finally, it is important to note that, although in the present study an analytical approach is followed, modelling of grass as a diffusive medium can also be done with a fine-mesh numerical scheme. In such a case, our method could be used for numerical inverse modelling to empirically obtain effective thermal properties of the vegetation and soil. The potential of this alternative route is briefly discussed at the end of this paper.
A short overview of the Haarweg measurement site and the used observations are presented in Sect. 2. These are followed by a description of the harmonic analysis method for the grasssoil system and the analysis procedure in Sect. 3. Section 4 gives the results and Sect. 5 gives the discussion of the current work. Finally, the summary and conclusions are presented in Sect. 6.

Haarweg Site and Observations
The observations used were obtained at the Haarweg measurement site located in the centre of the Netherlands (51 • 58 N, 5 • 38 W, altitude 7 m a.s.l.). This was the meteorological station operated by Wageningen University and Research from 1974-2012. The near surroundings are predominantly covered by perennial grassland (mostly Lolium perenne and Poa trivialis), which is mowed on a weekly basis during May-November and has an average height of 0.1 m. The soil predominantly consists of heavy river clay (Jacobs et al. 2008(Jacobs et al. , 2011, and is kept wet through seepage and experiences regular rainfall amounting to approximately 765 mm yr −1 (Jacobs et al. 2010). To keep the field from being waterlogged, it is superficially drained towards canals surrounding the field.
The incoming and outgoing shortwave (S in and S out ) and longwave (L in and L out ) radiative flux components are measured at 1.5 m above the surface with an aspirated pyranometer (model CM14, Kipp & Zonen B.V., Delft, the Netherlands) and pyrgeometer (model CG2, Kipp & Zonen B.V., Delft, the Netherlands), respectively, and averaged over 10 min. The vegetation temperature at the top of the grass T veg,0 is estimated from the incoming and outgoing longwave radiation via where = 0.99 is the emissivity of the grass and σ is the Stefan-Boltzmann constant. Due to a small sensitivity of the longwave radiometers to the shortwave radiation, a correction is applied to the radiative measurements (see Jacobs et al. 2008). The turbulent fluxes of heat and water vapour are measured using the eddy-covariance technique and averaged over 30 min. The observational heights are 3.44 m for the sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, Utah, USA) and 3.15 m for the gas analyser (LiCor7500, LI-COR, Lincoln, Nebraska, USA). Both are installed roughly in the centre of the meteorological site. The eddy-covariance data were processed through the fluxsoftware package EddyPro v5.1.1 (Fratini and Mauder 2014) from LI-COR Biosciences Inc. (Lincoln, Nebraska, USA). All necessary standard data-treatment and flux-correction procedures were included, such as axis rotation with the double-rotation method (see, e.g., Wilczak et al. 2001), raw data screening including spike removal (Vickers and Mahrt 1997), interval linear detrending, and low-pass filtering correction (Moncrieff et al. 1997). The fetch of undisturbed grass field extends to 300-350 m in all directions except for the north-east to east sector, which contains agricultural fields. As a result, wind directions between 020-090 • were excluded.
Soil temperatures are obtained (and averaged over 10 min) at depths of 0.05, 0.10, 0.20, 0.50, and 1.00 m using shielded Pt-100 elements. Small variations in these depths may occur throughout time due to soil compaction, frost events, and bioactivity. The soil heat flux is measured at a depth of approximately 0.05 m by a heat-flux plate (WS31, TNO, Delft, the Netherlands). The measured values are corrected for the shape of the instrument and for the discrepancy between the instrument conductivity and soil conductivity following Overgaard Mogensen (1970). As the correction factor also depends on the conductivity of the soil, the correction is included in our (fitting) procedure, where the soil parameters are determined (see Sect. 3.3). The patch of grass used for the soil temperatures and soil-heat-flux measurements is approximately 2 m × 2 m in size, and positioned within the larger meteorological station measuring approximately 120 m × 120 m.
For the year 2006, soil volumetric water content measurements are available at multiple depths in the soil (0, 0.025, 0.05, 0.10, 0.30, and 0.60 m). They are measured by a time-domain reflectometry system (TDR100 Reflectometer with S610 probes and CR23X datalogger, Campbell Scientific Inc., Logan, Utah, USA) at an interval time of 60 min. Sensors are placed in a triangle-formation spaced 8-8.5 m around the main logger to provide an areaaveraged value. Sensor measurements are converted in volumetric water content θ according to Topp et al. (1980). An overview of the Haarweg site, used observations, and naming scheme is shown in Fig. 2. A detailed description of the Haarweg site, instrumentation, and measurements, including the corrections, can be found in Jacobs et al. (2008). The patches for the soil temperature and heat-flux measurements are approximately 2 m × 2 m in size. Note that, the photograph is not fully representative for the investigated periods in this study, as the grass was still developing in the above photograph, but matured quickly after. The bare soil tile was not used in this study. During the summer, the grass height is mowed regularly to 0.10 m

Method: Harmonic Analysis
The standard soil harmonic method uses Fourier series to mathematically describe the soil temperatures and soil heat fluxes (see, e.g., Carslaw and Jaeger 1959;van Wijk and de Vries 1963). For a homogeneous, semi-infinite soil, an analytical solution over the entire depth of the soil is possible after determination of the soil diffusivity κ soil [m 2 s −1 ] and conductivity λ soil [W m −1 K −1 ]. In its simplest form, the analytical solution describes the evolution of the temperature T soil (z, t) in time t, the damping of its amplitude (the exponential part) and its phase shift (z D −1 k ) with depth z, when cosine-shaped temperature perturbations are imposed at the surface (z = 0) (Heusinkveld et al. 2004;Moene and van Dam 2014) whereT is the temperature at infinite depth, A k is the amplitude of the cosine wave at the surface, D k = √ 2κ soil /ω k is the damping depth, and ω k is the angular frequency. Extension to more realistic temperature signals containing multiple frequencies is then given by superposition of the individual frequency contributions, i.e., a Fourier series.
Here, we present an extension of the harmonic method for a two-layered system in which both layers are homogeneous and have their own thermal properties. The vegetation layer has small thickness and is placed on top of the semi-infinite soil layer. First, the mathematical solution describing the temperature and heat flux as a function of depth and time for an imposed periodic temperature forcing T veg,0 (t) of single frequency at the top of the vegetation is presented in Sect. 3.1. As above, this approach can then be extended to a temperature signal of any shape by applying the superposition principle (see Sect. 3.2). Finally, the procedure to determine the thermal properties of both layers is presented (see Sect. 3.3). From a historical perspective, similar two-layered solutions were presented earlier by Lachenbruch (1959), who studied inhomogeneous soils. However, the key difference is our application to the vegetation layer. We therefore derive and apply this solution to this new context.

Two-Layer Diffusion Model
A key assumption is that both the vegetation layer (indicated by j = veg) and the underlying soil layer ( j = soil) are homogeneous in composition, such that their thermal properties are constant with depth. Additionally, it is assumed that no horizontal gradients of temperature are present. Then, their respective temperatures T j (z, t) and the corresponding heat fluxes G j (z, t) are described by the one-dimensional diffusion equation of heat, in which λ j [W m −1 K −1 ] is the thermal conductivity in layer j, ρ j [kg m −3 ] the density, c j [J kg −1 K −1 ] the specific heat capacity, and κ j [m 2 s −1 ] the thermal diffusivity. The vertical coordinate z is defined positive downward into the ground with z = 0 placed at the atmosphere-vegetation interface. The interface between both layers, viz., the thickness of the vegetation, is situated at z = δ. Here, we opt for this coordinate system to simplify the mathematical expressions. The following boundary conditions are specified: where A and ω are the amplitude and the angular frequency of the top vegetation temperature T veg,0 . Note that, the long-term mean temperatures T j are assumed to be equal and are therefore subtracted. The conditions in Eqs. 4b and 4c specify the continuity of temperature and flow of heat across the interface in the absence of source terms. The last condition Eq. 4d specifies that all temperature perturbations vanish at infinite depth.
As the top vegetation temperature is given by a cosine, the following Ansatz is taken for the temperatures (single-frequency contribution) in both layers, in which A j (z) and φ j (z) are the depth-dependent amplitude and phase shift. For both arithmetic and computational convenience (see next Section), this Ansatz is rewritten in terms of complex wave components, where U j (z) is a complex amplitude and the asterisk denotes the complex conjugate. This complex amplitude is related to the original amplitude and phase shift according to A j = Re(U j ) 2 + Im(U j ) 2 and φ j = arctan(Im(U j )/Re(U j )). Substitution of this Ansatz into Eq. 3 results in the following ordinary differential equation for U j (and similarly for its complex conjugate) of which the general solution is where a j and b j are the weights, and β j = ω/(2κ j )(1 + i) is introduced as shorthand notation. Application of the boundary conditions results in the following relations The condition that T soil vanishes for infinite depth immediately results in setting b soil = 0, viz., the growing exponential (positive power) cannot be present. However, due to finite size of the vegetation layer, this is not the case for T veg and the growing exponential has to be retained. After rearrangement, the following solution is obtained The expression for T veg is only valid on the domain z ∈ [0; δ], whereas that of T soil is only valid for z ≥ δ. Using Eq. 10a, the heat flux at the atmosphere-vegetation interface (top of vegetation) corresponding to this single-frequency temperature signal becomes

Implementation of the Two-Layer Model
In this section, the implementation of the mathematical model is given, such that the model can be directly used for observational analysis. First, we note that our observed temperatures T j [n] or heat fluxes G j [n] consist of multiple frequencies, and are both discrete (averaged on interval t) and finite in length with N the number of samples. Therefore, the time series are decomposed using the fast Fourier transform (FFT) algorithm. For example, the vegetation temperature at the top of the vegetation T veg,0 is decomposed according to The frequency components for temperature at any desired position in the vegetation or soil are derived from those at the surface by use of the mathematical solution (Eq. 10). Here, we give the frequency coefficients for the vegetation temperature T veg [k](z) (valid for 0 ≤ z ≤ δ), and the soil temperatures The frequencydependent terms in the square brackets can be regarded as so-called transfer functions that transform the coefficients to other depths. The coefficients for the heat fluxes are formed by taking the derivates with respect to z and multiplying by −λ j . The modelled time series of temperature or heat flux follows from taking the inverse FFT. Furthermore, it follows from these transfer functions that, for example, the frequency coefficients at depths z 2 and z 1 in the soil are related via In other words, we regain the standard solution for the relation between two soil temperatures depending only on the difference in depth (see Eq. 2).

Optimization Procedure
We describe the procedure to apply the model introduced in Sect. 3.1 to obtain an estimated surface heat flux at the atmosphere-vegetation interface. As the thermal parameters (i.e., conductivity and diffusivity of both layers) are not known a priori, they need to be determined from the observed time series of temperature and soil heat flux. We assume the grass height to be fixed at 0.10 m. A relatively small change in grass height will result in a different set of vegetation parameters, with the overall qualitative behaviour remaining the same. The procedure consists of four general steps and requires two observed temperature series in the soil (T soil,5+δ and T soil,10+δ ), the observed heat flux in the soil (G soil,5+δ ) and the (radiative) top vegetation temperature (T veg,0 ). These steps are outlined below and are also schematically indicated in Fig. 2: 1. The thermal diffusivity of the soil κ soil is determined by optimization. Observations of T soil,5+δ are decomposed into its Fourier components according to Eq. 12. Using Eq. 10b, the Fourier components from T soil,10+δ are calculated from those of T soil,5+δ assuming a κ soil , and the inverse Fourier transform is applied to give a modelled T soil,10+δ . The value of κ soil is determined by minimizing the sum of squared differences between modelled and observed T soil,10+δ , i.e., nonlinear least squares. 2. The thermal conductance of the soil λ soil is determined. Now, the heat flux at z = 0.05 + δ m is calculated from the observations of T soil,5+δ and the optimized κ soil , corresponding to taking the derivative of Eq. 10b with respect to z. This calculation is performed in Fourier space, viz., the Fourier components G soil [k] are calculated from those of T soil [k] and the inverse FFT is used to obtain G soil,5+δ . The optimized λ soil results from minimizing the sum of squared differences between modelled and observed heat flux G soil,5+δ . 3. The thermal diffusivity and conductance of the vegetation layer κ veg and λ veg can subsequently be determined using the observed T veg,0 as input and T soil,5+δ as target series with the use of Eq. 10. These vegetation thermal parameters are optimized simultaneously as they are simultaneously present in the model (the factor m in, e.g., Eq. 10a). 4. After optimization of the thermal parameters, the heat flux at the atmosphere-vegetation interface can be calculated with Eq. 11 by superposition of the individual frequency components.
All Fourier transforms are performed on time series covering a full year. This means that all frequency components from the yearly scale to the 10-min scale are included. It is assumed that the temperature is quasi-periodic over this time period, such that a yearly averaged value exists. Any interannual trend is removed from the observations. Additionally, using a time series containing a full year increases the frequency resolution, i.e., the detectable step size in discrete frequency, and maintains the lowest observable frequency. The optimization procedure, however, is performed over relatively small time intervals of approximately 10 days. This is to minimize the effects of (unknown) changes in soil moisture, which strongly influences conductivity, and/or grass thickness. As a result, the optimized thermal properties are effective parameters for this time interval of ≈ 10 days.
Finally, the aforementioned yearly-average temperaturesT are subtracted from all observations before the optimization procedure due to depth-dependent biasses in these values. Such biases, although relatively small (≈ 1 K), may dominate the sum of squared differences. For interpretation of the results, the averageT (0.1 m) is added to all temperature time series after the optimization procedure.

Results
The procedure described in Sect. 3.3 is first applied to determine the thermal properties of the vegetation and soil layers. Subsequently, these thermal properties are used to calculate the conductive heat flux at the atmosphere-vegetation interface enabling the evaluation of the surface energy budget. Specifically, this is done for two 10-day periods in August-September 2005 and July 2006, corresponding to (late) summer. These two separate periods are chosen because they both contain relatively clear-sky and cloudy days, and they highlight the strengths and weaknesses of the current model. The averaged volumetric water content at 0.025-m depth in the soil was 0.11 m 3 m −3 during period 2. Unfortunately, no measurements of θ were available during period 1. An additional moist case during the spring of 2006 is briefly shown in Appendix 1. No suitable moist case in the summer period, viz., with relatively constant moisture and availability of turbulent flux measurements, was found. Figure 3 shows the time series of the temperature at 0.1-m depth as observed (in red) and modelled (in black) from the measurement at 0.05 m using Eqs. 10b and 12. For both periods, it is observed that the model accurately represents the observed time series of temperature assuming a constant diffusivity in the soil for each respective period. Apart from an accurate Where perhaps this could be expected by optimizing the diffusivity via a one-dimensional heat diffusion model, it is still not trivial that a model with a single fitted value suffices for the full 10-day time series. In other words: the close correspondence supports a simple diffusion concept to describe the soil, which is in reality a complex and composite medium. Similar arguments also apply to the vegetation layer, as will be discussed later on.

Soil Thermal Parameters
The diffusivities of the soil as found by the optimization procedure are approximately equal for both periods with the difference being < 5% (see Table 1). Our values for κ soil are close to that reported by Jacobs et al. (2008), who found κ soil = 3.2 × 10 −7 m 2 s −1 for a multi-year composite of the Haarweg site. It is well known that both the conductivity and heat capacity tend to increase with soil moisture content (Cosenza et al. 2003;Jacobs et al. 2008). As the diffusivity is the ratio of these, the influence of moisture on the diffusivity largely cancels out. Our results confirm that κ soil tends to be mostly constant in time.
At 0.10-m depth, the daily cycle and lower frequencies are clearly visible, whereas the higher frequencies are almost completely damped out (compare with Fig. 1). This is consistent with a damping (e-folding) depth of approximately 0.10 m for the 24-h cycle. This means that the amplitude of this frequency component has decreased by a factor e over this distance. Higher frequency components will have decreased by a larger fraction as they have smaller (a) (b) Fig. 4 As in Fig. 3, but for the soil heat flux at a depth of 0.05 m fitting one constant conductivity per period. The average standard deviations of the optimized soil heat flux, as given by the Monte Carlo method, are 1.6 W m −2 for period 1 and 1.2 W m −2 for period 2, respectively damping depths, D ∼ √ 1/ω. For example, variations at the hourly scale will have decreased by a factor of ≈ e 5 over the same distance. Note that the highest frequencies on time scales < 1 h have already been damped out by the vegetation layer (see Fig. 1).
With the optimized soil diffusivities, an analytical expression for the soil temperature T mod soil is available and the soil conductivities can be found. In the optimization procedure, this corresponds to minimizing the following cost function, in which (λ soil ) represents the instrument correction factor for the soil heat flux sensor following Overgaard Mogensen (1970). This correction factor is a function of the soil conductivity itself and thus has to be explicitly included in the procedure. The shape parameters of the heat flux plate can be found in Jacobs et al. (2008). The result of the optimization procedure is shown in Fig. 4 and the values of the conductivity are given in Table 1. The results are remarkably close, considering the fact that the relevant frequency-dependent information to fully model the heat flux at 0.05 m, i.e., information about phase and amplitude, is contained in just one temperature observation. However, one additional temperature observation and heat-flux observation are required to determine the soil thermal parameters. This indicates that the soil heat transport can be quite accurately modelled using Eq. 14 (or Eq. 2) assuming one constant κ soil and λ soil over the 10- (3.2 ± 0.8) · 10 −7 (0.90 ± 0.16) ( 8.0 ± 1.6) · 10 −7 (0.59 ± 0.08) 2 7-17 July 2006 (3.0 ± 0.8) · 10 −7 (0.52 ± 0.09) ( 12 ± 3) · 10 −7 (0.44 ± 0.06) The error bounds are one standard deviation around the optimal value and are determined using a Monte Carlo method (see Appendix 2) day intervals. Generally, the absolute value of the residuals is small, with the 90% percentiles being 3.3 W m −2 for period 1 and 2.0 W m −2 for period 2, respectively. A further improvement may be obtained by optimizing the depth of the sensors (which is seldom known exactly), or a better determination of these. For example, a small error in the phase is present for period 1. The diurnal component of the modelled signal consistently lags approximately 20 min with respect to the observations. This phase difference depends on the difference in depth of the sensors, when it is assumed that the optimized κ soil is correct. A tentative correction based on the diurnal frequency component suggests that the heat-flux sensor is positioned at a depth of z ≈ 0.042 m instead of 0.05 m during this period. Such a change will however also impact the optimized value of λ soil and further research is therefore needed.
Comparing both periods, it is found that the conductivity is significantly lower during period 2 than during period 1 (see Table 1). This likely indicates that the soil is wetter during period 1 than during period 2, which cannot be confirmed as no measurements of volumetric water content were available during period 1. It is however known that the conductivity increases with soil moisture, although the specific functional relation can be complicated (see, e.g., Cosenza et al. 2003). As we diagnose the conductivity for two relatively short periods, a full quantification of the effect of moisture is beyond the scope of the current study. Our diagnosed values for soil diffusivity and conductivity for this river clay area are in the range of values found from laboratory experiments on clay soils. A summary of laboratory values in Moene and van Dam (2014) gives diffusivity values ranging between 1.0 × 10 −7 m 2 s −1 for very dry clay to 4.5 × 10 −7 m 2 s −1 for wet clay (40% volumetric water content; compare also Abu-Hamdeh 2003). Likewise, the conductivity values range between 0.15 and 1.4 W m −1 K −1 for those conditions.

Vegetation Thermal Parameters
The optimized values of κ soil and λ soil enable the further determination of the diffusivity and conductivity of the vegetation layer corresponding to the third step in Sect. 3.3. As mentioned previously, κ veg and λ veg have to be determined simultaneously as they both appear in Eq. 10b. Large differences between the amplitude and phase of the input data T veg,0 and the modelled outcome temperature at 5-cm depth exist: T veg,0 may be as much as 10 K higher or lower than T soil,5+δ (see the black and blue lines in Fig. 1). A good representation therefore requires both sufficient damping of the (dominant) amplitudes and accurate phase shifts after the optimization. However, a direct optimization for both periods with all frequency components in the signal resulted in inaccurate representations of T soil,5+δ , which is likely the result of reaching a local optimum instead of a global one (not shown). To circumvent this issue, the model (step 3 in Sect. 3.3) is first only applied to the Fourier components corresponding to the diurnal variation (24-h cycle), such that an improved initial guess of the optimized κ veg and λ veg is obtained. Subsequently, the final values of κ veg and λ veg for the complete signal (again) containing all frequency contributions are optimized using the improved guess as initial value in the least-squares procedure.
The time series of T soil,5+δ as observed and modelled using T veg,0 as input are shown in Fig. 5. The corresponding values of κ veg and λ veg are given in Table 1. Reasonable agreement is obtained between the observed and modelled signals, which is quite remarkable if one realizes that a prediction of heat transfer over a discontinuous medium is made. The daily frequency component, which is most prominent, and the phase of the signal are accurately reproduced. Additionally, the model accurately reproduces the higher frequency components (on time scales < 24 h), which are necessary for intradaily 'wiggles' (see, e.g., 9 July 2006 in Fig. 5b). Therefore, it appears that the complexity of the vegetation-soil system can be modelled as a discontinuous medium of two homogeneous layers with their own thermal properties.
In spite of this, deviations are present between the observed and modelled temperature signal for both periods. This is most pronounced for period 1 (see Fig. 5a), in which the observed temperature is underestimated by ≈ 1-2 K during the beginning of the period, while this difference decreases towards the end of this period. The differences during period 2 are smaller in the beginning of the period (< 0.9 K) than near the end (< 1.6 K). These temporally The modelled temperature has been determined using T veg,0 as input (see Fig. 1) The average standard deviations of the modelled temperature, as given by the Monte Carlo method, are 0.023 K for period 1 and 0.026 K for period 2, respectively varying differences are caused by the lower frequency components. This is confirmed by repeating the procedure including only the higher frequency components (time scales < 48 h), which reduces the differences (not shown).
The values κ veg and λ veg are in the same order of magnitude as the corresponding parameters for the soil. Generally, the diffusivity of the vegetation is larger than that of the soil mainly due to its smaller density caused by the large air fraction, whereas the conductivity is smaller than that of the soil. This is to be expected due to the large fraction of air in the vegetation compared to the biomass. At the same time, moisture variations due to, e.g., rain and dewfall may have have a non-negligible effect on the thermal properties of the vegetation. This is not investigated here.

Estimated Surface Heat Flux and Energy Balance
After optimization of the thermal parameters of both vegetation and soil, the heat flux G veg,0 at the atmosphere-vegetation interface is calculated using Eq. 11. This is one of the main components of the surface energy budget alongside the total radiation Q n , turbulent sensible heat flux H and the turbulent latent heat flux L v E. These components of the SEB are shown in Fig. 6. Both periods show an alternation of relatively cloudy and clear-sky days and nights. During clear days, the net radiation typically exhibits a smooth behaviour in time and peaks around 500 W m −2 . On the other hand, large variations in time are present during the cloudy days. During clear nights, the radiation levels off to a near constant value at the end of the night. The latent heat flux is large during the day, but often negligibly small during the night. The turbulent sensible heat flux has smaller values compared to the latent heat flux for most observed days, indicating the ample presence of moisture to sustain evapotranspiration.
The modelled G veg,0 is of comparable magnitude to the sensible heat flux during the day. During the night, G veg,0 varies between zero and the net radiation in absolute value. The latter value being reached, when there is little turbulent activity. In such case, the radiative energy loss at the surface has to be mainly compensated by a supply from the vegetation-soil system. Although G veg,0 generally precedes T veg,0 and thus the outgoing longwave by ≈ 1.5 h, as determined by the maximum of the cross-correlation between these two signals (not shown), this is not necessarily the case for the other radiative components. Indeed, G veg,0 follows the variation in incoming radiative components at the 30-min to 1-h scale as a result of variation in cloudiness. An example is present during 27 August 2005, when two spikes occur in Q n (see Fig. 6a). As a result of the first sharp decrease in the incoming shortwave radiation, the heat flux at the atmosphere-vegetation interface momentarily reverses sign between these spikes and thus energy is released upward from the vegetation. The surface temperature T veg,0 then responds to this and temporarily lowers with respect to the average during that day. This illustrates the potential of a diffusive approach to heat transfer in contrast to, e.g., a bulk method between T veg,0 and T soil,5+δ , which will be unable to detect such spikes in G veg,0 as the associated frequency component is already damped out at a few centimetres in the soil (cf. Fig. 1). A skin-layer formulation will pick up such spikes, but only with a delay in time, as the skin-layer formulation changes in time with T veg,0 (see Sect. 4.4). Note that, the variations in incoming radiation at the 30-min to 1-h scale are also picked up by the turbulent fluxes. Figure 7 compares the available energy Q n − G veg,0 with the sum of the turbulent fluxes H + L v E in time. Here, the first quantity is averaged on 30-min intervals for the comparison. This figure shows a good correspondence between these both quantities. This correspondence is particularly strong during the increase (decrease) phase in the morning (evening) transition. Those are notoriously difficult to model with traditional approaches, such as calorimetric methods, which cause phase lagging (Oncley et al. 2007;EBEX). Discrepancies are mainly present either in the peak values at midday or in the values during night. For example, the turbulent fluxes tend to zero early at night on 29 August 2005 indicating negligible turbulent activity, while Q n − G veg,0 > 0 (compare with |G veg,0 | > |Q n |, Fig. 6a). This indicates that the 10-day optimized λ veg may be too large for this particular night. This situation is reversed during 3 September 2005 in which the relative difference during midday is larger than during night. Similar contrasts can also be found during period 2, for example, on 14 July 2006, where the correspondence at night is very strong, but worse the following day. Such discrepancies hint at a variation of κ veg and/or λ veg within the 10-day optimization period, which may be caused by variations in moisture.
Apart from such variation in the thermal vegetation properties, the turbulent fluxes as measured by the eddy-covariance technique may also be erroneous. Although a (small) minimum in H can occur at the beginning of the night (see, e.g., van der Linden et al. 2017), this minimum appears to be too strong during a number of nights in relation to the other SEBcomponents. No explanation is found for this particular behaviour. Additionally, erroneous behaviour of the sensible and latent heat fluxes can be present in the observed time series, for example, the sharp decrease in L v E during the night of 27-28 August 2005. Notwithstanding these differences, the comparison of the principal SEB-components confirms the potential of the modelled vegetation heat flux to provide an adequate closure of the SEB. Such a closure can generally not be expected when using, for example, the observed heat flux G soil,5+δ at the 0.05-m depth in the soil or the bulk method (see Sect. 4.4). Note that the optimization procedure for all thermal parameters does not depend on the turbulent fluxes H and L v E, and the SEB closure is therefore independently evaluated. A second view on the SEB closure is presented in Fig. 8, which shows a scatterplot of the turbulent fluxes against the available energy term. Generally, most observation points scatter around the 1 : 1 line, with some exceptions. The clusters with negative sum of turbulent fluxes and available energy close to zero, correspond with the aforementioned strong minimum in H .
Finally, we note that a rigorous study on closure of the SEB would require accounting for advection and additional storage terms, such as the photosynthesis energy flux, the air enthalpy change, and the moisture changes in both the air and vegetation (see, e.g., Jacobs et al. 2008;Meyers and Hollinger 2004;Mauder et al. 2020). Here, as explained, only the so-called crop enthalpy change, i.e., the change of energy in the vegetation due to heating, is automatically accounted for in our model by evaluation of vegetation heat flux at the top of the vegetation. A quantification of the unaccounted additional storage terms is beyond the scope of the current study. However, Jacobs et al. (2008) found that the inclusion of these additional storage terms may improve the SEB components at the Haarweg site by approximately 5% of which the energy change due to photosynthesis is the most dominant contribution (3%).

Comparison with Other Methods
The current results show the impact of accounting for the vegetation layer as a separate homogeneous medium ('sponge layer') with its own effective diffusivity, conductivity and heat capacity. However, it remains unclear how our heat flux performs in relation to other heat-flux approaches for use in the surface energy budget. Here, we briefly present such a comparison. Figure 9 shows the heat flux at the atmosphere-vegetation interface G veg,0 , the heat flux calculated with a skin layer G skin (see Viterbo and Beljaars 1995), the heat flux at the vegetation-soil interface G soil,δ and the observed soil heat flux at 0.05 m in the soil G obs soil,5+δ for a part of period 1. Here, G skin is evaluated using the observed temperature at the top of the vegetation and the (modelled) temperature at the top of the soil, according to where skin is an empirical transfer coefficient. Instead of an empirical determination, this value is calculated as skin = √ 2λ veg /δ ≈ 8.3 W m −2 K −1 (to be explained below). The current value is in range of values used for grass, i.e., 3-10 W m −2 K −1 (see, e.g., Duynkerke 1991;Steeneveld et al. 2006;ECMWF 2021).
Comparing these heat fluxes, large differences in amplitude are immediately observed between the heat fluxes evaluated at the top of the vegetation (G veg,0 and G skin ) and the heat fluxes evaluated either at the top of the soil (G soil,δ ) or measured in the soil (G obs soil,5+δ ). These differences amount up to ≈ 40 W m −2 between top of the vegetation and top of the soil, and ≈ 60 W m −2 between top of the vegetation and 0.05 m in the soil, respectively. These Fig. 9 Comparison of different heat fluxes during 28 August to 1 September 2005, which constitutes a 4-day part of period 1. This subperiod is taken to highlight the differences between the heat fluxes. G veg,0 and G soil,δ are modelled using the current method, G skin follows from Viterbo and Beljaars (1995) and G obs soil,5+δ is observed with the heat-flux plate differences in magnitude are caused by the heat storage in the respective layers themselves (compare with Fig. 1). Additionally, the damping of high-frequency components and the phase shift over depth (of the remaining frequency components) can also be clearly observed.
The comparison of G veg,0 to G skin is a bit more complicated (see black and dark blue lines). Globally, these two heat fluxes correspond well to each other both in terms of magnitude and temporal alignment. A closer inspection reveals, however, that G skin lags G veg,0 by approximately 10 to 50 min throughout the time series shown. As a result of the superposition of multiple frequencies, one precise value of the lag cannot be given. The occurrence of a lag follows directly from the mathematical structure of these heat fluxes. Whereas G skin is linearly related with T veg,0 (see Eq. 16), G veg,0 is related to the z-derivative of this temperature (see Eq. 11). Expressed in terms of the Fourier coefficients, this becomes in which B[k] contains the correction for the two-layer system (compare with the large brackets in Eq. 11). The term 1 + i corresponds to multiplying the amplitudes with √ 2 and shifting the signals ahead by one-eighth of a full period, i.e., 45 • . The effect of B[k] on the amplitude and phase shift is assumed minor here. The fact that G veg,0 thus precedes T veg,0 and G skin is a direct consequence of the use of Fourier's heat equation. Physically, this means that the temperature will only change some time after the flux of energy has occurred. Apart from the difference in phase shift between G veg,0 and G skin , a more subtle effect is present in the treatment of the amplitude. Assuming that T veg,0 and T soil,δ are accurately modelled by our model or can be accurately measured, a relation between frequency, κ veg , λ veg , the distance between both temperature evaluation points, and veg can be found. In such a case, veg may be chosen such that the amplitude response, not the phase response, between both formulations is equal for one single chosen frequency. In other words, this effectively makes veg a kind of tuning parameter by which the user can choose one frequency (e.g., the diurnal frequency), for which the amplitude response will be correct. However, in general, veg will be determined empirically or an educated guess will be made. In the above example (see Fig. 9), we made such an educated guess and took veg equal to √ 2λ veg δ −1 . Here, the factor δ −1 follows from a naive first-order approximation of the temperature derivative over the grass layer. The factor √ 2 follows from the amplitude change caused by 1 + i as described above. This value results in a high degree of correspondence between G skin and G veg,0 , especially in the diurnal amplitude. Implicitly, this particularly high degree of correspondence for the diurnal amplitude is caused by the fact that ω day /(2κ veg )δ −1 ≈ 3/2 for our grass layer during period 1. This results in a near optimal tuning of veg for the amplitude of the diurnal component, as described above. Note that, in general, large errors can be introduced, when using a naive first-order approximation for the derivative of an exponential function.
Notwithstanding all these considerations, we recognize that by our choice of veg the current G skin has very similar overall magnitude and temporal behaviour compared to G veg,0 . This implies that the average closure of the SEB on 24-h time scales or even shorter time scales will more-or-less perform equally to the SEB closure based on our G veg,0 . At the same time, when looking at semi-instantaneous closure or closure on short time scales of a few hours, the two methods will in principle differ, due to the inherent lack of phase shift in G skin . An example during which the closure on a very short time scales may be interesting is the rapid passage of clouds (see, e.g., the aforementioned spikes in Fig. 6a). However, it should also be noted that a true comparison during such very rapid events requires observations at even higher temporal resolution in all main SEB components to resolve the expected phase shift.
A final remark is made with respect to the use of the skin-layer formulation in numerical models. The skin-layer formulation does not allow for storage of energy in the vegetation layer, apart from the lack of phase shift itself. This also implies that during numerical integration too large fluxes are instantaneously inputted into (extracted from) the top soil layers, thereby inducing a too strong warming (cooling) over time.

Further Challenges
Despite the potential benefit of our model with frequency-dependent responses, we realize that our two-layer diffusion model is still a crude simplification of the realistic outdoor environment in which a number of assumptions have been made and physical processes are unaccounted for.
First, moisture is not explicitly taken into account in either the vegetation or the soil, while it is known to (strongly) impact conductivity and heat capacity (Jacobs et al. 2008;Moene and van Dam 2014). Whereas the relatively short optimization period of 10 days limits the impact of (long-term) variation of moisture in the deeper soil, short-term variation in the vegetation layer and top of the soil cannot be completely excluded. Apart from the resulting variation in the thermal properties itself, another effect from changes in moisture is potentially present. Moisture may be evaporated from within the vegetation, thereby extracting energy locally, which corresponds to sources (sinks) of energy inside the vegetation medium. The effect of this evaporation on the vegetation heat flux will be observed with some time delay, as the signal also has to diffuse through the vegetation. At the same time, it will only be observed as latent heat flux, when it leaves the vegetation layer and passes the sensor after another unknown time delay. Additionally, moisture being evaporated from the top soil may also be intercepted and deposited within the vegetation as "dew rise" thereby leading to an internal redistribution of energy.
An extension of the analytic model to account for moisture currently seems only partly possible. Temporal variation in the thermal parameters, e.g., λ veg , may possibly be related to temporal variation in moisture content, provided the change is relatively slow, such that moisture is a rather passive actor. In such case, the optimization procedure may be performed with time-dependent thermal parameters. This approach would require that detailed observations of moisture content are available throughout the vegetation and soil (not available in this study). Another possible extension to investigate is the inclusion of more layers in which every layer is assigned its separate thermal properties. Although this may make the simple transfer model impracticable, a numerical implementation is possible.
Second, a number of assumptions have been made to derive the model. One of these assumptions is the constant (and rigid) height δ of the vegetation layer. Also, the vegetation layer may respond to the ambient wind and the turbulent eddies may penetrate some distance into the vegetation. A similar infiltration may occur for the radiative components in which different fractions are absorbed at different heights in the vegetation (Duynkerke 1992;Ronda and Bosveld 2009). As a result, the exchanges of energy will not be solely restricted to the top of the vegetation. This becomes particularly important if the vegetation considered is higher than grassland and when it is more open. It remains to be assessed whether our current approach with modelling the vegetation as sponge layer is still a useful analogy for those cases, in comparison to traditional parameterizations.

Numerical Implication
Notwithstanding the challenges identified in the previous section, the current findings seem to have potential to improve the next generation of land-surface models (LSMs). The results show that a diffusion analogy offers an attractive alternative to bulk layer approaches and/or simple vegetation resistance laws in models, and an in-depth comparison between different approaches at multiple sites would be interesting for future work. In any case, care should be taken for the accurate representation of the highest frequency components and the day-night transitions.
An accurate numerical approximation of heat diffusion in the vegetation and soil layers requires that both layers are discretized in multiple vertical levels with sufficiently high vertical resolution. Here, the required vertical resolution z is dictated by the temporal resolution t of the physical processes to be resolved. The relation between these two resolutions follows from the damping depth of the temperature signal This damping depth is the distance over which the amplitude is attenuated by a factor of e 1 , i.e., reduced to ≈ 37%, and a phase shift of ≈ 58 • has occurred with respect to the original values. Therefore, it sets an upper bound on the required vertical step size z. In other words, an accurate representation of, for example, the effect of clouds on a 30-min time scale (cf. Fig. 6a) requires a resolution of ≈ 0.02 m assuming κ veg ≈ 10 −6 m 2 s −1 (see Table 1) at least for the first discretization levels. Note that this resolution requirement would be even more stringent for the soil ( z ≈ 0.008 m) as the diffusivity is a factor of 10 smaller. In preliminary work by our group (Aulbers 2021), this numerical route was explored. It was shown that the observational and analytical results could be accurately reproduced using a centimetre-scale resolution in both soil and vegetation (cf. Best 1998; Met Office model).
Such an accurate representation also allows for a realistic two-way coupling between the atmosphere and the surface, viz., a coupling in which they respond to each other on the relevant time scales. Indeed, Liu and Shao (2013) found that negative feedbacks exist in large-eddy simulations (LESs) between the atmosphere and the surface on the time scale of the large eddies (≈ 10 min), but only when a fine resolution is used. On the contrary, such a feedback is not allowed when the resolution is too coarse in which case the "land surface mainly acts as a quasi-stationary external forcing to the atmospheric turbulence" (Liu and Shao 2013). Similarly, Steeneveld et al. (2006) found that a high resolution in the soil is required to model the fast dynamics of the surface heat flux in single-column model simulations of stable conditions. These fast dynamics in turn are found to be important for the development of the stable boundary layer.
With respect to atmosphere-vegetation coupling using LES, a lot of work has been done on tall canopies, where turbulent length scales are large enough to be resolved. In Brunet (2020), an extensive overview is given of LES studies making a detailed investigation of the formation and the three-dimensional structures of characteristic canopy eddies. For short canopies like grass, attention to flow-canopy interaction has been limited, because of the difficulties in both observing and explicit modelling of processes. The current sponge approach may serve as an intermediate step of model complexity, where at least aspects of signal damping and phase lagging are incorporated.

Summary and Conclusions
In the present work, we explore an alternative view on heat transfer through grass by modelling it as a simple diffusion process. As such, the vegetation layer is treated as a homogeneous sponge layer with uniform physical properties like the thermal diffusivity and the thermal conductivity. Diffusion of temperature in both vegetation and the underlying soil is modelled via Fourier's heat equation, which is solved with harmonic methods. The method is applied to observational data from the Haarweg climatological station in the Netherlands, with recordings of multiple temperature sensors below the vegetation, and of the (sub)surface energy fluxes.
Our results indicate that soil temperatures can be modelled satisfactorily from knowledge of infrared temperature dynamics of the grass surface. The attractiveness of the diffusion analogy is that it is physically consistent with respect to the frequency-response of the system: it requires no specific tuning to a daily cycle (such as with bulk methods) and hence responds to all frequencies that are represented in the input data time series (temperature forcing). This means that the system responds to quick changes in, e.g., cloud cover and rapid daynight transitions, as well as to diurnal variations of surface temperature. Each frequency will behave differently with respect to its amplitude damping and phase shift with depth.
From knowledge of the temperature dynamics in the vegetation-soil continuum, the actual surface flux at the top of the atmosphere-vegetation interface can be inferred (which is different from the soil heat flux). By comparing this flux to the other components of the surface energy balance (radiative, sensible and latent heat fluxes), it is shown that a consistent picture is obtained with respect to closure of the budget. This is particularly true for the most challenging cases with rapidly varying cloud cover and for the day-night transition in which quick changes in radiation occur on a time scale of one to two hours. In those cases, a small error in phase would lead to large errors in the energy budget closure. Therefore, we conclude that the diffusion analogy to describe heat transport through grass forms a promising alternative to existing (bulk) methods.
Yet, we realize that the description is still a crude simplification of reality: both the grass and the soil layers are complex media with time-dependent and sometimes non-uniform physical properties. Therefore, future research would benefit from an intercomparison between different methods under varying climatological conditions. In particular, a thorough analysis of parameter variability, parameter estimation, and robustness of the modelled time series would be valuable. Similar to discussions on aerodynamic roughness, the question is, if vegetation properties like effective thermal diffusivity and conductivity can be tabulated with any reasonable generality. In that context, effects of temporal variation of the moisture conditions in the biomass call for attention.
Finally, the multi-frequency response of the surface temperatures and fluxes allows for studying a class of problems where the atmospheric boundary layer and the top-surface interact on sub-hourly time scales. As such, it would be interesting to couple the diffusive grass-soil layers to turbulence-resolving methods, such as large-eddy simulations, to study rapid feedbacks between the lower atmosphere and, for example, surface heterogeneities. Herein, alternative to the present analytical harmonic methods, diffusion through grass can equivalently be modelled with numerical approaches, provided that the discretization in the grass and soil continuum is fine enough to model the frequency range of interest.  Figure 10b shows the sum of net radiation and the estimated heat flux against the sum of the turbulent fluxes. The optimized thermal parameters for this moist period are: κ soil = (3.4±0.9)×10 −7 m 2 s −1 , λ soil = (1.4±0.2) W m −1 K −1 , κ veg = (2.3±0.5)×10 −7 m 2 s −1 , and λ veg = (0.41 ± 0.05) W m −1 K −1 .
The optimized value of the soil conductivity is almost three times larger than during period 2, which stresses the effect of moisture on this parameter. Yet, the diffusivity is only 14% larger, due to the fact that the soil heat capacity also increased, which largely cancels out the conductivity effect (see Sect. 4.1). Due to the above mentioned uncertainties, care should be taken in interpreting the vegetation thermal parameters.

Appendix 2: Uncertainty Analysis Using a Monte Carlo Method
Uncertainty values on the optimized parameters and the modelled time series of temperature and heat flux are determined using a Monte Carlo method. A Monte Carlo method means that the optimization procedure of Sect. 3.3 is repeated a number of C times under perturbed input variables and input time series. This method allows propagation of uncertainty or measurement errors in the input, represented by probability distributions around their mean values, to the modelled output. Standard deviations are taken as v /3, where v is the maximum uncertainty and/or measurement error. For illustration, a relative error of 0.667% in the radiative components results in an error in the vegetation temperature T veg,0 of ≈ 0.5 K It is assumed that the input variables, such as depths of the sensors, are normally distributed N (μ v , σ v ) around their mean values μ v with standard deviation σ v . Observations are assumed to have normally distributed (relative) errors with zero mean, i.e., N (0, σ v ). When the error is given as a relative value (in percents), the observed value is increased by a randomized percentage drawn from the normal distribution. An overview of the mean values and standard deviations on the inputs is given in Table 2.
Using these values, the optimization procedure is performed C = 2500 times. In every iteration, all variables and datapoints are perturbed according to their distributions. The uncertainty in the optimized thermal parameters is taken as the standard deviation of their distributions. The mean values and the default values (as given in the main text) are almost equal. The method also gives the associated spread of the modelled time series of temperature and vegetation heat flux. The average standard deviation of this spread is given in the figure captions.
Note that, in reality, measurement errors may be correlated in time, i.e., a slight bias (or drift) could be present in a sensor over a longer period of time [∼ O(days)]. By applying an uncorrelated perturbation per datapoint, additional variability is introduced, which makes our estimate here an upper bound of the correlated case.