Introducing lateral subsurface flow in permafrost conditions in a distributed
land surface scheme
Petra Koudelova1 and Toshio Koike1(1: Dept. of Civil Engineering, The University of Tokyo)
Petra Koudelova, Dept. of Civil Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
113-8656 Tokyo, Japan. e-mail: petra@hydra.t.u-tokyo.ac.jp
Abstract
This study introduces a quasi-3D land surface scheme, which accounts for both vertical and horizontal hydrologicalprocesses in regions with permafrost conditions. The attribute “quasi-3D” implies an explicit linkage between thevertical and horizontal componentss of the model. The model is developed by implementing the land surface scheme(LSS) SiB2, which incorporates a 1-D frozen soil parameterization (FSP), in a Distributed Hydrologic Model (DHM).
Comparative simulations, in which the uncoupled SiB2 and the quasi-3D model considering surface lateral flow areapplied to a small catchment in Tibet, demonstrate the impact of surface lateral flow on surface water and energybalance. A numerical experiment focusing on the effects of the permafrost conditions confirms the importance of lateralflow in hilly regions with frozen soil, and shows considerable differences in surface water and energy balance between
the upper parts of slopes and valleys.
Keyword: Tibet, land surface scheme, permafrost, distributed hydrological model

1. Introduction
The Tibetan Plateau is not totally flat; instead,undulating hills alternate with plains of various extents.
The Plateau is also a region with wide areas of frozen
soil, either permafrost (50%) or seasonally frozen soil
(44%) (Li and Koike, 2003). In regions with uneventerrain, the lateral flow is an important process
affecting the spatial distribution of the surface wetness,
which in turn essentially influences surface fluxes tothe atmosphere. Moreover, specific phenomena areassociated with frozen soil in hilly regions. Since thefrozen soil serves as an impervious table, liquid waterprovided by either the melting of the soil ice content or
the diffusion of soil moisture from upper layers
accumulates above this table. Consequently, theoverlying soil becomes saturated and thus saturated
flow is initiated along slopes. In addition, if the soil isfrozen up to the surface, surface flow occurs in the caseof rain or snowmelt. As a result, the upper parts of hills
are relatively dry, while the valleys are wet. Due to thelower liquid soil moisture content and hence the lower
heat capacity, the upper parts of the slopes thaw much
faster than the bottom ones. Accordingly, the active(unfrozen) layer at the bottom is quite shallow incomparison with the top, holding abundant liquid soil
moisture, which may produce a saturation excess anda ponding of water on the surface.
Because of its impact on surface fluxes, surface
wetness and its spatial distribution must be properly
represented in a land surface scheme (LSS). Regarding
the specific conditions in Tibet and the phenomena
described above, an LSS applied to this region should
account for the lateral transport of water as well as theprocesses associated with the freezing/thawing of soil.
However, common LSSs are one-dimensional (1-D)
models describing only vertical processes. In this study,
we introduce a quasi-3D land surface scheme, which
accounts for both vertical and horizontal hydrologicalprocesses, while regarding the frozen soil conditions.
The attribute “quasi-3D” implies an explicit linkage
between the vertical and horizontal components of themodel. The model is developed by implementing theland surface scheme SiB2 (Sellers et al., 1996)

incorporating a frozen soil parameterization (FSP) (Liand Koike, 2003) in a Distributed Hydrological Model(DHM).
We carry out two numerical experiments. In thefirst experiment, we investigate the effects of surface
lateral flow on the spatial distribution of the surface
layer soil moisture, surface temperature, and heat
fluxes as well as on the spatially averaged values ofthese variables over the computational domain. For thispurpose, we carry out comparative simulations, inwhich the coupled model and the uncoupled SiB2 areapplied to a small catchment in the upper reach of theNaqu basin. The second experiment focuses on theeffects of surface and subsurface lateral flows on a
local water budget in hilly permafrost regions. The
coupled model is applied to a single slope representing
typical topography in the Plateau.

