Simple Process-Led Algorithms for Simulating Habitats (SPLASH v.1.0): Robust Indices of Radiation, Evapotranspiration and Plant-Available Moisture

Bioclimatic indices for use in studies of ecosystem function, species distribution, and vegetation dynamics under changing climate scenarios depend on estimates of surface fluxes and other quantities, such as radiation, evapotranspiration and soil moisture, for which direct observations are sparse. These quantities can be derived indirectly from meteorological variables, such as near-surface air temperature, precipitation and cloudiness. Here we present a consolidated set of Simple Process-Led Algorithms for Simulating Habitats (SPLASH) allowing robust approximations of key quantities at ecologically relevant time 5 scales. We specify equations, derivations, simplifications and assumptions for the estimation of daily and monthly quantities of top-of-the-atmosphere solar radiation, net surface radiation, photosynthetic photon flux density, evapotranspiration (potential, equilibrium and actual), condensation, soil moisture, and runoff, based on analysis of their relationship to fundamental climatic drivers. The climatic drivers include a minimum of three meteorological inputs: precipitation, air temperature, and fraction of bright sunshine hours. Indices, such as the moisture index, the climatic water deficit, and the Priestley-Taylor coefficient, are 10 also defined. The SPLASH code is transcribed in C++, FORTRAN, Python, and R. One year of results are presented at the local and global scales to exemplify the spatiotemporal patterns of daily and monthly model outputs along with comparisons to other model results.


Introduction
Despite the existence of dense networks of meteorological monitoring stations around the world, plant ecophysiology and biogeography suffer from a lack of globally distributed observational data, especially those central to the estimation of ecosystem-level photosynthesis, including photosynthetic photon flux density and soil moisture.To overcome this deficiency, we present simple process-led algorithms for simulating habitats (SPLASH) for generating driving datasets for ecological and land-surface models (e.g., monthly carbon and water fluxes or seasonal plant functional trait distributions) from more readily available meteorological observations.
SPLASH is a continuation of the STASH (static shell) model, which was originally developed for modeling the climatic controls on plant species distributions at a regional scale (Sykes andPrentice, 1995, 1996;Sykes et al., 1996).The intention of STASH was to provide bioclimatic indices, reflecting the environment experienced by plants more closely than either standard summary variables such as mean annual temperature, or such constructions as "mean precipitation of the warmest quarter", while requiring only standard meteorological data as input.A key component of STASH was a simple, physically based soil-moisture accounting scheme, first developed by Cramer and Prentice (1988), which has been used inter alia in the original, highly cited BIOME model (Prentice et al., 1992), the general forest succession model (FORSKA) described by Prentice et al. (1993), and the Simple Diagnostic Biosphere Model (Knorr and Heimann, 1995).Despite the subsequent development of more complex dynamic global vegetation models (Cramer et al., 2001;Sitch et al., 2003;Woodward and Lomas, 2004;Quillet et al., 2010;Prentice and Cowling, 2013;Fisher et al., 2014) and land-surface models, the relatively simple algorithms in STASH continue to have many applications, including to new areas such as the distribution of plant functional traits (Harrison et al., 2010;Meng et al., 2015), assessment of climate-change impacts on specific biomes (Gallego-Sala and Prentice, 2012), large-scale water resource assessments (e.g., Ukkola et al., 2015), and simple first-principles modeling of primary production (Wang et al., 2014).The continuing utility of these algorithms owes much to their robustness, which in turn depends on the implicit assumption that vegetation functions predictably -so that, for example, evapotranspiration occurs at a potential rate under well-watered conditions, and is reduced as soil water is drawn down.STASH is thus unsuitable to answer questions like the effect of imposed vegetation changes on runoff, or modeling vegetationatmosphere feedbacks.Much more complex models that dynamically couple soil, vegetation, and atmospheric boundary layer processes exist for such applications; however, their complexity brings a burden in terms of lack of robustness and, potentially, large inter-model differences (Prentice et al., 2014).
Despite their long history of use, no single publication documents the algorithms of the STASH model.This work aims to fill that gap to allow for the continued development and use of these algorithms.As the new incarnation of STASH, SPLASH provides the same physically based soil-moisture accounting scheme with updated and corrected analytical expressions for the calculation of daily radiation, evapotranspiration, and soil moisture.Included in this documentation are the equation derivations, variable definitions, and information regarding model assumptions and limitations.One notable improvement is that we have discontinued the approximation of constant angular velocity in the orbit of Earth around the Sun.This version is thus suitable for palaeoclimate applications, whereby orbital precession (as well as changes in obliquity and eccentricity) influences the seasonal distribution of insolation.SPLASH also includes explicit consideration of elevation effects on biophysical quantities.
Key model outputs include daily insolation (incoming solar radiation at the top of the atmosphere) and net surface radiation (H o and H N , respectively); daily photosynthetic photon flux density (Q n ); daily condensation, soil moisture, and runoff (C n , W n , and RO); and daily equilibrium, potential, and actual evapotranspiration (E q n , E p n , and E a n ).Unlike the STASH model, SPLASH explicitly distinguishes potential and equilibrium evapotranspiration, recognizing that under well-watered conditions the excess of the former over the latter is a requirement for foliage to be cooler than the surrounding air, as has long been observed under high environmental temperatures (e.g., Linacre, 1967).
Input values of latitude, (rad), elevation, z (m), mean daily near-surface air temperature, T air ( C), and fractional hours of bright sunshine, S f (unitless), are used for calculating the daily quantities of net radiation and evapotranspiration.Daily precipitation, P n (mm d 1 ), is used for updating daily soil moisture.T air and P n may be derived from various sources, including the freely available daily-averaged air temperature and precipitation reanalysis data from the Water and Global Change (WATCH) program's meteorological forcing dataset (Weedon et al., 2014).Meteorological variables are also available in the Climatic Research Unit (CRU) gridded monthly time series datasets (Harris et al., 2014), which may be downscaled to daily quantities by means of quasi-daily methods (e.g., linear interpolation).Cloud cover fraction, for example the simulated quantities given in the CRU TS3.21 dataset, may be used to approximate S f .Penman's one-complement approximation based on the cloudiness fraction is regarded here as a sufficient estimate of S f (Penman, 1948).The piecewise linear method of Hulme et al. (1995) -an adaptation of the Doorenbos-Pruitt estimation procedure (Doorenbos and Pruitt, 1977) -as used in the development of the CRU cloudiness climatology (New et al., 1999) gives similar results.
We present SPLASH comprehensively re-coded in a modular framework to be readable, understandable, and repro-ducible.To facilitate varied application requirements (including computational speed), four versions of the code (C++, FORTRAN, Python, and R) are available in an online repository (see Sect. 6).The algorithms as presented here focus on application to individual site locations, but a natural extension is towards spatially distributed grid-based datasets.
In line with the intention of the original STASH algorithms, we also present bioclimatic indices at the monthly and annual timescales to exemplify the analytical applications of the SPLASH model outputs.

