HYDROLOGIC CYCLE AND HEAT TRANSPORT MODELINGAND ITS APPLICATION TO THE KANDA RIVER WATERSHED, TOKYO
By
Tsuyoshi Kinouchi
Department of Environment System Management, Fukushima University, Fukushima-Shi, Fukushima, Japan
and
Yangwen Jia
Department of Water Resources, Institute of Water Resources and Hydropower Research, Beijing, China

SYNOPSIS
A physically based distributed hydrological model (WEP model) was coupled with a newly developed model of heattransport in river water to simulate watershed-scale hydrological and heat transport processes and the resulting effects on
the stream environment.

The model is capable of simulating runoff from a highly urbanized watershed with a drainage

system of river channels and a sewer network, which enables appropriate prediction of river floods and inflows towastewater treatment plants. Such a model is designed to investigate factors affecting spatial and temporal variations of
the stream temperature by incorporating anthropogenic heat impacts and urban canopy processes.

By applying the modelto the Kanda River Watershed in Tokyo, Japan, stream flow and temperature as well as inflows to a WWTP were simulated
well.

The anthropogenic heat impact on stream temperature was found to be very significant in winter.
INTRODUCTION
Streams are thought to have important natural functions such as natural landscaping, conserving habitats, andcontrolling the micro-climate in urban environments.

Such functions can be effective when appropriate stream

temperatures and other quality levels are maintained. However, urban streams are vulnerable to thermal pollution due tovarious direct and indirect anthropogenic impacts such as the effluent from wastewater treatment plants, destruction of
riparian vegetation, river channelization, the runoff from paved surfaces and global warming.

Increased stream

temperature can decrease the environmental importance of urban streams and various ecological impacts have been
observed such as the increase in invasive species. Because of increasing urban water and energy waste due to population
growth and high energy consumption life-styles, the temperature and magnitude of urban wastewater has increased
significantly (Kinouchi, 1), which may have resulted in the long-term increase in stream temperature.

Kinouchi et al. (2)indicated that the stream temperature in winter and early spring in the central Tokyo increased at a rate of 0.11 –0.21 °C/year in river segments that had a considerable increase in wastewater heat input.
To conserve environmental values of urban streams, a framework needs to be provided for predicting stream
temperatures under various natural and anthropogenic controlling factors. Although there are many physically based
in-stream models that are capable of simulating the shading effects (Sinokrot and Stefan, 3), the impact of urbanization(LeBlanc et al., 4) and global warming (Tung et al., 5), few studies have modelled and validated stream temperature
variations considering hydrological processes in a watershed scale (Chen et al., 6; Haag and Luce, 7). However, thepredictive framework of urban stream temperature needs to include hydrological modeling as the stream temperature isvery dependent on the hydrological behavior of watersheds as well as on the riparian conditions and anthropogenic
heating/cooling factors. In addition, characteristic radiation environment in channels run through highly urbanized areas
requires special modeling of air-water heat fluxes. For this purpose, a physically based distributed hydrological model
developed by Jia et al. (8~10) was coupled with a newly developed model of heat transport in river water, and theperformance of the model was demonstrated by applying it to a highly urbanized watershed in Tokyo, Japan, whereanthropogenic impacts are significant in modifying natural hydrologic cycle and the water temperature.
MODEL
Modeling hydrological system with dual drainage networks
The simulation of hydrological processes in an urbanized watershed with dual drainage system is based on the WEPmodel, a physically based distributed hydrological model developed by Jia et al. (8~10). The main features of this model
include coupling simulation of hydrological processes and energy transfer processes, simulation of interactive impact of
surface water and groundwater, and the mosaic structure within a grid cell. The simulated hydrological processes include
evapotranspiration, infiltration, surface runoff, subsurface runoff, groundwater flow, overland flow, and river flow. The
simulated energy transfer processes include short-wave radiation, long-wave radiation, latent heat flux, sensible heat flux,
and soil heat flux.
In urbanized watersheds equipped with sewer networks, storm runoff flows down through gutters, lateral pipes andsewer pipes during periods of rainfall.

In a combined sewer system, the stormwater mixed with waste from houses andbuildings eventually flows into a wastewater treatment plant (WWTP).

Untreated wastewater in the sewer system

overflows to a river system through combined sewer overflow (CSO) outlets following heavy rainfalls.

In no-rainfall

periods, the waste water from houses and buildings flows into the WWTP and the treated wastewater is discharged into
streams.

To take these characteristics of dual drainage networks into account in the WEP model, flows in sewer pipes arerouted as an overland flow in underground spaces, instead of directly solving unsteady flows in a sewer network, to reduce
modeling complexity.

