III. Surface Layer Processes
The atmosphere and the Earth’s surface, whether land or ocean, exchange momentum, heat and moisture through a variety of processes. The exchange takes place in a relatively thin layer of the atmosphere called the surface layer (distinct portion of the planetary boundary layer - see section IV - nearest the Earth’s surface). The COLA GCM employs an empirical relation that approximates Monin-Obukhov similarity theory to determine the transfer or drag coefficients. The empirical relations are different for ocean and land surfaces.
A. Oceans
Latent and sensible heat fluxes from the surface of the world's oceans to the atmosphere are important energy sources for the atmosphere. Momentum flux between the atmosphere and oceans produces a wind stress on the ocean surface that drives most ocean currents. These are parameterized according to bulk aerodynamic schemes in which the flux is assumed to be proportional to the surface wind speed and the potential (moisture or temperature difference between ocean surface and air near the surface). The coefficient of proportionally or drag coefficient is determined by an empirically determined functional fit. Over the ocean surface the functional fit is similar to that shown in Sato et al. (1989b):
Cm and Ch are multiplied by the magnitude of the wind velocity to yield the drag coefficients for momentum and heat/moisture transfer. The functions Fm and Fh are determined as follows. First, set
In the stable region: Ri > 0
where
For the unstable region: Ri < 0
where
and
where
In the above formulae, zr is the reference height, z0 is the roughness length and Ri is the bulk Richardson number. This approach, which differs from the Monin and Obukhov formulation, was adopted to improve computational efficiency.
Sea ice points are treated differently. The temperature of the sea ice is computed by solving the following equation:
Where the subscript, i, refers to ice and the subscript, 1, refers to the lowest s level of the model. Fluxes and values of model variables from the previous time step are used to estimate the new values, the drag coefficients are recomputed, and the new sea ice temperature is then computed.
B. Land
The land surface of the Earth is composed of a variety of different plants, soil and geographical formations, all of which exchange mass, momentum and heat with the atmosphere through several processes in varying degrees. Many studies have shown the differences in various fluxes between the land surface and the atmosphere depend upon whether the fluxes are measured over bare soil or over vegetated areas. Further refinement of the dependence of the fluxes on vegetation type has also been achieved. The fluxes which are most significantly dependent upon the vegetation type are the radiative flux (e.g., Charney et al., 1977), the latent heat flux (e.g., Shukla and Mintz, 1982), and the momentum flux (e.g., Sud and Smith, 1985). It has been shown that seasonal simulation of the atmospheric circulation near the surface depends heavily upon the nature of the surface fluxes, particularly of latent heat (Sato et al., 1989a). Other studies have shown that the latent heat flux has a direct bearing upon the accuracy of the short to medium range forecast of surface temperature and humidity (Shukla et al., 1990). As a result of these studies, the COLA GCM includes an explicit formulation of the land surface vegetation and its exchanges with the at
Simplified SiB(SSIB) Model:
In SSIB, the following processes are taken into account:
As shown by Sato et al.(1989a), using a realistic treatment of atmosphere-biosphere interaction leads to:
The vegetation-soil layer affects the radiative transfer at the surface, the partitioning of surface energy into sensible heat flux and latent heat flux and the transfer of momentum between land and atmosphere. SSiB is intended to realistically simulate the biophysical exchange processes. The biophysical controls on these exchange processes are modeled in a mutually consistent manner by representing the effects of vegetation explicitly. Each COLA GCM land grid point is assigned one of 12 vegetation types (or permanent ice cover), whose physical and morphological properties are specified and are dependent upon the time of year (Dorman and Sellers, 1989).
The development of SSIB was based upon observational data from the Amazon rain forest (Xue et al., 1991). A number of observational data sets from different sites and vegetation types have been used to validate and evaluate this model, including data from: the Anglo-Brazilian Amazonian Climate Observation Study (ABRACOS) - a field experiment over an Amazon deforestation site (Xue et. al., 1996); the Sahelian Energy Balance experiment (SEBEX) - a field experiment over a semi-arid region (Xue and Allen, 1995); the First ISLSCP field experiment (FIFE) in Kansas (F. Chen et al., 1995); the series of Russian hydrological measurements (Robock et al., 1995, Schlosser, 1995); the Cabauw experiment on a grassland site in Netherlands (T. Chen et al., 1995b); and the Hapex-Mobilhy field experiment - a series of observations in a cultivated (soya) field mixed with coniferous forest (Xue et al., 1996b). Most of the following description of SSiB is taken from Xue et al. (1996a and 1996b).
SSiB has three soil layers and one canopy layer, and eight prognostic variables: soil wetness in the three soil layers; temperatures at the canopy, ground surface and deep soil layers; water stored on the canopy, and water (liquid or frozen) stored on the ground (Fig. 2). The governing equation for canopy temperature Tc is based on the energy conservation equation:
where Cc, Rnc, Hc , and l Ec are the heat capacity of the canopy, and the net radiation, sensible heat, and latent heat fluxes at the canopy level, respectively. In the radiative transfer submodel, the optical and geometric properties of the leaves and stems, and the optical properties of the soil affect the surface albedo and the attenuation of photosynthetically active radiation (PAR) down through the canopy. The surface albedo is modeled to have a diurnal variation with a minimum at local noon. Since the variation of the albedo with solar angle is quite regular, a quadratic equation is used to describe the albedo and its diurnal variations. The equation was fitted to the results of two-stream radiative transfer model calculations. For a specific vegetation type, the albedo a is also a function of the solar zenith angle q and snow cover S:
The coefficients of the quadratic equation (a,b,c,d, and e) depend on vegetation type.
The force-restore method is used to predict the time variation of the soil temperature Tgs within a shallow layer where the diurnal cycle of temperature predominates:
where t is the day length, Cgs is the effective heat capacity of soil, Td is the temperature for deep soil, and Rngs, Hgs, and l Egs are net radiation, sensible heat, and latent heat fluxes at the ground, respectively. The equation for deep soil temperature, Td , representative of the soil layer dominated by the annual cycle is:
The depths at which Tgs and Td are computed are unspecified, and actually vary in space and time as a function of soil properties and soil moisture, through variations in heat capacity.
The governing equation for the canopy interception water store Mc is based on water mass conservation:
where Pc, Dc, and Ewc are the precipitation rate, water drainage rate, and evaporation rate from the wetted portions of the vegetation canopy, respectively. The canopy has a maximum instantaneous water store of (1/10 LAI) mm liquid water equivalent. The governing equation for snow depth Ms is
where Ps, Sm, and Wf are snowfall, snow melt, and the freezing of water, respectively.
The water that accumulates on the ground and canopy surfaces may be in either liquid or solid form (snow). As precipitation nears the ground in the model, it is assumed to instantly equilibrate with the air at the reference temperature, Tr. If Tr is below freezing, then the precipitation is assumed to be snow. If the temperature of the vegetation canopy or the ground surface is below freezing, water storage on the ground is assumed to be frozen and accumulates. Snow is assumed to accumulate entirely on the ground. Once snow has accumulated, if Tr rises above freezing, then the snow begins to melt. If the temperature at either one of the two subsurface model levels is above freezing when snow is melting, then the meltwater infiltrates into the ground. If the surface layer soil moisture becomes saturated, then meltwater runs off. Finally, if snow is melting but both subsurface temperatures are below freezing, then meltwater is assumed to run off. Snow accumulation on the ground thus has a direct effect on the GCM surface hydrologic balance. Snow increases the surface albedo substantially in the model. Snow also significantly alters the partitioning of the surface latent and sensible heat fluxes in favor of the latent heat flux. The effects of sublimation, cloud physics, and blowing snow are neglected in the model. The initial value of snow depth is specified based on observations (see section VIIA).
In the three soil layers, water movement is described by a finite-difference approximation to the diffusion equations:
where q 1, q 2, q 3, D1, D2, and D3 are the volumetric soil water content and soil thickness of the top, middle, and lower soil layers, respectively. The term I is the infiltration where I = Dc - Ru (Ru is the instantaneous surface runoff). The variable Edc is the transpiration rate, and Egs is the evaporation from bare soil. The coefficients bi (i=1,2,3) are fraction factors which depend on the root distribution. Qi,j is the transfer of water between the ith and jth layers and is defined to be positive upward as:
where k is the hydraulic conductivity, y is the soil water potential, and z is the thickness between two soil layers. The soil water potential in SSiB is taken from the empirical relationship of Clapp and Homberger (1978):
where y s is the volumetric soil water content at saturation and the B parameter is an empirical constant dependent on the soil type. The drainage of water out of the bottom layer is:
The first term of the right hand side is contributed by gravity only, with no diffusive transport occurring, as modeled by Sellers et al. (1986). The variable k is the mean slope angle and set to larger than 3° and k3 is the hydraulic conductivity of the third layer, and Qb is the base flow runoff and is proportional to the soil wetness in the lowest soil layer, Qb = Kc * (q 3/q s). q 3 is the volumetric soil water content in the third soil layer. This term was suggested by Liston et al. (1994) to account for the GCM sub-grid scale variation of soil moisture. The empirical constant Kc is derived from large river basins. The hydraulic conductivity at layer i, ki, is
where ks is the hydraulic conductivity at saturation. The logarithm of k decreases linearly as B increases. Hydraulic conductivity changes more dramatically with B in the dry soil than in the wet soil; and has a larger variation with soil moisture when B is large.
The parameterization of the stomatal resistance, rc, in SSiB was based on the work of Jarvis (1976). Three stress terms are included in this scheme which describe the dependence of stomatal resistance on atmospheric temperature, soil water potential, and vapor pressure deficit. Unlike the traditional approach in which the leaf water potential is used to calculated the stomatal resistance (Federer, 1979), soil moisture is used to control stomatal resistance in SSiB. Using leaf water potential makes the computation very complex and a large number of parameters are required. Evidence (Blackman and Davies; 1985, Wetzel and Chang, 1987) indicates that the stomatal response to water supply is not controlled by the plant's internal water potential. Instead, it appears that roots "sense" the soil moisture supply and directly transmit chemical messages (cytokinin) to the guard cells to keep stomata open. Thus, stomatal closure occurs in direct response to soil moisture and is not delayed while a plant's internal water supply is squandered. The empirical formula for the adjustment factor f(y ) for soil water potential is:
where c2 depends on the vegetation type and c1 is a constant which is obtained using the wilting point. The stomata completely close at the wilting point in the model. The variable c2 is a slope factor; a large value of c2 means that f(y ) changes from 0 to 1 very fast when soil water content varies from the wilting point to the point stomata start to close. Note that in Table 1 of Xue et al. (1991) the values of c1 and c2 should be interchanged.
The equation for transpiration from the canopy is:
where q(Tc) and qa are the saturated specific humidity at the canopy temperature and specific humidity of the canopy air space, respectively. The term rb is the bulk boundary layer resistance, rc is the stomatal resistance, and wc is the wetness fraction of the canopy. Bare soil evaporation in the model is given by:
where q(Tgs) is the saturated specific humidity at the surface temperature Tgs, and Vg is the vegetation cover. The resistance to the transfer of water vapor from the upper soil layer to the canopy air space includes aerodynamic resistance rd and soil surface resistance rsurf (Fig. 1). The results of Camillo and Gurney (1986) were used to fit a simple relationship between soil surface resistance and soil moisture in the first layer:
The relative humidity of the air at the soil surface is:
As usual, g is the acceleration due to gravity and Rv is the gas constant for water vapor. Similarity theory was used to calculate the aerodynamic resistance from the canopy to the reference height. Based upon the equations of Paulson (1970) and Businger et al. (1971), a linear relationship between the Richardson number and the aerodynamic resistance was developed. The equation for momentum flux transfer is:
where U* is the friction velocity, Ur is the wind velocity at the reference height, Cun is the neutral momentum transfer coefficient, and Cu is the non-neutral coefficient. Louis (1979) developed a parameterization scheme in which the total aerodynamic resistance including both neutral and non-neutral parts is a function of Richardson number and surface roughness. The accuracy of this parameterization strongly depends on the surface roughness length. In a biosphere model, the range of values of surface roughness length can by very large. Parameterizing the total resistance may cause large errors. In SSIB, we only parameterize the non-neutral part; Cun, which is dependent on vegetation and soil properties, is given by:
where k0 is von K´ rm´ n’s constant, zm is the reference height, d is the displacement height, and z0 is the roughness height. For the non-neutral part,
and
where u is the wind speed and q is the potential temperature. Xue et al. (1991) fitted an exponential equation to relate the bulk Richardson number to surface resistance that reproduces the results of similarity theory. Eddy flux in each grid point in a GCM is assumed to be the average of the influence of a large number of sub-grid-scale eddies. Sud and Smith (1985) found that the relationship between aerodynamic resistance and ensemble mean Richardson number had a tendency to change from exponential to linear. The linear parameterization described here has produced very good flux simulations in many cases (e.g., Shao et al. 1994). The equations for heat flux transfer are:
for latent heat flux, E, and
for sensible heat flux, F. Variables Tm, qm, Ta, and qa are the temperatures and specific humidity at the reference height and in the canopy air space, respectively. Parameter CTN is the neutral heat transfer coefficient and ra is the aerodynamic resistance. Coefficient CTN is given by
where z2 is the height of the canopy and g3 is the ratio of actual ra to ra calculated using a log-linear wind profile assumption. It is held constant at 0.75. The depth of the transition layer above the canopy is zr. Above this transition layer, the log-linear assumption is valid. Following Sellers et al. (1989), we use:
For the non-neutral heat transfer coefficient, we have:
for 0 > Ri > -10, where
For 0.16 > Ri > 0, we have:
