Advantages of analytically computing the ground heat flux in land surface models

It is generally accepted that the ground heat flux accounts for a significant fraction of the surface energy balance. In land surface models, the ground heat flux is typically estimated through a numerical solution of the heat conduction equation. Recent research has shown that this approach introduces errors in the estimation of the energy balance. In this paper, we calibrate a land surface model using a numerical solution of the heat conduction equation with four different vertical spatial resolutions. It is found that the thermal conductivity is the most sensitive parameter to the spatial resolution. More importantly, the thermal conductivity values are directly related to the spatial resolution, thus rendering any physical interpretation of this value irrelevant. The numerical solution is then replaced by an analytical solution. The results of the numerical and analytical solutions are identical when fine spatial and temporal resolutions are used. However, when using resolutions that are typical of land surface models, significant differences are found. When using the analytical solution, the ground heat flux is directly calculated without calculating the soil temperature profile. The calculation of the temperature at each node in the soil profile is thus no longer required, unless the model contains parameters that depend on the soil temperature, which in this study is not the case. The calibration is repeated, and thermal conductivity values independent of the vertical spatial resolution are obtained. The main conclusion of this study is that care must be taken when interpreting land surface model results that have been obtained using numerical ground heat flux estimates. The use of exact analytical solutions, when available, is recommended.


Introduction
An accurate estimate of the surface energy balance is very important for climate modeling and numerical weather prediction (Orth and Seneviratne, 2014).Of the three components of the net radiation (the latent, sensible and ground heat fluxes), the latent, and sensible heat fluxes provide a direct coupling of the surface energy balance to the atmosphere.For this reason, and also because typically the amplitude of the ground heat flux is smaller than the turbulent fluxes, it can be argued that climate and land surface modelers have paid more attention to an accurate estimation of these fluxes than to the ground heat flux.
However, Gentine et al. (2011) showed that the ground heat flux acts as a high-pass filter because of the strong contrast in the soil-and air-heat capacities and thermal conductivities.Because numerical solutions of the heat conduction equation can miss high-frequency fluctuations, errors in the estimation of the surface energy balance may arise.Gentine et al. (2012) showed that both models and measurements indeed miss these high-frequency fluctuations, and suggested a correction method.Wang and Bou-Zeid (2012) also noted the difficulty of accurately measuring ground heat fluxes, and used an analytical solution of the heat diffusion equation to correct ground heat flux measurements.The problem of ground heat flux estimation errors when using spatially discrete data was also shown by Lunati et al. (2012), who derived an analytical expression for the energy residual, using a cosine boundary condition at the surface.The impact of the ground heat flux parameterization on the energy balance was also studied by Russell et al. (2015).They concluded that the methods that allowed for the most variation in inputs between time steps outperformed the use of diurnal or constant input values.Other studies that intercompared ground heat Published by Copernicus Publications on behalf of the European Geosciences Union.
The problem with the numerical estimation of the ground heat flux in land surface models is that their vertical spatial resolution is too coarse to accurately estimate the soil temperature gradient.This gradient can be very steep near the soil surface, and errors in its estimation are compensated for by adopting fictitious values for the soil thermal parameters (the thermal conductivity and heat capacity).The use of analytical solutions of the heat conduction equation can be expected to partially solve this problem.
A few attempts have been undertaken to derive analytical solutions of the heat conduction equation that can easily be implemented in land surface models.A number of solutions can be found in Carslaw and Jaeger (1959).Shao et al. (1998) solved the equation analytically and compared the solution to temperature observations.Gao et al. ( 2003) derived an analytical solution as well, in order to determine thermal conductivity values.Cichota et al. (2004) compared analytical and numerical solutions for specific conditions.In an application for large-scale modeling, Bennett et al. (2008) used an analytical solution to estimate the global ground heat flux.Another example can be found in Nunez et al. (2010), where an analytical solution to model the ground heat flux was used.Wang et al. (2012) derived an analytical solution using a sine wave as a boundary condition, and Hu et al. (2016) used the Fourier transform to derive an analytical solution.What these studies have in common is that the solutions have been derived for specific initial and boundary conditions.These solutions also assume vertical homogeneity in the soil thermal parameters, which is very rarely the case.This makes it very difficult to apply them in land surface models, where an analytical expression for the boundary conditions is impossible to determine.
This paper focuses on the estimation of ground heat fluxes and soil thermal properties using a land surface model.It is first examined whether or not calibrated soil thermal properties are independent of the vertical spatial resolution of the model, if the heat conduction equation is solved numerically.An analytical solution of the heat conduction equation is then derived, with temporally varying boundary conditions, which can be applied in the model.Using this analytical solution instead of the numerical approximation, the dependence of the obtained soil thermal properties on the model spatial resolution is then investigated.