Methodology
The implementation of the soil-moisture accounting scheme follows the steps outlined by Cramer and Prentice (1988), where daily soil moisture, W n (mm), is calculated based on the previous day's moisture content, W n 1 , incremented by daily precipitation, P n (mm d 1 ), and condensation, C n (mm d 1 ), and reduced by daily actual evapotranspiration, E a n (mm d 1 ), and runoff, RO (mm): where P n is a model input, C n is estimated based on the daily negative net radiation, E a n is the analytical integral of the minimum of the instantaneous evaporative supply and demand rates over a single day, and RO is the amount of soil moisture in excess of the holding capacity.An initial condition of W n is assumed between zero and the maximum soil-moisture capacity, W m (mm), for a given location and is equilibrated over an entire year by successive model iterations (i.e., model spin-up).Under steady-state conditions, the SPLASH model preserves the water balance, such that To solve the simple "bucket model" represented by Eq. (1), the following steps are taken at the daily timescale: calculate the radiation terms, estimate the condensation, estimate the evaporative supply, estimate the evaporative demand, calculate the actual evapotranspiration, and update the daily soil moisture.Daily quantities may be aggregated into monthly and annual totals and used in moisture index calculations.

Top-of-the-atmosphere solar radiation
The calculation of C n and E a n begins with modeling the extraterrestrial solar radiation flux, I o (W m 2 ).The equation for I o may be expressed as the product of three terms (Duffie and Beckman, 2013): where I sc (W m 2 ) is the solar constant, d r (unitless) is the distance factor, and cos ✓ z (unitless) is the inclination factor.
Values for I sc may be found in the literature (e.g., Thekaekara and Drummond, 1971;Willson, 1997;Dewitte et al., 2004;Fröhlich, 2006;Kopp and Lean, 2011); a constant for I sc is given in Table 2.The distance factor, d r , accounts for additional variability in I o that reaches the Earth.This variability is due to the relative change in distance between Earth and the Sun caused by the eccentricity of Earth's elliptical orbit, e (unitless), and is calculated as follows (Berger et al., 1993): where ⌫ (rad) is Earth's true anomaly.True anomaly is the measure of Earth's location around the Sun relative to its position when it is closest to the Sun (perihelion).
The last term, cos ✓ z , attenuates I o to account for the Sun's height above the horizon (measured relative to the zenith angle, ✓ z ), accounting for the off-vertical tilt of Earth's rotational axis, " (i.e., obliquity).The inclination factor is calculated as follows (Duffie and Beckman, 2013): where (rad) is the latitude, (rad) is the declination angle, and h (rad) is the hour angle, measuring the angular displacement of the Sun east or west of solar noon ( ⇡  h  ⇡ ).Declination is the angle between Earth's equator and the Sun at solar noon (h = 0), varying from +" at the June solstice to " at the December solstice; the changing declination is responsible for the change in seasons.For the purposes of ecological modeling, may be assumed constant throughout a single day.See, e.g., Woolf (1968) for the precise geometric equation representing : = arcsin (sin sin ") , where (rad) is Earth's true longitude (i.e., the heliocentric longitude relative to Earth's position at the vernal equinox) and " (rad) is obliquity (i.e., the slowly varying tilt of Earth's axis).Several other methods are widely used for the estimation of for a given day of the year (e.g., Cooper, 1969;Spencer, 1971;Swift Jr., 1976) but are not recommended because they do not account for the change in Earth's orbital velocity with respect to the distance between Earth and the Sun, while Eq. ( 5) does.The relationship between true longitude, , and true anomaly, ⌫, is by the angle of the perihelion with respect to the vernal equinox, e !(rad) (Berger, 1978): While the three orbital parameters (i.e., e, ", and e !) exhibit long-term variability (on the order of tens of thousands of years), they may be treated as constants for a given epoch (e.g., e = 0.0167, " = 23.44 , and e != 283.0for 2000 CE), or they may be calculated using the methods of Berger (1978) or Berger and Loutre (1991) for palaeoclimate studies.Berger (1978) presents a simple algorithm to estimate for a given day of the year (see Appendix A).The daily top-of-the-atmosphere solar radiation, H o (J m 2 ), may be calculated as twice the integral of I o measured between solar noon and the sunset angle, h s , assuming that all angles related to Earth on its orbit are constant over a whole day: where r u = (sin sin ) and r v = (cos cos ), both unitless.
The sunset angle can be calculated as the hour angle when the solar radiation flux reaches the horizon (i.e., when I o = 0) and can found by substituting Eq. ( 4) into Eq.( 2), setting I o equal to zero, and solving for h: To account for the undefined negative fluxes produced by Eq. ( 2) for h h s and h  h s , I o should be set equal to zero during these nighttime hours.To account for the occurrences of polar day (i.e., no sunset) and polar night (i.e., no sunrise), h s should be limited to ⇡ when r u /r v 1 and zero when r u /r v  1.

Net surface radiation
The net surface radiation, H N (J m 2 ), is the integral of the net surface radiation flux received at the land surface, I N (W m 2 ), which is classically defined as the difference between the net incoming shortwave radiation flux, I SW (W m 2 ) and the net outgoing long-wave radiation flux, I LW (W m 2 ): The calculation of I SW is based on the reduction in I o due to atmospheric transmittivity, ⌧ (unitless), and surface shortwave albedo, sw (unitless): A constant value for sw is given in Table 2. Atmospheric transmittivity may be expressed as a function of elevation (to account for attenuation caused by the mass of the atmosphere) and cloudiness (to account for atmospheric turbidity).At higher elevations, there is less atmosphere through which shortwave radiation must travel before reaching the surface.To account for this, Allen (1996) presents an equation based on the regression of Beer's radiation extinction function at elevations below 3000 m with an average sun angle of 45 , which can be expressed as follows:
H N can be decomposed into its net positive, H + N (J m 2 ), and net negative, H N (J m 2 ), components (i.e., H N = H + N + H N ).Assuming I LW is constant throughout the day and making substitutions for I o into Eq.( 10), an expression for H + N may be derived as twice the integral of I N between solar noon (i.e., h = 0) and the net surface radiation flux cross-over hour angle, h n (rad): where Here, h n is the hour angle when I SW equals I LW and can be found by setting I N = 0 in Eq. ( 9) and solving for h, following the same substitutions as used for h s in Eq. ( 8), and may be expressed as follows: To account for the occurrences when the net surface radiation flux does not cross the zero datum, h n should be limited to ⇡ when (I LW r w r u )/(r w r v )  1 (i.e., net surface radiation flux is always positive) and zero when (I LW r w r u )/(r w r v ) 1 (i.e., net surface radiation flux is always negative).
Complementary to H + N , H N may be calculated as twice the integral of I N between h n and solar midnight, defined by the piecewise function of I N between h n and h s and I LW between h s and solar midnight (i.e., h = ⇡ ), given as follows (note that H N is a negative quantity): Figure 1 shows an example of a half-day I N curve used in the integrals defined in Eqs. ( 14) and ( 16).I N , which is at its peak at solar noon, crosses zero at h n and reaches a minimum at h s .After sunset (i.e., h > h s ), when I SW is zero, under the positive net radiation curve (solid gray line), above the zero line (solid black line), and between the vertical lines of solar noon and h n .H N is represented as twice the integral below the zero line and above the negative net radiation curve (dashed gray line).

Photosynthetically active radiation
The daily photosynthetically active radiation in units of photon flux density, Q n (mol m 2 d 1 ), is calculated based on the number of quanta received (moles of photons) within the visible light spectrum, which also corresponds to the action spectrum of photosynthesis (Monteith and Unsworth, 1990): where vis (unitless) is the visible light albedo and fFEC (µmol J 1 ) is the flux-to-energy conversion factor (Ge et al., 2011).This factor takes into account both the portion of visible light within the total solar spectrum (approximately 50 % (Stanhill and Fuchs, 1977)), and the mean number of quanta in the visible light energy band (approximately 4.6 µmol J 1 (McCree, 1972)).The 1⇥10 6 converts the units of Q n from µmol m 2 d 1 to mol m 2 d 1 .Values for vis and fFEC are given in Table 2.

Condensation
The daily condensation, C n , may be expressed as the waterequivalent of the absolute value of negative net radiation, H N : where E con (m 3 J 1 ) is the water-to-energy conversion factor that relates the energy released or required for a unit volume of water to evaporate or condense at a given temperature and pressure, based on a simplification of the Priestley and Taylor (1972) potential evapotranspiration for a horizontally uniform saturated surface, which may be expressed as follows: where s (Pa K 1 ) is the slope of the saturation vapor pressure-temperature curve, L v (J kg 1 ) is the latent heat of vaporization of water, ⇢ w (kg m 3 ) is the density of water, and (Pa K 1 ) is the psychrometric constant.Standard values may be assumed for certain parameters (e.g., L v = 2.5⇥10 6 J kg 1 ; ⇢ w = 1⇥10 3 kg m 3 ; = 65 Pa K 1 ); however, equations for the temperature dependence of s and L v (e.g., Allen et al., 1998;Henderson-Sellers, 1984) and the temperature and pressure dependence of ⇢ w and (e.g., Kell, 1975;Chen et al., 1977;Allen et al., 1998;Tsilingiris, 2008) are available (see Appendix B).
The barometric formula may be used to estimate the atmospheric pressure, P atm (Pa), at a given elevation, z (m), when observations are not available.Assuming a linear decrease in temperature with height, which is a reasonable approximation within the troposphere (i.e., for z < 1.10 ⇥ 10 4 m), the following equation may be used (Berberan-Santos et al., 1997): where P o (Pa) is the base pressure, T o (K) is the base temperature, z (m) is the elevation above mean sea level, L (K m 1 ) is the mean adiabatic lapse rate of the troposphere, g (m s 2 ) is the standard gravity, M a (kg mol 1 ) is the molecular weight of dry air, and R (J mol 1 K 1 ) is the universal gas constant.Values for the constants used in Eq. ( 20) are given in Table 2.

Evaporative supply
The evaporative supply rate, S w (mm h 1 ), is assumed to be constant over the day and can be estimated based on a linear proportion of the previous day's soil moisture, W n 1 (Federer, 1982): where S c (mm h 1 ) is the supply rate constant (i.e., maximum rate of evaporation) and W m (mm) is the maximum soil-moisture capacity.Constant values for S c and W m are given in Table 2.

Evaporative demand
The evaporative demand rate, D p (mm h 1 ), is set equal to the potential evapotranspiration rate, E p (mm h 1 ), as defined by Priestley and Taylor (1972).E p usually exceeds the equilibrium evapotranspiration rate, E q (mm h 1 ), due to the entrainment of dry air in the convective boundary layer above an evaporating surface (Raupach, 2000(Raupach, , 2001)).E p is related to E q by the Priestley-Taylor coefficient, which may be defined as 1 plus an entrainment factor, !(Lhomme, 1997): The constant value used for ! is given in Table 2.The calculation of E q is based on the energy-water equivalence of I N , ignoring the soil heat flux (Lhomme, 1997): where 3.6 ⇥ 10 6 converts the units of E q from m s 1 to mm h 1 .Note that E q is defined only for positive values (i.e., E q = 0 for I N < 0).The Priestley-Taylor potential evapotranspiration is preferred in this context to the general Penman-Monteith equation for actual evapotranspiration (Penman, 1948;Monteith, 1965), which requires knowledge of stomatal and aerodynamic conductances, or to any of the "reference evapotranspiration" formulae (Allen et al., 1998) that specifically relate to agricultural crops.Daily equilibrium evapotranspiration, E q n (mm d 1 ), is based on the integration of Eq. ( 23) for values of positive I N , or simply the energy-water equivalence of H + N : where 1 ⇥ 10 3 converts E q n from m d 1 to mm d 1 .The daily demand, which is equal to the daily potential evapotranspiration, E p n (mm d 1 ), may be calculated from E q n , as in Eq. ( 22): (25)

Actual evapotranspiration
The calculation of daily actual evapotranspiration, E a n (mm d 1 ), is based on the daily integration of the actual evapotranspiration rate, E a (mm h 1 ), which may be defined as the minimum of the evaporative supply and demand rates (Federer, 1982): where S w (mm h 1 ) is the evaporative supply rate, defined in Eq. ( 21), and D p (mm h 1 ) is the evaporative demand rate, defined in Eq. ( 22).The analytical solution to E a n may be expressed analogous to the methodology used for solving H o and H N and is defined as twice the integral of E a between solar noon and h n , which comprises two curves: S w for 0  h  h i and D p for h i  h  h n , where h i (rad) is the hour angle corresponding to the intersection of S w and D p (i.e., when S w = D p ): which may be expressed as follows: where r x = 3.6⇥10 6 (1+!) E con (mm m 2 W 1 h 1 ).The intersection hour angle, h i , is defined by setting Eq. ( 21) equal to Eq. ( 22) and solving for h: To account for the occurrences when supply is in excess of demand during the entire day, h i should be limited to zero when cos h i 1.For occurrences when supply limits demand during the entire day, h i should be limited to ⇡ when cos h i  1. Figure 2 shows an example of the half-day evaporative supply and demand rate curves.D p (dashed red line) is at a maximum at solar noon and decreases down to zero at h n , www.geosci-model-dev.net/10/689/2017/Geosci.Model Dev., 10, 689-708, 2017 while S w (dotted blue line) is constant throughout the day.
The point where S w equals D p is denoted by the vertical bar at h i .E a (solid gray line), limited by supply during most of the day, follows the S w line between solar noon and h i .During the time between h i and h n , E a is no longer limited by supply and follows the D p curve.After h n , both D p and E a are zero.E a n is represented by twice the area above the zero line (horizontal solid black line), below the E a line, and between the vertical bars of solar noon and h n .

Runoff
The calculation of daily runoff, RO, is based on the excess of daily soil moisture without runoff compared to the holding capacity, W m , and is given by the following equation: where W n ⇤ (mm) is the daily soil moisture without runoff (i.e., Eq. 1 where RO = 0).

Soil moisture
With analytical expressions for C n , E a n , and RO (i.e., Eqs. ( 18), (27b), and (29), respectively), W n may now be calculated by Eq. ( 1).Once W n is calculated, the following limits are checked: The calculation of RO in Eq. ( 29) should prevent W n from being greater than W m , thus satisfying the upper limit of Eq. ( 30).The limiting effect of S w on E a n , through Eqs. ( 27) and ( 28), should, in most cases, prevent W n from falling below zero and satisfy the lower limit of Eq. ( 30); however, due to the assumption that S w is constant throughout the day, there is the possibility that E a n + RO may exceed W n 1 +P n +C n , resulting in negative W n .In these rare cases, in order to maintain the mass balance of the bucket model presented in Eq. ( 1), E a n is reduced by an amount equal to the magnitude of the negative soil moisture.

Bioclimatic indices
One application of the SPLASH model is for the estimation of the surface fluxes required for the calculation of bioclimatic indices.Typically described at longer timescales (e.g., monthly or annually), the daily SPLASH fluxes can be integrated to monthly and annual totals: where X is a model output parameter at a given day (X d ), month (X m ), or year (X a ) and N is the total number of days to sum over for a given month (N m ) or a given year (N a ).
The following sections describe three common bioclimatic indices.

Moisture index
There exists a long history that includes several variants of the moisture index, MI, also commonly referred to as the aridity index, AI, or moisture ratio, MR (Thornthwaite, 1948;Budyko, 1961).A current definition describes MI as the ratio of annual precipitation to annual potential evapotranspiration (Middleton and Thomas, 1997), given as follows: where P a (mm a 1 ) is the annual precipitation and E p a (mm a 1 ) is the annual potential evapotranspiration as calculated by Eq. ( 31); P a and E p a may be substituted with their multi-year means (i.e., P a and E p a ) if available.Values less than 1 are indicative of annual moisture deficit.

Climatic water deficit
The climatic water deficit, 1E, defined as the difference between the evaporative demand (i.e., potential evapotranspiration) and the actual evapotranspiration, has been shown to be a biologically meaningful measure of climate as it pertains to both the magnitude and length of drought stress experienced by plants (Stephenson, 1998).At the monthly timescale, this index is calculated as follows: where 1E m (mm mo 1 ) is the monthly climatic water deficit, E p m (mm mo 1 ) is the monthly potential evapotranspiration and E a m (mm mo 1 ) is the monthly actual evapotranspiration.E p m and E a m are the monthly totals of E p n and E a n , respectively, calculated by Eq. ( 31).Values of 1E may also be computed at the annual timescale.

Priestley-Taylor coefficient
The Priestley-Taylor coefficient, ↵, is the ratio of actual evapotranspiration to equilibrium evapotranspiration, which represents the fraction of plant-available surface moisture (Priestley and Taylor, 1972;Sykes et al., 1996;Gallego-Sala et al., 2010).At the monthly timescale, this is defined as follows: where ↵ m is the monthly Priestley-Taylor coefficient, E a m is the monthly actual evapotranspiration and E q m (mm mo 1 ) is the monthly equilibrium evapotranspiration.Due to the entrainment factor, ↵ m may vary between zero (i.e., no moisture) and 1 + !(i.e., unlimited moisture).Values of ↵ may also be computed at the annual timescale.

Results
The methodology described in Sect. 2 was translated into computer application code (C++, FORTRAN, Python and R).The following sections describe the year-long SPLASH simulation results (2000 CE) at the local and global scales along with comparisons with other model results.

Local temporal trends and bioclimatic indices
The SPLASH model was run at six locations across North America (see Fig. 3), representing six distinct climate regions across latitudinal and elevational gradients.A total of 10 years (i.e., 1991-2000) of monthly CRU TS3.23 data (i.e., precipitation, air temperature, and cloudiness fraction) were extracted from the 0.5 ⇥ 0.5 pixel located over each site.Air temperature and cloudiness fraction were assumed constant and monthly precipitation was divided equally across each day of the month.Fractional sunshine hours were calculated as the one-complement of the cloudiness fraction.Orbital parameters (for paleoclimatology studies) were assumed constant and calculated for the 2000 CE epoch based on the methods of Berger (1978).Model constants were assigned as per Table 2.
The first year of the simulation (i.e., 1991) was iterated (approximately twice) until the daily soil moisture, initialized at zero, reached equilibrium, after which the model was spun-up for eight years (i.e., 1992-1999).The results, shown in Figs. 4 and 5, are for the year 2000.Accompanying the daily SPLASH results in Fig. 4, shown in red, are daily surface fluxes based on the three-layer variable infiltration capacity (VIC) model, extracted from the 1/16 pixel over each of the six sites from the datasets provided by Livneh et al. (2015).
Figure 4a shows the daily results for a tundra region over Banff National Park in Alberta, Canada, with a mean annual temperature of 4 C and annual precipitation of 986 mm.The SPLASH H N depicts a bell-shaped curve characteristic for the Northern Hemisphere.During the spring and summer months, SPLASH H N is higher than the VIC results, which exhibit a lower H N during the first half of the year.The SPLASH W n remains saturated throughout the year at a level between the second and third layers of VIC.SPLASH and VIC E p n are similar in magnitude throughout the year with SPLASH E a n following E p n all year.Figure 4b shows the daily results for a continental warm summer region over the Adirondack region of New York with a mean annual temperature of 5 C and annual precipitation of 1080 mm.There is agreement between the SPLASH and VIC H N and E p n throughout the year.The SPLASH W n remains saturated throughout most of the year with a dry-down period during mid-to late summer and a recovery period during the autumn.
Figure 4c shows the daily results for a temperate region with dry summers over the Bay Area of California with a mean annual temperature of 14 C and annual precipitation of 594 mm.During the dry summer months, SPLASH H N is slightly higher than VIC, during which time the SPLASH W n is depleted, causing moisture-limited E a n to occur.Before which, during the winter and early spring, the SPLASH W n is saturated and E a n follows E p n .Figure 4d shows the daily results for a hot arid desert region in southwestern Arizona with a mean annual temperature of 23 C and annual precipitation of 39 mm.Over the entire year, SPLASH H N is higher than VIC, with the largest differences occurring during the summer months.The SPLASH and VIC W n are both consistently low throughout the year.This water limitation is expressed in the low www.geosci-model-dev.net/10/689/2017/Geosci.Model Dev., 10, 689-708, 2017 Figure 4e shows the daily results for an equatorial monsoonal region near the southern tip of Florida with a mean annual temperature of 24 C and annual precipitation of 1500 mm.There is agreement between SPLASH and VIC H N throughout the year; however, SPLASH W n is higher than all three layers of VIC except for a few days following a large rain event in October.During the drier winter, there is a slight moisture limitation shown in the SPLASH E a n .Throughout the year, SPLASH E p n is slightly higher than VIC.
Figure 4f shows the daily results for a cold arid steppe region in San Luis Potosí, Mexico with a mean annual temperature of 18 C and annual precipitation of 346 mm.During the winter, SPLASH H N is slightly higher than VIC.The SPLASH W n remains low throughout the year at a level between the first and second layers of VIC.The moisture limitation results in a lower SPLASH E a n throughout the year.The SPLASH and VIC E p n agree during the year.Figure 5 shows the SPLASH monthly integrated evapotranspiration results (E p m in solid black, E q m in dotted black, and E a m in dashed red) along with two monthly bioclimatic indices: 1E m and ↵ m .For both the tundra and continental climate sites (Figs.5a and b, respectively), E a m is equivalent to E p m , which results in constant indices for 1E m (i.e., 0 mm) and ↵ m (i.e., 1.26).At the annual timescale, 1E a and ↵ a for these two sites are the same as their monthly values and MI is greater than one, suggesting that these are water-available sites.
Figure 5c shows the monthly SPLASH results for the temperate region with dry summers.Similar to the daily results (i.e., Fig. 4c), during the dry summer, E a m falls below the E p m and E q m curves.This results in a positive 1E m and a drop in ↵ m during the summer months.At the annual timescale, 1E a is 619 mm, which is slightly higher than the annual precipitation (i.e., 594 mm), and both ↵ a and MI are less than 1 (i.e., 0.633 and 0.477, respectively), suggesting that this is a water-limited site.
The hot arid desert region presents a more extreme case as shown in Fig. 5d, where E a m is constantly below both E p m and E q m .This results in a positive bell-shaped 1E m curve and a shallow bowl-shaped ↵ m curve.At the annual timescale, 1E a is 1450 mm, which is significantly higher than the annual precipitation (i.e., 39 mm).Also, both ↵ a and MI are significantly less than 1 (i.e., 0.236 and 0.0219, respectively), suggesting that this is a critically water-limited site.
In contrast, at the equatorial monsoonal site, shown in Fig. 5e, E a m closely follows the E p m curve, which results in an almost zero 1E m and a nearly constant 1.26 ↵ m .At the annual timescale, 1E a is 29 mm, ↵ a is 1.24, and MI is 0.985, which all suggest that this site is not water limited.
Similar to the hot arid desert, at the high elevation of the cold arid steppe, shown in Fig. 5f, E a m is constantly below both E p m and E q m .Unlike the hot arid desert site, the cold arid steppe site is at a lower latitude, which results in a flatter H N curve (as shown in Fig. 4f) that leads to a more constant E p m curve.At the annual timescale, 1E a is 905 mm, which is greater than the annual precipitation (i.e., 346 mm).Both ↵ a and MI are less than 1 (i.e., 0.482 and 0.236, respectively), which suggests that this is a water-limited site.( C), and monthly cloudiness fraction.Monthly precipitation was converted to daily precipitation by dividing the rainfall equally amongst the days in the month.Fractional sunshine hours were calculated based on the one-complement of cloudiness fraction and assumed constant over the month.Mean daily air temperature was also assumed constant over each day of the month.Half-degree land-surface elevation (m above mean sea level) was provided by CRU TS3.22 (Harris et al., 2014).Once again, orbital parameters were assumed constant over the year and calculated for the 2000 CE epoch based on the methods of Berger (1978), and model constants were assigned as per Table 2.

Global simulation of spatial patterns
The SPLASH simulations were driven by the data described above, one pixel at a time, starting each pixel with an empty bucket and terminating when a steady state of soil moisture was reached between the first and last day of the year.Following the spin-up to equilibrate the soil moisture, the model was driven once again to produce daily simulations of net radiation and soil moisture.Figure 7b and e show the SPLASH results of the mean daily relative soil moisture (%) for the months of June and December, respectively.An ice sheet was imposed over Greenland (i.e., no soil moisture).For comparison, the National Center for Environmental Prediction (NCEP) Climate Prediction Center (CPC) Version 2 mean soil moisture (van den Dool et al., 2003) for June and December of the same year are plotted in Fig. 7a and d, respectively.The relative soil moisture in both datasets is computed as the ratio of millimeters of soil moisture over the total bucket size (i.e., 760 mm in NCEP CPC and 150 mm in SPLASH).
Overall, the SPLASH model simulates soil-moisture patterns similar to the NCEP CPC model results.The differences between the SPLASH and NCEP CPC model results are highlighted in Fig. 7c and f.Once again, regions in red indicate where the SPLASH model results are higher than the NCEP CPC model results and regions in gray indicate areas where the SPLASH model results are lower.In contrast to the NCEP CPC soil moisture, the SPLASH model produces a relatively full bucket across wet vegetated regions.The lower relative fullness of the NCEP CPC bucket may be contributed to its significantly larger bucket size.Despite the differing magnitudes of soil moisture, the spatial distributions of soil moisture show consistently drier regions in both simulations at both time periods, especially across mid-northern latitudes (e.g., eastern North America, northern Africa, and central Asia).Seasonal shifts in soil moisture from June to December are also consistently shown (e.g., southern transition in Africa, eastern transition in South America, and northern transition in Australia).There are discrepancies in the spatial distribution of soil moisture across the high-latitude regions (especially Russia).The predominantly saturated conditions in the SPLASH simulations across Russia for December (Fig. 7e) may actually be representative of an increasing snow pack, which could account for these differences.

Discussion
The results presented in Sect. 4 are intended to illustrate the dynamic patterns and trends in the SPLASH model outputs across regions and seasons for a single year under a steady state.The SPLASH model results are promising despite the model's simplifications and limited climatic drivers.At the local scale, the comparison between SPLASH and VIC across climate and elevation gradients (i.e., Fig. 4) shows relatively good agreement for E p n .There are some discrepancies between H N , especially at the high-latitude, high-elevation tundra site (i.e., Fig. 4a) and at the low-elevation hot arid desert site (i.e., Fig. 4d), where the SPLASH simulations were higher than VIC for portions to all of the year.These discrepancies are likely due to local deviations from the glob-ally averaged surface albedo.This is especially true when there is snow cover, as SPLASH does not model snowpacks.Soil moisture also showed relatively good agreement, except at the equatorial monsoonal site (i.e., Fig. 4e), where the SPLASH simulation was consistently higher than VIC.This discrepancy may be due to the assumed constant maximum soil-moisture holding capacity.Furthermore, at the global scale, the SPLASH model reasonably captures the latitudinal gradation of net surface radiation flux (where surface emission assumptions are valid) compared to the CERES EBAF results (i.e., Fig. 6) and produces similar spatial patterns of soil moisture, albeit at different magnitudes, compared to the NCEP CPC soil-moisture results (i.e., Fig. 7).
While the methodology presented in Sect. 2 makes numerous assumptions and simplifications (e.g., saturationexcess runoff generation, invariant soil properties, and constant global parameterization), it provides a simple and robust framework for the estimation of radiation components, evapotranspiration, and plant-available moisture requiring only standard meteorological measurements as input.The SPLASH model currently only produces saturation-excess runoff.For more realistic runoff generation, other water balance models allow runoff to occur when the bucket is less than full, for example the empirical relationship of runoff to the weighted relative soil moisture in the simple water balance model (Orth et al., 2013).Regarding the bucket size, in principle, W m in Eq. ( 21) could be formulated as a property of soil type (as was done, for example, in the original BIOME model); there are some objections to doing so.While W m has a standard definition in agronomy (i.e., the difference between field capacity and wilting point), the wilting point in reality depends on plant properties.Also, the effective "bucket size" depends on rooting behavior, which is highly adaptable to the soil wetness profile.The absolute value of daily soil moisture will be influenced by the bucket size (as shown in Fig. 7) and can have an impact on the local hydrology (e.g., Fig. 4e); however, plant-available moisture indexes, such as ↵ (i.e., the ratio of supply-limited to non-supply-limited evapotranspiration), have commonly been found to be relatively insensitive to the bucket size.Regarding localized effects, the standard values presented in Table 2 are representative of reasonable global means; however, it is recommended that local parameterization (e.g., shortwave albedo) be used if and when data are available.
Over the years, a common misconception has developed regarding the calculation of daily actual evapotranspiration (as defined by Federer, 1982), whereby the integration of Eq. ( 26) is mistakenly interpreted as follows: where D (mm d 1 ) is the total daily demand, given by Eq. ( 25), and S (mm d 1 ) is the total daily supply over the hours of positive net radiation, which may be given by the following: where h n is the net radiation cross-over angle, given by Eq. ( 15), and the constant coefficient converts the units of radians to hours.As shown in Fig. 2, E a n is a piecewise function consisting of two curves overlaid throughout the course of a single day that must be accounted for simultaneously; however, even in some recent model developments, E a n is calculated using Eq. ( 35), including the equilibrium terrestrial biosphere models BIOME3 and BIOME4 (Haxeltine and Prentice, 1996;Kaplan, 2001) and the Lund-Potsdam-Jena dynamic global vegetation model (Sitch et al., 2003).Only under specific circumstances will Eq.( 35) produce correct results.It is the intention of this work to provide a simple analytical solution that correctly accounts for the integration of Eq. ( 26), which has been provided in the form of Eq. (27b).

Code availability
The code, in four programming languages (C++, FOR-TRAN, Python, and R), is available on an online repository under the GNU Lesser General Public License (https: //bitbucket.org/labprentice/splash).The repository includes the present release (v1.0) and working development of the code (with makefiles where appropriate), example data, and the user manual.All four versions of the code underwent and passed a set of consistency checks to ensure similar results were produced under the same input conditions.The following describes the requirements for compiling and executing SPLASH v.1.0.
For the C++ version, the code was successfully compiled and executed using the GNU C++ compiler (g++ v.4.8.2) provided by the GNU Compiler Collection (Free Software Foundation, Inc., 2016).It utilizes the C numerics library (cmath), input/output operations library (cstdio), and the standard general utilities library (cstdlib), and it references the vector container and string type.
For the FORTRAN version, the code was successfully compiled and executed using the PGI Fortran compiler (pgf95 v.16.1-0) provided by The Portland Group -PGI Compilers and Groups (NVIDIA Corporation, 2016) and the GNU Fortran compiler (gfortran v.4.8.4) provided by the GNU Compiler Collection (Free Software Foundation, Inc., 2016).
For the Python version, the code was successfully compiled and executed using Python 2.7 and Python 3.5 interpreters (Python Software Foundation, 2016).It requires the installation of third-party packages, including NumPy (v.1.10.4 by NumPy Developers, 2016) and SciPy (v.0.17.0 by SciPy Developers, 2016) and utilizes the basic date-and time-type (datetime), logging facility (logging), Unix-style pathname pattern extension (glob), and miscellaneous operating system interface (os) modules.
For the R version, the code was successfully compiled and executed using R-3.2.3 "Wooden Christmas-Tree" (The R Foundation for Statistical Computing, 2015).Constants for M a and M v are given in Table 2.The temperature dependence of C p may be assumed to be negligible (e.g., C p = 1.013 ⇥ 10 3 J kg 1 K 1 ) or calculated by the following (Tsilingiris, 2008): for T air between 0-100 C. The coefficients of T air are given in Table B1.

Figure 1 .
Figure 1.Example of the net radiation flux curve between the hours of solar noon (i.e., h = 0) and solar midnight (i.e., h = ⇡ ).The I N curve is equal to the difference between the net incoming shortwave radiation flux, I SW (solid red line), and the net outgoing long-wave radiation flux, I LW (dotted blue line).Positive I N , shown decreasing from solar noon to zero at the cross-over hour angle, h n , is denoted with a solid gray line, while negative I N , shown decreasing from zero to I LW between h n and the sunset hour angle, h s , and constant between h s and solar midnight, is denoted with a dashed gray line.The solid black horizontal line marks the datum of zero radiation.

Figure 2 .
Figure 2. Example of actual evapotranspiration curve between the hours of solar noon (i.e., h = 0) and solar midnight (i.e., h = ⇡).The evaporative demand, D p (dashed red line), is at a maximum at solar noon and zero at the cross-over hour angle, h n .The evaporative supply, S w (dotted blue line), is constant throughout the day.The point where supply is equal to demand denotes the intersection hour angle, h i .Actual evapotranspiration (solid gray line) is defined as the minimum of S w and D p throughout the day.

Figure 4 .
Figure 4. Daily simulations of net radiation, H N , soil moisture, W n , and evapotranspiration, E n , for the six climate regions defined in Fig. 3: (a) tundra, (b) continental with warm summers, (c) temperate with dry summers, (d) hot arid desert, (e) equatorial monsoon, and (f) cold arid steppe.Black lines represent SPLASH modeled net radiation, soil moisture, and evapotranspiration (potential in solid black and actual in dashed black).Red lines represent VIC three-layer modeled surface fluxes from Livneh et al. (2015) for net radiation, soil moisture (layer 1 in solid red, layer 2 in dashed red, and layer 3 in dotted red), and potential evapotranspiration.Days of the year are represented along the x axis.Data are for 1 year (2000 CE).

For
Figure 5. Monthly SPLASH results of evapotranspiration, E m (potential in solid black, actual in dashed red, and equilibrium in dotted black), climatic water deficit, 1E m , and Priestley-Taylor coefficient, ↵ m , for the six climate regions defined in Fig. 3: (a) tundra, (b) continental with warm summers, (c) temperate with dry summers, (d) hot arid desert, (e) equatorial monsoon, and (f) cold arid steppe.Months of the year are represented along the x axis.Results are of 1 year (2000 CE).
Figure 6b and e show the SPLASH results of the mean daily net surface radiation flux (MJ m 2 ) for the months of June and December, respectively.For comparison, the Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) average all-sky surface net total flux for June and December 2000 are plotted in Fig. 6a

Figure 6 .
Figure 6.Global mean net downward surface radiation flux (MJ m 2 ) for June (left) and December (right) 2000 CE from CERES EBAF (top), from the SPLASH model (middle), and showing the differences between SPLASH and CERES EBAF (bottom).(a) CERES EBAF June 2000; (b) SPLASH June 2000; (c) difference between SPLASH and CERES EBAF June 2000; (d) CERES EBAF December 2000; (e) SPLASH December 2000; and (f) difference between SPLASH and CERES EBAF December 2000.Color bars on the right are linear interpolations of results in megajoules per square meter.Differences in red indicate higher SPLASH model results.

Figure 7 .
Figure 7. Global mean relative soil moisture (%) for June (left) and December (right) 2000 CE from NCEP CPC (top), from the SPLASH model (middle), and showing the difference between SPLASH and NCEP CPC (bottom).(a) NCEP CPC June 2000; (b) SPLASH June 2000; (c) difference between SPLASH and NCEP CPC June 2000; (d) NCEP CPC December 2000; (e) SPLASH December 2000; and (f) difference between SPLASH and NCEP CPC December 2000.Color bars on the right are linear interpolations of results in units of relative soil moisture (%).The relative soil moisture is based on the total bucket size (i.e., 760 mm for NCEP CPC and 150 mm for SPLASH).Differences in red indicate higher SPLASH model results.

Table 2 .
Constants and Standard Values.

Table B1 .
Coefficients of T air .