GEOtop 2 . 0 : simulating the combined energy and water balance at and below the land surface accounting for soil freezing , snow cover and terrain effects

GEOtop is a fine-scale grid-based simulator that represents the heat and water budgets at and below the soil surface. It describes the three-dimensional water flow in the soil and the energy exchange with the atmosphere, considering the radiative and turbulent fluxes. Furthermore, it reproduces the highly non-linear interactions between the water and energy balance during soil freezing and thawing, and simulates the temporal evolution of the water and energy budgets in the snow cover and their effect on soil temperature. Here, we present the core components of GEOtop 2.0 and demonstrate its functioning. Based on a synthetic simulation, we show that the interaction of processes represented in GEOtop 2.0 can result in phenomena that are significant and relevant for applications involving permafrost and seasonally frozen soils, both in high altitude and latitude regions.


Introduction
Frozen soil and snow cover interact in various ways with hydrology, climate, ecosystems and with human infrastructures.These natural systems are complex and characterised by many non-linear processes that operate and interact over different scales.Their mathematical representation and quantification is gaining in importance, especially in the light of global climate change.This importance derives on one hand from the requirement to study more and more complex systems, and, on the other hand, this representation of more complex systems can inform decisions about their simplification (Freeze and Harlan, 1969).In fact, the systems of equations required for representing such environments are often simplified by excluding processes that are considered less important for the problems addressed.Such an a priori exclusion, however, may not be quantitatively justified and mostly dictated by the need for mathematical tractability.Estimating the error inherent in model simplifications is therefore desirable for weighing the costs and benefits of differing options.
There is a great diversity of models (understood here as mathematical representations of one or more processes) and simulators (computer programs, usually comprising implementations of several models to represent a natural system) to simulate cold-region processes.For example, models applied in permafrost environments are normally: (i) models applied at a local, regional, and continental scale that integrate a one-dimensional form of heat and water flow equation with phase change and predict the evolution of the depth of thaw; (ii) hydrological models commonly applied at a large scale that predict river discharge without accurately describing the coupling with the soil energy balance; and (iii) models that very accurately describe and couple water and energy subsurface flow in frozen soil, but do not consider the S. Endrizzi et al.: Simulating the combined energy and water balance with soil freezing, snow and terrain heat flux exchanged with the atmosphere.Table 1 provides an overview of some common frozen soil and permafrost models.
GEOtop 2.0 is a process-based model developed from the blueprint GEOtop 0.75 described in Rigon et al. (2006) and Bertoldi et al. (2006).It was built on the idea that the combination of terrain effects, energy and water balance produce unique results for different meteorological forcings, which makes difficult the a priori exclusion of any of the processes.GEOtop covers the full spectrum of hydrological fluxes, from the energy balance in the complex terrain to snow and vegetation.Therefore, it makes possible the modelling of the interactions between several hydrological, cryospheric, ecological and geomorphological processes in an interdisciplinary research framework.For example, it has been applied to model geomorphological processes like landslide triggering (Simoni et al., 2008), ecohydrological processes (Bertoldi et al., 2010;Della Chiesa et al., 2014;Kunstmann et al., 2013) and peat hydrology (Lewis et al., 2012).However, this paper mainly focuses on the aspects related to the cryosphere.In this context, permafrost research provides an intersection of phenomena related to frozen soil, the flow of water, and the snow pack.
As a novel combination of processes represented, this paper presents GEOtop 2.0, an improved version of the opensource software GEOtop, which now simulates the energy and water balance at and below the land surface, soil freezing, snow cover dynamics, and terrain effects.It is a research tool for studying, for example, the hydrological and thermal phenomena at locations that differ in soil types and topography to specific climatic forcings.Output consists of variables such as temperature, water and ice contents, or of integrated variables such as stream discharge.The software operates in pointwise and distributed modes and can be flexibly controlled, because all relevant parameters that govern e.g.discretisation, input/output or numerics can be set via keywords.
GEOtop describes the evolution in time of temperature and water content in the soil and snow cover and is driven by meteorological forcings.This is accomplished by solving the heat and water flow equations with boundary conditions accounting for the interactions with the atmosphere at the surface in terms of energy and water fluxes.The solution of the equations is obtained numerically in the soil domain and snow cover.
GEOtop 2.0 is significantly different from GEOtop 0.75.It includes a fully three-dimensional description of the Richards equation, whereas in the previous version the equation was only solved in the vertical direction and the lateral flow was parameterised, in a similar way as in large-scale land surface models.In the new version, a multilayer snow cover and the surface energy balance are fully integrated in the heat equation for the soil, which is solved with a rigorous numerical method based on Kelley (2003), while in the previous version, snow cover was described with a bulk method (Zanotti et al., 2004) and the surface energy balance, though complete in its components and accommodating complex terrain, was not numerically coupled to the soil heat equation.In GEOtop 2.0 (hereafter GEOtop), soil freezing and thawing are represented, meteorological forcings are distributed, and channel routing is described as overland flow with the shallow water equation neglecting the inertia.The description of vegetation with a double-layer surface scheme in order to more accurately represent the heat and vapour exchanges of vegetation with the soil surface and the atmosphere has also been included in GEOtop and is described in Endrizzi and Marsh (2010).The code of GEOtop is publicly available in the terms described in Appendix A.
The core components of GEOtop are here presented.The description will particularly consider the soil volumetric system and the equations to be solved, the interaction with the atmosphere, the effects of complex terrain, the numerics, the representation of the snow cover, and the distribution of the meteorological data.It is shown that the simulator produces plausible results in its major components.In addition, a simulation experiment is presented in order to demonstrate that the combination of terrain effects, energy and water balance produce unique results for different meteorological forcings, making an a priori exclusion of any of these processes difficult and providing one important rationale for developing and using a simulator such as GEOtop.

Volumetric system
The volumetric system consists of a soil volume of a userspecified uniform depth (typically a few metres to hundreds of metres) and is discretised in several layers parallel to the surface.Close to the surface the layers are usually prescribed thinner than at depth, as the gradients of temperature and water content resulting from the interaction with the atmosphere are stronger.The surface can be additionally Table 1.Overview of existing simulators dealing with soil-freezing processes.Abbreviations: in the snow column: dd = degree day factor, ly = one-layer energy balance, mly = multilayer energy balance; in the soil energy and water balance columns: 0d = simplified, 1d = onedimensional (vertical) solution of the heat or Richards equation, 3d = three-dimensional solution; in the column about the energy exchange with atmosphere: wiT = with complex topography, woT = without complex topography.discretisedspatially using a regular square grid.Therefore, the elementary units, here referred to as "cells", are given by the volumes resulting from the intersection of layers parallel to the surface with the columns defined in a direction normal to the surface.The heat and water flow equations are not fully coupled numerically, but they are linked in a time-lagged manner (e.g.Panday and Huyakorn, 2004).This method allows keeping the complexity of the numerics moderate, while the equations are solved reasonably fast.