Site and data description
The data used in this study have been acquired in the framework of the AgriSAR (AGRIcultural bio/geophysical retrieval from frequent repeat pass SAR (Synthetic Aperture Radar) and optical imaging) 2006 campaign, for which the test site was located in Mecklenburg-Vorpommern in northeastern Germany, approximately 150 km North of Berlin.Pauwels et al. (2008).

Model description
For the purpose of this study, the water and energy balance model developed in Scheerlinck et al. (2009) and applied in Pauwels and De Lannoy (2011) was used.Only a short description will be provided in this section, and for a full model description we refer to Scheerlinck et al. (2009).
The model couples three physical equations.The movement of soil water in the unsaturated zone is modeled using a numerical solution of the Richards equation (Richards, 1931), which results in the profile of the pressure head (ψ; m).This equation requires the evapotranspiration as a boundary condition, which is calculated through an iteration of the surface energy balance, resulting in the surface skin temperature (Shuttleworth, 1992).This skin temperature is then used as a boundary condition for a numerical solution to the heat conduction equation, which results in the soil temperature profile (T ; K).
Table 1 lists the eleven parameters that need to be estimated.All parameters are constant in time.The soil parameters are also assumed to be homogeneous throughout the profile.We acknowledge the fact that the model represents a very strong simplification of the physical reality.However, Scheerlinck et al. (2009) and Pauwels and De Lannoy (2011) obtained excellent results for the test site with this model.For this reason, the model is deemed sufficiently realistic.
The model is applied with four different uniform vertical spatial resolutions, namely, 0.01, 0.025, 0.05, and 0.1 m.
4 The parameter estimation algorithm: particle swarm optimization

General description
The parameter estimation algorithm used in this paper, particle swarm optimization (PSO), is based on the collective behavior of individuals in decentralized, self-organizing systems.These systems are created through a population of individuals that interact locally with each other and with the community.These interactions lead to global behavior, which can result in the achievement of certain objectives.Examples of such systems in nature are abundant: ant colonies, swarms of birds, and schools of fish (Kennedy and Eberhart, 1995).For a complete description of the algorithm we refer to Scheerlinck et al. (2009).

Application in this study
In order to estimate the model parameters, observations of the net radiation (R n ), latent heat flux (LE), sensible heat flux (H ), ground heat flux (G), and soil moisture at 5 (θ 1 ), 9 (θ 2 ), 15 (θ 3 ), and 25 cm (θ 4 ) are used.The energy balance variables are in W m −2 , and the soil moisture values are dimensionless.A global objective function is calculated: The RMSE (root mean square error) values for each variable are calculated as where x o and x s are the observed and simulated values, respectively, and σ x is the standard deviation of the variable.
The global objective function is then minimized through the application of PSO; 36 iterations are performed, ensuring convergence of the algorithm, and the method is repeated 24 times.In order to ensure an analysis of the most optimal parameter values, the parameter sets corresponding to the eight lowest objective function values are retained for further analysis.More specifically, the average parameter and objective function values over these eight repetitions are examined.

Model calibration using the numerically calculated ground heat flux
The model simulations resulting in the lowest objective function values for the different spatial resolutions will be analyzed in this section.Figure 1 shows the comparison between the modeled energy balance components and the observations, for the spatial resolution of 0.01 m.Table 2 shows the statistics of the linear regressions between the observed and simulated energy balance components for the four spatial resolutions.Figure 2 shows the comparison between the modeled and the observed soil moisture values for the four spatial resolutions.From these figures and tables, it can be concluded that the model adequately simulated the water and energy balance processes, and that the results are very similar for the four different resolutions.
Table 3 shows the parameter values obtained from the model calibration with a spatial resolution of 0.01, 0.025, 0.05, and 0.1 m.In order to determine which parameters are significantly effected by the model spatial resolution, we applied a t test to the slopes from the linear regressions between the spatial resolution (x axis) and the parameter values (y axis), at the 95 % confidence level.In other words, it is tested whether or not there is a significant linear trend between z (the model resolution) and the parameter values.The objective function value has been found to change significantly with the resolution, as well as λ, ψ c , K s , f , and κ.All the other parameters are not significantly dependent on the spatial resolution.Of all the parameters affected by the resolution, the parameter that shows the largest variation in values is the thermal conductivity κ, with the value at z = 0.1 being more than 4 times the value at z = 0.01 m.No other parameter shows such a dramatic variation.
This result can be explained by the independence of the albedo (α) and the parameters determining the roughness lengths and zero-plane displacement height (f v , f h , and f d ) on the model spatial resolution.Since these parameters are similar for the different resolutions, and the model is calibrated using R n and the three heat fluxes, the resulting skin temperatures are similar for the different resolutions.A similar skin temperature leads to similar sensible and ground heat fluxes, as is proven in Table 2.The ground heat flux is defined as the gradient in the soil temperature multiplied by the thermal conductivity.Since the skin temperatures and the ground heat fluxes (the latter used in the calibration) are equal, but the model spatial resolution is different, the same ground heat flux can only be obtained with a different thermal conductivity.Because the finest spatial resolution will result in the steepest gradient of the soil temperature, it can be expected that this will result in the lowest thermal conductivity value.Table 3 shows that this is indeed the case.This renders the physical interpretation of the thermal conductivity values impossible.In order to solve the issue related to the dependence on the grid resolution in the use of a numerical solution of the heat conduction equation, we propose the use of an analytical solution.First, the steady-state temperature profile for a constant temperature at the bottom (T b,0 ) and top (T u,0 ) of the profile is calculated.The depths of the top and bottom of the profile are denoted as z u and z b , respectively.This solution is then used as the initial condition for the same equation, now with different bottom (T b,1 ) and top (T u,1 ) temperatures as boundary conditions.It should be clarified that the time t starts at zero for this new solution.The temperature profile at time t is then calculated, and used as the initial condition for the same equation, again with different temperatures (T b,2 and T u,2 ) as boundary conditions, and time starting at zero.The solution at time t is then again used as the initial condition for the same equation with different boundary conditions (T b,3 and T u,3 ) and time starting at zero, and so on.
Using this recursive methodology, the temperature profile for the Mth time step can be written as where T u,M and T b,M are the top and bottom temperatures (boundary conditions) for the Mth time step, with M > 0, and The surface ground heat flux then becomes (5) Appendix A describes the details of the derivation of this solution as well as a methodology to apply these equations in a computationally efficient manner.
It should be noted that with this analytical solution it is no longer necessary to calculate the soil temperature profile in order to calculate the ground heat flux.In the original model formulation, the heat conduction equation needed to be solved numerically, using the surface skin temperature as a boundary condition, so the temperature of the first soil layer could be calculated, and the ground heat flux could be computed.However, Eq. ( 5) shows that this first layer temperature is no longer a variable in the calculation of the ground heat flux.More complex land surface models may contain parameters that depend on the soil temperature profile, and thus would require the application of Eq. (3).However, in this model, this is not required.

Comparison to the numerical solution
A synthetic test case is used to intercompare the analytical and the numerical solutions.Equations (3) and ( 5) are applied to a soil column of 1 m depth, with a thermal conductivity of 0.5 W m −1 K −1 , and a heat capacity of 2.5 × 10 6 J m −3 K −1 .In a first application, the spatial resolution is 0.01 m, and the temporal resolution is 60 s.The numerical solution thus needs to be applied 60 times per 1 h time step in this case.
In a second application, a spatial resolution of 0.05 and a temporal resolution of 1 h (3600 s) are used.The air temperature values from the AgriSAR data set are used as boundary conditions at the top of the profile, and the bottom boundary conditions are assumed to linearly increase from 3 to 13 • throughout the simulation.
Figure 3 shows the comparison of the temperature profiles obtained from the numerical and analytical solutions, for the test case with the fine spatial and temporal resolutions.An excellent agreement between both methods can be seen.Figure 4 shows the same comparison for the coarse resolutions.In many of these profiles, especially when sharp changes of temperature occur, a strong deviation of the numerical solution from the analytical solution can be observed.Since these coarse resolutions correspond to values typically used in land surface models, this leads to the conclusion that care must be taken when interpreting ground heat flux simulations from these models.This is demonstrated in Fig. 5, in which the ground heat fluxes from both solutions are compared.For the fine resolutions, a very good agreement is obtained, but relatively strong differences between both methods are found when coarse spatial and temporal resolutions are used.