Overland flow routing is carried out by means of the kinematic wave method with equivalentroughness parameters deduced from land use data, while a single value is given for the flow routing in a sewer network.
The magnitude of CSO depends on the structure and dimension of the CSO outlet and the water level within the CSO
outlet. As this model does not predict water levels in the CSO outlet, we alternatively modeled flows at the CSO outlet bysetting the maximum capacity of stormwater inflows to the WWTP as follows:
j

j

j

j

Qof = Qin − Qlm ， Qlm = Vlm CA j / CA(1)

where Qof j = the overflow to the river from the CSO outlet j, Qin j = the inflow to a computational mesh where the CSO
outlet j is located, Qlm j = the outflow to a downstream computational mesh from the mesh where the CSO outlet j islocated, Vlm = the maximum inflow capacity at the WWTP, CA = upstream catchment area of the WWTP, CAj = thecatchment area at CSO outlet j. Vlm is given by Vlm = βVTP, where VTP is the treatment capacity of the WWTP, β is aconstant. In practice, common sewer design allows inflow three times larger than the designed magnitude, whichmeans β = 3. However, β could be different from the design setting in reality, thus β is a calibration factor in this study.
Sewerage from houses and buildings during periods without rainfall is routed to compute inflows to the WWTP.
The flow routing in a sewer network considers the seepage of groundwater and soil moisture into sewer pipes. The soil
moisture seepage is modeled so that it is proportional to the unsaturated conductivity in the ground, and the seepage ofthe groundwater depends on the level difference of groundwater table and the water level in a pipe.
Modeling heat transport in streams
A basic equation for heat transport in an open channel is written as

Wh
Jh
Js
∂T
∂( AT ) ∂ (QT ) 1 ∂++=+( AD ) +ρ w c pwV ρ w c pw H ρ w c pw H
∂x
A∂t
A∂x
A ∂x

(2)

where T = the water temperature [°C]; A = the cross-sectional area [m2]; Q = the flow rate [m3/s]; D = the longitudinal
dispersion coefficient; Wh = the net heat load from surface flow, subsurface flow, groundwater flow, treated wastewater andcombined sewer overflow [J/s]; Jh = the air-water heat flux [W/m2]; Js = the heat flux from the stream bed to the stream
water [W/m2]; V = the water volume of stream element [m3]; H = the water depth [m]; ρw = the density of water [kg/m3];
cpw = the specific heat of water [J/kg/°C]; x = the flow direction [m]; t = the time [s]. The net heat load Wh from sources iscomputed as

Wh = ρ w c pw (∑ Q ps T ps + ∑ Qnps Tnps )(3)

where Qps and Tps are the magnitude and temperature of point source (wastewater effluents from WWTPs and CSO),
respectively; Qnps and Tnps represent the magnitude and temperature of non-point source (surface flow, subsurface flow andgroundwater flow), respectively.

Thermal runoff may cause a rise of stream temperature to a lethal level during rainfall

periods in hot seasons (Herb et al, 11). However, since we mainly focus on stream temperature variations during the low
flow periods, the temperature of the surface runoff is assumed to be the same as the air temperature. The observed
volume and temperature of the wastewater effluents were used, though it can be computed from Kinouchi (1). The
temperature of subsurface and groundwater flows is assumed to be the same as the ground temperature at a specified depth.
The ground temperature Tgd(z, t) at depth z is computed from equation (4) (Carslaw and Jaeger, 12) with the ground surface
temperature given by equation (5).

∑ A exp(− z nω 2a )sin (nωt − z nω 2a + ε )6

T gd ( z , t ) = A0 +

n

n =1

n

(4)

6

T gd (0, t ) = a 0 +∑ [a cos(nωt ) + b sin(nωt )]
n

n

(5)

n =1

2

An = an + bn2

(6)

where sin(εn) = an/An; A0 = a0; ω = 2π /τ; τ = the period; a = the thermal diffusivity; an and bn are the Fourier coefficients.
The air-water heat flux Jh is computed as**

J h = RN − J c − J e = I s + J s − J c − J e(7)where RN = the net radiation at the water surface; Jc = the sensible heat flux; Je = the latent heat flux; Is* = the net solar
shortwave radiation at the water surface; Js* = the net atmospheric longwave radiation at the water surface. The net
radiation RN is associated with many factors such as the solar zenith angle, the river transversal and cross-sectional
dimensions, riparian vegetation and buildings. Here, we take such factors into account, except for the shading effect ofbuildings.
The net solar shortwave radiation Is*, Iw* and the net atmospheric longwave radiation Js*, Jw* are given as

I s = S DIR ( w − L)(1 − α s ) / w + S DIFψ s (1 − α s ) + S DIR Lα w (1 − ψ s )(1 − α s ) /(2h)*+ S DIFψ w (1 − ψ s )α w (1 − α s )I w = S DIR L(1 − α w ) /( 2h) + S DIFψ w (1 − α w ) + S DIR ( w − L)ψ wα s (1 − α w ) / w + S DIFψ w (1 − α w )*+ S DIR L(1 − 2ψ w )α w (1 − α w ) /( 2h) + S DIFψ w (1 − 2ψ w )α w (1 − α w )J s = ψ sε s L ↓ +(1 − ψ s )ε sε wσTw4 − ε sσT 4
*

(9)

(10)

J w = ψ wε w L ↓ +ψ wε wε sσT 4 + (1 − 2ψ w )ε w σTw4 − ε wσTw4
*

(8)

2

(11)

where subscripts s and w denote the water surface and the river side wall, respectively, SDIR and SDIF are downward direct
and diffuse solar radiations on a horizontal surface at the top of a stream channel canopy, L ↓ = the downward
atmospheric longwave radiation; ψs = the sky-view factor for the water surface; ψw = the sky-view factor for one side wall;
Tw = the surface temperature of the river side wall; h = the channel depth; w = the channel width; L = the shadow length; εs
and εw are the emissivity of the water surface and the river side wall, respectively, αs and αw are the albedo of the water
surface and the river side wall, respectively, σ = the Stefan–Boltzmann constant.
The sky-view factors of the water surface ψs and the river wall ψw are given from Masson (13) as[]ψ s = (h w)2 + 1 − h w

1
2

1/2

ψ w = {h w + 1 − [(h w )2 + 1]1 / 2 } / (h w )

(12)

(13)

The shadow length L is computed as (Kusaka et al., 14)⎧h tan θ z sin θ nL=⎨⎩w(L < w)(L ≥ w)(14)

where θz = the solar zenith angle; θn = the difference between the solar azimuth angle and the stream channel orientation.
The riparian vegetation and structures crossing the stream modify the shortwave and longwave radiative fluxes. We
modelled the attenuation of the solar radiation by the riparian vegetation as follows:

S DIR = f s SDIR

(15)

S DIF = f s SDIF

(16)

f s = 1 − λI veg f LAI(17)

f LAI = ( LAI − Lmin ) ( Lmax − Lmin )(18)

where fs = an attenuation factor for the riparian vegetation, SDIR = direct solar radiation at the above-canopy level; SDIF =diffuse solar radiation at the above-canopy level; Iveg = a riparian vegetation index (0 : no vegetation, 1 : vegetation along
one side, 2 : vegetation along both sides); LAI = the leaf area index of the riparian vegetation; Lmax = maximum value of
LAI; Lmin = minimum value of LAI; λ = a calibrated constant to reflect the shading effect by riparian vegetation.
The sensible heat flux is expressed as (Chapra et al., 15)

J c = c1 f (u )(T − Ta )(19)

where c1 = 0.47 ; f (u ) = 19 + 0.95u 2 ; Ta = the air temperature (°C); u = the wind speed at a height of 7 m (m/s). The
incoming radiations at the above-canopy level (SDIR, SDIF and L↓) and the latent heat flux in equation (7) are computed
by Jia et al. (8).
The heat flux from the river bed to the stream water (Js) is calculated from equations (20) through (22).

J s = ρ sed c ps

2α sed(Tsed − T )
H sed

(20)

dTsed
= k h (T − Tsed )
dt

(21)

k h = 2α sed H sed

(22)

2

where Tsed = the temperature of the river bed layer; Hsed = the effective thickness of the river bed layer; ρsed = the density of
the river bed layer; cps = the specific heat of the river bed layer; αsed = the thermal diffusivity of the river bed layer.

The surface temperature of the river side wall (Tw) is calculated by a similar formulation using equations (23) through
(27).
*

*

J w + I w = H C + GC

(23)

dTg

= k hg (Tw − Tg )(24)

H C = c1 f (u )(Tw − Ta )

(25)

dt

GC = ρ g c pg

2α g
Hg

k hg = 2α g H g(26)(Tw − Tg )

2

(27)

where HC = the sensible heat flux from the river side wall; GC = the heat flux from the river side wall into the ground; Hg =the thickness of the river side wall; Tg = the internal temperature of the side wall; ρg = the density of the side wall; cpg = thespecific heat of the side wall; αg = the thermal diffusivity of the side wall.
A discrete form of equation (2) is written by equation (28) after Chapra et al. (15). Equation (28) is implicitly
solved by means of the Gaussian elimination method.
D′D′
Q
dTi Qi −1=
Ti −1 − i Ti + i −1 (Ti −1 − Ti ) + i (Ti +1 − Ti )
Vi
Vi
Vi
dt
Vi
−

Wh ,i

J h ,i(28)

J s ,i

1 dAi
Ti +
+
+
Ai dt
ρ w c pwVi ρ w c pw H i ρ w c pw H i

Di′ = Di ⋅ Ai /(Δxi + Δxi +1 ) / 2(29)

where subscript i denotes river element i, Δx = the length of river element i.
APPLICATION
Study area
We applied the hydrological model coupled with the heat energy transport model to the Kanda River watershedlocated in the heart of Tokyo, Japan (Figure 1).

Tributaries run in four open channels and one covered channel in

upstream sub-watersheds, and eventually merges to the mainstream.
the stream (Figure 1).

The Momozono River is a covered channel all along

The channels have a trapezoidal cross-sectional shape in most of their reaches, some of which havea small ditch on the bottom of the channel. The Kanda River watershed has a combined sewer system and the watershed
can be divided into five sewer catchments. Sewage generated in the catchment bounded by a dashed line (Figure 1) flows

Figure 1. Outline of the Kanda Watershed
into two WWTPs within the watershed (the Ochiai WWTP and the Nakano WWTP), where treated effluents are discharged
into the Kanda River. In all other catchments, the sewage is transported to WWTPs located outside the Kanda River
watershed through the sewer pipes across the watershed boundary.
rainfall with relatively light intensity.

Stormwater flows in the same route during periods ofDuring periods of heavy rainfall, the CSO in all sewer catchments directly

discharges into the Kanda River system from the CSO outlets located within the watershed.
Input data and parameters
The main input data required to run the model are listed in Table 1.
prepared using ArcGIS.

Spatial distribution of correspondent data was

Stream network data was reproduced on the GIS platform from the digital maps by the

Geographical Survey Institute.

The watershed was divided into 9 sub-watersheds.

River cross-sectional area wasapproximated as a trapezoidal shape. The shape data was provided by the maps obtained from the Tokyo Metropolitan
Government. The orientation of the stream channel was calculated for each river element by utilizing a tool of ArcGIS.
Categorized data of riparian vegetation and structures crossing channels, which were prepared from aerial photographs, areattributed to individual computational element of river channels. The actual sewer network data was not used. Instead, weTable 1. List of main input data used for the coupled model simulation
Category

Item

Land surface conditions

Watershed boundary
Topography
Land use

Meteorological and hydrologicalconditions

Soil type
Precipitation (19 stations)
Relative humidity (1 station)

River network and hydraulicdimensions

Air temperature (3 stations)Wind velocity (3 stations)
Sunshine duration (3 stations)Stream channel
River profiles

Anthropogenic conditions

Riparian vegetation
Crossing structure
Population
Sewer catchment boundary

SourceNational digital information*
Digital map* (50m by 50m DEM)
Fine digital information* (13 types ofland use, 100m by 100m raster data)
Geologic map
Hourly data from the Tokyo
Metropolitan Government
Hourly data from AMeDAS by theJapan Meteorological Agency
Digital Map (1:25,000)
Maps obtained from Tokyo
Metropolitan Government
Aerial photographs
Aerial photographs
National census data by the Ministry ofInternal Affairs and Communications
Maps obtained from the Tokyo
Metropolitan Government*Products of the Geographic Survey Institute
Table 2. Parameters used for the heat transport model
parameters

River bed

River side wall

Thickness Hsed, Hg (m)Density ρsed, ρg (g/cm3)
Specific heat cps, cpg (cal/g/K)
Thermal diffusivity αsed, αg (cm2/s)
Albedo αs, αw
Emissivity εs, εw

0.1
2.2
0.21
0.008
0.08*
0.95*

0.3
2.4
0.21
0.02
0.3
0.90

*Values for water surface

assumed that the flows in the sewer pipes followed the terrain of the watershed, which was found to be a reasonableassumption.
There are three categories of parameters used in the WEP model (Jia et al., 8~10): (1) parameters of land surface andriver channel system; (2) parameters of vegetation; and (3) parameters of soil and aquifer. Parameters for the heattransport model are additionally specified following Chapra et al. (15) and Oke (16) as listed in Table 2.

All ofparameters were initially specified from literatures or estimated according to land cover information, observation data, andsome parameters are selected for model calibration.

The model calibration was performed for both stream flows andtemperatures at locations where observed values are obtained. For flow simulations, some flood hydrographs and low flows
were checked to calibrate parameters such as the Manning’s roughness of the river flow, β in equation (1) and the hydraulicconductivity of the riverbed. Parameter β used to give Vlm in equation (1) was set to 5 following the calibration. Calibrationof parameters for the stream temperature including λ in equation (17) was carried out to reproduce its seasonal and diurnal
variations.

Ratios of the impervious area, four types of pervious areas and water surfaces for each computational mesh

(100m by 100m) were given by reclassifying the fine land use data using the parameters in Table 3, which were deduced
referring to previous studies (Jia et al., 9 and 10).
The spatial distribution data of the topsoil was generated from geological maps.

The topsoil in the watershed isTable 3. Land use classification parameters
Land use classification used in the WEP model
Water
Land use type ofthe fine digital surface ratio
information

Paddy ratio
in PA

Pervious area (PA)
Short
Tall
vegetation
vegetation
ratio in PA ratio in PA

Bare soil
ratio in PA

Impervious area (IA)Urbancanopy ratein IA

forest

0

0

0.8

0.2

0

0.27

0

paddy

0

1

0

0

0

0

0

other farmland

0

0

0.2

0.8

0

0.21

0

under preparing

0

0

0.2

0.3

0.5

0.23

0

vacant land

0

0

0.2

0.3

0.5

0.27

0

industrial use

0

0

0.2

0.3

0.5

0.75

0.36

housing area

0

0

0.2

0.3

0.5

0.75

0.36

housing area

0

0

0.2

0.3

0.5

0.75

0.36

commercial area

0

0

0.2

0.3

0.5

0.80

0.36

road

0

0

0

1

0

0.86

0

park

0

0

0.3

0.5

0.2

0.16

0.36

public works

0

0

0.2

0.3

0.5

0.45

0.36

river/lake

1

0

0

0

0

0.14

0

mainly the Kanto loam, except for the alluvial soil in riparian zones (Figure 1). The physical properties for each soil typesuch as soil porosity, soil saturated hydraulic conductivity and soil moisture–suction relation curve were given from Jia etal. (9 and 10). Values of Manning’s roughness coefficient for river flow and overland flow were set to 0.015 and 0.01,
respectively. The air-water heat flux in the Momozono river, a fully covered channel, was computed by setting the direct
and diffuse solar radiations as zero and the downward longwave radiation from the cover was assumed to be the same as

L ↓.
Model application
The coupled model was run for the simulation from January 1 to December 31 of 2005 in a time step of 1 hour(except for overland flow and river flow routing with a time step of 10 minutes) and a grid cell size of 100m×100m.
computational domain includes 9062 grid cells.
longitudinal length of 390m - 450m.

The

River channels were divided into computational elements with the

The headwater of the stream was maintained by the pumped groundwater. Thus,

the upstream boundary conditions were set by the constant stream flow and variable water temperatures deduced from theground temperature.

The temperature of the base flow from the groundwater was likewise set at an equivalenttemperature of the ground at the depth of the river bed.
Results
Streamflow and inflow to WWTPs
The simulated results of low streamflows were compared with those observed at Kotobuki and Chitose (Figure 1,

Table 4. Comparisons of simulated stream flow with observed during no-rainfall periods
Location

PeriodObserved average flow

Simulated average flow

3

Kotobuki

2005/1/20-21

0.52 m /s

0.55 m3/s

Chitose2005/1/18-19

0.07 m3/s

0.10 m3/s

70
Observed
Simulated

30
20
10

7/10
0:00

7/10
3:00

20

7/9
18:00

8/16
6:00

90
80
70
60
50
40
30
20
10
0
8/15
18:00

Observed
Simulated

40

Flow rate (m3/s)

Flow rate (m3/s)

30

7/10
6:00

30
20
10

8/15
21:00

8/16
0:00

8/16
3:00

40

7/9
21:00

7/10
0:00

7/10
3:00

7/10
6:00

Observed
Simulated

8/15
21:00

8/16
0:00

8/16
3:00

8/16
6:00

60

Observed
Simulated

30
20
10

Observed
Simulated

50
Flow rate (m3/s)

Flow rate (m3/s)

40

0
7/9
21:00

50

0
8/15
18:00

50

10

0
7/9
18:00

Observed
Simulated

60
3

40

Flow rate (m /s)

Flow rate (m3/s)

50

40
30
20
10

0
8/25
6:00

8/25
12:00

8/25
18:00

8/26
0:00

0
8/25
6:00

8/26
6:00

200
150
100
50

Flow rate (m3/s)

Observed
Simulated

250
Flow rate (m3/s)

8/25
18:00

8/26
0:00

8/26
6:00

400

300

Observed
Simulated

350
300
250
200
150
100
50
0

0
9/4
12:00

8/25
12:00

9/4
18:00

9/5
0:00

9/5
6:00

9/5
12:00

9/4
12:00

9/4
18:00

9/5
0:00

9/5
6:00

9/5
12:00

Figure 2. Comparison of simulated flood discharges with observed (left: Honan, right: Asahi)
Table 4). Most of the river reaches are covered by a concrete bottom slab and side walls, which prevent water interaction
between river flows and groundwater. As a result, the streamflow could be quite low during periods without rainfall.
This situation was reproduced well by the simulation.