Heat equation
As explained in Appendix B, the equation representing the energy balance in a soil volume subject to phase change is where U ph is the volumetric internal energy of soil (J m −3 ) subject to phase change, t (s) time, ∇• the divergence operator, G the heat conduction flux (W m −2 ), S en the energy sink term (W m −3 ), S w the mass sink term (s −1 ), L f (J kg −1 ) the latent heat of fusion, ρ w the density of liquid water in soil (kg m −3 ), T ( • C) the soil temperature and T ref ( • C) the reference temperature at which the internal energy is calculated.Writing G according to the Fourier's law and considering Eq. (B9), Eq. (1) becomes where λ T is the thermal conductivity (W m −1 K −1 ) and u f is defined in Eq. (B9).Equation ( 2) is numerically solved one-dimensionally neglecting the lateral gradients, and is integrated assigning the heat fluxes at the upper and lower boundaries of the domain.The upper boundary is given by the interface with the atmosphere or snowpack.At the lower boundary an energy flux is prescribed.This can be assigned externally as a parameter and, depending on the conditions and depth of the soil column, this depends on terrain geometry and transient effects that often overprint the deep geothermal heat flow locally (Gruber et al., 2004).The sink term S en can also be assigned externally.
Since the total mass of water is kept constant and is given by the resolution of the water balance equation in time lagged manner, the unknown of Eq. ( 2) is T .However, the equation determines the mass that changes phase.Since ice has a lower density than liquid water, freezing would lead to unrealistically large gauge pressures that cannot be converted into an expansion of the soil matrix, due to the lack of a mechanical model.Therefore, similarly to Dall'Amico et al. (2011a), a rigid soil scheme is assumed, which implies that no volume expansion during freezing is allowed, and the densities of ice and liquid water are equal, and set to 1000 kg m −3 .
The expression of dU ph , defined in Eq. (B12), implies a proper description of soil freezing and thawing processes.Water phase change from liquid to solid state in the soil is not an isothermal process like in free-surface water (e.g.Wettlaufer and Worster, 2006).Rather, phase change occurs over a range of temperatures.Several authors (e.g.Spaans and Baker, 1996) have defined fixed relations between unfrozen water content and temperature, referred to as "freezing soil characteristic curve".This is a simplification, since more complex behaviours have been observed (Koopmans and Miller, 1966).The ice volumetric content can then be calculated as the difference of total and unfrozen water contents.The definition of the freezing soil characteristic curve allows expressing dU ph in the following way: where θ w (-) is the volumetric fraction of liquid water in the soil, dθ ph w its variation due to phase change, C the volumetric heat capacity (J m −3 K −1 ) defined in Appendix B, c i and c w (J kg −1 K −1 ) are the specific thermal capacity of ice and liquid water, respectively, and C a is referred to as apparent heat capacity, function of temperature, defined as

Thermal conductivity
The thermal conductivity (λ T ) in Eq. ( 2) is a combination of the thermal conductivities of each component of the soil multiphase mixture (λ sp for soil particles, λ i for ice, λ w for liquid water, and λ a for air).It is, therefore, a non-linear function of temperature, since the proportion of liquid water and ice contents depends on temperature.While a simple additive mixing law is exact for the heat capacity, the behaviour of a multiphase mixture concerning the thermal conductivity is much more complex.Several non-linear mixing laws have been proposed (e.g. de Vries, 1963;Johansen, 1975;Balland and Arp, 2005).GEOtop uses the one proposed by Cosenza et al. (2003), which was derived in analogy to the dielectric permittivity, namely where θ sp (-) is the soil porosity, θ i (-) the volumetric fraction of ice, and θ a (-) the volumetric fraction of air and gaseous components.For continuity it is

Soil freezing characteristic curve
Dall'Amico et al. (2011a) derived the soil freezing characteristic curve from the soil water retention curve using the Van Genuchten parameters (Van Genuchten, 1980).They assumed a rigid soil scheme and use the "freezing = drying" assumption (Miller, 1965), which implies that: (i) the freezing (thawing) water is like evaporating (condensing) water; (ii) the ice pressure is equal to the air pressure; (iii) the water and ice content in the soil are related to the soil water retention curve.However, the assumption that ice is always at the air pressure may be restrictive in permafrost modelling, since ice pressure at depth may be significantly high.Nevertheless, this can be extended.Instead of assuming that ice is at the air pressure, it can be more generally supposed that liquid water is not subjected to external pressures, which, on the other hand, are completely supported by the soil matrix and the ice.This entails that, when pore water is subjected to an external pressure (e.g.hydrostatic) it starts to freeze, and the liquid water is completely unloaded of this pressure once the first ice is formed.Therefore, liquid water pressure would unrealistically undergo a pressure jump when freezing starts.
However, since the focus of GEOtop is to simulate soil temperature and moisture dynamics (and not ice pressure), this is deemed reasonable.

Water flow equation
As explained in Appendix B, the system of equations representing the water balance in the soil is where dθ ph w (-) is the fraction of liquid water content in soil subject to phase change, dθ fl w (-) is the fraction of liquid water content transferred by water flux, ρ i the density of ice (kg m −3 ), with θ i (-) the fraction of ice in soil and J w (m s −1 ) the flux of liquid water.This equation describes the water flow occurring below the soil surface (subsurface flow) and is normally referred to as the variably saturated Richards equation.According to Darcy's law, J w can be written as where K (m s −1 ) is the hydraulic conductivity, ψ (m) the liquid water gauge pressure head and z f (m) the elevation head, i.e. the elevation above a reference level.When ψ is positive, water pressure is higher than the atmospheric pressure, and soil is saturated.When ψ is negative, soil is unsaturated.According to (Dall'Amico et al., 2011a), in variably saturated conditions it is where ψ T (T ) is the temperature-dependent soil matric potential determining the contribution of freezing below the melting temperature T * ( • C), and ψ w0 (m) is the matric potential corresponding to the total water content.The function ψ T (T ) is defined in Dall'Amico et al. (2011a).As described in Sect.2.1.2,it is assumed that when soil is freezing, the external pressure is completely carried by the ice.When soil is unsaturated (ψ < 0), the water content θ w is calculated by means of the soil water retention curve according to the Van Genuchten (1980) model.When soil is saturated (ψ ≥ 0), θ w should always be equal to the saturated value.However, a biunique relation between θ w and ψ is needed in order to numerically solve the equation.Therefore, the concept of specific storativity (S s in m −1 ) is used (Ray, 1996), which is defined as the volume of water added to storage, per unit volume and per unit rise in pressure head.Therefore, it is where θ r (-) is the residual water content, and α (m −1 ), n (-) and m (-) are the soil-specific parameters in the model of Van Genuchten (1980), normally referred to as Van Genuchten parameters.
Defining H (m) as the sum of the pressure and potential heads: the second part of Eq. ( 7), combined with Eq. ( 8), becomes Since temperature and ice content are kept constant (given by the resolution of the energy balance equation in timelagged manner), the matric potential ψ is just a function of ψ w0 (defined in Eq. 9), which is eventually the only unknown of Eq. ( 12).
Analogous to the case of water content controlled only by drying processes, the hydraulic conductivity K is dependent on the soil matric potential ψ associated with liquid water (Mualem, 1976).However, the presence of ice may significantly reduce the hydraulic conductivity due to the apparent pore blockage effect exerted by ice: this is accounted for by further reducing the hydraulic conductivity by an impedance factor smaller than 1 and equal to 10 −ωq (Hansson et al., 2004;Kurylyk and Watanabe, 2013), where ω is a coefficient and q is the ice fractional content given by θ i /(θ s − θ r ).
Equation ( 12) is solved in a fully three-dimensional way in order to describe the two gradients of H in the direction parallel and normal to the surface.Once the soil becomes saturated as a result of a precipitation or melting snow, normal gradients may become very small in comparison to those in the parallel direction, which, in turn, are responsible for the routing of water through the soil.

Overland flow
The surface (or overland) water flow must also be considered to consistently describe the water balance in the soil and the runoff mechanisms.This process is described with the approximation proposed by Gottardi and Venutelli (1993), who extended to the surface flow the validity of Darcy's law, which, strictly speaking, would not be valid with the flow being turbulent.Using the water conservation and Darcy's law for the overland flow, the surface water balance can be written as where ψ| z=0 and z f | z=0 are respectively the liquid water pressure head and the elevation head at the soil surface, K sur (m s −1 ) the conductance, in analogy to K defined in Eq. ( 8), and P e (m s −1 ) the effective precipitation per unit horizontal surface that reaches the soil surface, including snowmelt flow and deducting evaporation from the soil.The variable ψ| z=0 cannot be negative in this equation and is written in place of the water depth above the surface.Following Gottardi and Venutelli (1993) the conductance is where s (m) is the length along the direction of maximum local slope, c s the surface roughness coefficient (m 1−γ s −1 ) and γ an exponent between 0 and 1 that varies according to the formulation of c s .For example, in the formulation of Manning it is c s = n −1 r and γ = 2 3 , where n r is the Manning coefficient.In the formulation of Chezy it is c s = C r and γ = 1 2 , where C r is the Chezy coefficient).Equation ( 13) actually works as a boundary condition at the soil surface for Eq. ( 12).

Numerics
In order to reduce the complexity of the numerical method, Eqs. ( 2) and ( 12) are linked in a time-lagged manner, instead of solving them in a fully coupled way.Both equations have the same form, which can be generalised as where χ is the unknown function of space and time, F a nonlinear function of the unknown (corresponding to the internal energy content for the heat equation and the total water content for the water flow equation), S the sink term and κ a conductivity function of the unknown.
All the derivatives are discretised as finite differences.Therefore, the following relation is obtained: where the equation is written for the generic ith cell; n represents the previous time step, at which the solution is known, n + 1 is the next time step, at which the solution is unknown.t is the time step, j is the index of the M adjacent cells with which the ith cell can exchange fluxes, m represents a time instant between n and n+1, κ ij the conductivity between the cell i and j , D ij the distance between the centres of the cells i and j , S i the sink terms, and G i the residual that is to be minimised for finding a solution.Equation ( 16) is a system of N equations, and the second term of the left-hand side is the sum of the fluxes exchanged with the neighbouring cells.The variables at the instant m are represented with a linear combination between the instant n and n+1.If m = n the method is fully explicit and unstable, if m = n + 1 2 the method has a second-order precision but might not be always stable, and if m = n+1 the method has a first-order precision but is unconditionally stable.Since there are more concerns on stability than precision, the latter is the chosen method.
A solution of Eq. ( 16) is sought with a special Newton-Raphson method, with the following sequence (Kelley, 2003): where χ is the vector χ i that appears in Eq. ( 16), d denotes the Newton direction, and λ d is a scalar smaller than or equal to 1 referred to as path length, found with a line search method like the Armijo rule (Armijo, 1966).The quantity λ d d (χ n ) is also referred to as the Newton step.The Newton direction is obtained solving the following linear system: where G is the vector G i that appears in Eq. ( 16) and G (χ ) denotes the Jacobian matrix G (χ ) ij = ∂G i (χ ) /(∂χ j ).
If Eq. ( 15) is solved neglecting the lateral gradients, the number of adjacent cells that are actually considered is maximum 2 (i.e the cell below and above).Therefore, the matrix G (χ n ) is tridiagonal and symmetric, and then invertible with simple direct methods (El-Mikkawy and Karawia, 2006).On the other hand, if Eq. ( 15) is solved fully threedimensionally M can be up to 6, and, therefore, G (χ n ) is a symmetric and sparse matrix.Its inversion is a more complex problem (Niessner and Reichert, 1983).In this case, the linear system in Eq. ( 18) is solved approximatively with an iterative method, the BiCGSTAB Krylov linear solver (Van Der Vorst, 1992).This iterative process becomes an inner iteration, nested in the outer iteration defined in Eq. ( 17).

Energy exchange with the atmosphere
The heat flux exchanged with the atmosphere (S), hereafter referred to as "surface heat flux", is given by the sum of net shortwave (solar) radiation (SW), net longwave radiation (LW), and turbulent fluxes of sensible (H ) and latent heat (LE), namely The surface heat flux is dependent on the temperature of the surface, which is, in turn, the unknown of the equation.In addition, the latent heat flux also depends on the soil moisture at the surface, which is a further coupling term to the water flow equation.All the fluxes in Eq. ( 19) are positive if they are directed towards the surface.The following section discusses how the components of the surface heat flux are calculated in the simple case of horizontal flat terrain.Then, in Sect. 4 the case with complex terrain is presented.Depending on the input data available, radiation components can be either assigned directly at input or calculated by the model.

Shortwave radiation
The net shortwave radiation appearing in Eq. ( 19) is a balance given by the incoming radiation SW in from the atmosphere and the reflected radiation SW out , which is given by SW in multiplied by the broadband albedo.
Incoming shortwave radiation on a flat ground surface is the result of the top-of-atmosphere (SW toa ) shortwave radiation, and atmosphere and cloud transmissivities (respectively τ a and τ c ): While SW toa can easily be expressed with analytical formulae depending on solar height and azimuth (e.g.Iqbal, 1983), the transmissivities are more complex and uncertain to calculate.Their calculation is fully described in Appendix D1.
Albedo is treated differently according to whether the ground surface is snow free or snow covered.In the former case the albedo varies linearly with the liquid water contents of the top soil layer, while, in the latter, the formulation of Dickinson et al. (1993) is used.This formulation (i) accounts for the decrease of the snow reflectance with the time from the last significant snowfall, (ii) partitions the spectrum into visible and near-infrared components and considers different coefficients, (iii) considers an increase of albedo at lower sun angles as a result of the Mie scattering properties of snow grains (Hock, 2003).In addition, snow albedo is decreased for shallow snowpack since a significant portion of incoming shortwave radiation is actually absorbed by the soil surface (Tarboton and Luce, 1996).

Longwave radiation
The net longwave radiation in Eq. ( 19) is a balance of the component LW in coming from the atmosphere and LW out emitted by the surface.Differently from the shortwave radiation, the two components are independent and calculated separately.
The incoming longwave radiation at the surface is the integrated result of the radiation emitted at different levels in the atmosphere with different temperatures and gas concentrations.Clear-sky radiation is calculated with one of the several empirical formulations present in the literature (e.g.Brutsaert, 1975a;Satterlund, 1979;Idso, 1981;Konzelmann et al., 1994;Dilley and O'Brien, 1997), which in general apply the Stefan-Boltzmann law using the air temperature measured at the surface (T a in K) with an effective atmosphere emissivity a (-) dependent on air temperature T a and water vapour pressure e a (bar), namely where σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ).The relations in the literature differ on the expression of a (e.g.Brutsaert, 1975a;Idso, 1981;Konzelmann et al., 1994;Prata, 1996;Dilley and O'Brien, 1997).In cloudy skies the emissivity of the atmosphere is increased to a value c , which may be significantly higher than a .There is a rather high uncertainty in the formulations for c , which is evident in the spread of the results if different empirical formulations are used (e.g.Gubler et al., 2012;Flerchinger et al., 2009).In particular, the problem is more evident in the case of intermediate fractional cloud covers.Most formulations for c use information on the cloud-covered fraction of the sky (Deardorff, 1978), which is mostly obtained with visual observation.As an alternative, Crawford and Duchon (1998) proposed a direct relation between the shortwave radiation cloud transmissivity (τ c ) and c : which provides a linear interpolation between the clear-sky emissivity value and the black-body emissivity in the ideal case of cloud cover completely obscuring the ground.This relation is used in GEOtop since it provides a direct estimation of c from incoming shortwave radiation without using the cloud-cover fraction of the sky.However, τ c is primarily affected by the cloud cover around the solar disc and not much by cloud cover far from it.Sicart et al. (2006) showed that this bias is eased if relative long time averages (such as daily) of τ c are taken.In GEOtop this average time is set as a parameter, but normally ranges between 4 h and 1 day.Longer average intervals reduce the cloudiness directionality bias, while shorter ones allow an estimation of the day evolution of the cloud cover.This feature is particularly useful for the estimation of τ c during the night, which is estimated with a linear interpolation between the value before sunset and after sunrise (Gubler et al., 2012).
The outgoing longwave radiation (LW out ) emitted by the surface can also be calculated with the Stefan-Boltzmann law: where T sur is the temperature of the surface (K) and s is the emissivity of the surface.The latter is set as a parameter if the surface is snow covered.Otherwise, similarly to albedo, it is interpolated between a dry and wet value according to the water content of the first soil layer (Snyder et al., 1998).

Turbulent fluxes
The turbulent fluxes of sensible (H ) and latent heat (LE) are calculated with the flux-gradient relationship (Brutsaert, 1975b;Panoksky and Dutton, 1984;Garratt, 1992): where ρ a is the air density [kg m −3 ], c p the specific heat at constant pressure (J kg −1 K −1 ), w s the wind speed (m s −1 ), L e the specific heat of vaporisation (J kg −1 ), Q s the saturated specific humidity (kg kg −1 ) at the surface, Q a the specific humidity of the air, and r a the aerodynamic resistance (-).The α YP and β YP coefficients take into account the soil resistance to evaporation, and only depend on the liquid water pressure close to the soil surface.They are calculated according to the parameterisation of Ye and Pielke (1993), which considers evaporation as the sum of the proper evaporation from the surface and diffusion of water vapour in soil pores at greater depths.The aerodynamical resistance is obtained applying the Monin-Obukhov similarity theory (Monin and Obukhov, 1954), which requires that known values of wind speed, air temperature and specific humidity are available at least at two different heights above the surface.Known values at only one height above the surface are sufficient if it is assumed that just above the surface (properly at zero height above the surface): (i) the value of air temperature is equal to the value of soil temperature at the surface (this assumption also leads to the boundary condition non-linearity), (ii) the specific humidity is equal to α YP Q s , and (iii) wind speed is zero.

Complex terrain
Complex terrain significantly complicates the representation of the surface heat flux with respect to the ideal flat terrain case.This section shows how the effects of complex terrain are taken into account in the calculation of the components of the surface heat flux.

Shortwave radiation
Incoming shortwave radiation is always partitioned into two components: a direct component that comes from the direction of the sun, and a diffuse component assumed to be isotropic.Since incoming shortwave radiation is often measured as global (i.e.sum of direct and diffuse components), it becomes important to differentiate in its direct and diffuse portions since the two components react differently to complex terrain.Erbs et al. (1982) provided an empirical expression relating k t , the ratio of the hourly diffuse radiation to the hourly global radiation, to the ratio of the hourly global radiation to the hourly radiation at the top of the atmosphere (namely τ a • τ c ).So far all the radiation components are relative to flat surfaces.In complex terrain: (i) the direct component is obtained by multiplying the direct component for the flat surface by the ratio cos(θ n ) cos(θ v ) , where θ n is the angle between the normal to the surface and the direction of the sun, and θ v the angle between the vertical and the direction of the sun (solar angle deviation); (ii) the direct component may be shaded by the surrounding topography (cast shadow) or it can happen that the angle θ n be larger than 90 degrees (self shadow); (iii) there is also a component of incoming shortwave resulting from reflections from the surrounding terrain, which are assumed to be isotropic in the terrain view angle; and (iv) the diffuse radiation coming from the sky, assumed isotropic, has to be reduced according to the visible sky angle.
The incoming shortwave radiation is calculated as follows: -Clear-sky global radiation on a flat terrain surface is obtained from Eq. ( 20).
-Cloud transmittance is obtained either from Eq. (D5) or as the ratio of measured global radiation to clear-sky global radiation (both on a flat terrain surface).
-The diffuse and direct portions of incoming shortwave radiation on a flat ground surface are calculated according to the formula of Erbs et al. (1982).The diffuse radiation obtained is referred to as hemispheric diffuse radiation, because it considers the sky unobstructed.
-Direct radiation is corrected according to topography, accounting for shadowing and solar incidence angle.
-Diffuse radiation is calculated as two components: (i) one coming from the atmosphere calculated multiplying the hemispheric diffuse radiation by the sky view factor (V f ), which is a topographical parameter that accounts for the portion of the sky that is actually seen from a pixel and varies from 0 (sky not visible) to 1 (flat terrain case, where the sky view angle is the entire hemisphere); and (ii) a second component resulting from the shortwave radiation reflected from the terrain seen from the pixel (surrounding terrain).This component can be calculated with complex algorithms that account for all pixels visible from the pixel of interest p (Helbig et al., 2009).This is not performed in GEOtop, which, instead, calculates this radiation component in the following way: either it is considered that the terrain surrounding a generic pixel p has the same outgoing shortwave radiation as p (that is SW out,p ), and, so, the radiation from the surrounding terrain is given by (1 − V f ) • SW out,p (one-dimensional approximation), or the average of outgoing shortwave radiation in a certain area is taken (SW out,av ) and the radiation from the surrounding terrain is considered as

Longwave radiation
As for diffuse shortwave radiation, in complex terrain incoming longwave radiation comes from both the atmosphere and the surrounding terrain.The former is given by the incoming longwave radiation calculated as in the flat case multiplied by V f .The component emitted from the surrounding terrain is calculated with the same two methods shown for shortwave radiation.

Turbulent fluxes
The turbulent exchange in complex terrain has been observed to significantly deviate from the Monin-Obukhov similarity theory, which is built from the premise that terrain is flat and infinitely homogeneous (e.g.De la Casiniere, 1974).In particular, a maximum wind speed is often observed near the surface as a result of wind gravity flows, whereas according to the theory the wind profile should be logarithmic (in neutral atmosphere) with a small deviation due to temperature gradients (Halberstam and Schieldge, 1981;Meesters et al., 1997;Wagnon et al., 1999).Including these effects in a model like GEOtop would require solving the Navier-Stokes equations for the wind field, which is beyond the purposes of the model.However, it has been also observed (e.g.Denby and Greuell, 2000) that if the measurements of wind, temperature and relative humidity are performed as close as possible to the surface, the conditions are actually closer to the assumptions of the Monin-Obukhov similarity theory.This justifies the application of the theory in GEOtop also for complex terrain.Since meteorological variables are measured only at a limited number of locations, a statistical distribution method is required to assign meteorological forcings (wind, temperature and relative humidity) for the surface energy balance to the whole surface.This issue is described in Appendix C.

Snow cover
The snow cover plays an important role as it buffers the energy and mass exchanges between the atmosphere and soil.Important processes related to the snow cover dynamics include snow warming and cooling, melting and refreezing, water percolation, accumulation due to snow precipitation, avalanches, deposition of wind blown snow, erosion due to wind and densification due to snow metamorphism.
The system of equations for snow is similar to the set of Eqs. ( 1) and ( 7) used for the soil matrix, as snow may also be considered a porous material.However, snow has the following peculiarities: (a) the snow volume of control is ephemeral, i.e. it may disappear as a result of melting; (b) the rigidity to the structure is given by the ice grains; (c) the porosity of snow φ s is variable and depends on the ice volumetric content θ is ; (d) the control volume is not fixed, but is subject to variations due to accumulation, compaction and melting processes; (e) the capillarity effects are in general not significant (Jordan, 1991).
In GEOtop snow is computed solving in sequence: (i) the heat equation, (ii) snow metamorphism, (iii) water percolation, (iv) accumulation.The simplified effects of avalanches and blowing snow are also considered in a simplified way, respectively following Gruber (2007) and Pomeroy et al. (1993), but they will not be dealt with in this paper.

Snow volumetric system
Let us call V s the control volume (m 3 ) of snow: the sum of the liquid, solid and gaseous contents in the volume must equal 1 for continuity: where θ ws (-) is the liquid water volumetric fraction and θ as (-) the volumetric fraction of air and gaseous components.
Considering that the rigidity to the structure is given by the ice grains, snow porosity may be calculated as the available pores excluding the ice grains, i.e.
Snow porosity is not constant as it depends on the snow metamorphism.

Heat equation
Following the approach used for the soil, the heat equation for snow becomes similar to Eq. (B7), the internal energy of the snow being where T s ( • C) is the snow temperature.The heat equation is solved neglecting the lateral gradients with the same numerical method used to solve the heat equation in the soil.The boundary condition at the interface with the atmosphere is given by the surface heat flux, as described in the previous section.At the interface with the soil surface the heat exchange is given by the conduction heat flux, which is dependent on the temperature gradient in proximity to the interface.This last exchange flux actually couples the heat equations in the soil and snow, which have to be solved together in a system.
A freezing characteristic curve is defined also for snow in order to derive an expression for the apparent heat capacity.However, the phase change in the snow takes place virtually at 0 • C, since the pores are large enough that temperature depressions due to capillarity effects are not significant (Jordan, 1991).The definition of a freezing characteristic curve has mostly a numerical reason, therefore the curve must approach as much as possible the unit step function or Heaviside function, which has value zero when temperature is negative (Celsius) and value one when temperature is positive, but must preserve continuity.A simple relation relating temperature and the ratio between liquid water content and total water content in the snow is used (Jordan, 1991): where a ( • C −1 ) is a constant.The higher value this constant is set to, the closer to a step function the curve is, but, at the same time, the more difficult the numerical resolution is.Jordan et al. (1999) set the constant to 10 2 .However, values up to 10 5 can be assigned.
The heat flux at the soil-snow interface is calculated defining an effective thermal conductivity at the interface, and the temperature gradient calculated with the temperatures of the lowest snow layer and the top soil layer.The effective thermal conductivity is sought in a similar way as the effective thermal conductivities at the interfaces between soil layers with Eq. ( 5) proposed by Cosenza et al. (2003).

Metamorphism
GEOtop describes the densification that the newly fallen snow rapidly undergoes (destructive metamorphism) as well as the slow compaction process as a result of the snow weight (overburden), using the empirical formulae of Anderson (1976), improved by Jordan (1991) and Jordan et al. (1999).The constructive metamorphism leading to new shapes of the snow crystals, like hoar layers, is not represented.
The equation describing snow densification is applied separately for each snow layer and is written as (Anderson, 1976) 1 where D is the thickness of the snow layer, and C 1 and C 2 are respectively the total fractional compaction rate (s −1 ) as a result of destructive metamorphism and overburden.Equation (30) is integrated with the same time step t used for the heat equation, assuming that the snow layer has a thickness D 0 (known) at the beginning of the time step and a thickness D 1 at the end.The integration is carried out as follows: which gives The formulations of the compaction rates are reported in Appendix D2.

Water percolation
The equation governing the water balance in each snow layer is the second part of Eq. ( 7), which, considering a 1-D discretisation and integrating along the vertical direction, results in where J up ws (m s −1 ) is the incoming flux from above, J dw ws the outgoing downward flux and S ws (m s −1 ) the liquid water sink term integrated in the snow layer.
Following Colbeck (1972), the flux J ws occurs as soon as θ ws reaches a certain threshold value accounting for the capillary retention.This threshold is set as a snow porosity fraction, namely S r • φ s , where S r is referred to as irreducible water saturation and normally ranges between 4 and 7 % (Colbeck, 1972).The flux J ws is calculated according to Darcy's law but neglecting the capillarity effects and using only the gravimetric gradient: where δ ( • ) is the local slope angle and K s (m s −1 ) is the hydraulic conductivity for snow.K s is calculated according to the model of Brooks and Corey (1964), as in Jordan (1991): where K s,max (m s −1 ) the maximum hydraulic conductivity for the snow, and S e (-) is the effective saturation, which is given by The maximum hydraulic conductivity may significantly vary with snow properties and aging; however, a constant value of 5 × 10 −3 m s −1 is used, which is typical for an isothermal snowpack according to Shimizu (1970).
The water flux in an isothermal snowpack is usually very fast, as a result of the high porosity of snow.On the other hand, in a non-isothermal snowpack water may percolate into a cold snow layer and be there refrozen.This has the net effect to slow down the percolation process.Therefore, the hydraulic and thermal control on water percolation have significantly different timescales.The water flow is then calculated in an uncoupled way in the following steps: (i) the heat Eq. ( 1) for the snow cover is solved assuming that the liquid water does not move, obtaining a "static" solution; (ii) the incoming water flow from above is added (also rain for the upper layer); (iii) the internal energy content of the inflowing water (in terms of latent heat) is added to the internal energy of the snow, assuming that there is an instantaneous energy mixing that may lead to the partial or complete refreezing of the liquid water; (iv) finally the outgoing downward water flux is calculated with the new state variables.

Accumulation
In most cases, precipitation data are given only as total precipitation.The most common method for splitting this into rain and snow (e.g.US Army Corps of Engineers, 1956;Auer, 1974) is to define two air temperature thresholds: a higher value above which precipitation is only rain, and a lower value below which precipitation is only snow (Kienzle, 2008).Garen and Marks (2005) proposed to use dew temperature instead, as the temperature interval in which both rain and snow precipitation are coexistent is in this case much smaller, so that just one threshold value can be defined.In GEOtop both methods are available where the threshold temperatures are set as parameters (for air temperature normally around −1 and 3 • C) and in between a linear interpolation is performed.
Fresh snow density depends on grain size and crystal type, as they affect the way fresh snow is deposited.Smaller grains with simpler shape pack more efficiently and lead to denser snow.These effects, however, are indirectly described parameterising the density as a function of air temperature and wind, which are easier to measure and have been correlated to fresh snow density in several studies (e.g.McGurk et al., 1988).The following formula proposed by Jordan et al. (1999) is used, which incorporates both temperature and wind effects: where ρ ns is the density of the new snow (kg m −3 ).

Discretisation
The snowpack, according to the thermal gradients, may be roughly classified into three regions: an upper, middle and bottom portion.In the upper and bottom regions the vertical gradients are often high, as a result of the interactions with the atmosphere and the underlying soil, respectively.On the other hand, in the middle region the vertical gradients are weaker.The snow discretisation in GEOtop is done in order to accurately describe the thermal gradients in the snowpack and avoid the allocation of unnecessary memory.The total number of layers, in fact, depends on the mass of snow present, whereas the distribution of layers in the snowpack privileges the upper and bottom zones.
The details of the snow discretisation scheme are reported in Appendix D3.

Testing GEOtop
A full model evaluation and validation is not performed here.Since the number of possible test cases is nearly infinite, a complete validation for a complex simulator like GEOtop can never be claimed.The actual testing is performed in each application study where the simulator is used.Some studies have been already published (e.g.Bertoldi et al., 2010;Della Chiesa et al., 2014;Endrizzi and Marsh, 2010;Endrizzi et al., 2011;Gubler et al., 2013;Fiddes et al., 2013), others are ongoing or planned for the near future.Table 2 shows a list of the major model components or modelled processes that have been tested or require further evaluation.We demonstrate here that GEOtop produces plausible results in its major components in order to add credibility to the results of the subsequent experiment.

Time-lagged vs. coupling solution
The error due to linking the heat and water flow equations in a time-lagged manner instead of full numerical coupling is here evaluated experimentally.For this purpose, we initialise a 0.3 m deep column of sandy loam with a uniform pressure head of −0.1 m and a linear temperature profile between 0 • C at the top and −1 • C at the bottom.No energy or water is exchanged with the outside.Running this simulation over time, the temperature of the soil column will approach a uniform value and pressure gradients due to freezing will redistribute water.Finally, to estimate the magnitude of error inherent in using a linked solution, we compare results based on time steps of 1 s and 1 h.The temperature and total water content (i.e.liquid and frozen) as well as deviations between both simulations are shown in Fig. 1.As the deviations in temperature are negligible and those in liquid water are of the order of ±6 %, the linked solution is deemed acceptable in this context.

Frozen soil scheme
The performance of the frozen soil model has been evaluated by Dall'Amico et al. (2011a) by comparison with the analytical solution to the Stefan problem (Lunardini, 1981) and with the experimental results of Hansson et al. (2004).

Snow model
The performance of GEOtop with respect to snow has been evaluated using the data published by Morin et al. (2012) at Col de Porte, a mountain pass in the French Alps at 1326 m a.s.l.near Grenoble.The model was driven by hourly measured near-surface meteorological data: air temperature, relative humidity, downwelling shortwave and longwave radiation, wind speed, air pressure and precipitation provided by the Meteo France weather station located at the pass.The simulation was performed with standard parameters without improving the results by trial and error or fitting.Similar to the evaluation of the snow simulator CROCUS (Brun et al., 1992) by Vionnet et al. (2012), we quantify performance based on daily averages of snow depth and water equivalent using the bias and the root mean squared deviation (RMSD) for the months December to May during the years 2001 to 2011. Figure 2  GEOtop.inpts.ColDePorte).For snow water equivalent, we obtain a RMSD of 37.1 mm (39.7, 37.0 mm) and a bias of −3.2 mm (−17.3,−2.3 mm).For snow depth, the RMSD is 0.15 m (0.11, 0.13 m) and the bias is 0.07 m (−0.01, 0.08 m).
For comparison, the range of values reported for two different versions of Crocus (Vionnet et al., 2012) are given in brackets.Based on visual inspection, soil temperatures below the snow cover are represented reasonably well.

Complex topography and temperature
In order to demonstrate the performance with respect to topography, we used near-surface temperatures measured in steep bedrock that accumulates nearly no snow.The sites chosen for this demonstration are two contrasting sites: "Jungfrau ridge south" and "Jungfrau ridge north" (PER-MOS, 2009), positioned at a horizontal distance of about 20 m -the first is sun-exposed and the latter is shaded, resulting in a mean annual temperature difference of nearly 8 • C. Both sites are located in proximity to Jungfraujoch, a mountain pass in the Swiss Bernese Alps at 3470 m a.s.l.. Horizon shading and sky view factor are parameterised based on fish-eye photography (Gruber et al., 2003) Endrizzi (2009), Endrizzi and Marsh (2010) and here Distributed snow cover modelling Endrizzi (2009), Buri (2013) and ongoing Soil energy and water balance in one-dimensional simulations Gubler et al. (2013), Fiddes et al. (2013) and here Soil energy balance in complex terrain Bertoldi et al. (2010) and Kunstmann et al. (2013) Interactions between soil andvegetation Della Chiesa et al. (2014) Interactions between snow cover and vegetation Endrizzi and Marsh (2010) Interactions between soil freezing and lateral water flow Endrizzi et al. (2011), Endrizzi and Gruber (2012) and here Water balance in complex terrain Bertoldi et al. (2014) and Brenner (2014) Runoff production planned distance of about 1.2 km to the rock temperature measurements and is about 150 m lower in elevation.Air temperature measurements have been extrapolated using the standard lapse rate of 0.0065 The results are reported in Fig. 3 and the parameter file used is available in the Supplement (GEOtop.inpts.JungfrauJoch).
Based on the comparison with around 50 000 hourly measurements at each location (south/north), we obtain a bias of −0.63/ − 0.68 • C and a RMSD of 4.37/2.01• C.

Sensitivity study
The work of Gubler et al. (2013) demonstrates the robustness of GEOtop as it is built on more than two million highly diverse simulations that converged.This study further revealed that the sensitivity of ground temperatures to changes in temporal of spatial discretisation are moderate, and that equilibria independent of initial conditions can be reached reliably.