Model calibration using the analytically calculated ground heat flux
Figure 6 shows the comparison between the modeled energy balance terms and the observations, for the simulations with a spatial resolution of 0.01 m, and for the parameter set corresponding to the lowest objective function value.A similar model performance as obtained with the numerical solution of the heat conduction equation is achieved.Table 4 shows the results of the linear regressions between the modeled energy balance terms and the observations, again obtained using the parameter values corresponding to the lowest objective function value, for the four different spatial resolutions.
Comparing Tables 4 to 2 leads to the conclusion that the energy balance terms are simulated practically identically when the model uses the numerical or the analytical solution of the heat conduction equation.
Figure 7 shows the comparison between the simulated soil moisture values and the observations for the same parameter sets.The comparison of Figs. 7 and 2 shows that the analytical and numerical solutions of the heat conduction equation lead to very similar simulations of the soil water balance as well.
Table 3 shows the average parameter values from the eight retained PSO results.A slope t test for the linear regressions between the spatial resolutions (x axis) and the parameter values (y axis), at the 95 % confidence level, showed that the objective function value changes significantly with the resolution, as well as λ, ψ c , K s , and f .The only difference with the results obtained with the numerical solution of the heat conduction equation is that κ is no longer dependent on    the spatial resolution.This could be expected, because the expression for the ground heat flux is not dependent on the temperature of the first soil layer, and thus on the spatial resolution.Furthermore, the value for the heat capacity (C) is now less variable than when the ground heat flux was calculated numerically.More specifically, the standard deviation has been reduced from 557 943 to 278 984 J m −3 K −1 .A pooled variance t test with 95 % confidence was applied to the parameter values obtained with the analytical and numerical solutions, to investigate which parameters are significantly different.This test showed that the only parameters that are significantly different are the objective function value and the heat capacity for all spatial resolutions, and the thermal conductivity for a spatial resolution of 0.01 and 0.1 m.Interestingly, the objective function values obtained with the analytical solution are slightly higher than with the numerical solution.However, since this objective function is composed of eight rescaled RMSE values, the resulting model performance is very similar, as is shown in Table 4 and Fig. 7. Since the thermal conductivity values no longer depend on the spatial resolution, on which they depended when the numerical solution was used, it can be expected that the obtained thermal conductivity values are different.Because of the different manner of solving the heat conduction equation, it can also be expected that the obtained soil-heat capacity values are different as well.The key conclusion from these simulations is that the overall model performance is independent of the type of calculation of the ground heat flux (analytically or numerically), but that the results of the model calibration are more robust (i.e., independent of the spatial resolution) if an analytical solution of the ground heat flux equation is used.