2. Model description
The full descriptions of the SiB2, the FSP
incorporated in the SiB2, and the DHM adopted in the
coupled model are omitted in this paper, since they areavailable in referred literature: (Sellers et al., 1996) for
SiB2, (Li and Koike, 2003) for the FSP, and (Dutta etal., 2000) for the DHM. Here, we introduce
modifications of individual parts, new features and thecoupling scheme.
In the coupled model, the SiB2 incorporating the
FSP is embedded into the framework of the DHM andsolves all of the vertical processes for each grid cell
individually using the meteorological forcing data. It
considers vertical interactions between the surface andthe overlying atmosphere, while neglecting thehorizontal interactions between the neighboring gridcells. Among other prognostic variables, the SiB2
estimates surface runoff as a result of infiltration andsaturation excesses, and the thickness of the saturated
zone above the frozen table as well as the soil watercontent available for lateral subsurface saturated flow.
These quantities are then treated by the saturated zone
component and the surface flow component of theDHM. The updated values of the saturated zone andthe surface water storage are used as an initial state in
the next time step in the SiB2.
Since the SiB2 is a vertical scheme, the only
possible water input is precipitation. In order tointroduce the effect of surface lateral flow on thesurface wetness and fluxes, a simple procedure
analogous to the precipitation-interception-infiltrationrunoff scheme was introduced to treat the surface
lateral flow contribution (lateral inflow) in the SiB2.
The procedure is implemented at the very beginning of
the SiB2 calculation. Firstly, the lateral inflow fills thesurface water interception store. Then, the lateral
inflow, lowered by interception, is partitioned into theinfiltration rate and the infiltration excess contributing
to the surface runoff. It is assumed that the surface
water is transported by means of small rills occupyinga fraction of the whole of the grid cell area.
Once the surface runoff contribution has been
determined, some calculations are carried out toaccount for possible phase changes. Prior to the lateral
inflow, the ground surface is at a specified temperature
Tg and can be covered with snow. An energy balancecalculation is performed to calculate a single newtemperature for the ground surface, the surface runoff
contribution, and the frozen or melted water produced
as a result of phase changes. The patchy snow
treatment included in the SiB2 is applied when thelateral inflow comes to a grid cell with non-zerosnow/ice storage. As a result, some portion of thesurface flow may refreeze on the surface if the surface
temperature is less than Tf, or some snow melts due tothe lateral inflow. When the treatment of the lateral
inflow has proceeded, the handling of precipitation
starts and the next parts follow as in the original SiB2.
The FSP in the SiB2 is based on an approximationof Stefan solution and predicts frozen/thawed depth
and phase changes of soil water content over time. The
three-soil-layer structure in the SiB2 is maintained in
the original FSP (Li and Koike, 2003), but the
governing equations of water balance and surface heatbalance are modified to involve the soil
freezing/thawing process. In the contrast with theoriginal SiB2 (Sellers et al., 1996), the presented modelaccounts for vertically heterogeneous soil, i.e. different
porosities and hydraulic conductivities in the threelayers. The resolution of the three SiB2 soil layers is,however, too coarse for the prediction of a saturated
zone above the frozen table. Therefore, we introduce a
multi-layer soil model in this study. The calculation ofliquid water and ice contents is kept unchanged. The
solution of vertical unsaturated flow is expanded forthe multi-layer structure.
The depth of the saturated zone is determined after
the unsaturated flow calculation. Firstly, the layers
above the frozen table are checked for saturation,
starting from the bottom layer. If there is a continuous
saturated zone comprising at least one layer just above
the frozen table, it represents the saturated zone. If thebottom layer is not saturated, the concept of Ishidaira
et al. (1998) is adopted. It is assumed that the bottom
layer is wetter than the upper ones. If this assumption isvalid, the thickness of the saturated zone is determinedfrom the equation:

D θ l ,b = D ′θ s ,b + (D − D ′)θ l ,b −1 ,(1)where D is the thickness of the bottom layer, D’ is the
thickness of the saturated zone, θl,b is the averagevolumetric liquid content in the bottom layer b, θl,b-1 isthe average volumetric liquid content in the layer b-1,
and θs,b is the saturated water content in the bottom
layer. If the condition θl,b > θl,b-1 is not fulfilled or if theactive layer consists of only the uppermost layer, thevalue of θl,b-1 in Eq. (1) is replaced with 0.9θl,b.
The 2-D subsurface saturated flow scheme is basedon the groundwater component of the DHM, which
employs a non-steady Boussinesq equation. . In thecase of frozen soil, the impervious bed moves up or
down over time according to the freezing/thawing
process. Moreover, the thickness of the saturated zone
may vary greatly along a slope and over time.
Accordingly, the depth of the impervious bed and thesaturated zone are calculated in every time step in the
SiB2. Interactions between the saturated zone and theoverlying soil are neglected in the 2-D subsurface flowrouting. The moisture content available for the routing
is determined by subtracting the residual liquid content
from the saturated value. The obtained value represents
the aquifer storage coefficient and is calculated in eachtime step. The growth/drop of the water head due to thesaturated flow is converted into an increase/decrease in
the liquid moisture content in the affected layers. The
subsurface flow calculation is followed by the 2-D
surface flow component and the river flow component
in the DHM.

3. Numerical experiments
3.1. Surface lateral flow effects
In this experiment, we intend to demonstrate theeffects of the surface lateral flow on the land surface
wetness, and hence on the surface fluxes. We perform
simulations by both the proposed quasi 3-D model andthe uncoupled SiB2 and mutually compare thetemporally and spatially averaged values of heat fluxes,
surface layer soil moisture and surface temperature.
This experiment focuses on surface flow, thus
subsurface lateral flow is not considered in the quasi
3-D model.
The small catchment Mardi (350 km2) in the upper
reach of the Naqu basin (Nu-Jiang river) was chosenbecause of its relative homogeneity in soil conditionsand vegetation. The Mardi catchment is locatedbetween N 32°22’ and N 32°44’, and between E 91°38’
and E 91°55’. Its altitude varies from 4800 m asl up to
5350 m asl, with the most frequent values around 4850
– 5100 m asl. The vegetation type is short grass
combined with bare soil and the soil is classified as a
sandy-loam type. The river channel network is quite
dense, but consists of many shallow streams.
Two simulations were carried out from July 4th toAugust 4th with a time step of 1 hour. The first half ofthe simulation period was rather dry. The rainy phasestarted on July 20th and continued until the end of thesimulation. The summer period was chosen in order toinsure a sufficient amount of rainfall for surface runoff
generation (monsoon season). The forcing data wereassumed to be uniform over the whole of the catchment.

Tab. 1: Statistical parameters of 2D distributed results for the wet period.

Sensible heat (Whm-1)
Sim. 1
Sim. 2Average
14200
17100
St. deviation
350
3Rel. deviation
3%
0.02 %

Wet period

Latent heat (Whm-1)
Sim. 1
Sim. 2
30500
25900
580
5
2%
0.02 %

Time cumulative values
Sensible heat flux

19
16
13

31

0.30

10.0

28

0.27

9.6

-3

25

3

TS2 only
Simulation 2

22
19
16
13

dry wetPeriodGround surface temperature

GST (°C)

22

Surface soil moisture

VWC (m m )25

TS2-IISDHM
Simulation 1

LHF (kWh m-2)

SHF (kWh m-2)

28

Temperature (°C)
Sim. 1
Sim. 2
8.9
9.4
0.06
0.002
0.6 %
0.02 %Catchment average values
Time average values

Latent heat flux

31

Soil moisture (m3m-3)Sim. 1
Sim. 2
0.290
0.278
0.003
0.0001
1%
0.04%

0.24
0.21
0.18

8.8
8.4

0.15
dry wetPeriod9.2

8.0
dry wetPerioddry wetPeriodFig. 1: Comparison of catchment averages: total or averaged values over the dry and wet periods. The black
bar represents Simulation 1 (quasi 3-D model); the white bar stands for Simulation 2 (uncoupled SiB2).
Simulation 1 was performed by employing the quasi
3-D model, while Simulation 2 was carried out by
applying the uncoupled SiB2 to each grid cell within
the catchment. All of the parameters and input data
were identical for both simulations and thus the results
reveal the effects of surface lateral flow only.
A profound spatial variability in the surface
wetness and fluxes is caused by the surface flow acting
in two ways: (1) the direct effect on surface fluxes
through the temporary surface water storage, and (2)the impact on the surface layer soil moisture throughinfiltration from the surface flow and in turn on heat
fluxes. As a result, flat areas located below slopes are
extremely wet with high latent heat and low sensible
heat flux, while the tops of hills and declining areas are
relatively dry, producing lower latent heat and higher
sensible heat fluxes. Such a variability is provided by
the quasi 3-D model, while the uncoupled SiB2
produces a uniform spatial distribution of theinvestigated quantities, which is demonstrated in Table
1, which lists some of the statistic parameters. Because
the mesoscale spatial variability of surface fluxes
stimulates atmospheric convective activity, such anactivity may be missed in mesoscale atmosphericmodels when coupled with LSSs disregarding lateralflow.
Fig. 1 compares the catchment-averaged values ofheat fluxes, surface layer soil moisture and surface
temperature from Simulations 1 and 2. The values areeither summed (fluxes) or averaged (moisture andtemperature) over the dry period and the wet period.
There is no difference between Simulation 1 andSimulation 2 during the dry period because surface
runoff is not generated. Nevertheless, when surface
flow occurs during the wet period, evident differences
appear between Simulations 1 and 2. In Simulation 1,
the surface runoff generated by the SiB2 isconsequently treated by the surface flow component,
i.e. the total water is conserved within the catchment.
However, in Simulation 2, the surface runoff is lost andthus the total water within the catchment isunderestimated. The higher surface and subsurface
water contents in Simulation 1 lead to a higher latent
heat flux and lower sensible heat flux and surface
temperature in comparison with Simulation 2.

3.2. Lateral flow in permafrost conditions
The main aim of this experiment is to show theeffects of surface and subsurface lateral flows on a
local water budget in hilly permafrost regions.
Therefore, we set up a simple catchment represented bya single slope inclining in one direction only. The
catchment follows a real slope in the Naqu basin. The
model is run using the forcing data provided by the
CEOP Tibet 2002 observation. The simulation period
is about 40 days starting on April 1st, when the soil isfrozen up to the surface along the whole of the slope.
Due to the lack of observations of soil moisture andtemperature, we apply hypothetical initial conditions.
We set up two scenarios. One scenario (C1) starts withhomogeneous initial conditions along the slope
represented by the water content at saturation (ice +liquid). The second scenario (C2) starts withheterogeneous initial conditions represented bydifferent soil water contents along the slope: saturation
at the bottom of the slope and gradually decreasing up
the slope. The forcing data and all of the parameters
associated with soil and vegetation are homogeneousover the catchment and are identical for both scenarios.
The results are shown in Figures 2 – 4.
The C1 scenario reveals a strong impact of thelateral flow on the surface wetness. At the bottom, thesoil near the surface stays saturated, while the top of
the slope is gradually drying up (Fig. 2). Consequently,

-1

1
0

Soil moisture - C2 scenario

top - surface
top - root

5/6

2
5/11

bottom - surface
bottom - root

Fig. 2 Volumetric water content (liquid + ice) in
the surface layer (full lines) and the root zone
(dash-dotted lines) at the top (bold lines) and at
the bottom (thin lines) of the slope. The thicknessof the surface layer is 4 cm and of the root zone is
16 cm. The top chart – homogeneous initialconditions and the bottom chart – heterogeneous.
thawed depth (m)Active layer evolution
0.0
0.2
0.4
0.6
0.8
1.0
4/1

top - C2
top - C1
bottom - C2bottom - C1
4/6 4/11 4/16 4/21 4/26 5/1

150
100
50
0

5/6 5/11