Simulation experiment
In order to demonstrate the relevance of the modelling approach, GEOtop has been run in a catchment made up of two hillslopes forming a convergent topography (Fig. 4).Different simulations have been set up that differ in (i) topography, and (ii) model configuration with respect to the water balance.It is considered to have a gravel soil with a saturated hydraulic conductivity of 0.002 m s −1 , thermal conductivity (of the soil matrix) of 2.5 W m −1 K −1 , θ r = 0.057, θ sp = 0.487 and the following Van Genuchten parameters: α = 2.0 m −1 and n = 1.8.The soil has been discretised with 80 layers: the first 4 layers starting from the surface have thicknesses ranging from 0.01 to 0.08 m, in consideration of the high vertical gradients of temperature and water pressure, the lowest 15 layers have thicknesses ranging from 0.2 to 0.5 m whereas the remaining layers have thickness of 0.1 m.Overall, the soil domain reaches a depth of 10.5 m.The surface has been represented with square pixels of a dimension of 20 m.
Topographical differences have been created by varying the inclination angle of the lateral slopes (β in Fig. 4) from 5 • to 20 • .The longitudinal slope (from point 3 to 4 in Fig. 4) of 5 • and the average elevation (3000 m) are kept constant.The two topographies are hereafter referred to as 5 • and 20 • topography, respectively.
The following water balance configurations, which are here referred to for simplicity as 3-D, 1-D, and 0-D, are considered.3-D means that the full three-dimensional variably saturated Richards equation with the surface flow is used.In 1-D, the Richards equation is solved only in one dimension, i.e. in the vertical, and no lateral subsurface drainage is considered.Surface flow occurs as lateral flow only on the surface.In 0-D, the soil water balance is not solved, and no infiltration is described.Therefore, the total water content always remains at its initial value, but soil water will undergo freezing and thawing.
Six simulations have been then performed: 5 and 20 • 0-D.They have been run for the hydrological year 2001-2002 (from 1 October 2001 to 30 September 2002) using meteorological data measured at the station of Davos, Switzerland, located at 1595 m a.s.l. and operated by Meteoswiss.Air temperature was extrapolated using the standard lapse rate of 0.0065 K m −1 .This year was chosen because it is the most similar to the average year in the 1981-2010 period, if the cumulated winter precipitation (from October to May) and the average air temperature in the summer (from June to September) are considered.The first quantity represents a proxy for snow precipitation and the second one approximates summer warming.Respective values are 452 mm and 10.6 • C for the hydrological year 2001-2002, and 479 mm and 10.6 • C for the average year in the period 1981-2010 period.
The system was initialised considering an absence of snow cover, uniform soil temperature of −1 • C, and soil deeper than 1 m saturated, and hydrostatic pressure profile in both the saturated and unsaturated portions (prolonging the hydrostatic pressure profile also for negative pressures).However, in order to reduce the influence of the arbitrariness of the initial condition, a spin-up simulation has been performed for 100 years.At the end of the spin-up the catchment is still completely underlain by permafrost, with active layer depths ranging from 1 to 7 m.Since the response time of permafrost soil is extremely slow, a longer spin-up simulation would probably enhance the spatial variability of the active layer depth and completely thaw permafrost in some part of the catchments.However, a 100-year spin-up is considered acceptable in this context, since we are mainly interested in the approach evaluation.
The set of simulations are evaluated with the following steps: (i) annual average variables in the single points shown in Fig. 4 are compared, (ii) annual average distributed variables are compared and their spatial distribution in the catchment discussed, and (iii) the time evolution of some variables for a selected point is shown.Figure 5 shows the mean annual temperature at a 4 m depth, the active layer depth (considered as the lowest depth of thaw reached during the summer), and the annual average water table depth (the average was calculated only when there is actually a water table) for the six points for the six simulations.Points 1 and 2 are both on the south slope at short distance apart and have very similar temperature values, but the different topographies and water balance descriptions give significantly different values of temperatures, which range between +0.3 • C (permafrost thawed) and slightly negative values (permafrost still present at a 4 m depth).The active layer and water table depth show significant differences in the simulations (of the order of 1-2 m), with slightly deeper values for Point 2, which is downstream.Point 3 and 4 are located in the channel portion and also exhibit similar values.The results of the six simulations give temperature ranges between slightly positive and slightly negative values and differences of the order of 1-2 m in active layer and water table depths.Points 5 and 6 are located on the north slope and are significantly colder than the other points, with temperatures ranging from −1 • C and slightly negative values, and active layer and water table depths ranging from 0.3 to 1 m.It is therefore important to notice that both different topographies and different hypotheses on soil water balance have a significant effect on the results.
Figure 6 shows the spatial distribution of the active layer depth and the depth of thaw.If the water balance is not solved (0-D simulation), the spatial variability of the active layer is given by aspect and elevation only.In the north slope the values are rather homogeneous around 1 m, while in the south slope they are significantly dependent on altitude, ranging from 1 m at the top to 3 m at the bottom for the 5 • topography, and from 4 to 7 m for the 20 • topography.If the water balance is solved 1-D, the effect of altitude and aspect is attenuated, and the spatial variability reduced.If a full 3-D water balance is considered, a clear relation between active layer depth and water table can be recognised.Where water table is shallower and, therefore, soil is wetter, the active layer is deeper.In the  given by the six simulations that have been performed.5 • topography, the water table at the lowest elevations of the south slope is at the surface.For mid elevations water table depth reaches a maximum value of about 2 m, and then for higher elevations it decreases to values of about 1-1.5 m.This is probably a consequence of the strong dependance of active layer on elevation.However, the active layer depth has a larger spatial variability than in the 0-D water balance case, since it reaches values up to 5 m at the bottom of the south slope (instead of 3 m) and up to 2.5 m at the bottom of the north slope (instead of 1 m).This is probably a result of the interplay between soil moisture and freezing-thawing energetics, which is added to the spatial variability induced by slope and aspect.The 20 • topography has similar features, but the lateral drainage in this steeper topography adds active layer spatial variability in a lesser degree than in the 5 • topography.Figure 7 shows the temporal evolution of temperature, and liquid, solid and total water content in soil and snow for point 1 shown in Fig. 4. All charts also report the lower and upper boundaries of the thawed layer, and the water table.Temperature is in general negative and close to 0 • C except at the surface.The soil is completely frozen from mid-December to mid-June, and the thawed portion reaches the maximum depth in mid-September, which is kept until almost mid-November, when the freezing front from above, starting in October, gets significantly low.The charts of liquid water, ice and total water content show that three regions, from top to bottom, can roughly be distinguished: (i) a region relatively rich in total water at the top (hereafter referred to as "wet region"), (ii) a drier region in the lower portion of the active layer ("dry region"), and (iii) the "undisturbed frozen region" below.The wet region starts to develop in the autumn when the active layer starts freezing and is initially bordered below by the freezing front.It gets then wetter as a result of early snow melting and, probably, water transfer from the slightly warmer wet region below.Liquid water content is low (< 0.1) and ice content is high (> 0.3) when the wet region is frozen, and temperature is significantly low (around −2 • C).When the wet region thaws in the summer it gets even wetter as a result of infiltration.However, thawing takes place at a relatively slow rate in this region, as a result of the high ice content and, therefore, the high energy content required to melt it.The dry region has slowly lost the initial water content over many years due to the lateral drainage during the summer.When the thawing front reaches this region in late summer, the thawed soil depth sharply increases because the energy required to thaw the soil is relatively low.At the same time a large amount of liquid water is drained from the wet region above, which virtually disappears.This amount of liquid water is gradually drained horizontally, and it eventually accumulates in a thin layer in the lower part of the region.This layer gets thinner as the season progresses, but it does not immediately disappear when the region freezes in early winter, as a result of relatively high (though negative) temperature.The undisturbed region below was never affected by summer thawing, and the total water content is still dependent on the initial condition.Even though a spinup simulation has been performed repeating for 100 years the meteorological data corresponding to the hydrological year 2001-2002, it is not guaranteed that an equilibrium state has been reached, and, therefore, a longer spin-up simulation could entail different conditions.