Simulated flood discharges during four major storm events in 2005

were compared with hydrographs at Asahi and Honan obtained from observed rating curves and measured water levels.
Simulated hydrographs relatively well reproduced observed ones at Asahi for three events (July 9-10, August 15-16 andInflow (m3/s)

20

25

Observed
Simulated

15
10
5

7/9
6:00

7/9
18:00

7/10
6:00

7/10
18:00

7/11
6:00

Observed
Simulated

0
8/12
12:00
25

8/13
12:00

8/14
12:00

8/15
12:00

8/16
12:00

Observed
Simulated

20
15

15

10

10

5

5
0
8/24
12:00

10

Inflow (m3/s)

Inflow (m3/s)

20

15

5

0
7/8
18:00
25

Observed
Simulated20
Inflow (m3/s)

25

0
8/25
0:00

8/25
12:00

8/26
0:00

8/26
12:00

8/27
0:00

9/3
0:00

9/4
0:00

9/5
0:00

9/6
0:00

9/7
0:00

9/8
0:00

Figure 3. Comparison of simulated inflows to the Ochiai WWTP during four rainfall events

August 25-26) and at Honan for two events (July 9-10 and August 25-26) as shown in Figure 2. Most abrupt changes ofdischarge were observed during the rainfall event of August 25-26, in which the water level rose about 1.5m within 40
minutes. Similar features are simulated for this event even though hourly rainfall data are used. Simulated hydrographs
for rainfall periods from September 4 to 5, which was caused by the approach of the typhoon No.14, overestimated thepeak discharge and the total volume at both locations. During this period, flooding occurred along the stream channels
and the runoff was temporally stored in detention reservoirs located upstream of the watershed, while the model does notdeal with these processes. Figure 3 shows the correspondent inflows to the Ochiai WWTP. Diurnal variations andabrupt changes of the inflow were simulated well, although the simulated base inflows underestimated measured ones
probably due to the setting of the seepage modeling. Our observations showed that the seepage accounts for about 20% of
the total inflow to the WWTPs in Tokyo (Nakayama et al., 17).

Therefore, further verifications may be needed for themodeling of soil-water seepage from the ground surrounding the sewer pipes.
Water temperature
The simulated stream temperatures are compared with those observed (Figure 4). Modeled stream temperatures
matched well with those observed at Goryo throughout the year. Modeled stream temperatures at Manki show larger
amplitudes of the diurnal variation than those observed in the fall and winter seasons. Akebono, where we started
temperature measurement from November 2005, exhibited flatter temperature variations than those observed at Goryo andManki. This is mainly because the wastewater effluents from two WWTPs comprise the main portion of the streamflow
at Akebono and thus the diurnal temperature variation is dominated by the effluent temperature that has small diurnal
amplitudes. In addition, our monitoring detected sudden descends of the stream temperature at Akebono, followed by
the soon recovery to the previous temperature level. For example, the phenomenon occurred in the morning of the 23rdof November, probably due to the temporal reduction in wastewater effluents to 0.8m3/s from its daily average of 2.1m3/s.
These abrupt temperature changes at Akebono were also well captured by the simulation (Figure 4).

Stream20
15
10

Observed

5

Simulated

Stream temperature (°C)

0
3/25
35

30

Stream temperature (°C)

Goryo

25

Stream temperature (°C)

Stream temperature (°C)

30

30

3/30

4/4

4/9

4/14

4/19

4/24

4/29

5/4

5/9

5/14

5/19

5/24

9/10

9/15

9/20

9/25

Goryo

25
20
15

Observed

10

Simulated

5
7/27

25

8/1

8/6

8/11

8/16

8/21

8/26

8/31

9/5

Manki

20
15

Observed

10

Simulated

5
0
9/27

10/2

10/7

10/12 10/17 10/22 10/27

11/1

11/6

11/11 11/16 11/21 11/26

12/1

30
25

AkebonoObserved

Simulated

Simulated (No wastewater effluent)

20
15
10
5
0
11/21

11/22

11/23