Fig. 3 Evolution of the depth of the active layerover time at the top of the slope (bold lines) andthe bottom of the slope (thin lines). Full lines
represent the C2 scenario and the dash-dottedlines show the C1 scenario.
thawing occurs faster at the top than at the bottom(Fig. 3). In addition, the wet bottom produces a higher
latent heat flux and a very low sensible heat flux
(Fig. 4). The extremely low sensible heat flux during
the first 20 days of the simulation is caused by thefrozen surface covered with an ice film. The C2
scenario demonstrates the effect of the initialconditions prior to thawing. Because there is only alittle soil ice content at the top of the slope, there is no
moisture available for lateral flow after thawing. In
addition, only little rain falls during the simulated
period and thus the soil moisture content is lowering
even at the bottom of the slope. The results for thebottom in C2 are very similar to the results at the top in
C1 because the initial state at the bottom in C2 is thesame as at the top in C1 and because of the lack ofwater supply in C2 (little rain, little moisture in theupper part of the slope).

4. Summary
Regarding the specific conditions in the Tibetan
Plateau, which make the lateral transport of water an
important factor in determination of the local water andenergy budget, we proposed the quasi 3-D land surface
model accounting for both the vertical and horizontal-2
SHF (Wm )01

4/6 4/11 4/16 4/21 4/26 5/1

Latent heat flux - daily average

200-2
LHF (Wm )2

Sensible heat flux - daily average

-1

0.5
0.4
0.3
0.2
0.1
0.0
4/1

Precipitation (mmh )Precipitation (mmh )3 -3
vwc (m m )3 -3
vwc (m m )Soil moisture - C1 scenario
0.5
0.4
0.3
0.2
0.1
0.0

200
150
100
50
0
-50
4/1

4/6 4/11 4/16 4/21 4/26 5/1 5/6 5/11
top - C2
bottom - C2
top - C1
bottom - C1

Fig. 4 Daily averages of latent (top chart) andsensible (bottom chart) heat fluxes. The lines
have the same meaning as in Figure 2.
hydrological processes. The numerical experiments
showed that:(1) The lateral flow in the quasi 3-D model causes a
profound spatial variability in the surface wetness,
temperature and heat fluxes. Consequently, theheterogeneity of the surface characteristics may
enhance the convective activity in an atmosphericmodel.
(2) The bulk response of the catchment to theatmosphere is significantly affected by the treatmentof surface runoff in the quasi 3-D model, which, incontrast to the 1-D SiB2, assures conservation ofwater balance within the catchment.
(3) The thawing/freezing of soil in hilly regions isstrongly influenced by the lateral flow. The spatially
non-uniform evolution of the active layer, together
with the lateral flow cause a considerable spatial
heterogeneity of surface wetness and surface fluxes.
These results are qualitatively consistent withobservations suggesting a realistic performance of thepresented model. However, it remains to be validated
quantitatively against reliable data before couplingwith an atmospheric model.

References
Dutta, D., S. Herath, and K. Musiake, Flood inundation
stimulation in a river basin using a physically based
distributed hydrologic model, Hydrological Processes,
14, 497-519, 2000.
Ishidaira H., T. Koike, M. Lu, and N. Hirose, Developmentof a 2-D soil model for heat and water transfer inpermafrost regions, Annual J. of Hydraulic Eng., JSCE,
42, 133-138, 1998, in Japanese.
Koike, T., and P. Koudelova: Introducing the lateral surface
flow into a land surface scheme, Proc. of 2003 IUGG
conference, Sapporo, A116, 2003.
Li, S-X., and T. Koike, Frozen soil parameterization in SiB2
and its validation with GAME-Tibet observations, Cold
Region Science and Technology, 36, 165-182, 2003.
Sellers, P. J., D. A. Randall, G. J. Collatz, J. A. Berry, C. B.
Field, D. A. Dazlich, C. Zhang, G. D. Collelo, and L.
Bounoua, A revised land surface parameterization (SiB2)
for atmospheric GCMs. Part I: Model Formulation, J.
Climate, 9, 676-705, 1996.