Conclusions
A water and energy balance model, using a numerical solution of the heat conduction equation, has been calibrated against energy balance and soil moisture observations, for four different vertical spatial resolutions (0.01, 0.025, 0.05, and 0.1 m).It has been found that a number of parameters are dependent on this resolution, with the soil thermal conductivity values showing the largest dependence.An analytical solution of the heat conduction equation has then been derived, allowing for the bottom and top boundary conditions (i.e., the bottom and surface skin temperatures) to vary over time.Using this analytical solution has the advantage that the soil temperature profile no longer needs to be computed.For fine spatial and temporal resolutions the analytical and numerical solutions cannot be distinguished, but different results are obtained for resolutions typically used in land surface models.When the ground heat flux is calculated using this analytical solution, and the model is calibrated, the obtained soil thermal conductivity is no longer dependent on the model spatial resolution.Furthermore, the variability in the obtained soil-heat capacity is also strongly reduced.The results in this paper indicate that a similar model performance is obtained when the ground heat flux is calculated analytically or numerically.However, the calibration is more robust, and the parameter values more physically interpretable, if the analytical solution is used.One must thus be careful when using numerical solutions of the heat conduc-tion equation in land surface models, and preference should be given to the use of analytical solutions.The solution derived in this paper does not allow for temporally varying soil thermal properties, and ongoing research is focusing on the derivation of an analytical solution that is straightforward to apply in land surface models in these conditions.
Appendix A: Analytical solution A1 The governing equation A solution of the heat conduction-convection equation is derived first, since this equation is analytically more straightforward to solve than the heat conduction equation because of the easier inversion from the Laplace domain.Furthermore, this general solution can be used for purposes outside the scope of this paper.The limit case with zero convection is then calculated.The governing equation is where t is time (s), z is the depth positive upwards (m), T is the soil temperature (K), C is the volumetric heat capacity of the soil (J m −3 K −1 ), κ is the soil thermal conductivity (W m −1 K −1 ), and v is the water velocity (m s −1 ; positive upwards).We assume the parameters are uniform throughout the profile, so the equation becomes A2 Steady-state solution In order to obtain a realistic initial condition, we will calculate the steady-state solution.For example, the profile at the end of a very long, hot day.The equation becomes We will calculate the poles of this equation (the values for y for which the denominator is zero), and calculate the residuals of each pole.The analytical solution is then simply the sum of the residuals (Brutsaert, 1994).Following this theorem, writing F (y) as P (y)/T (y), we can write

Figure 1 .
Figure 1.Comparison between the modeled and the observed energy balance terms for the simulation with z = 0.01 m and a numerical solution of the heat conduction equation.

Figure 3 .
Figure 3.Comparison between soil temperature profiles obtained using the numerical (+) and analytical (solid line) solutions for fine spatial and temporal resolutions.Both solutions are very similar and therefore difficult to distinguish.

Figure 5 .
Figure 5.Comparison of the resulting ground heat fluxes from the fine and coarse spatial and temporal resolutions to the analytical solution.

Figure 6 .
Figure 6.Comparison between the modeled and the observed energy balance terms for the simulation where z = 0.01 m and an analytical solution of the heat conduction equation (Eq.A31).

Figure 7 .
Figure 7.Comparison between the modeled and the observed soil moisture values for the four different spatial resolutions, and an analytical solution of the heat conduction equation.
conditions T = T b,0 , z = z b T = T u,0 , z = z u .(A4)The solution of this equation isT = T b,0 + T b,0 − T ufirst new boundary conditionsWe will use a constant time steps t.For the first new boundary conditions, we will use the steady-state profile as the initial condition.The boundary conditions areT = T b,1 , z = z b T = T u,1 , z = z u .(A6)We will solve this equation through a Laplace transform.We will denote the transform of T (z, t) as F (z, y), with y the Laplace variable.The differential equation becomesCyF − C T b,0 + T u,0 − T b,this differential equation is F = T b,1 − T b,0 e b(z−z b ) y sinh b 2 + Cy κ (z − z u ) sinh b 2 + Cy κ (z b − z u ) + T u,1 − T u,0 e b(z−z u ) y sinh b 2 + Cy κ (z − z b ) sinh b 2 + Cy κ (z u − z b ) has poles for y equal to zero and for the hyperbolic sine equal to zero.For y equal to zero this simply becomesT 1 (z) = T b,1 − T b,0 e b(z−z b ) sinh (b (z − z u )) sinh (b (z b − z u )) + T u,1 − T u,0 e b(z−z u ) sinh (b (z − z b )) sinh (b (z u − z b ))+ T u,0 − T b,0 the hyperbolic sine, we defineb 2 + Cy κ (z u − z b ) = j x,(A14) Hydrol.Earth Syst.Sci., 20, 4689-4706, 2016 www.hydrol-earth-syst-sci.net/20/4689/2016/

Table 1 .
The model parameters that need to be estimated.

Table 2 .
Results of the linear regressions between the energy balance observations (x axis) and the simulations (y axis) for the simulations with a numerical solution of the heat conduction equation.Units are W m −2 .

Table 3 .
Averages of the eight lowest objective function values and the corresponding parameter values.

Table 4 .
Results of the linear regressions between the energy balance observations (x axis) and the simulations (y axis) for the simulations with an analytical solution of the heat conduction equation.Units are W m −2 .