11/24

11/25

11/26

11/27

11/28

11/29

11/30

12/1

Figure 4. Comparison of simulated stream temperatures with observed
The bottom figure includes results for the hypothetical case of no wastewater effluent.
temperatures simulated by assuming no wastewater inputs were far below the temperatures simulated with wastewater
inputs (Figure 4, bottom), and the maximum difference of stream temperature was found to be approximately 12°C. Thissuggests that there is a strong influence of wastewater effluents on the stream temperature in the Kanda River, andpotential impacts on the habitat for fish and invertebrates and algae growth.
CONCLUSIONS
A framework was presented for the stream temperature prediction under distributed anthropogenic and natural heatsources in the urbanized watershed. The predictive model developed in this study is based on a model of heat transport inriver water incorporated into a physically based distributed hydrological model (WEP model). This coupled model cansimulate dynamic hydrological responses of an urbanized watershed with a drainage system of river channels and a
combined sewer network, which enables appropriate prediction of inflows to wastewater treatment plants and streamflows.
Stream temperature prediction considers the influence of heat input by wastewater effluents, groundwater seepage, surface
runoff as well as the atmospheric radiations, and the heat budget within a river canopy. This model was applied to the
Kanda River watershed in Tokyo, and showed that it could reproduce well low flows and flood runoff, including inflows tothe WWTP. Overall magnitudes of the simulated stream temperature agreed well with those of observed. Impacts ofanthropogenic heat input from WWTPs were found to be very significant in winter.

Abrupt changes of stream

temperature in a time scale of several hours were attributed to the fluctuation of the volume of wastewater effluents.
Further improvements are necessary to more accurately simulate observed seasonal and diurnal variations in the stream
temperature throughout the year.
ACKNOWLEDGMENTS
This research project was financially supported by CREST (Core Research for Evolutional Science and Technology)of the Japan Science and Technology Corporation as a part of the research project ‘Water and Energy Forcing due toUrbanization in Land-Atmosphere-Coastal System’ organized by Prof. Manabu Kanda (Tokyo Institute of Technology).
REFERENCES
1. Kinouchi, T.: Impact of long-term water and energy consumption in Tokyo on wastewater effluent: implication forthermal degradation of urban streams, Hydrological Processes, Vol.21, pp.1207 - 1216, 2007.
2. Kinouchi, T., Yagi, H. and Miyamoto, M.: Increase in stream temperature related to anthropogenic heat input fromurban wastewater, Journal of Hydrology, Vol.335, pp.78-88, 2007.
3. Sinokrot, B.A. and Stefan, H.G.: Stream temperature dynamics: measurement and modeling, Water Resources Research,
Vol.29, No.7, pp.2299–2312, 1993.
4. LeBlanc, R.T., Brown, R.D. and FitzGibbon, J.E.: Modeling the effects of land use change on the water temperature in
unregulated urban streams, Journal of Environmental Management, Vol.49, pp.445–469, 1997.
5. Tung C., Lee T. and Yang Y.: Modelling climate change impacts on stream temperature of Formosan landlockedsalmon habitat, Hydrological Processes, Vol.20, pp.1629–1649, 2006.
6. Chen, Y.D., McCutcheon, S.C., Norton, D.J. and Nutter, W.L., Stream temperature simulation of forested riparian
areas: II. Model application, Journal of Environmental Engineering, ASCE, Vol.124, No.4, pp.316–328, 1998.
7. Haag, I. and Luce, A., The integrated water balance and water temperature model LARSIM-WT, Hydrological
Processes, Vol.22, pp.1046-1056, 2008.
8. Jia, Y., Ni, G., Kawahara, Y. and Suetsugi, T.: Development of WEP model and its application to an urban watershed,
Hydrological Processes, Vol.15, No.11, pp.2175-2194, 2001.
9. Jia, Y., Ni, G., Yoshitani, J., Kawahara, Y. and Kinouchi, T.: Coupling simulation of water and energy budgets andanalysis of urban development impact, Jour. of Hydrologic Eng., ASCE, Vol.7, No.4, pp.302-311, 2002.
10. Jia, Y., Kinouchi, T. and Yoshitani, J.: Distributed hydrologic modeling in a partially urbanized agricultural watershedusing water and energy transfer process model. Jour. of Hydrologic Eng., ASCE, Vol. 10. No. 4, pp.253-263, 2005.
11. Herb, W. R., Janke, B., Mohseni, O. and Stefan, H. G.: Thermal pollution of streams by runoff from paved surfaces,
Hydrological Process, Vol.22, pp.987 – 999, 2008.
12. Carslaw, H. S. and Jaeger, J.C.: Conduction of heat in solids 2nd edition, Oxford Univ. Press, London, 1959.
13. Masson, V.: A physically-based scheme for the urban energy budget in atmospheric models, Boundary-Layer
Meteorology, Vol.94, No.3, pp.357-397, 2000.
14. Kusaka, H., Kondo, H., Kikegawa, Y. and Kimura, F.: A simple single-layer urban canopy model for atmosphericmodels: comparison with multi-layer and slab models, Boundary-Layer Meteorology, Vol.101, No.3, pp.329-358, 2001.