Conclusions
GEOtop 2.0 describes the energy balance in the soil and snow taking into account the interactions with the atmosphere and solves a fully three-dimensional form of the water flow equation.It uses a simplified, but physically consistent parameterisation of the soil water retention curve to describe soil freezing characteristics in saturated and unsaturated conditions.In its numerical implementation, GEOtop 2.0 uses sophisticated integration methods, which allow convergence even in cases where parameters have nearly discontinuous behaviour, and results in proper conservation of mass and energy.This allows the investigation of complex hydrologic phenomena in cold regions for which no compound parameterisation may exist.
The model has shown a consistent physical realism, and the comparison with snow cover data has resulted in a reasonable agreement.An experiment with varying model structure, carried out by differing the treatment of water transport in the soil, has highlighted that significant differences of temperature, water fluxes and water table depths can result from this.Furthermore, the spatial differentiation of these results in response to topography highlights the complex nature of the phenomena investigated and represented.
GEOtop 2.0 represents a wider range of processes in the water and the energy budgets at fine scales than most other simulators.It thus allows studying their interactions without introducing "ad hoc" solutions that may compromise the representation of complexity.In this paper, we have described the details and functioning of fine-scale hydrology and its interaction with frozen soil and snow cover.However, the model could be applied for a much wider range of environments and scientific issues.
where all the symbols are described in Table D3.The variation of the water content θ w and internal energy U may be divided into the component due to "phase change" (superscript "ph") and the component due to water "flow" (superscript "fl"): According to this assumption, Eq. (B1), after some rearrangements, becomes Equalizing both members of Eq. (B5) to a common value, say zero, eventually one obtains an equivalent system for the water balance equation: The components "fl" and "ph" of the internal energy may be derived starting from the definition of the internal energy U (Dall' Amico et al., 2011a): where C := ρ sp c sp (1−θ sp )+ρ i c i θ i +ρ w c w θ w is the volumetric heat capacity (J m −3 K −1 ).Differentiating Eq. (B7) one obtains where From the first equation in (B6) the variation of the ice content may be related to the variation of θ ph w : Substituting Eqs.(B3) and (B10) into Eq.(B8), after some calculations, it is obtained that dU where one can define The flux J is the heat advected by flowing water and equals Substituting Eqs.(B4) and (B13) into Eq.(B2) and considering that ∇ From the second equation of Eq. (B6), Eventually Eq. (B15) becomes

Appendix C: Distribution of meteorological data
In the GEOtop code a complete set of routines is finalised at the spatial interpolations of meteorological variables.Air temperature is distributed according to Liston and Elder (2006a).All the measurements at different elevations are converted into values corresponding to a unique reference elevation according to a spatially constant lapse rate, which can however vary in time.The obtained values are spatially interpolated with the geostatistical method of Barnes (1964).
The elevation correction, given by the lapse rate multiplied by the difference between actual elevation and reference elevation, is then applied on the interpolated temperature field, which is related to the reference elevation.Relative humidity is converted into dew temperature, and the same interpolation method used for air temperature is used with a dew temperature lapse rate, which is normally smaller than the air temperature lapse rate (Liston and Elder, 2006a).Precipitation is also distributed with the same method.An adjustment factor considering the dependence of precipitation on the elevation, in analogy to lapse rate for air temperature, is also included and calculated according to Thornton et al. (1997).
Wind speed is an important factor affecting the turbulent fluxes of sensible and latent heat.In order to describe the effect of topography on the surface energy balance it is important to consider in the model a topographically dependent wind field.A full resolution of the fluid dynamic equations would be too computationally heavy for GEOtop.The wind field is instead parameterised using topography (Liston and Elder, 2006a).In particular, the parameterisation is implemented correcting the wind speed with factors depending on slope and curvature of the surface: where w s is the wind speed modulus, w s0 the wind speed modulus ideally unaffected by topography, S the slope of the curve given by the intersection of the surface and by a vertical plane oriented in the direction of the wind, C v the curvature of the same curve, and χ s and χ c are calibration parameters.
Wind direction is corrected according to Ryan (1977) in order to represent wind skirting round topographic obstacles.
Recently GEOtop has been also enabled to exploit Me-teoIO (Bavay and Egger, 2014), a library developed by the Snow and Avalanche Research Institute of Davos (Switzerland) aimed at caching, filtering, resampling and spatially interpolating meteorological variables.
increased by 50 %.For densities larger than the cutoff density, α s decreases very rapidly, becoming one-tenth at a density of 150 kg m −3 and one-hundredth at 200 kg m −3 .Furthermore, the compaction decreases with snow temperature, and at a temperature of about −17 • C it is two times as much as the value at a temperature of 0 • C. The parameter C 2 is related to the ratio between the weight of the overlying snow column P s (N m −2 ) calculated at the centre of the considered snow layer and the snow viscosity η (N s m −2 ): where with ρ sn := ρ i θ is + ρ w θ ws (kg m −3 ) the snow density.Densification takes place also when melting occurs.When ice melts in a specific snow layer, it is considered that the thickness of the layer reduces proportionally to the ice content.This makes sense, as the snow depth is a result of the ice grain structure.This leads to an increase in density, since the same total water content will occupy a smaller volume.Densification occurs also when liquid water starts to refreeze and the new ice will fill the empty pores.On the other hand, a snow layer may be subject to a density decrease as a result of the percolation process, because the total water volume in the layer will decrease, but not its volume.All these processes are also included in the model.
Another snow densification process is related to wind loads.When the wind speed is higher than a threshold (set to the wind value at which snow starts to be drifted), the wind load is considered as an additional overburden.This describes the snow packing at the surface, which leads to a progressive resistance to being drifted (Liston and Elder, 2006b).

D3 Snow discretisation
The mass of snow per unit area (kg m −2 ) present in a snow layer l is referred to as M l , and the total mass of snow per unit area present in the snowpack is named M tot .GEOtop requires four parameters to set the snow layering scheme: -M * up and M * dw (kg m −2 ), which are the maximum mass per unit area of the snowpack in the upper and bottom regions, respectively; At each time step the layers are re-organised (in number, thickness, mass content and internal energy) according to the evolution of the snowpack.In particular, as outlined in Table D2, three processes may occur: -Layer splitting: if the mass of the top snow layer as a result of new snow accumulation exceeds the maximum allowed mass (M l > M * l ), then it is split into two layers in such a way that the lower new layer keeps a mass equivalent to M * l and the new surface layer has the remainder of the mass.In the case that the total mass of the upper region exceeds the threshold M * up , then the lower layer is pushed to the middle or lower region.
-Layers merging: two adjacent layers may be merged into one layer that will have the sum of ice and liquid water content of the layers prior to merger, and temperature resulting from the energy content given by the sum of the energy contents of the layers prior to merger.Layer merging happens in the following cases: (i) if the number of layers of the middle zone N mid exceed N * mid , then two adjacent layers of the middle zone are merged.The choice falls on the two adjacent layers that have the smallest combined mass: this allows keeping the layers in the middle region of similar snow mass content and prevents from excessively smoothing the snow vertical profile; (ii) if a snow layer completely loses its ice mass in a one-time-step integration of the heat equation, the snow layer that would disappear is merged with the underlying layer and then the heat equation is reintegrated.
et al.: Simulating the combined energy and water balance with soil freezing, snow and terrain

Figure 1 .
Figure 1.Evolution of temperature (a) and total water content (c) for simulations with a time step of 1 s.The difference to the solutions with a time step of 1 h are also shown for temperature (b) and total water content (d) based on subtracting the hourly solution from the 1 s solution.

Figure 3 .
Figure 3. Temperature in complex topography: comparison of measured (red line) and simulated (black line) results at Jungfrau ridge south (left) and north (right) during year 2005/2006.Measured data are taken from PERMOS (2009).

Figure 4 .
Figure 4. Synthetic catchment and location of points analysed.The converging topography has a sun-exposed and a more shaded side.It is varied with respect to the inclination angle of its hillslopes β.Channel inclination, in the direction from point 3 to 4, is always 5 • .

Figure 5 .
Figure5.Mean annual ground temperature at a depth of 4 m, active layer depth, and water table depth for the six points shown in Fig.4given by the six simulations that have been performed.

Figure 6 .
Figure 6.Distributed results: active layer depth for the simulations 5 • 0-D (1), 20 • 0-D (2), 5 • 1-D (3), 20 • 1-D (4), active layer depth (5) and mean annual water table depth (in metres) (6) for the simulation 5 • 3-D, and active layer depth (7) and mean annual water table depth (8) for the simulation 20 • 3-D.All depths and the scale are in metres, and the elevations in metres above sea level.Consistently with Fig.4, north is towards the left, and the left and right slopes have respectively a south and north aspect.

Figure 7 .
Figure 7. Time evolution of temperature, and liquid, solid, and total water content for soil and snow at Point 1, as shown in Fig. 4, for the simulation for the 5 • topography 3-D.The blue line indicates the water table, while the lower and upper borders of the thawed layer are shown in red.Negative depths correspond to soil, and positive depths to the snow cover.

Table 2 .
List of the major model components or modelled processes, and corresponding publication where they are evaluated, even partially, if available.

Table D1 .
Definition of snow discretisation regions.

Table D2 .
Available processes in a layer and triggering condition.As described inTable D1, the number of regions used to describe the snowpack depends on M tot : (i) if M tot < M * up , only the upper region is used, which extends throughout the whole snowpack; (ii) if M tot is larger than M * up , but smaller than M * up + M * dw , also the bottom region is created with mass M tot − M * up ; (iii) if M tot > M * up + M * dw , then also the middle region is defined.

Table D3 .
Table of symbols for soil balance and surface and subsurface water flow equations.

2857, 2014 2852 S. Endrizzi et al.: Simulating the combined energy and water balance with soil freezing, snow and terrainTable D4 .
Table of symbols for the atmospheric variables.

Table D5 .
Table of symbols for the snow.