15. Chapra, S.C., Pelletier, G.J. and Tao, H.: QUAL2K: A modeling framework for simulating river and stream water
quality, version 2.04: documentation and users manual, Civil and Env. Eng. Dept., Tufts University, Medford, MA.,2006.
16. Oke, T.: Boundary Layer Climates 2nd edition, Routledge, New York. 1987.
17. Nakayama, Y., Kanda, M. and Kinouchi, T.: The investigation of the water and heat transfer in urban sewage system
based on the temperature observation at the wastewater treatment plants, Journal of Japan Society of Hydrology andWater Resources, Vol.20, No.1, pp.25-33, 2007. (In Japanese)

APPENDIX − NOTATION

The following symbols are used in this paper.
a= thermal diffusivity;

a n, b n= the Fourier coefficients;

cpg= specific heat of the side wall;cps= specific heat of the river bed layer;

cpw= specific heat of water;fs= an attenuation factor for the riparian vegetation;h= channel depth;t= time;u= wind speed at a height of 7 m;w

= channel width;x= flow direction;z= depth from the ground surface;A= cross-sectional area;CA= upstream catchment area of the WWTP;

CAj= catchment area at CSO outlet j;D= longitudinal dispersion coefficient;GC= heat flux from the river side wall into the ground;H= water depth;HC= sensible heat flux from the river side wall;Hg= thickness of the river side wall;

Hsed= effective thickness of the river bed layer;Is*= net solar shortwave radiation at the water surface;

Iw*= net solar shortwave radiation at the river side wall;Jc= sensible heat flux;Je= latent heat flux;Jh= air-water heat flux;Js= heat flux from the stream bed to the stream water;Js*= net atmospheric longwave radiation at the water surface;

Jw*= net atmospheric longwave radiation at the river side wall;Iveg= riparian vegetation index (0 : no vegetation, 1 : vegetation along one side, 2 : vegetation along both sides);L= shadow length;

LAI= leaf area index of the riparian vegetation;

Lmax= maximum value of LAI;

Lmin

= minimum value of LAI;

L↓

= downward atmospheric longwave radiation;Q= flow rate;

Qin j= inflow to a computational mesh where the CSO outlet j is located;Qlm j= outflow to a downstream computational mesh from the mesh where the CSO outlet j is located;

Qnps= magnitude of non-point source (surface flow, subsurface flow and groundwater flow);

Qof j= overflow to the river from the CSO outlet j;

Qps= magnitude of point source (wastewater effluents from WWTPs and CSO);

RN= net radiation at the water surface;

SDIF= diffuse solar radiation at the above-canopy level;

SDIF= downward diffuse solar radiation on a horizontal surface at the top of a stream channel canopy;SDIR= direct solar radiation at the above-canopy level;

SDIR= downward direct solar radiation on a horizontal surface at the top of a stream channel canopy;T= water temperature;Ta= air temperature;Tg= internal temperature of the side wall;

Tgd(z, t) = ground temperature at depth z and time t;
Tnps= temperature of non-point source (surface flow, subsurface flow and groundwater flow);

Tps= temperature of point source (wastewater effluents from WWTPs and CSO);

Tsed= temperature of the river bed layer;Tw= surface temperature of the river side wall;V= water volume of stream element;Vlm= maximum inflow capacity at the WWTP;

VTP= treatment capacity of the WWTP;

Wh= net heat load from hydrological sources;

αg= thermal diffusivity of the side wall;

αs= albedo of the water surface;

αsed= thermal diffusivity of the river bed layer;

αw= albedo of the river side wall;β= a calibrated constant to compute Vlm;εs= emissivity of the water surface;εw= emissivity of the river side wall;

λ= a calibrated constant to reflect the shading effect by riparian vegetation;

θn= difference between the solar azimuth angle and the stream channel orientation;

θz= solar zenith angle;

ρg= density of the side wall;ρsed= density of the river bed layer;

ρw= density of water;σ= the Stefan–Boltzmann constant;ψs= sky-view factor for the water surface;ψw= sky-view factor for one side wall;Δx= length of a river element.