Description of MATSIRO6

MATSIRO6 document writing team∗April 13, 2021

∗ MATSIRO6 document writing team (in alphabetical order)

Authors: Qiang Guo, Kanon Kino, Shuping Li, Tomoko Nitta, and Akira Takeshima
Supervisors/Editors: Tomoko Nitta, Yukihiko Onuma, Yusuke Satoh, Tatsuo
Suzuki, Kumiko Takata, Natsuki Yoshida, and Kei Yoshimura

Remarks

How to cite
MATSIRO6 Document Writing Team (2021), Description of MATSIRO6, CCSR Report No. 66, Division of Climate System Research, Atmosphere and Ocean ResearchInstitute, The University of Tokyo, https://doi.org/10.15083/0002000181

Publisher Information
Division of Climate System Research,
Atmosphere and Ocean Research Institute,
The University of Tokyo
5-1-5, Kashiwanoha, Kashiwa-shi,
Chiba, 277-8568, Japan
Phone : 04-7136-4371 Fax : 04-7136-4375
Contact: Kei Yoshimura (kei@iis.u-tokyo.ac.jp)© 2021 MATSIRO6 Document Writing Team1

Acknowledgment
This document describes about the Sixth version of Minimal Advanced Treatments of Surface
Interaction and RunOﬀ (MATSIRO) model, which is implemented in MIROC6 climate model.
The original document was written in 2001 by Dr. Seita Emori in Japanese about the very ﬁrst
version of MATSIRO, which was implemented in CCSR/NIES AGCM5. There has been no
update since then, even though there are signiﬁcant updates in the model itself. Furthermore,thinking about Integrated Land Simulator or ILS, the new framework of land model, which
adopts MATSIRO as a part of physical model components, it is highly needed to make eﬃcient
ways to update the document by community eﬀort. These are the main motivations of making
this document.
To update the description, we organized MATSIRO6 document writing team in the summerof 2020. Furthermore, it was decided to utilize Github version management platform to maintainthe documentation by the team. The team consists of two groups; the authors and the supervisors/editors. All of the authors are doctor course students of Global Hydrology Group of the
University of Tokyo, who use or are interested in using MATSIRO6, and the supervisors/editors
are researchers, who contributed to development of MATSIRO6 and/or its previous versions.
The original Japanese document written by Dr. Emori is in LaTeX format and there wasEnglish translated document done by Dr. Hanasaki and Dr. Hijioka in MS Word format. It
was needed to make the document into Markdown format from these two ﬁles. This part was
mainly done by Ms. Kanon Kino. Then the team was divided into several groups, and eachgroup became in charge of each section. These groups and corresponding sections are as follows:• Mr. Akira Takeshima supervised by Dr. Yukihiko Onuma for Snow (Ch.8)• Mr. Shuping Li supervised by Dr. Natsuki Yoshida for Runoﬀ (Ch.9)• Ms. Kanon Kino supervised by Dr. Tatsuo Suzuki for Lake (Ch.11)• Mr. Qiang Guo supervised by Dr. Yusuke Satoh for Wetland (Ch.12)
Later, in addition to the sections above, it was realized that some sections needed moreamendments and help, so that following additional assignments were conducted.
• Dr. Tomoko Nitta for Prognostic Variables (Ch.1.2)• Mr. Shuping Li for Tile scheme (Ch.13)• Mr. Qiang Guo for Physical Formulation and Process in Lake (Ch.11.4)• Dr. Tomoko Nitta for Lake-river Coupling (Ch.11.5)
It took about ﬁve months to draft, and the whole draft was reviewed by all supervisors/editors. At last, the rest of the sections were amended and improved by all.
The team thanks the support of “Integrated Research Program for Advancing Climate Models (TOUGOU Program)” from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan. Finally, the team would sincerely express our respect to Dr. Seita
Emori, the ﬁrst developer of the model and the author of the original version of this document.
Kei Yoshimura
lead member of the teamMarch 31, 2021

2

Revision History
2021/3/31 First version completed.
2021/4/13 Add Remarks

3

Contents
1 Introduction
1.1 Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Prognostic variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Output data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Input data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 External parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
8
8
10
10
11

2 Vegetation type parameters

14

3 Radiation parameters
3.1 Calculation of ground surface (forest ﬂoor) albedo . . . . . . . . . . . . . . . . .
3.2 Calculation of canopy albedo and transmissivity . . . . . . . . . . . . . . . . . .
3.3 Calculation of surface radiation ﬂux, etc. . . . . . . . . . . . . . . . . . . . . . . .

14
14
15
18

4 Turbulence parameters (bulk coeﬃcient)
4.1 Calculation of roughness with respect to momentum and heat . . . . . . . . . . .
4.2 Calculation of bulk coeﬃcient with respect to momentum and heat . . . . . . . .
4.3 Calculation of bulk coeﬃcient with respect to vapor . . . . . . . . . . . . . . . .

18
18
20
21

5 Stomatal resistance
5.1 Calculation of soil moisture stress factor . . . . . . . . . . . . . . . . . . . . . . .
5.2 Calculation of amount of photosynthesis . . . . . . . . . . . . . . . . . . . . . . .
5.3 Calculation of stomatal resistance (2) . . . . . . . . . . . . . . . . . . . . . . . .
5.4 Calculation of ground surface evaporation resistance . . . . . . . . . . . . . . . .

21
21
22
24
25

6 Surface energy balance
6.1 Calculation of surface turbulent ﬂuxes . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Calculation of heat conduction ﬂuxes . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 Solution of energy balance at ground surface and canopy . . . . . . . . . . . . . .
6.3.1 Energy balance at ground surface and canopy . . . . . . . . . . . . . . . .
6.3.2 Case 1: When there is no melting at the ground surface . . . . . . . . . .
6.3.3 Case 2: When there is melting at the ground surface . . . . . . . . . . . .
6.3.4 Conditions for solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.5 Updating of ground surface and canopy temperatures . . . . . . . . . . .
6.3.6 Updating of ﬂux values . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25
25
28
28
28
29
30
30
31
31

7 Canopy Water Balance
7.1 Diagnosis of canopy water phase . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Prognosis of canopy water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 Evaporation (sublimation) of canopy water . . . . . . . . . . . . . . . . .
7.2.2 Interception of precipitation by the canopy . . . . . . . . . . . . . . . . .
7.2.3 Dripping of the canopy water . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.4 Updating and melting of canopy water . . . . . . . . . . . . . . . . . . . .
7.3 Fluxes given to the soil, snow, and runoﬀ process . . . . . . . . . . . . . . . . . .

32
32
33
33
33
34
35
35

4

8 Snow
8.1 Diagnosis of snow cover fraction . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.1.1 Case 1: When OPT SSNOWD is active . . . . . . . . . . . . . . . . . . .
8.1.2 Case 2: When OPT SSNOWD is inactive . . . . . . . . . . . . . . . . . .
8.2 Vertical division of snow layers . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3 Calculation of snow water equivalent . . . . . . . . . . . . . . . . . . . . . . . . .
8.3.1 Sublimation of snow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3.2 Snowmelt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3.3 Freeze of snowmelt water and rainfall in snow . . . . . . . . . . . . . . . .
8.3.4 Snowfall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3.5 Redivision of snow layer and rediagnosis of temperature . . . . . . . . . .
8.4 Calculation of snow heat conduction . . . . . . . . . . . . . . . . . . . . . . . . .
8.4.1 Snow heat conduction equations . . . . . . . . . . . . . . . . . . . . . . .
8.4.2 Case 1: When snowmelt does not occur in the uppermost layer . . . . . .
8.4.3 Case 2: When snowmelt occurs in the uppermost layer . . . . . . . . . . .
8.5 Fluxes given to the soil or the runoﬀ process . . . . . . . . . . . . . . . . . . . .
8.6 Glacier formation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.7 Dust in snow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.7.1 Dust fall on the snow cover . . . . . . . . . . . . . . . . . . . . . . . . . .
8.7.2 Redistribution of dust . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.8 Albedo of snow and ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.8.1 Albedo of snow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.8.2 Albedo of ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36
38
38
39
39
40
40
41
42
43
43
44
44
45
46
46
47
47
47
48
49
49
50

9 Runoﬀ
9.1 Outline of TOPMODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.2 Application of TOPMODEL assuming simpliﬁed topography . . . . . . . . . . .
9.3 Calculation of runoﬀ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.3.1 Estimation of grid average water table depth . . . . . . . . . . . . . . . .
9.3.2 Calculation of groundwater runoﬀ . . . . . . . . . . . . . . . . . . . . . .
9.3.3 Calculation of surface runoﬀ . . . . . . . . . . . . . . . . . . . . . . . . . .
9.4 Water ﬂux given to soil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50
53
54
55
56
56
57
58

10 Soil
10.1 Calculation of soil heat conduction . . . . . . . . . . . . . . . . . . . . . . . . . .
10.1.1 Soil heat conduction equations . . . . . . . . . . . . . . . . . . . . . . . .
10.1.2 Solution of heat conduction equations . . . . . . . . . . . . . . . . . . . .
10.2 Calculation of soil moisture movement . . . . . . . . . . . . . . . . . . . . . . . .
10.2.1 Soil moisture movement equations . . . . . . . . . . . . . . . . . . . . . .
10.2.2 Solution of soil moisture movement equations . . . . . . . . . . . . . . . .
10.3 Phase change of soil moisture . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.3.1 Ice sheet process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58
58
58
59
60
60
61
62
63

11 Lake11.1 Calculation of lake surface conditions . . . . . . . . . . . . . . . . . . . . . . . . .
11.1.1 Calculation of lake surface albedo . . . . . . . . . . . . . . . . . . . . . . .
11.1.2 Lake surface roughness . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.2 Solution of energy balance at lake surface . . . . . . . . . . . . . . . . . . . . . .
11.3 Calculation of lake ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.3.1 Calculation of heat ﬂux and growth rate . . . . . . . . . . . . . . . . . . .

63
63
66
66
67
69
70

5

11.3.2 Sublimation and freshwater ﬂux for lake . . . . . . . . . . . . . . . . . . .
11.3.3 Updating lake ice fraction . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.3.4 Growth and Melting of lake ice and snow . . . . . . . . . . . . . . . . . .
11.4 Physical formulation and process . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.4.1 Setting the vertical diﬀusion coeﬃcients . . . . . . . . . . . . . . . . . . .
11.4.2 Estimate the diﬀusion terms of the tracer equations . . . . . . . . . . . .
11.4.3 Time integration of the tracer equations . . . . . . . . . . . . . . . . . . .
11.4.4 The vertical convection . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11.5 Lake river coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73
75
76
81
81
82
83
86
87

12 Wetland
12.1 Outline of wetland scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12.2 Inﬂow and outﬂow of the wetland . . . . . . . . . . . . . . . . . . . . . . . . . . .
12.3 Storage of the surface runoﬀ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12.4 Water input of soil surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88
88
88
89
89

13 Tile scheme

90

References

93

6

1

Introduction

Minimal Advanced Treatments of Surface Interaction and RunOﬀ (MATSIRO) is a land surface parameterization originally formulated in the early 2000’s (Takata et al., 2003) to produce
terrestrial boundary condition of the atmospheric general circulation model developed by theCenter for Climate System Research at the University of Tokyo and the National Institute forEnvironmental Studies (CCSR/NIES AGCM), as well as to other global climate models. It was
originally designed to be primarily used for integration of climate simulations such as those involving long time scales from one month to several hundred years coupled with the atmosphericmodel at grid resolutions of tens of kilometers or more. The main objective in its development
was to represent all of the important water and energy exchange processes between land and atmosphere as fully and accurately as possible (i.e., advanced treatment) in such time and spatial
scales, while modeling them as simply as possible (i.e., minimal treatment) so as to allow theresults to be easily interpreted.
At the very early stage (2000’s), when the ﬁrst description was written by Emori (Emori,2001), MATSIRO was developed based on the land surface submodel of CCSR/NIES AGCM5.4g
coupled with the parameterization for a vegetated surface (canopy) by (Watanabe, 1994), whileat the same time improving certain processes such as those related to snow and runoﬀ. Subsequently, with modiﬁcations in the structure of AGCM, changes were made dealing with ﬂux
couplers and parallel processing so as to make it compatible with AGCM5.6. With regard to
the physiological processes of vegetation, a Jarvis-type function was initially used for stomatal
resistance. Later, however, the Farquhar-type photosynthesis scheme, which served as a de facto
standard in the world due to the progress of studies on climate-ecosystem interactions at thattime, was ported from SiB2 code.
Since then, MATSIRO has been used in CCSR/NIES/FRCGC AGCM5.7, which is adopted
by an atmosphere-ocean coupled model MIROC3 (K-1 Model Developers, 2004), and MIROC5(Watanabe et al., 2010) with only minor model updates and bug ﬁxes. Instead, during thesetimes, MATSIRO had been intensively used with so-called oﬄine setting (i.e., without couplingwith the atmospheric model, but forced with given surface atmospheric forcings) with original
updates and improvements of some processes. Such eﬀorts generate some branches of the model,
called MAT-GW [water table depth dynamics; Koirala et al. (2014)], HiGW-MAT [coupled withanthropogenic activity and groundwater; Pokhrel et al. (2015)], IsoMATSIRO [stable water
isotopes included; Yoshimura et al. (2006)]. MATSIRO also has been validated and comparedwith other land surface parameterizations under some international protocols like PILPS, GSWP,
and ISIMIP. By these attempts, it was revealed that, apart from its original design, MATSIRO
was able to be used for higher resolution (up to about 1 kilometer) and shorter time scales(subdaily), so that MATSIRO is now used for short-term ﬂood prediction over Japan with a
new hydro-dynamical model, CaMa-Flood (Yamazaki et al., 2011), for example (Ma et al., 2021).
Here, with development of MIROC6 (Tatebe et al., 2019), MATSIRO adopted signiﬁcantmodel updates, so we decided to update its description for the ﬁrst time since 2001. We named
this version of MATSIRO as MATSIRO6. Major and minor updates in this document are as
follows:• Snow cover fraction diagnosis (Chapter 8)• Runoﬀ representation (Chapter 9)• Lake process coupled with river (Chapter 11)• Snow-fed wetland process (Chapter 12)• Tiling with lake, cropland, and natural vegetation (Chapter 13)
It should be noted that even though river routing process is also implemented in MIROC6
7

and other versions, this document does not include description of river process.

1.1

Structure

MATSIRO consists of ﬂux calculation section and land integration section. In the ﬂux calculation
section, the calculations are conducted separately for snow-covered and snow-free potions. For
each snow-free portion (l = 1) and snow-covered portion (l = 2), the subroutines for various
processes are called, the ﬂuxes are calculated, and the ground surface temperature and canopy
temperature are updated. Speciﬁcally, the following subroutines are called in the order shown
below:(a) MATLAI: vegetation type parameter (LAI, vegetation height) set(b) MATRAD: calculation of radiation parameters (albedo, vegetation transmissivity, etc.)(c) MATBLK: calculation of turbulence parameters (bulk coeﬃcients)(d) MATRST: calculation of stomatal resistance, bare soil surface evaporation resistance, etc.
(e) MATFLX: calculation of surface ﬂux(f) MATGHC: calculation of ground heat conduction(g) MATSHB: solution of surface heat balance
Then, ﬂuxes from lake surface are calculated separately for ice-covered and ice-free portions.
The following subroutines are called in the order shown below:(h) LAKEBC: calculation of lake surface conditions (albedo, roughness, etc.)(i) SFCFLX: calculation of surface ﬂux(j) RADSFC: calculation of radiation ﬂux (downward and upward shortwave radiation)(k) LAKEHB: solution of energy balance at lake surface
In the land integration section, the subroutines for various processes are called and land
surface prognostic variables are updated. Speciﬁcally, the following subroutines are called in theorder shown below:(l) MATCNW: calculation of canopy water balance(m) MATSNW: calculation of snow water equivalent, snow temperature, and snow albedo(n) MATROF: calculation of runoﬀ(o) MATGND: calculation of soil temperature, soil moisture, and frozen soil
Finally, the lake modules are called and the related prognostic variables are updated.
(p) SETSCNV: calculation of convergence of shortwave radiation(q) LAKEIC: calculation of lake ice(r) LAKEPO: calculation of lake water and temperature(s) PUTDEFF: sending a lake water deﬁcit for the lakes whose water levels are below thelower limit

1.2

Prognostic variables

MATSIRO has the following prognostic variables for each tile:8

Table 1.1: Prognostic variables
Variable

Description

Units

Ts(l) (l = 1, 2)
Tc(l) (l = 1, 2)
Tg(k) (k = 1, . . . , Kg )w(k) (k = 1, . . . , Kg )wi(k) (k = 1, . . . , Kg )wc
SnTSn(k) (k = 1, . . . , KSn )αSn(b) (b = 1, 2, 3)
ASn
µ
Dm
Rsn
Isn
AL
Thist
ρd
ρd
αice(b)Ground surface temperature
Canopy temperature
Soil temperature
Soil moisture content
Frozen soil moisture content
Water content on the canopy
Snow water equivalent
Snow temperature
Snow albedo
Snow cover fraction
Accumulated snowfall
Accumulated snowmelt
Reset ﬂag for snow accumulation or ablation season
Snow index
Land cover fraction
Long-term mean temperature
Dust density in snow
Dust density in snow (mass)
Ice albedo

[K]
[K]
[K]
[m3 /m3 ][m3 /m3 ][m][kg/m2 ]
[K]
[−]
[−]
[kg/m2 ][kg/m2 ]
[−]
[−]
[−]
[K]
[ppmw][ppmw][−]

where l = 1, 2 denotes snow-free and snow-covered portions, respectively; k is the verticallayer number of the soil or snow (the uppermost layer is 1, with the number increasing as thelayer becomes deeper); Kg is the number of soil layers; KSn is the number of snow layers; andb = 1, 2, 3 denotes the bands of visible, near infrared, and infrared wavelengths, respectively.
As the default setting, the soil has six layers whose thicknesses are deﬁned by the depth
boundaries of 5, 20, 75, 100, 200, and 1000 cm from the surface. The deﬁnition points of soil
temperature, soil moisture, and frozen soil moisture are the same. The maximum number ofsnow layer is three by default, while the number is diagnosed from snow water equivalent. As
the default setting, the maximum number is three layers.
The ground surface temperature and canopy temperature are so-called surface temperatures
whose heat capacity is zero; however, they take the form of prognostic variables. (The current
calculation method depends on the values of the preceding step because the stability, etc. assessed by the values of the preceding step are used. If the stability, etc. were to be assessedby updated values and calculation iterated to the point of convergence, perfect diagnostic variables would be obtained that would not depend on the values of the preceding step.) The othervariables are all prognostic variables that always require the values of the preceding step.
The ground surface temperature and canopy temperature are updated in the ﬂux calculation
section. All of the other variables (original prognostic variables) are updated in the land surface
integration section.
The lake part has the following prognostic variables:
Table 1.2: Prognostic variables for the lake scheme
Variable

Description

Units

T

Lake temperature

[◦ C]

9

1.3

Variable

Description

Units

S
TI
AI
hI
hS
h

Lake salinity
Lake ice surface temperature
Lake ice concentration
Mean lake ice thickness over ice-covered part of a grid
Mean snow depth over lake ice
Lake level

[PSU]
[◦ C]
[−]
[cm]
[cm]
[cm]

Output data

The following variables are output from the ﬂux calculation section:
Table 1.3: Output data from the ﬂux calculation section
Variable

Description

Units

τx
τy
H
E
RS↑
↑
RL
αs(b) (b = 1, 2, 3)
TsR
Fg(1/2)
FSn(1/2)Et(i,j) (i = 1, 2; j = 1, 2, 3)∆Fconv
Froot(k) (k = 1, . . . , Kg )
LAI
ASnc

Surface eastward wind stress
Surface northward wind stressSensible heat ﬂux
Latent heat ﬂux
Upward shortwave radiation ﬂux
Upward longwave radiation ﬂux
Surface albedo
Surface radiation temperature
Surface heat transfer ﬂux
Heat transfer ﬂux for snow surface
Evapotranspiration
surface energy convergence
Root sucking ﬂux
leaf area index
Canopy freezing area ratio[N/m2 ]
[N/m2 ]
[W/m2 ]
[kg/m2 /s][W/m2 ]
[W/m2 ][−]
[K]
[W/m2 ]
[W/m2 ]
[kg/m2 /s][W/m2 ]
[kg/m2 /s][m2 /m2 ][−]where i = 1, 2 denotes liquid and solid evapotranspiration, respectively; and j = 1, 2, 3
denotes evaporation from the bare soil surface (forest ﬂoor), transpiration, and canopy water
evaporation, respectively. Other indexes are the same as described earlier.
The following variable is output from the land surface integration section:
Table 1.4: Output data
Variable

Description

Units

Ro

runoﬀ

[kg/m2 /s]

Runoﬀ is used as an input variable for the river channel network model.

1.4

Input data

The following variables are input in the ﬂux calculation section:

10

Table 1.5: Input data for land ﬂux calculation
Variable

Description

Units

ua
va
Ta
qa
Pa
Ps
↓
R(d,b)(d = 1, 2; b = 1, 2, 3)

Atmospheric 1st layer eastward wind
Atmospheric 1st layer northward wind
Atmospheric 1st layer temperature
Atmospheric 1st layer speciﬁc humidity
Atmospheric 1st layer pressure
Surface pressure
Surface downward radiation ﬂux

[m/s]
[m/s]
[K]
[kg/kg][Pa][Pa][W/m2 ]

where d = 1, 2 denotes direct and diﬀuse, respectively; and b = 1, 2, 3 denotes the bands ofvisible, near infrared, and infrared wavelengths, respectively.
The following variables are input in the land surface integration section:
Table 1.6: Input data for land integration section

1.5

Variable

Description

Units

P rc
P rl
PSnc
PSnl
Fg(1/2)
FSn(1/2)Et(i,j) (i = 1, 2; j = 1, 2, 3)∆Fconv
Froot(k) (k = 1, . . . , Kg )
LAI
ASnc
Ddust
DBC

Convective rainfall ﬂux
Large-scale Rainfall ﬂux
Convective snowfall ﬂux
Large-scale snowfall ﬂux
Surface heat transfer ﬂux
Heat transfer ﬂux for snow surface
Evapotranspiration
Surface energy convergence
Root water uptake
Leaf area index
Canopy freezing area ratio
Dust deposition ﬂux
Black carbon deposition ﬂux[kg/m2 /s][kg/m2 /s][kg/m2 /s][kg/m2 /s][W/m2 ]
[W/m2 ]
[kg/m2 /s][W/m2 ]
[kg/m2 /s][m2 /m2 ][−][kg/m2 /s][kg/m2 /s]External parameters

The external parameters necessary for the execution of MATSIRO are broadly divided into two
types: parameters whose values for each grid cell are given by horizontal distribution (map),and parameters whose values are given by land cover type or soil type tables. The land covertypes and soil types are the parameters given by map, and through this, each parameter givenby table is allocated to individual grid cells; that is,parameter given by a map:
ϕ(i, j)

(1.1)

ψ(I),

(1.2)

parameter given by a table:

11I = IL (i, j)(1.3)

I = IS (i, j)

(1.4)

or

where (i, j) are indexes of the grid horizontal location, IL is the land use type, and IS is the soil
type.
The types of external parameters ϕ given by map are as follows:
Table 1.7: External parameters given by maps
Variable

Description

Temporal resolution

Units

IL
IS
LAI0
α0(b) (b = 1, 2, 3)tan βsσzLand cover type
Soil Type
Leaf Area Index (LAI)
Ground surface (forest ﬂoor) albedo
Tangent of the mean surface slope
elevation standard deviation

constant
constant
every month
constant
constant
constant[−][−][m2 /m2 ]
[−]
[−]
[m]

The types of external parameters ψ given by table for each land cover type are as follows:

12

13

Description
vegetation height
Height of the bottom of the canopy
Leaf albedo
Leaf transmissivity
Percentage of root presence
Momentum exchange coeﬃcient between the individual leaves and the atmosphere
Heat exchange coeﬃcient between individual leaves and the atmosphere
Vegetation coverage
Rubisco reaction capacity
An -gs slope of the relationshipAn -gs relationship intercepts
Photosynthetic eﬃciency per photon
Coupling factor between wc and we
Coupling factor between wp and ws
Respiratory coeﬃcient
Critical temperature of high temperature suppression
Critical temperature of cryogenic suppression

Variableh0
hB0
rf (b) (b=1,2)
tf (b) (b=1,2)
froot(k) (k = 1, . . . , Kg )
cd
ch
fV
Vmax
m
b
ϵ3 , ϵ4
θce
θps
fd
s2
s4

Table 1.8: External parameters given by tables

[m]
[m]
[−]
[−]
[−]
[−]
[−]
[−]
[m/s]
[−]
[m/s]
[m/s/moll]
[−]
[−]
[−]
[K]
[K]

Units

The types of external parameters ψ given by table for each soil type are as follows:
Table 1.9: External parameters given by tables

2

Variable

Description

Unitscg(k) (k = 1, . . . , Kg )kg(k) (k = 1, . . . , Kg )wsat(k) (k = 1, . . . , Kg )Ks(k) (k = 1, . . . , Kg )ψs(k) (k = 1, . . . , Kg )b(k) (k = 1, . . . , Kg )Speciﬁc heat of soil
Thermal conductivity of soil
Soil porosity
Saturated permeability of soil
Soil saturation moisture potential
Index of soil moisture potential curve[J/m3 ]
[W/m/K][m3 /m3 ]
[m/s]
[m]
[−]

Vegetation type parameters

The leaf area index (LAI), vegetation height, etc. are set as vegetation type parameters.
Seasonally changing horizontal distributions are loaded as external parameters for LAI, andthe values according to land use type are loaded as external parameters for the heights of the
canopy top and bottom. When there is snow, only the vegetation above the level of the snow
depth is taken into consideration and the type parameters are corrected as follows:h = max(h0 − DSn , 0)

(2.1)

hB = max(hB0 − DSn , 0)h − hB
LAI = LAI0
h0 − hB0

(2.2)
(2.3)

where h is the height of the canopy top (vegetation height), hB is the height of the canopy
bottom (height of clear length), LAI is the leaf area index, and h0 , hB0 , and LAI0 are therespective values when there is no snow. DSn is the snow depth. LAI is approximated on theassumption that it is uniformly distributed vertically between the canopy top and bottom.
In terms of coding, next the mean values of the snow-free and snow-covered portions are
solved by weighting with the snow-covered ratio (ASn ), etc., as follows:h = ASn h + (1 − ASn )h0(2.4)

However, because the snow-free portion and snow-covered portion are respectively calculated,
it should be noted that ASn takes the value of either 0 (snow-free portion) or 1 (snow-covered
portion), so no mixing of values occurs (similar cases are also seen later).

3

Radiation parameters

Next, the radiation parameters (albedo, vegetation transmissivity , etc.) are calculated.

3.1

Calculation of ground surface (forest ﬂoor) albedo

The horizontal distributions of the ground surface (forest ﬂoor) albedo b = 1, 2 are loaded asexternal parameters, with b = 1, 2 denoting the wavelength bands of visible and near infrared,
respectively. The infrared ground surface albedo (α0(3) ) is set to a ﬁxed value (horizontal
distributions can also be prepared if desired).
14

With regard to the ice-sheet portion and snow-covered portion, the dependence of the incidence angle of albedo is considered by the following function form:
α0(d,b) = α̂0(b) + (1 − α̂0(b) ) · 0.4(1 − cos ψin(d) )5(3.1)

where b = 1, 2 are wavelength bands; d = 1, 2 are direct and diﬀuse, respectively; and α̂0(b) isthe value of albedo when the incidence angle is 0 (from directly overhead). The cosine of theincidence angle cos ψin(d) is expressed ascos ψin(1) = cos ζ,

cos ψin(2) = cos 50◦(3.2)

for direct insolation and diﬀuse radiation, respectively, where ζ is the solar zenith angle.
With regard to regions other than the ice-sheet portion and snow-covered portion, the zenith
angle dependence is not taken into consideration for the albedo of the ground surface (forestﬂoor) and the same values are given to direct insolation and diﬀuse radiation; that is,α0(d,b) = α0(b)(d = 1, 2; b = 1, 2)(3.3)

Moreover, in the case of the infrared wavelength only diﬀuse radiation needs to be considered.
The value of the infrared albedo for all ground surfaces is given independently of the zenith angle,
as follows:
α0(2,3) = α0(3)

3.2

(3.4)

Calculation of canopy albedo and transmissivity

The calculation of canopy albedo and transmissivity is based on the calculation of radiation
within a canopy layer proposed by Watanabe and Ohtani (1995).
Considering the canopy as vertically uniform and making use of several assumptions for
simpliﬁcation, the transfer equations of insolation within the canopy and the boundary condition
are expressed as
dSd↓
dL
dSr↓
dL
dSr↑
dL
Sd↓ (0)= −F sec ζSd↓

(3.5)

= −F (1 − tf (b) )df Sr↓ + F tf (b) sec ζSd↓ + F rf (b) df Sr↑(3.6)= F (1 − tf (b) )df Sr↑ − F rf (b) (df Sr↓ + sec ζSd↓ )(3.7)= Sdtop(3.8)

Sr↓ (0) = Srtop(3.9)

Sr↑ (LAI)=

α0(1,b) Sd↓ (LAI) + α0(2,b) Sr↓ (LAI)(3.10)

where Sd↓ is the downward direct insolation; Sr↑ and Sr↓ are the upward and downward diﬀuse
radiation, respectively; L is the leaf area cumulatively added downward from the canopy top; df
is the diﬀusivity factor (= sec 53◦ ), rf (b) ; tf (b) are the leaf albedo and transmissivity, respectively(the same value is used for diﬀuse radiation and direct insolation); and F is a factor denoting
the direction of the leaves with respect to the radiation. Here, the distribution of the directionof the leaves is assumed to be random (F = 0.5) for simplicity.

15

These can be solved analytically, giving the following solutions:Sd↓ (L) = Sdtop exp(−F · L · sec ζ)Sr↓ (L) =Sr↑ (L)=(3.11)

C1 eaL + C2 e−aL + C3 Sd↓ (L)A1 C1 eaL + A2 C2 e−aL + C4 Sd↓ (L)(3.12)(3.13)

where
a = F df [(1 − tf (b) )2 − rf2 (b) ]1/2
A1 =
A2 =

(3.14){1 − tf (b) + [(1 − tf (b) ) − rf2 (b) ]1/2 }/rf (b){1 − tf (b) − [(1 − tf (b) )2 − rf2 (b) ]1/2 }/rf (b)2A3 = (A1 − α0(2,b) )eaLAI− (A2 − α0(2,b) )e−aLAI(3.15)(3.16)(3.17)

C1 = {−(A2 − α0(2,b) )e−aLAI (Srtop − C3 Sdtop )+ [C3 α0(2,b) + α0(1,b) − C4 ]Sd↓ (LAI)}/A3(3.18){(A1 − α0(2,b) )e(Srtop − C3 Sdtop )[C3 α0(2,b) + α0(1,b) − C4 ]Sd↓ (LAI)}/A3(3.19)

C2 =
−

aLAI

C3 =

sec ζ[tf (b) sec ζ + df tf (b) (1 − tf (b) ) + df rf2 (b) ]d2f [(1 − tf (b) )2 − rf2 (b) ] − sec2 ζ(3.20)C4 =rf (b) (df − sec ζ) sec ζ2df [(1 − tf (b) )2 − rf2 (b) ] − sec2 ζ(3.21)

Albedo αs at the canopy top is expressed asSr↑ (0) = αs(1,b) Sd↓ (0) + αs(2,b) Sr↓ (0)(3.22)

therefore,
αs(2,b) = {A2 (A1 − α0(2,b) )eaLAI − A1 (A2 − α0(2,b) )e−aLAI }/A3

(3.23)

αs(1,b) = −C3 αs(2,b) + C4 + (A1 − A2 )(C3 α0(2,b) + α0(1,b) − C4 )e−F ·LAI·sec ζ /A3 (3.24)
are obtained.
If the canopy transmissivity (Tc ) (speciﬁcally, the ratio of incident insolation absorbed by
the forest ﬂoor to the incident insolation of the canopy top) is deﬁned by
Tc(2,b) = {(1 − A2 )(A1 − α0(2,b) ) − (1 − A1 )(A2 − α0(2,b) )}/A3(3.25)

Tc(1,b) = −C3 Tc(2,b)+ {(C3 α0(2,b) + α0(1,b) − C4 )((1 − A1 )eaLAI − (1 − A2 )e−aLAI )/A3 + C3 − C4 + 1}e−F ·LAI·sec ζ(3.26)
the following are obtained:
Sd↓ (LAI) + Sr↓ (LAI) − Sr↑ (LAI) = Tc(1,b) Sd↓ (0) + Tc(2,b) Sr↓ (0)(3.27)

The above calculations are performed for b = 1, 2 (visible and near infrared), respectively.
16

The leaf albedo rf and transmissivity tf are loaded as external parameters for each land
cover type; however, the following two modiﬁcations are made before these parameters are usedin the above calculations.
1. Snow (ice) eﬀect on leaf surface
When the canopy temperature does not exceed 0 °C, the canopy water is regarded as snow(ice). In this case, using the snow albedo (αSn(b) ) and canopy water (wc ), the following assumptions are made:rf (b) = (1 − fcwet )rf (b) + fcwet αSn(b)(3.28)

fcwet = wc /wc,cap(3.29)

where wc,cap is the canopy water capacity. With regard to transmissivity, the following assumption is made for convenience so that the absorptivity 1 − rf − tf ) does not take a negativevalue:
tf (b) = (1 − fcwet )tf (b) + fcwet tSn(b) ,tSn(b) = min(0.5(1 − αSn(b) ), tf (b) )(3.30)

When the canopy water is in the liquid state, the resultant change in the leaf radiation
parameter is ignored. Moreover, although the cases of snowfall after interception by the canopy(snow cover) and of frozen water in the canopy (ice) can be considered, with each having diﬀerent
radiation characteristics, the same albedo of snow on the forest ﬂoor is used for all cases here.
2. Eﬀect considering the directions of reﬂection and transmission
In solving the above equations, all reﬂected light is assumed to return in the direction ofincidence. However, if, for example, the diﬀusion of only a portion of the reﬂected light in the
direction of incidence is taken into consideration, the leaf radiation parameters can be replaced
as follows (Watanabe, personal communication):rf (b) = 0.75rf (b) + 0.25tf (b)(3.31)

tf (b) = 0.75tf (b) + 0.25rf (b)(3.32)

The above calculations are performed for b = 1, 2 (visible and near infrared), respectively.
In addition, in consideration of the uneven distribution of vegetation (such as savanna) inparts of the grid cells, before calculating the albedo, etc., the LAI (taking the original LAI as
the grid mean value) of the vegetation-covered portion is calculated as follows:
LAI = LAI/fV(3.33)↓and this is used for the calculation for albedo, etc. mentioned above. (R(d,b)) is the vegetationcovered ratio in the grid cell. After the albedo, etc. are calculated, the area-weighted mean ofthe vegetation-covered portion and non- vegetation-covered portion are obtained asαs(d,b) = fV αs(d,b) + (1 − fV )α0(d,b)(3.34)

Tc(d,b) = fV Tc(d,b) + (1 − fV )(1 − α0(d,b) )(3.35)

17

3.3

Calculation of surface radiation ﬂux, etc.

↓
Using the surface downward radiation ﬂux (R(d,b)) and albedo calculated above, the following
radiation ﬂuxes are calculated:

RS↓ =
RS↑ =
↓
RL

=

RSgnd =

2 X
2X
b=1 d=1
2 X
2X↓
R(d,b)(3.36)↓
αs(d,b) R(d,b)(3.37)

b=1 d=1↓
R(2,3)
2 X
2X(3.38)↓Ts(d,b) R(d,b)(3.39)

b=1 d=1

P AR =

2
X↓
R(d,1)(3.40)d=1↓
iswhere RS↓ and RS↑ are the downward and upward shortwave radiation ﬂux, respectively; RL
gnd
the downward longwave ﬂux; RS is the shortwave ﬂux absorbed by the forest ﬂoor; and P AR
is the downward photosynthesis active radiation (PAR) ﬂux.
The canopy transmissivity of shortwave and longwave radiation, and the emissivity of longwave radiation, are then calculated as follows:

TcS = RSgnd /(RS↓ − RS↑ )(3.41)

TcL = exp(−F · LAI · df )(3.42)

ϵ = 1 − αs(2,3)

4(3.43)

Turbulence parameters (bulk coeﬃcient)

Next, the turbulence parameter (bulk coeﬃcient) is calculated.

4.1

Calculation of roughness with respect to momentum and heat

The calculation of roughness is based on Watanabe (1994). In that study, using the results of amultilayer canopy model by Kondo and Watanabe (1992) as a function form for the roughnessof a bulk model best ﬁtting those results, Watanabe (1994) proposed the following:

18




0.45

h − d −1z0s −1/0.45
exp(−2A+ )ln= 1 − exp(−A+ ) + − ln
(4.1)
z0
h
!−1
P 2

1
h−d
P1
=
(4.2)
ln †
− ln(zT s /h) P1 + A+ exp(A+ )
zT
"
#0.9

 

 0 1/0.9
Ch − d −1h − d −1Tlnln
= CT∞ 1 − exp(−P3 A+ ) +exp(−P4 A+ ) (4.3)
z0
zT
CT∞
h − d = h[1 − exp(−A+ )]/A+
cd LAIA+ =2k 2
1h−d h−d
ln †
= ln
0
z0
CT
z

(4.4)
(4.5)
(4.6)

T

−1 + (1 + 8FT )1/2
2  z 0.1
zT s 0.22
Ts= 0.0115
exp 5
h
h

 z 0.35 Ts= 0.55 exp −0.58
h

CT∞ =

(4.7)

P1

(4.8)

P2

P3 = [FT + 0.084 exp(−15FT )]0.15(4.10)

2FT1.1(4.11)= ch /cd(4.12)

P4 =
FT

(4.9)

where z0 and zT are the roughness of the overall canopy with respect to momentum and heat,
respectively; z0 s and zT s are the roughness of the ground surface (forest ﬂoor) with respect tomomentum and heat, respectively; cd and ch are the exchange coeﬃcient between an individualleaf and the atmosphere with respect to momentum and heat, respectively; h is the vegetation
height; d is the zero-plane displacement; and LAI is LAI. zT† is the roughness with respect toheat when assuming no transfer of heat to or from the leaf surface, and is used when solving thecoeﬃcient of heat transfer from the forest ﬂoor.
z0s and zT s are given as external data for each land cover type. Their values (z0s = 0.05
m, zT s = 0.005 m) are ﬁxed as standards. However, the following modiﬁcations are made withrespect to the snow-covered portion:z0s = max(fSn z0s , z0Sn )(4.13)

zT s = max(fSn z0s , zT Sn )(4.14)

fSn = 1 − DSn /z0s

(4.15)

where DSn , z0Sn and zT Sn are the roughness of the snow-covered portion with respect to momentum and heat, respectively.
cd and ch are parameters determined by the leaf shape, and are given as external data foreach land cover type.

19

4.2

Calculation of bulk coeﬃcient with respect to momentum and heat

After Watanabe (1994), the bulk coeﬃcient is also calculated using Monin-Obukhov similarity
as

CM
CH
CHs

−2

za − d+ Ψm (ζ)= k 2 ln
z0

−1 
−1
za − dza − d= k 2 ln+ Ψm (ζ)
ln+ Ψh (ζ)
z0
zT
#−1

−1 "
za − dza − d2= k lnln+ Ψm (ζg )+ Ψh (ζg )z0z†(4.16)(4.17)(4.18)T

CHc = CH − CHs(4.19)

where CM and CH are the bulk coeﬃcients of the overall canopy (leaf surface + forest ﬂoor)
with respect to momentum and heat, respectively; CHs is the bulk coeﬃcient of the groundsurface (forest ﬂoor) ﬂux with respect to heat; CHc is the bulk coeﬃcient of the canopy (leafsurface) ﬂux with respect to heat; Ψm and Ψh are Monin-Obukhov shear functions with respect
to momentum and heat, respectively; and za is the reference height of the atmosphere (heightof the troposphere). Using the Monin-Obukhov lengths ζ and ζg related to the overall canopy
and ground surface (forest ﬂoor), respectively, L and Lg are respectively expressed as:za − dLza − dζg =Lg
ζ=(4.20)(4.21)and the Monin-Obukhov lengths are expressed as:
3/2

L =

Θ0 CM |Va |2
kg(CHs (Ts − Ta ) + CHc (Tc − Ta ))Ls =Θ0 CM |Va |2
kgCHs (Ts − Ta )(4.22)

3/2(4.23)

where Θ0 =300K; |Va | is the absolute value of the surface wind speed; k is the Karman constant;g is the gravitational acceleration; and Ta ,T c$ and Ts are the temperature of the troposphere,
canopy (leaf surface), and ground surface (forest ﬂoor), respectively.
Since the bulk coeﬃcient is necessary for calculation of the Monin-Obukhov length, andthe Monin-Obukhov length is necessary for calculation of the bulk coeﬃcient, the calculation isiterated (twice as a standard) with a neutral bulk coeﬃcient as the initial value.
Prior to this calculation, the snow depth in the snow-covered portion is added to the zeroplane displacement. However, the upper limit is set so that the zero-plane displacement does
not exceed the value of za :d = min(d + DSn , fmax · za )As a standard, fmax is set at 0.5.

20

(4.24)

4.3

Calculation of bulk coeﬃcient with respect to vapor

This calculation is performed after the calculation of stomatal resistance, described later.
When the stomatal resistance (rst ) and ground surface evaporation resistance (rsoil ) havebeen solved, the bulk coeﬃcient with respect to vapor is solved as:−1(CHc |Va |)−1 + rst /LAI−1
CEs |Va | = (CHs |Va |)−1 + rsoil
CEc |Va | =(4.25)(4.26)(Previously, this parameter was solved by converting stomatal resistance, etc. into a decreaseof the exchange coeﬃcient via roughness. However, since this approach seems to be problematic,
a simpler method had been adopted in its place.)
In addition, when there is no stomatal resistance, etc. (such as evaporation from wet surfaces), the same value as for the bulk coeﬃcient of heat is used for the bulk coeﬃcient of vapor.

5

Stomatal resistanceFor the calculation of stomatal resistance, a photosynthesis-stomatal model based on Farquhar
et al. (1980), Ball (1988), and Collatz et al. (1990; 1991; 1992) is used. The code of SiB2
(Sellers et al., 1996b) is used virtually unchanged, with the exception of the method for solving
the resistance of the overall canopy. A Jarvis-type empirical equation could be used instead;however, the explanation of this point is omitted here.

5.1

Calculation of soil moisture stress factor

Soil moisture stress with respect to transpiration is solved. By solving the soil moisture stressfactor in each soil layer, and weighting with the root distribution in each layer, the stress factorof the overall soil is calculated.
Referring to SiB2 (Sellers et al., 1996b), the soil moisture stress in each layer is evaluatedby the following equation:
fw(k) = [1 + exp(0.02(ψcr − ψk ))]−1(k = 1, . . . , Kg )(5.1)

The stress factor of the overall soil is then obtained by

fw =

Kg
X

fw(k) froot(k)

(5.2)

k=1

where froot(k) is the root distribution fraction in each layer, which is an external parameter for
PKgeach land cover type.
k=1 froot(k) = 1.
Furthermore, the weight of transpiration distributed to the root uptake ﬂux in each layer isexpressed as

frootup(k) = fw(k) froot(k) /fw
Note thatPKg

k=1 frootup(k) = 1 here.

21(k = 1, . . . , Kg )

(5.3)

5.2

Calculation of amount of photosynthesis

The amount of photosynthesis is calculated after SiB2 (Sellers et al., 1996b).
The amount of photosynthesis is considered to be regulated by the following three upper
limits:A ≤ min(wc , we , ws )(5.4)

where wc is the upper limit set by the eﬃciency of photosynthesis enzymes (Rubisco), and weis the upper limit set by photosynthetically active radiation. ws is the upper
limit of theeﬃciency of use of photosynthate (sink) in the case of C3 vegetation, or the upper limit set by
CO2 concentration in the case of C4 vegetation (Collatz et al., 1991, 1992).
The respective magnitudes are estimated as follows:


wc =




Vm

ci − Γ∗ci + Kc (1 + O2 /KO )(in case of C3 vegetation)

Vm

(5.5)

(in case of C4 vegetation)
(5.6)



we =




P AR · ϵ3P AR · ϵ4ci − Γ∗ci + 2Γ∗
(in case of C3 vegetation)(5.7)(in case of C4 vegetation)
(5.8)


ws =

Vm /2
Vm ci /5(in case of C3 vegetation)(in case of C4 vegetation)(5.9)(5.10)

where Vm is the Rubisco reaction capacity, ci is the partial pressure of CO2 in the stoma, O2 isthe partial pressure of oxygen in the stoma, and P ARis the photosynthetically active radiation(PAR). is the CO2 compensation point, which is expressed by Γ∗ is the compensation point of
CO2 and is represented by Γ∗ = 0.5O2 /S Kc , KO , and S are functions of temperature, whosefunction form is shown later. ϵ3 and ϵ4 are constants determined by the vegetation type.
In order to express a smooth transition between the diﬀerent upper limits, 5.4 is actually
solved as
βce wp2 − wp (wc + we ) + wc we = 0(5.11)

βps A2 − A(wp + ws ) + wp ws = 0(5.12)

The amount of net photosynthesis An can be obtained when solving the two equations inorder while selecting the smaller of the two solutions for each equation. βce , βps are constants
determined by the vegetation type. Note that when β = 1, coincidence is achieved with a simple
minimum-value operation.

22

When the amount of photosynthesis has been solved, (An ) is solved as
An = A − Rd(5.13)

where Rd is the amount of respiration, expressed as
Rd = fd Vm(5.14)

Here, fd is a constant determined by the vegetation type.
Vm , etc. depend on the temperature and soil moisture, as follows (note that although thetemperature dependence diﬀers according to the term in which Vm appears, the value is expressedby the same Vm ):Vm = Vmax fT (Tc )fw(5.15)

Kc = 30 × 2.1(5.16)

QT

KO = 30000 × 1.2

QT(5.17)

S = 2600 × 0.57

QT

(5.18)


2.1QT /{1 + exp[s1 (Tc − s2 )]}



QT


 1.8Q /{1 + exp[s3 (s4 − Tc )]}
2.1 T /{1 + exp[s1 (Tc − s2 )]}/{1 + exp[s3 (s4 − Tc )]}fT (Tc ) =(when wc and we are for C4 )



Q

T
 1.8

 QT
2 /{1 + exp[s5 (Tc − s6 )]}QT(when wc and we for C3 )(when ws is for C3 )(5.19)(when ws is for C4 )(when Rd )= (Tc − 298)/10(5.20)

where Vmax , s1 , . . . , s6 are constants determined by the vegetation type.
If Vmax , P AR, ci , Tc and fw are given by the above, the amount of photosynthesis in an
individual leaf can be calculated. In reality, these values can be considered to be distributed
unevenly even in the same canopy; however, ci , Tc , and fw are approximated here as being thesame for all leaves, whereas vertical distribution is taken into consideration in the case of Vmax
and P AR. P AR is greater at the top of the canopy, and the lower the position in the canopy,
the more it is attenuated. Vmax is also considered to conform with this property of P AR and tohave a similar distribution.
The average vertical distribution ofP AR (and therefore the vertical distribution of Vmax ) isexpressed as
P AR(L) = P ARtop exp(−fatn aL)(5.21)

where L is the leaf area added cumulatively from the canopy top, P ARtop is P AR at the canopy
top, a is the attenuation coeﬃcient deﬁned in 3.15, and fatn is a constant for adjustment. Using
this, the factor favr expressing the average value of P AR, is deﬁned as follows:
Z LAI.
P AR(L)dL

favr =(LAI · P ARtop ) =023

1 − exp(−fatn aL)
fatn a(5.22)

Since each of the terms An (wc , ws , we , Rd ) is proportional to Vmax or P AR, based on theassumption that the vertical distributions of Vmax and P AR are proportional, by multiplying
An , which was solved using the values of Vmax and P AR at the top of the canopy, by favr , theaverage amount of photosynthesis of leaf (An ) can be solved:
An = favr An(5.23)

This parameter is expressed as An hereafter.

5.3

Calculation of stomatal resistance (2)The net photosynthesis (An ) and stomatal conductance (gs ) are related by the semiempiricalequation of Ball (1988) as follows:
gs = mAn
hs + bfw
cs(5.24)

where cs is the molar fraction of CO2 (number of mol of CO2 per 1 mol of air) at the leaf surface,
fw is the soil moisture stress factor, and m and b are constants determined by the vegetation
type.
hs is the relative humidity at the leaf surface and is deﬁned as
hs = es /ei(5.25)

where es is the molar fraction of vapor at the leaf surface, ei is the molar fraction of vapor
in the stoma, and ei = e∗ (Tc ) is the mole fraction of water vapor in the stomata. e∗ denotes themolar fraction of saturated vapor.
Assuming that the vapor ﬂux from the inside of the stoma to the leaf surface is equal tothe vapor ﬂux from the leaf surface to the atmosphere (i.e., that there is no convergence anddivergence of vapor at the leaf surface),gs (ei − es ) = gl (es − ea )(5.26)

es = (gl ea + gs ei )/(gl + gs )(5.27)

from which we obtain

where ea is the molar fraction of vapor in the atmosphere and gl is the conductance from theleaf surface to the atmosphere. gl is expressed by gl = CHc |Va |/LAI using the bulk coeﬃcient.
Similarly, assuming that there is no convergence and divergence of CO2 at the leaf surface,
An = gl (ca − cs )/1.4 = gs (cs − ci )/1.6(5.28)

cs = ca − 1.4An /gl(5.29)

ci = cs − 1.6An /gs(5.30)

from which we obtain

where ca and ci are the molar fractions of CO2 in the atmosphere and in the stoma, respectively.
The numerical values 1.4 and 1.6 are constants that appear due to the diﬀerence in the diﬀusion24

coeﬃcients of vapor and CO2 .
If we order the equations by substituting 5.24, the following equation is obtained for gs :Hgs2 + (Hgl − ei − Hbfw )gs − gl (Hbfw + ea ) = 0(5.31)

H = (ei cs )/(mAn )(5.32)

However, since

5.30 is used for cs .
Among the two solutions of (5.31), the larger one is the signiﬁcant solution. From the above,
if An is known, gs can be solved; however, when solving gs , ci is used. ci can be solved by (5.30)
if gs is solved. That is, An is necessary in order to solve gs , whereas ci , namely gs , is necessary
in order to solve An . Iterative calculation is therefore required.
The algorithm for the iterative calculation is ported from SiB2, which uses the method ofquickening the convergence by iterating six times and putting the errors in decreasing order toestimate the next solution.
Lastly, using stomatal conductance, the stomatal resistance is expressed as
rst = 1/gst

5.4(5.33)

Calculation of ground surface evaporation resistanceThe ground surface evaporation resistance (rsoil ) and relative humidity of the uppermost soil
layer (hsoil ) are calculated as follows:
rsoil = a1 (1 − W(1) )/(a2 + W(1) )ψ(1) ghsoil = exp
Rair Tg(1)(5.34)(5.35)

where W(1) = w(1) /wsat(1) is the degree of saturation of the uppermost soil layer, ψ1 is the
moisture potential of the uppermost soil layer, g is the gravitational acceleration, Rair is thegas constant of the air, and Tg(1) is the temperature of the uppermost soil layer. a1 and a2 areconstants, with a1 = 800, a2 = 0.2. as standard values.

6

Surface energy balance

6.1

Calculation of surface turbulent ﬂuxes

The turbulent ﬂuxes at the ground surface are solved by bulk formulae as follows. Then, bysolving the surface energy balance, the ground surface temperature (Ts ) and canopy temperature(Tc ) are updated, and the surface ﬂux values with respect to those values are also updated. The
solutions obtained here are temporary values. In order to solve the energy balance by linearizing
with respect to Ts and Tc , the diﬀerential with respect to Ts and Tc of each ﬂux is calculated
beforehand.
• Momentum ﬂux

25τx = −ρCM |Va |ua(6.1)

τy = −ρCM |Va |va(6.2)

where τx and τy are the momentum ﬂuxes (surface stress) of the zonal and meridional directions,
respectively.
• Sensible heat ﬂux

Hs = cp ρCHs |Va |(Ts − (Ps /Pa )κ Ta )(6.3)

Hc = cp ρCHc |Va |(Tc − (Ps /Pa ) Ta )

(6.4)

κ

∂Hs /∂Ts = cp ρCHs |Va |(6.5)∂Hc /∂Tc = cp ρCHc |Va |(6.6)

where Hs and Hc are the sensible heat ﬂux from the ground surface (forest ﬂoor) and canopy(leaf surface), respectively; κ = Rair /cp and Rair are the gas constants of air; and cp is the
speciﬁc heat of air.
• Bare soil surface (forest ﬂoor) evaporation ﬂux∗
g
Et(1,1) = (1 − ASn )(1 − fice ) · ρC
Es |Va |(hsoil q (Ts ) − qa )

(6.7)

∗
g
Et(2,1) = (1 − ASn )fice · ρC
Es |Va |(hsoil q (Ts ) − qa )

(6.8)

∗
g
∂Et(1,1) /∂Ts = (1 − ASn )(1 − fice ) · ρC
Es |Va |hsoil · dq /dT |Ts

(6.9)

g
∂Et(2,1) /∂Ts = (1 − ASn )fice · ρC
Es |Va |hsoil · dq /dT |Ts(6.10)∗

where Et(1,1) and Et(2,1) are the water evaporation and ice sublimation ﬂuxes at the baresoil surface, respectively; q ∗ (Ts ) is the saturation speciﬁc humidity at the ground surface temperature; hsoil is the relative humidity at the soil surface layer; ASn is the snow-covered ratio;
and fice is the ratio of ice in the uppermost soil layer, expressed as

fice = wi(1) /w(1)(6.11)

Since the snow-free portion and snow-covered portion are calculated separately, it should
be noted that ASn takes the value of either 0 (snow-free portion) or 1 (snow-covered portion).
When the ﬂux is downward (i.e., dew formation), there is no soil moisture resistance; therefore,
the bulk coeﬃcient is taken as:

g
C
Es =



CEs (hsoil q ∗ (Ts ) − qa > 0)
CHs (hsoil q ∗ (Ts ) − qa ≤ 0)

• Transpiration ﬂux

26(6.12)

∗g
Et(1,2) = (1 − fcwet ) · ρC
Ec |Va |(q (Tc ) − qa )(6.13)

Et(2,2) = 0(6.14)∗g∂Et(1,2) /∂Tc = (1 − fcwet ) · ρC
Ec |Va | · dq /dT |Tc(6.15)∂Et(2,2) /∂Tc = 0(6.16)

where Et(1,2) and Et(2,2) are transpiration of water and ice, respectively; and Et(2,2) is always
0. fcwet = wc /wc,cap is the wet fraction of the canopy. When the ﬂux is downward, which isconsidered to be dew formation on the dry part of the leaf, the bulk coeﬃcient is taken as:CEc (q ∗ (Tc ) − qa > 0)
CHc (q ∗ (Tc ) − qa ≤ 0)

(6.17)

Et(1,3) = fcwet · ρCHc |Va |(q ∗ (Tc ) − qa )(6.18)

Et(2,3) = 0(6.19)gC
Ec =• Canopy evaporation ﬂux
When Tc ≥ 0 ◦ C:∗

∂Et(1,3) ∂Tc = fcwet · ρCHc |Va | · dq /dT |Tc(6.20)∂Et(2,3) ∂Tc = 0(6.21)

when Tc < 0 ◦ In case of C:

Et(1,3) = 0(6.22)∗

Et(2,3) = fcwet · ρCHc |Va |(q (Tc ) − qa )∂Et(1,3) ∂Tc = 0(6.23)(6.24)∗∂Et(2,3) ∂Tc = fcwet · ρCHc |Va | · dq /dT |Tc(6.25)

where Et(1,3) and Et(2,3) are the evaporation of water and the sublimation of ice at the canopy
surface, respectively.
• Snow sublimation ﬂux

ESn = ASn · ρCHs |Va |(q ∗ (Ts ) − qa )∗∂ESn /∂Ts = ASn · ρCHs |Va | · dq /dT |Ts(6.26)(6.27)

where ESn is the snow sublimation ﬂux. Since the snow-free portion and snow-covered portionare calculated separately, it should also be noted here that ASn takes the value of either 0(snow-free portion) or 1 (snow-covered portion).

27

6.2

Calculation of heat conduction ﬂuxes

The heat conduction ﬂuxes in the snow-free and snow-covered portions are calculated. Similarly
to the turbulent ﬂuxes, when the energy balance is solved later and the surface temperature is
updated, the heat conduction ﬂux values are updated with respect to that value.
In addition, it should also be noted here that since the snow-free portion and snow-covered
portion are calculated separately, ASn takes the value of either 0 (snow-free portion) or 1 (snowcovered portion).
• Heat conduction ﬂux in the snow-free portion
Fg(1/2) = (1 − ASn ) · kg(1/2) /∆zg(1/2) (Tg(1) − Ts )∂Fg(1/2) /∂Ts = −(1 − ASn ) · kg(1/2) /∆zg(1/2)(6.28)(6.29)

where Fg(1/2) is the heat conduction ﬂux, kg(1/2) is the soil heat conductivity, ∆zg(1/2)
is the thickness from the temperature deﬁnition point of the uppermost soil layer to theground surface, and Tg(1) is the temperature of the uppermost soil layer.
• Heat conduction ﬂux in the snow-covered portion
FSn(1/2) = ASn · kSn(1/2) /∆zSn(1/2) (TSn(1) − Ts )∂FSn(1/2) /∂Ts = −ASn · kSn(1/2) /∆zSn(1/2)(6.30)(6.31)

where FSn(1/2) is the heat conduction ﬂux, kSn(1/2) is the snow heat conductivity, ∆zSn(1/2)
is the thickness from the temperature deﬁnition point of the uppermost snow layer to theground surface, and TSn(1) is the temperature of the uppermost snow layer.

6.3

Solution of energy balance at ground surface and canopy

The energy balance is solved for two cases: (1) when there is no melting at the ground surface,and (2) when there is melting at the ground surface. In case (2), the solution is obtained byﬁxing the ground surface temperature (Ts ) at 0 °C, and the energy available for use in melting
is diagnosed from the energy balance. Snowmelt on vegetation is treated by correction later on;
therefore, that case is not solved separately here. Moreover, the case of the snow completely
melting within the time steps is also treated by correction later on.
6.3.1

Energy balance at ground surface and canopy

The energy divergence at the ground surface (forest ﬂoor) is∆Fs = Hs + Rsnet + lEt(1,1) + ls (Et(2,1) + ESn ) − Fg(1/2) − FSn(1/2)(6.32)

where l and ls are the latent heat of evaporation and sublimation, respectively; and Rsnet is the
net radiation divergence at the ground surface, given by↓
Rsnet = −(RS↓ − RS↑ )TcS − ϵRL
TcL + ϵσTs4 − ϵσTc4 (1 − TcL )where σ is the Stefan-Boltzmann constant.
The energy divergence at the canopy (leaf surface) is
28(6.33)

∆Fc = Hc + Rcnet + l(Et(1,2) + Et(1,3) ) + ls (Et(2,2) + Et(2,3) )(6.34)

where Rcnet is the net radiation divergence at the canopy, given by↓
Rcnet = −(RS↓ − RS↑ )(1 − TcS ) − ϵRL(1 − TcL ) + (2ϵσTc4 − ϵσTs4 )(1 − TcL )6.3.2(6.35)

Case 1: When there is no melting at the ground surface

When there is no melting at the ground surface, ∆Fs = ∆Fc = 0 are solved so that Ts and Tc
holds true for the energy balance at the ground surface and canopy.
The energy balance equation linearizing each term with respect to Ts and Tc can be expressed
as


∆Fs
∆Fc

current


=

∆Fs
∆Fc

past


+

∂∆Fs /∂Ts ∂∆Fs /∂Tc∂∆Fc /∂Ts ∂∆Fc /∂Tc



∆Ts
∆Tc




=

0
0


(6.36)The part with pst on the right-hand side is where the ﬂuxes calculated 6.35 using the valuesof Ts and Tc obtained in the previous step.
The diﬀerential terms are as follows:
∂∆Fs
∂Ts
∂∆Fs
∂Tc
∂∆Fc
∂Ts
∂∆Fc
∂Tc

=
=
=
=



∂Et(1,1)∂Fg(1/2) ∂FSn(1/2)∂Et(2,1) ∂ESn
∂Hs ∂Rsnet
+
+l
+ ls
+
−
−
(6.37)
∂Ts
∂Ts
∂Ts
∂Ts
∂Ts
∂Ts
∂Ts
∂Rsnet(6.38)∂Tc∂Rcnet(6.39)∂Ts
∂Et(1,2) ∂Et(1,3)∂Et(2,2) ∂Et(2,3)∂Hc ∂Rcnet
+
+l
+
+
+ ls
(6.40)
∂Tc
∂Tc
∂Tc
∂Tc
∂Tc
∂Tc

where
∂Rsnet∂Ts∂Rsnet∂Tc∂Rcnet∂Ts∂Rcnet∂Tc= ϵ4σTs3(6.41)= −(1 − TcL )ϵ4σTc3(6.42)= −(1 − TcL )ϵ4σTs3(6.43)= 2(1 − TcL )ϵ4σTc3(6.44)

Using the above equations, Eq. (140) is solved for Ts and Tc .

29

6.3.3

Case 2: When there is melting at the ground surface

When either there is snow on the ground surface or the land cover type is ice sheet, and also
the ground surface temperature solved in case 1, Tscurrent = Tspast + ∆Ts , is higher than 0 °C,
melting at the ground surface occurs. When there is melting at the ground surface, the groundsurface temperature is ﬁxed at 0 °C. That is:∆Ts = ∆Tsmelt = Tmelt − Tspast(6.45)

where Tmelt is the melting point (0 °C) of ice.
With Tc known, ∆Ts is solved by the following equation similarly to 6.36:

.
∂∆Fc
∂∆Fc
past
melt
∆Tc = −∆Fc
−
∆Ts
∂Ts
∂Tc

(6.46)

Thus, ∆Ts and ∆Tc are determined, and the energy convergence at the ground surface to be
used for melting is solved by the following equation:∆Fconv = −∆Fscurrent = −∆Fspast −
6.3.4

∂∆Fs
∂∆Fs
∆Tsmelt −
∆Tc
∂Ts
∂Tc

(6.47)

Conditions for solutions

Several conditions are set for the solution of the ground surface energy balance. After solving
the energy balance, if any of the conditions are not followed, the ﬂux that has contravened theconditions is ﬁxed at the limit value that satisﬁes the conditions, and the energy balance issolved again.
1. Vapor in the troposphere should not be excessively removed.
Due to the instability of temporal calculations, it is possible that large downward latent heatis produced. The conditions are set so that even in such a case, the vapor in the tropospherefrom the surface is not completely removed; that is,> −qa (Ps − Pa )/(g∆t)
Etcurrent(i,j)(i = 1, 2; j = 1, 2, 3)

currentESn> −qa (Ps − Pa )/(g∆t)(6.48)(6.49)

where g is the gravitational acceleration and ∆t denotes the time steps of the atmospheric model.
For the values of Et etc. to be used for judgment, the updated ﬂux values (current) with respect
to the values of Ts and Tc that have been updated so as to satisfy the energy balance are used.
The same applies to all of the other conditions listed below. Updating of the ﬂux values is
described later.
2. Soil moisture should not take a negative value.
Soil moisture should not take a negative value due to transpiration; that is,30

X

Etcurrent<(1,2)

ρw wk ∆zg(k) /∆tL

(6.50)

k∈rootzone

where ρw is the water density and ∆tL denotes the time steps of the land surface model.
3. Canopy water should not take a negative value.
Canopy water should not take a negative value due to evaporation; that is,

Etcurrent< ρw wc /∆tL(i,3)(i = 1, 2)(6.51)

4. The snow water equivalent should not take a negative value.
The snow water equivalent should not take a negative value due to sublimation of snow; thatis,current
ESn< Sn/∆tL

6.3.5(6.52)

Updating of ground surface and canopy temperatures

The ground surface temperature and canopy temperature are updated as follows:

Tscurrent = Tspast + ∆Ts(6.53)

Tccurrent = Tcpast + ∆Tc(6.54)

Based on the updated canopy temperature, the canopy water is diagnosed in advance as
being either liquid or solid. This information is used when treating freezing and melting of the
canopy water, as follows:ASnc =

0 (Tc ≥ Tmelt )1 (Tc < Tmelt )(6.55)

where ASnc is the frozen fraction on the canopy.
6.3.6

Updating of ﬂux values

The ﬂux values are updated with respect to the updated values of Ts and Tc . When F denotes
any given ﬂux, updating of the values is performed as follows:
F current = F past +

∂F
∂F
∆Ts +
∆Tc
∂Ts
∂Tc

(6.56)

Using the updated ﬂux values, the ﬂuxes output into the atmosphere, etc. are calculated as31

follows:
H = Hs + Hc
3 X
2
X
E =
Et(i,j) + ESn(6.57)(6.58)

j=1 i=1↑↓
RL= TcL ϵσTs4 + (1 − TcL )ϵσTc4 + (1 − ϵ)RL

(6.59)↑
TsR = (RL/σ)1/4(6.60)

where TsR is the radiation temperature at the ground surface.
The root uptake ﬂux in each soil layer is then calculated as follows:
Froot(k) = frootup(k) Et(1,2)(k = 1, . . . , Kg )(6.61)

where Froot(k) is the root uptake ﬂux and frootup(k) is the weighting for distribution of thetranspiration to the root uptake ﬂux in each layer.

7

Canopy Water Balance

The canopy water balance is calculated.

7.1

Diagnosis of canopy water phase

With regard to canopy water, the liquid phase (intercepted rainfall, dew formation, and frozen
water content that has melted) and solid phase (intercepted snow, icing, and liquid water content
that has frozen) are considered separately and the coexistence of the two phases is allowed. The
only prognostic variable is the water content (wc ) encompassing both the liquid and solid phases,
and depending on whether the canopy temperature (Tc ) is higher or lower than Tmelt = 0◦ C,
it is diagnosed as liquid or solid, respectively. The reason why the liquid and solid phases cancoexist is that Tc in snow-covered and snow-free portions is calculated separately. That is, thefrozen fraction on the canopy (ASnc ) is deﬁned (in actuality, it is obtained as a result of spatial
averaging by the coupler) as follows:


 0

(1 − ASn )ASnc =
A
 Sn


1

(Tc(1) ≥ Tmelt , Tc(2) ≥ Tmelt )(Tc(1) < Tmelt , Tc(2) ≥ Tmelt )(Tc(1) ≥ Tmelt , Tc(2) < Tmelt )(Tc(1) < Tmelt , Tc(2) < Tmelt )(7.1)

where wcl = wc (1−ASnc ) and wci = wc ASnc are the liquid and solid water content of the canopy,
respectively.
+1
For ASnc , the value updated in the ﬂux calculation section AτSnc
is given by the coupler,τ
but the value of the previous step ASnc is stored in MATCNW. τ denotes the time steps. This issolved from the initial values of Tc and Sn at the time of initiating the calculation, and therefore
does not become a new prognostic variable.

32

7.2

Prognosis of canopy water

The prognostic equations for the canopy water in the liquid and solid phases are given respectively as
τ +1
τ
wcl
− wcl
= PIl − El − Dl + Mc∆tLτwτ +1 − wciρw ci= PIi − Ei − Di − Mc
∆tL

ρw

(7.2)
(7.3)

where PIl and PIi are the precipitation interception in the respective cases, El and Ei are theevaporation (sublimation), Dl , Di are the dripping, and Mc is the melting. Note that here, theτ and w τ are deﬁned using Aτ
values before the updated wcl
ci
Snc before it is updated, as follows:

7.2.1

τ
wcl
= wcτ (1 − AτSnc )

(7.4)

τ
wci
= wcτ AτSnc(7.5)

Evaporation (sublimation) of canopy waterFirst, by subtracting the evaporation (sublimation), the canopy water is partially updated asfollows. The evaporation (sublimation) has already been solved in the ﬂux calculation section.
∗
τ
wcl
= wcl
− El ∆tL /ρw

(7.6)

∗
τ
wci
= wci
− Ei ∆tL /ρw(7.7)

El = Et(1,3)(7.8)

Ei = Et(2,3)(7.9)

Then, if either wcl or wci become negative in value, it is supplemented by the other until thevalue returns to 0, and the melting (negative value in the case of frozen water) that is assumed
to be produced is then inserted in Mc .
7.2.2

Interception of precipitation by the canopy

The precipitation interception and dripping are considered by separating the places of convective
precipitation and nonconvective precipitation. The fraction of the convective precipitation area
(Ac ) is assumed to be uniform (0.1 as a standard value). Stratiform precipitation is also assumed
to be uniform.

PIlc = fint (P rc /Ac + P rl )(7.10)

PIlnc
c
PIi
nc
PIi

= fint P rl

(7.11)= fint (PSnc /Ac + PSnl )(7.12)= fint PSnl

(7.13)

33

c denote the interception in the convective precipitation area, and P nc andwhere PIlc and PIiIlnc
PIi denote the interception in the nonconvective precipitation area. fint is the interception
eﬃciency, and is simply given by
fint =LAI (LAI < 1)
1
(LAI ≥ 1)(7.14)

By adding the intercepted precipitation, the canopy water is further partially updated asfollows:
c∗
∗
wcl
= wcl
+ PIlc ∆tL /ρw
nc∗
wcl
c∗
wci
nc∗
wci

7.2.3

=
=
=

∗
wcl
+ PIlnc ∆tL /ρw
∗
c
wci
+ PIi
∆tL /ρw∗nc
wci + PIi ∆tL /ρw(7.15)(7.16)(7.17)(7.18)Dripping of the canopy water

For dripping, dripping due to the canopy water capacity being exceeded and natural dripping
due to gravity are considered, as follows:c∗c∗)− wc,cap , 0) + Dg (wcl
Dlc = max(wcl

Dlnc
Dic
Dinc

=
=
=

nc∗
nc∗
)
− wc,cap , 0) + Dg (wcl
max(wclc∗c∗max(wci − wc,cap , 0) + Dg (wci )nc∗
nc∗
max(wci− wc,cap , 0) + Dg (wci)(7.19)(7.20)(7.21)(7.22)

where the canopy water capacity (wc,cap ) is, from the water capacity per unit leaf area (wc max )and LAI, assumed to be
Wc,cap = Wc max LAI(7.23)

Wc max is set at 0.2 mm as a standard value, and the same value is used with respect to theliquid and solid phases.
The natural dripping due to gravity Dg is, after Rutter et al. (1975), assumed to be
Dg (wc ) = D1 exp(D2 wc )(7.24)

D1 = 1.14 × 10−11 and D2 = 3.7 × 103 are standard values, and the same values are usedwith respect to the liquid and solid phases.
By subtracting the dripping, the values are updated as follows:
c∗∗
c∗
c
wcl
= wcl
− DIl
∆tL /ρw
nc∗∗
wcl
c∗∗
wci
nc∗∗
wci

=
=
=

nc∗
nc
wcl
− DIl
∆tL /ρwc∗c
wci − DIi ∆tL /ρw
nc∗
nc
wci
− DIi
∆tL /ρw

34(7.25)(7.26)(7.27)(7.28)

7.2.4

Updating and melting of canopy water

Moreover, by taking the average of the convective precipitation area and nonconvective precipitation area, the canopy water can be updated as follows:
∗∗
c∗∗
nc∗∗
wcl
= Ac wcl+ (1 − Ac )wcl∗∗
wci
wcτ +1

=
=

c∗∗
nc∗∗
Ac wci+ (1 − Ac )wci
∗∗
∗∗
+ wci
wcl

(7.29)(7.30)(7.31)

However, if updating of the frozen fraction (ASnc ) is considered,
τ +1
+1
wcl
= wcτ +1 (1 − AτSnc)(7.32)

τ +1
wci(7.33)=+1wcτ +1 AτSnc

The melting Mc is therefore diagnosed asτ +1∗∗
Mc = −ρw (wci
− wci
)/∆tL

(7.34)

When the melting is produced during evaporation, that portion is added.
Here, the canopy temperature should be changed due to the latent heat of melting; however,
it is impossible because we are ignoring the heat capacity of the canopy. Moreover, although it
would be advantageous to change the temperature of the surrounding atmosphere, this is also
not possible in view of the need for agreement with the calculation in the land surface integration
section. Hence, for convenience, in order to conserve the energy of the system, the latent heatof melting is given as the heat ﬂux to the soil (or snow).

7.3

Fluxes given to the soil, snow, and runoﬀ process

The water ﬂux Fw given to the snow or the runoﬀ process after interception by the canopy
is respectively expressed with respect to the convective precipitation area and nonconvective
precipitation area, and the liquid and solid phases, as follows:c
Fwl= (1 − fint )(P rc /Ac + P rl ) + Dlc(7.35)=(1 − fint )P rl + Dlnc

(7.36)=(1 − fint )(PSnc /Ac + PSnl ) + Dic(1 − fint )PSnl + Dinc

(7.37)

nc
Fwl
c
Fwi
nc
Fwi

=

(7.38)For the calculation of runoﬀ, convective rainfall and stratiform rainfall are given separately,while snowfall is consolidated because separation is not necessary, as follows:cnc
P rc∗ = Ac(Fwl
− Fwl
)

P rl∗
∗
PSn

(7.39)=

nc
Fwl(7.40)=cnc
Ac Fwl+ (1 − Ac )Fwl

(7.41)∗ are the convective precipitation, the stratiform precipitation, andwhere P rc∗ , P rl∗ , and PSn
the snowfall after interception by the canopy, respectively.

35

The energy ﬂux correction portion for the soil or the snow is∆Fc,conv = −lm Mc(7.42)

where lm is the latent heat of melting.

8

Snow

The snow cover fraction, snow water equivalent, snow temperature and snow albedo are calculated here.
Most of the processes are included in SUBROUTINE MATSNW in matsnw.F, but the ice
albedo is calculated in SUBROUTINE ICEALB in matice.F.
The following tables present attributions of key variables in SUBROUTINE MATSNW.

36

37

Long nameAmount of accumulated snow
Snow temperature of the kth layer
Snow albedo for band b
Snow fraction
Accumulated snow
Accumulated snowmelt depth
Water input after canopy interception derived from convective precipitation (input) and input to ground (output)
Water input after canopy interception derived from layered precipitation (input) and input to ground (output)
Density of dust and black carbon in the kth snow layer calculated using Dw
Density of dust and black carbon in the kth snow layer calculated using Dm

VariableSnTSn(k) (k = 1, 2, 3)
αb (b = 1, 2, 3)
ASn
WSn
MSn
Prc
Prl
ρdw (k) (k = 1, 2, 3)
ρdm (k) (k = 1, 2, 3)

Table 8.1: Modiﬁed
Name in
GLSNW
GLTSN
GLASN
GLRSN
GLSDA
GLSDM
WINPC
WINPL
CDST
CDSTM

Unit
kg/m2
K
kg/m2
kg/m2
kg/m2 /s
kg/m2 /s
ppmv
ppmv

Table 8.2: Output
Variable

Description

Unit

Name in code

Rogl

Glacier formationkg/m2 /s

GGLACR

Table 8.3: Input
Variable

Description

Unit

Name in code

PSn
ESn
FSn(1/2)
Tg(k)
wg(k)
wc
ASnc
Dw
Dm
zsd
Thist
-

SnowfallSnow sublimation
Snow surface heat ﬂux
Soil temperature of the kth layer
Soil moisture
Canopy water
Canopy snow ratio
Weighted mean of the ﬂuxes of dust and black carbon
Flux of dust and black carbon
Standard deviation of topography
Annual mean temperature over the latest 30 years
Index of the surface condition
Soil typekg/m2 /s
kg/m2 /s
W/m2T
m3 /m3m
kg/m2 /s
kg/m2 /s
m
K
-

SNFAL
SNSUB
SNFLXS
GLG
GLW
GLWC
SNRATCDSTFAL
DSTFLMGRZSD
T2HISTILSFCILSOIL

8.1

Diagnosis of snow cover fraction

MATSIRO has two ways of calculation of the snow cover fraction, and the user can switch themwith the option OPT SSNOWD.
8.1.1

Case 1: When OPT SSNOWD is active

The snow cover fraction is diagnosed in the SUBROUTINE SSNOWD DRV, a driver of a Subgrid
SNOW Distribution (SSNOWD) submodel developed by Liston (2004), with a physically based
parameterization of sub-grid snow distribution considering various factors such as diﬀerences in
topography, the time of snowfall or snow melting, etc (Nitta et al., 2014; Tatebe et al., 2019).
The snow cover fraction is formulated for accumulation and ablation seasons separately. At
the grids where snow does not disappear completely, the ﬂag of transition from ablation season
to accumulation season is raised. It occurs at the time step next of 1 August in the northernhemisphere and 1 February in the southern hemisphere. The transition is judged based on thesnow amount, accumulated snow, net snow gain and snowfall. This ﬂag is included in the restartvariables similarly to the state variables.
For the accumulation season, snowfall occures uniformly and the snow cover fraction isassumed to be unity in the grid cell. For the ablation season, the snow cover fraction decreases
based on the sub-grid distribution of the snow water equivalent. Under the assumption of
uniform melt depth WSnm , the sum of snow-free and snow-covered fraction equals unity:
Z WSnmZ ∞f (WSn )dWSn +0f (WSn )dWSn = 1,
WSnm

38

(8.1)

where WSn is the snow water equivalent depth and f (WSn ) is the probability distribution function of snow water equivalent depth within the grid cell. The snow depth distribution withineach grid cell is assumed to follow a lognormal distribution:
" 
 #
1
1 ln(WSn ) − λ 2√ exp −f (WSn ) =,2ζWSn ζ 2π(8.2)

where λ = ln(WSn ) − 12 ζ 2 and ζ 2 = ln(1 + CV 2 ).
Here WSn is the accumulated snow and CV is the coeﬃcient of variation. With the default
settings, CV is diagnosed from the coldness index Thist , standard deviation of the subgrid
topography zsd and vegetation type that is a proxy for surface winds. For Thist , the annually
averaged temperature over the latest 30 years using the time relaxation method of Krinneret al. (2005), in which the timescale parameter is set to 16 years, is applied. The temperaturethreshold for a category diagnosis is set to 0 and 10 ◦ C. When Thist ≥ 10, CV takes constantvalue of 0.0.6. In other two cases of 0 ≤ Thist < 10 and Thist < 0, it is determined from zsd andvegetation type.
The snow amount Sn is given by
Z WSnmSn(WSnm ) =Z ∞0[f (WSn )]dWSn +0(WSn − WSnm )[f (WSn )]dWSn ,

(8.3)

WSnm

and this equation is rewritten to






zWSnm − ξ
zWSnm
ζ2
1
1
√
√
erfc
Sn(WSnm ) = exp λ +− WSnm erfc
,
2
2
2
2
2

(8.4)

√
)−λ
, and zWSnm is the value of z when WSn = WSnm and erfc iswhere ξ = (1 − 2)z, z = ln(WSnζ
the complementary error function. WSnm is calculated from this equation and the snow amountSn using Newton-Raphson methods (in SUBROUTINE SSNOWD ITR in ssnowd.F).
Then, the snow cover fraction ASn (WSnm ) is calculated by
ASn (WSnm ) = 1 −Z WSnm
0

8.1.2

z

1
WSnm
f (WSn )dD = erfc
.
2
2

(8.5)

Case 2: When OPT SSNOWD is inactive

The snow cover fraction is diagnosed in SUBROUTINE SNWRAT. The snow cover fraction isformulated as a function of the snow amount Sn:
Sn(WSnm ) = min(p
Sn/Snc ),(8.6)

where Snc is 100 kg/m2 as a standard.

8.2

Vertical division of snow layers

In order to express the vertical distribution of the snow temperature, when the snow waterequivalent is large, the snow is divided into multiple layers and the temperature is deﬁned ineach layer. The number of snow layers can be varied, with the number of layers increasing as

39

the snow water equivalent becomes larger. A minimum of one layer and a maximum of threelayers are set as a standard.
This process is treated in SUBROUTINE SNWCUT in matsnw.F.
The number of layers and the mass of each layer are determined uniquely by the snow waterequivalent. Consequently, the mass of each layer does not become a new prognostic variable.
f (k) (k = 1, 2, 3)) is determined as follows (k = 1
As a standard, the mass of each layer (∆Sn
is the uppermost layer):f < 20)
f

(Sn
 Sn
f
f < 40)f (20 ≤ Sn
∆Sn(1) =
0.5Sn

 20
f
(Sn ≥ 40)f < 20)

0
(Sn


 f
f < 60)Sn − ∆Sn(1) (20 ≤ Sn,
f (2) =
∆Sn
f
f

(60
≤
Sn
<
100)
0.5(
Sn
−
20)


 40
f ≥ 100)(Sn(f < 60)0(Snf (3) =∆Snf
f ≥ 60)Sn − (∆Sn(1) + ∆Sn(2) ) (Sn(8.7)

f = Sn/ASn .
where Snf is the snow water equivalent in the snowSn is the grid-mean snow water equivalent, and Sn
f (k) ) is also the value of the snow-covered
covered portion. Note that the mass of each layer (∆Sn2
portion, not the grid-mean value. The unit is kg/m .
From the above, it can be clearly seen that the number of snow layers (KSn ) is as follows,
as a standard:


0



1
KSn =

2



3

8.3

f = 0)(Snf < 20)(0 < Snf < 60)(20 ≤ Snf ≥ 60)(Sn(8.8)

Calculation of snow water equivalent

The prognostic equation of the snow water equivalent is given bySnτ +1 − Snτ∗= PSn− ESn − MSn + F rSn
∆tL

(8.9)

∗ is the snowfall ﬂux after interception by the canopy, E
where PSn
Sn is the sublimation ﬂux,
MSn is the snowmelt, and F rSn is the refreeze of snowmelt or the freeze of rainfall.

8.3.1

Sublimation of snow

First, by subtracting the sublimation, the snow water equivalent is updated:Sn∗ = Snτ − ESn ∆t,

40(8.10)

f ∗(1) = ∆Snf τ(1) − ESn /ASn ∆tL∆Sn(8.11)

The asterisk indicates that the variable is under updating in the time step.
In a case where the sublimation is larger than the snow water equivalent in the uppermost
snow layer, the remaining amount is subtracted from the layer below. If the amount in thesecond layer is insuﬃcient for such subtraction, the remaining amount is subtracted from thelayer below that.
8.3.2

Snowmelt

Next, the snow heat conduction is calculated to solve the snowmelt. The method of calculating
the snow heat conduction is described later. The updated snow temperature incorporating the∗heat conduction is assumed to be TSn(k). When the snow temperature is calculated and thetemperature of the uppermost snow layer becomes higher than Tmelt = 0◦ C, the calculation ofsnow heat conduction is performed again with the ﬁxed snow temperature of the uppermost
snow layer at Tmelt . In this case, the energy convergence ∆Feconv in the uppermost snow layer
is calculated. This is not the grid-mean value but the value of the snow-covered portion. The
snowmelt in the uppermost snow layer isfSn(1) = min(∆Feconv /lm , ∆Snf ∗(1) /∆t).
M(8.12)

With regard to the second snow layer and below, if the estimated snow temperature is higher
than Tmelt , it is adjusted to Tmelt and the desidual energy from the adjustment is applied to thesnowmelt. That is, it is assumed to be∗∗TSn(k)= Tmelt .

(8.13)∗f ∗ /∆t,∆Feconv = (TSn
− Tmelt )cpi ∆Sn(k)(k)(8.14)∆Feconv is newly deﬁned bywhere cpi is the speciﬁc heat of snow (ice), and the snowmelt is solved as in Eq. 8.12.
By subtracting the snowmelt, the mass of each layer is updated:∗∗∗f (k) = ∆Snf (k) − M
fSn .
∆Sn
(k)

(8.15)

During these calculations, when a certain layer is fully melted, the remaining amount of∆Feconv is given to the layer below to raise the temperature in that layer; that is,∗
fSn ,∆Feconv= ∆Feconv − lm M
(k)(8.16)∗∗∗∗f ∗(k+1) )∆t.
TSn
= TSn
/(cpi ∆Sn
+ ∆Feconv
(k+1)
(k+1)

(8.17)∗
When all of the snow is melted, ∆Feconv
is given to the soil.
The snowmelt of the overall snow is the sum of the snowmelt in each layer (note, however,41

that it is the grid-mean value):

MSn =

K
Sn
X

fSn(k) ASn
M(8.18)

k=1

By subtracting the snowmelt, the snow water equivalent is updated:Sn∗∗ = Sn∗ − MSn ∆t.
8.3.3(8.19)

Freeze of snowmelt water and rainfall in snow

The freeze of snowmelt water and rainfall in the snow is calculated next. With regard to thesnowmelt water, consideration is given to the eﬀect of the liquid water produced by the snowmelt
in the upper layer refreezing in the lower layer. The retention of liquid water content in theaccumulated snow is not considered, and the entire amount is treated whether it has frozen in
the snow or percolated under the snow.
The liquid water ﬂux at the snow upper boundary in the snow-covered portion is
FewSn(1) = P rc∗ + P rl∗ + MSn /ASn .

(8.20)

Here, the melted portion in the second layer of the snow and below is also assumed to havepercolated from the snow upper boundary (in actuality, snowmelt in the second layer or below
rarely occurs).
It is reasonable to assume the temperature of the snowmelt water as 0 ◦ C, and the temperature of rainfall on the snow is also assumed to be 0 ◦ C for convenience. The temperature of thesnow increases due to the latent heat of the freezing of water; however, when the temperatureof the snow in a certain layer is increased to 0 ◦ C, any additional water is assumed to be unable
to freeze and to percolate to the layer below. In addition, an upper limit is set on the ratio ofwater that can be frozen compared with the mass of snow in the layer. The amount of freeze in
frSn(k) is solved bya given layer F∗∗

frSn = min Few
F
Sn(k) ,
(k)

∗∗ )
f (k)
cpi (Tmelt − TSn
∆Sn
(k)

lm

∆t

∗∗

, fF max

f (k)
∆Sn
∆t

!
,

(8.21)

where FewSn(k) is the liquid water ﬂux ﬂowing from the top of the kth layer of snow cover. fF max
is assumed to be 0.1 as a standard value.
The snow temperature change is updated by
∗∗∗
TSn
=
(k)

frSn ∆t + cpi (T ∗∗ ∆Sn
f ∗∗ + Tmelt F
frSn ∆t)
lm F
(k)
Sn(k)
(k)
(k)
f ∗∗
f
cpi (∆Sn(k) + F r Sn(k) ∆t),(8.22)and the mass is updated as follows:
∗∗∗

∗∗

f +F
frSn ∆t.
f (k) = ∆Sn
∆Sn
(k)
(k)

(8.23)

The amount of freeze in the overall snow is the sum of the amounts of freeze in each layer42

(note, however, that it is the grid-mean value):

F rSn =K
Sn
X

frSn ASn .
F(k)(8.24)

k=1

By adding the amount of freeze, the snow water equivalent is partially updated as follows:Sn∗∗∗ = Sn∗∗ + F rSn ∆t.

(8.25)The liquid water that has percolated from the snow to the lower boundary is given to the
soil.
8.3.4

SnowfallLastly, by adding the snowfall after interception by the canopy, the ﬁnally updated snow waterequivalent is obtained:∗Snτ +1 = Sn∗∗∗ + PSn∆t.

(8.26)

However, when the temperature of the uppermost soil layer is 0 ◦ C or more, the snowfall isassumed to melt on the ground. In this case, the energy of the latent heat of melting is takenfrom the soil.
When snow is produced by snowfall in a grid where no snow was formerly present, the snowcovered ratio (ASn ) is newly diagnosed by Eq. 8.5 and the snow temperature (TSn(1) ) is assumed
to be equal to the temperature of the uppermost soil layer.
The snowfall is added to the mass of the uppermost layer:
τ +1

∗∗∗

∗
f (k) = ∆Snf (k) + PSn∆Sn∆t/ASn .

8.3.5(8.27)

Redivision of snow layer and rediagnosis of temperature

When the snow water equivalent is updated, the snow-covered ratio is rediagnosed as describedin the section 8.1 and the mass of each layer is reallocated as described in the section 8.2. The
temperature in each redivided layer is rediagnosed so that the energy is conserved as follows:

new
TSn
=
(k)

old
KSn

X


old
f old Aold f(lold ∈knew ) TSn(l)∆Sn(l) Sn

.

new

new
f
(∆Sn
(k) ASn ).

(8.28)

l=1

It should be noted that the variables with the index “old” and “new” are those before andafter redivision, respectively. f(lold ∈knew ) is the ratio of the mass of the kth layer after redivision
to the mass of the lth layer before redivision.

43

8.4
8.4.1

Calculation of snow heat conduction
Snow heat conduction equations

The prognostic equation of the snow temperature due to snow heat conduction is as follows:f (k)cpi ∆Sn∗τTSn(k)− TSn(k)∆t= FeSn(k+1/2) − FeSn(k−1/2)(k = 1, . . . , KSn )(8.29)

with the heat conduction ﬂux FeSn given by(FSn(1/2) − ∆Fconv(1) )/ASn − ∆Fc,convTSn(k+1) − TSn(k)kSn(k+1/2)e∆zSn(k+1/2)
FSn(k+1/2) =
TSn(B) − TSn(k)kSn(k+1/2)∆zSn(k+1/2)(k = 0)(k = 1, ..., KSn − 1),

(8.30)(k = KSn )where kSn(k+1/2) is the snow heat conductivity, assigned the ﬁxed value of 0.3 W/m/K as astandard. The subscript k + 1/2 of the ﬂux represents the ﬂux from the (k + 1)th snow layer tothe upper one. ∆zSn(k+1/2) is the thickness of each snow layer, deﬁned by
f (1) /ρSn

0.5∆Sn
f (k) + ∆Snf (k+1) )/ρSn∆zSn(k+1/2) = 0.5(∆Snf (K ) /ρSn
0.5∆SnSn(k = 1)(k = 2, ..., KSn − 1) ,(8.31)(k = KSn )where ρSn is the snow density, assigned the ﬁxed value of 300kg/m3 as a standard. The snow
density and heat conductivity are considered to change over time due to compaction and changes
in properties (aging), but the eﬀect of such changes is not considered here.
In Eq. 8.30, the snow upper boundary ﬂux FeSn(1/2) is given using three energy variables: theheat conduction ﬂux from the snow to the ground surface solved in the ground surface energybalance FSn(1/2) , the ground surface energy convergence produced when the ground surface
temperature is solved by the snowmelt condition ∆Fconv , and the energy correction produced
when a change has occurred in the phase of the canopy water ∆Fc,conv . ∆Fconv is assumed to
be given only to the snow-covered portion, while ∆Fc,conv is given uniformly to the grid cells.
Since the sign of the ﬂux is taken as upward positive, the convergence has a negative sign.
In the equation for the snow lower boundary ﬂux FeSn(KSn +1/2) , TSn(B) is the temperature of
the snow lower boundary (the boundary surface of the snow and the soil). However, since theﬂux from the uppermost soil layer to the snow lower boundary is
Feg(1/2) = kg(1/2)

Tg(1) − TSn(B)∆zg(1/2).

(8.32)

There is assumed to be no convergence at the snow lower boundary, and TSn(B) is solved byputting
FeSn(KSn +1/2) = Feg(1/2) .

44(8.33)

When this is substituted into Eq. 8.30, the following is obtained:"∆zg(1/2) ∆zSn(KSn +1/2)
FeSn(KSn +1/2) =+kg(1/2)
kSn(KSn +1/2)
8.4.2

#−1
(Tg(1) − TSn(KSn ) ).

(8.34)

Case 1: When snowmelt does not occur in the uppermost layer

The implicit method is used to treat the temperature from the uppermost snow layer to thelowest snow layer, as follows:

∗
τ
FeSn
= FeSn+(k+1/2)(k+1/2)∂ FeSn(k+1/2)∂TSn(k)∆TSn(k) +∂ FeSn(k+1/2)∂TSn(k+1)∆TSn(k+1) ,(FSn(1/2) − ∆Fconv )/ASn − ∆Fc,conv(k = 0)





k


τ
τ
 Sn(k+1/2) (TSn(k+1)− TSn(k))(k = 1, ..., KSn − 1)τ∆zeSn(k+1/2),
FSn(k+1/2) =

#−1
"


∆zg(1/2) ∆zSn(KSn +1/2)τ+(Tg(1) − TSn) (k = KSn ) k(KSn )
kSn(KSn +1/2)
g(1/2)(8.35)

kSn(k+1/2)−(k = 1, ..., KSn − 1) ∆zSn(k+1/2)∂ FeSn(k+1/2)
"
#−1
=
,

∂TSn(k)∆zSn(KSn +1/2)∆z
g(1/2)
(k = KSn )+− kkSn(KSn +1/2
g(1/2)
(k = 0)0∂ FeSn(k+1/2)
kSn(k+1/2)=∂TSn(k+1)(k = 1, ..., KSn − 1) ∆z
Snk+1/2)and Eq. 8.29 is treated asf (k)cpi Sn∆TSn(k)t∗∗= FeSn− FeSn(k+1/2)(k−1/2)
τ= FeSn+(k+1/2)∂ FeSn(k+1/2)∂TSn(k)τ− FeSn−(k−1/2)∆TSn(k) +∂ FeSn(k−1/2)∂TSn(k−1)∂ FeSn(k+1/2)∂TSn(k+1)∆TSn(k−1) −∆TSn(k+1)∂ FeSn(k−1/2)∂TSn(k)(8.36)∆TSn(k)and solved by the LU factorization method as ∆TSn(k) (k = 1, ..., KSn ) simultaneous equations
with respect to KSn . At this juncture, it should be noted that the ﬂux at the snow upperboundary is ﬁxed as the boundary condition, and the snow lower boundary ﬂux is treatedexplicitly with regard to the temperature of the uppermost soil layer, which is the boundary
condition of the snow lower boundary. The snow temperature is updated by
∗
τ
TSn
= TSn
+ ∆TSn(k)
(k)
(k)

45

(8.37)

8.4.3

Case 2: When snowmelt occurs in the uppermost layer

When the temperature of the uppermost snow layer solved in case 1 is higher than 0 ◦ C, snowmelt
occurs in the uppermost snow layer. In this case, the temperature of the uppermost snow layer
is ﬁxed at 0 ◦ C. The ﬂux from the second snow layer to the uppermost snow layer is thenexpressed as
∗
Fe3/2
=

kSn(3/2)∆zSn(3/2)τ(TSn− Tmelt ) +(2)∂ FeSn(3/2)∂TSn(2)∆TSn(2)(8.38)and solved similarly to case 1 (when there is only one snow layer, the snow temperature issimilarly ﬁxed in the ﬂux from the soil to the snow).
The energy convergence used for melting in the uppermost snow layer is given by:∗f (1)∆Feconv = (Fe3/2− Fe1/2 ) − cpi Sn∗
Tmelt − TSn
(1)

∆t

.

(8.39)

Even if the temperature of the second snow layer and below is higher than Tmelt , the calculation is not iterated and the snowmelt is corrected accordingly.

8.5

Fluxes given to the soil or the runoﬀ process

The heat ﬂux given to the soil through the snow process is∗∗∗∆Fconv= ASn (∆Feconv− FeSnKSn ) − lm PSn,melt,(8.40)∗
where ∆Feconv
is the energy convergence remaining when all of the snow has melted, FeSnKSn∗
is the heat conduction ﬂux at the lowest snow layer, and PSn,meltis the snowfall that melts
immediately when it reaches the ground, deﬁned as∗
PSn,melt
=

(
0
∗
PSn

(Tg(1) ≤ Tmelt )(Tg(1) > Tmelt )(8.41)

Since the energy of the snow-free portion is given to the soil as it is, the energy correction
term due to the phase change of the canopy water is as follows:∗∆Fc,conv= (1 − ASn )∆Fc,conv .

(8.42)

The water ﬂux given to the runoﬀ process through the snow process is then expressed as
P rc∗∗ = (1 − ASn )P rc∗ ,∗∗+ PSn,melt,
P rl∗∗ = (1 − ASn )P rl∗ + ASn FewSn(8.43)∗
where FewSn
is the ﬂux of the rainfall or snowmelt water that has percolated through the lowest
snow layer.

46

8.6

Glacier formation

In this case, the maximum value is set for the snow water equivalent, and the portion exceeding
the maximum value is considered to become glacier runoﬀ:
Rogl = max(Sn − Snmax , 0)/∆t,
Sn = Sn − Rogl ∆t,f (K ) = ∆Sn
f (K ) − Rogl /ASn ∆t,∆SnSn(8.44)

Sn

where Rogl is the glacier runoﬀ. The mass of this portion is subtracted from the lowest snow
layer. Snmax is uniformly assigned the value of 1000kg/m2 as a standard.

8.7

Dust in snow

The amount of dust on the snow cover and in the snow layers are calculated in SUBROUTINEDSTCUT in matsnw.F.
8.7.1

Dust fall on the snow cover

The ﬂux of dust and black carbon (DBC), the light-absorpting particles, is calculated in SUBROUTINE MATSIRO in matdrv.F when the option OPT SNWALB is active.
Flux of DBC Dm and weighted-mean ﬂux of DBC Dw is obtained by
Dm = Ddust + DBC ,γdust Ddust + γBC DBC
Dw =,γdust + γBC(8.45)

where Ddust and DBC is the ﬂux of dust and black carbon, respectively. γc (c = dust, BC) isdeﬁned by
γc = γc,vis ωvis + γc,nir ωnir + γc,if r ωif r ,(8.46)

where γc,b is the absorption coeﬃcient for the band b, and ωb is the radiation weight of the bandb. Three bands of wavelength, visible (vis), near infrared (nir) and infrared (ifr) are consideredin MATSIRO, and they take the following values as a standard.
γdust,vis = 6.36777 × 102 , γdust,nir = 3.34617 × 102 , γdust,if r = 8.62054 × 102 ,γBC,vis = 7.43200 × 104 ,ωvis = 0.46,

γBC,nir = 2.93200 × 104 ,ωnir = 0.36,

γBC,if r = 2.47174 × 103 ,(8.47)

ωif r = 0.18

The amount of DBC on the uppermost snow layer is updated as
Mdτm+1(1) = Mdτm (1) + Dm ,τ
Mdτw+1(1) = Mdw (1) + Dw ,where Mdm (k) and Mdw (k) are the amount of DBC on the kth snow layer.

47(8.48)

8.7.2

Redistribution of dust

The amount of DBC in each layer is calculated in SUBROUTINE DSTCUT based on the resultsof snow layer recutting (SUBROUTINE SNWCUT). Note that this subroutine is applied for both
M dm and M dw , so they are represented by M d. Also, in this section, ρd(k) represents ρdm (k)and ρdw (k) , the density of dust and black carbon calculated using Mdm and Mdw , respectively.
The snow mass of kth layer after updating of snow mass but before snow layer recuttingτ +1/2∆Sn(k)
is calculated inτ +1/2∆Sn(k)+1= ∆Snτ(k) AτSn /AτSn(k = 1, 2, 3),(8.49)

where τ and τ + 1 represent before updating of snow mass and after recutting of snow layer,
respectively.
τ +1/2+1
When ∆Snτ(1)> ∆Sn(1) , the amount of DBC in the 1st layer increases due to increasein the snow mass in this layer. This is calculated asτ +1/2
τ +1/2
τ +1/2
τ +1


ρ
∆Sn
+
ρ
∆Sn
−
∆Sn
−
∆Sn
d(2)
d(3)

(2)
(1)
(1)
(2)






τ +1/2
τ +1/2+1∆Snτ(1)− ∆Sn(1)> ∆Sn(2)
τ +1
τ


Md(1)
− Md(1)=,(8.50)
τ +1/2

τ +1

ρ
∆Sn
−
∆Sn
 d(2)
(1)
(1)






τ +1/2
τ +1/2+1∆Snτ(1)− ∆Sn(1)≤ ∆Sn(2)where ρd(k) is the density of DBC in the kth layer.
τ +1/2+1
When ∆Snτ(1)≤ ∆Sn(1), the amount of DBC in the 1st layer decreases, and thus
τ +1/2
τ +1
+1
τ
Md(1)
− Md(1)= −ρd(1) ∆Sn(1) − ∆Snτ(1).

(8.51)

τ +1
+
−
τ
Md(1)
− Md(1)= Md(1)− Md(1),(8.52)

It leads to

where




τ +1/2+,0,Sn=ρd(2) min max Sn+
Md(1)
(1)
(2)




τ +1/2+ ρd(3) max max Sn+
,
0
−
Sn
,
0
,
(1)
(2)


−
Md(1)
= ρd(1) max −Sn+
,
0
,
(1)

τ +1/2τ +1
Sn+(1) = Sn(1) − Sn(1).

(8.53)(8.54)(8.55)

The change in the amount of DBC in the 3rd layer is determined similarly, and thus in the48

2nd layer it is calculated as follows:
τ +1
−
+
−
+
τ
Md(2)
− Md(2)= Md(1)− Md(1)+ Md(3)− Md(3).

(8.56)

Finally, the density of DBC on the kth snow layer ρd(k) is updated byρd(k) = Md(k) /∆Sn(k) × 106 .

8.8
8.8.1

(8.57)

Albedo of snow and ice
Albedo of snow

The albedo of the snow is calculated in SUBROUTINE SNWALB in matsnw.F.
The albedo of the snow is large in fresh snow, but becomes smaller with the passage of time
due to compaction and changes in properties as well as soilage. In order to take these eﬀects
into consideration, the albedo of the snow is treated as a prognostic variable.
The nondimensional age of snow at the time step of τ , Aτg , is formulated in

Aτg =

falb
,
1 − falb(8.58)

where
falb = min


τ −α
αvis
vis,new, 0.999 .
αvis,old − αvis,new(8.59)

αbτ is the albedo of the snow for band b at the time step of τ . Three bands of wavelength,
visible (vis), near infrared (nir) and infrared (ifr) are considered in MATSIRO, and here thefactors for visible band are used. αb,new is the albedo of newly fallen snow for band b and αb,oldis that of old snow. In default, αvis,new , αnir,new , αif r,new , αvis,old , αnir,old and αif r,old are set to
0.9, 0.7, 0.01, 0.65 (or 0.4 if the options OPT SNWALB is active), 0.2 and 0.1, respectively.
The age of snow at the next time step τ + 1 is, after Yang et al. (1997), assumed to be givenby the following equation:
10Aτg +1 = Aτg + (fage + fage
+ rdirt )∆tL /τage ,(8.60)

where



1
1
,
−
fage = exp fageT
Tmelt TSn(1)(8.61)

fageT = 5000, τage = 1 × 10 s, Tmelt = 273.15 K.
6

TSn(1) is the temperature of the ﬁrst layer of snow.
rdirt represents the eﬀect of DBC. When the option OPT SNWALB is inactive,(
rdirt =

rdirt,c (over continental ice)rdirt,0 (elsewhere),(8.62)

where rdirt,c = 0.01 and rdirt,0 = 0.3. When this option is active, the density of DBC is considered
49

as
(
rdirt =

min(rdirt,c + rdirt,s ρd(1) , 1000) (over continental ice)min(rdirt,0 + rdirt,s ρd(1) , 1000) (elsewhere),(8.63)

where rdirt,s is the DBC factor for slope with a constant value of 0.1 and ρd(1) is the density ofweighted DBC on the 1st snow layer.
Using this, the albedo of the snow at the time step of τ + 1, αbτ +1 , is solved byτ +1
αbτ +1 = αb,new+Aτg +1
1 + Aτg +1(αb,old − αb,new ),(8.64)

When snowfall has occurred, the albedo is updated to the value of the fresh snow in accordance with the snowfall: ∗
PSn ∆tL, 1 (αb,new − αbτ +1 ).
(8.65)
αbτ +1 = αbτ +1 + min∆Snc∆Snc is the snow water equivalent necessary for the albedo to fully return to the value of the
fresh snow.
8.8.2

Albedo of ice

The albedo of the ice sheet, αb,surf , is calculated in ENTRY ICEALB in matice.F.
This is expressed in a following function of the water content above the ice according to
Bougamont et al. (2005):
 w

surf
αb,surf = αb,wet − (αb,wet − αb,ice ) exp − ∗ ,w(8.66)

where αb,ice is the land ice albedo without surface water, αb,wet is the one with surface water,
wsurf is the thisness of surﬁce water and w∗ is the characteristic scale for surﬁcial water. b represents the three bands of wavelength, visible (vis), nearinfrared (nir) and infrared (ifr), similar
to ice albedo. In default, αvis,ice , αnir,ice and αif r,ice are set to 0.5, 0.3 and 0.05, respectively,and αb,wet is set to 0.15 for all bands.
where αb,ice is the land ice albedo without surface water, αb,wet is the one with surface water,
wsurf is the thisness of surﬁce water and w∗ is the characteristic scale for surﬁcial water. b represents the three bands of wavelength, visible (vis), nearinfrared (nir) and infrared (ifr), similar
to ice albedo. In default, αvis,ice , αnir,ice and αif r,ice are set to 0.5, 0.3 and 0.05, respectively,and αb,wet is set to 0.15 for all bands.

9

Runoﬀ

The surface runoﬀ and groundwater runoﬀ are solved using a simpliﬁed TOPMODEL (Beven andKirkby, 1979). The calculation of runoﬀ are solved in SUBROUTINE: MATROF in matrof.F,
and the related variables and parameters are introduced as follows:• Output variables

50

Variable

Description

Code

Units

Ro
Ro(k)
Pr

Total runoﬀs
Runoﬀ ﬂux from the k th soil layer
Water ﬂux given to the soil through the runoﬀ process

RUNOFF
RUNOFL
WINPT

kg/m2 /s
kg/m2 /s
kg/m2 /s

• Input variables

Variable

Description

Code

Units

P rc
P rl
Msn
Tg(k)
w(k)
θi(k)
(Rogl )∆t
Xs-

Convective rainfall
Nonconvective rainfall
Snow melt
Soil temperature
Soil moisture
Soil ice
Glacier formation
Time step
Surface condition
Soil type

WINPC
WINPL
SNMLT
GLG
GLW
GLFRS
GLACR
DELT
ILSFC
ILSOILkg/m2 /s
kg/m2 /s
(kg/m2 /s)
Km3 /m3m3 /m3(kg/m2 /s)s-• Internal work variables

Variable

Description

Code

Units

Ros
Roi
Roo
Rob
Rs
wsat(k)
Ks(k)
K0
Ks0
ψ(k)
z(x)
Asat
zf
-

Saturation excess runoﬀ
Inﬁltration excess runoﬀ
Surface storage overﬂow
Base runoﬀ
Surface runoﬀ (CMIP5)
Saturated soil moisture
Saturated hydrological conductivity from each layer
Saturation hydraulic conductivity at the ground surface
Saturation hydraulic conductivity at the depth of 2m
Matric potentiald(PSI)/d(W)Water table depth
Surface saturation fraction
Frost table level
Water table level
Flag for saturation

RUNOFS
RUNOFI
RUNOFORUNOFB
SRUNOF
GWS
DFWS
DFWST
GPSI
DPDW
WTABD
ASSAT
KFTAB
KWTAB
ISAT

kg/m2 /s
kg/m2 /s
kg/m2 /s
kg/m2 /s
kg/m2 /s
m3 /m3
m/s
m/s
m/s
m

• Internal parameters

51

m
m
m
-

52

Tangent of mean surface slope
Mean length of surface slope
Inﬂow rate into surface tank
Outﬂow rate from surface tank
Critical water table depth
Surface water storage
Maximum water table depth
Epsilon for soil moisture
Convective storm area for runoﬀ
Eﬀect of macropore on base-runoﬀ
Ground depth for groundwater runoﬀ calculatation

tanβs
Ls
α
τ

S
Acum
Ep
z”

1
f

Description

PARAMETERGRTANSGRLENSFRACSRTAUSSWTCRITSFCSTRWTBMAXEPSGWFRACCR
EFFMP
-

Code
0.333
0.001
0.001
0.1
1
2

Initial values
m
m
m2
m

Units

9.1

Outline of TOPMODELIn TOPMODEL, the horizontal distribution of a water table along the slope in a catchment
basin is considered. The downward groundwater ﬂow at a certain point on the slope is assumed
to be equal to the accumulated groundwater recharge in the upper part of the slope above thatpoint (quasi-equilibrium assumption). Then, the groundwater ﬂow must be greater in the lower
part of the slope. Under another assumption described later, for the groundwater ﬂow to begreater, the water table needs to be shallow. Thus, the distribution is derived such that thelower the slope, the shallower the water table. When the mean water table is shallower than acertain level, the water table rises to the ground surface at an area lower than a certain point in
the slope to form a saturated area. In this way, TOPMODEL is characterized by the mean watertable, the size of the saturated area, and the groundwater ﬂow velocity, which are importantconcepts for estimating the runoﬀ, being physically connected in a coherent manner.
TOPMODEL contains the following major assumptions:
1. The soil saturation hydraulic conductivity attenuates exponentially as the soil depth increase;
2. The gradient of the water table is considered to be locally the same as that of the slope;
3. The downward groundwater ﬂow at a certain point on the slope is equal to the accumulated
groundwater recharge in the upper slope above that point.
The usage of the symbols below is in accordance with the usual practice in descriptions of
TOPMODEL (Sivapalan et al., 1987; Stieglitz et al., 1997).
The ﬁrst assumption can be expressed as

Ks (z) = K0 exp(−f z)(9.1)

where Ks (z) is the soil saturation hydraulic conductivity at depth z, K0 is the saturation
hydraulic conductivity at the ground surface which diﬀers among diﬀerent soil types, and f is
the attenuation coeﬃcient.
When the depth of the water table at a certain point i is designated as zi , the downward
groundwater ﬂux on the slope at that point qi is
Z Z
qi =

Ks (z)dz · tan β =ziK0tan β[exp(−f zi ) − exp(−f Z)]
f

(9.2)

where β is the gradient of the slope because of the second assumption. Z is the depth of
the impervious surface. Normally, however, the term exp(−f Z) can be ignored because Z is
assumed to be suﬃciently large compared with f1 . Moreover, the slope direction soil moisture
ﬂux in the unsaturated zone above the water table is small and thus also ignored.
If the groundwater recharge rate R is assumed to be horizontally uniform, with the thirdassumption, the above equation is expressed as:

aR =

K0
tan β exp(−f zi )
f

(9.3)

where a is the total upstream area (per unit contour line length at point i with respect topoint i.

53

When this is solved for zi , the following is obtained:1zi = − ln
f



f aR
K0 tan β(9.4)

The averaged water table depth z in domain A is described below, while the grid averaged
water table depth will be mentioned in section 9.3.1:
1
z=
A

Z

1
Λ≡
A

zi dA = −Λ −
A



Z
ln
A

fa
K0 tan β

1
ln R
f

(9.5)


dA

(9.6)

Therefore, the recharge rate R can then be expressed as a function of the mean water table
depth z as follows:

R = exp(−f z − Λ)

(9.7)

Under the third assumption, this is all of the groundwater runoﬀ discharged from domain A.
Next, if R is substituted into 9.4 , the following relationship of zi and z is obtained:
 


1
fa
zi = z −ln−Λf
K0 tan β(9.8)

The domain that satisﬁes zi ≤ 0 is the surface saturated area.

9.2

Application of TOPMODEL assuming simpliﬁed topography

Normally, when TOPMODEL is used, detailed topographical data on the target area is required.
Here, however, the average shape of the slope in a grid cell is roughly estimated from the average
inclination and the standard deviation of the altitude in the grid (this estimation method is
temporary at this stage, and further study is required).
The topography in the grid cell is represented by the slope with uniform gradient βs and the
distance from the ridge to valley Ls .
Ls is estimated using the standard deviation of altitude (σz ) as follows:
√
2 3σz
Ls =
tan βs(9.9)√
where 2 3σz is the altitude diﬀerence between the ridge and valley in serrate topographysuch that the standard deviation of altitude is σz .
The x-axis is taken from the ridge toward the valley on the horizontal surface. Then, thetotal upstream area at point x is x, and 9.4 becomes54

1
z(x) = − ln
f



f xR
K0 tan βs(9.10)

Using this, from 9.5 the grid average water table depth is

1
z=
Ls

Z Ls
0



 
f Ls R1ln−1z(x)dx = −f
K0 tan βs(9.11)

from 9.7 the groundwater recharge rate isR=

K0 tan βsexp(1 − f z)
f Ls

(9.12)

and from 9.8, the relationship between the water table depth at point x and the grid averagewater table depth is1
z(x) = z −
f



x
ln
+1
Ls

(9.13)

If z(x) ≤ 0 is solved for x, the following are obtained:x ≥ x0(9.14)x0 = Ls exp(f z − 1)

(9.15)

Therefore, the fraction of the saturated area is solved as

Asat =

Ls − x0= 1 − exp(f z − 1)Ls

(9.16)

However, when Asat ≥ 0 and z > f1 , no saturated area exists.

9.3

Calculation of runoﬀ

Four types of runoﬀ mechanisms are considered, and the total of the runoﬀs Ro by each mechanism is assumed to be the total runoﬀ from the grid cell:

Ro = Ros + Roi + Roo + Rob(9.17)

where Ros is the saturation excess runoﬀ (Dunne runoﬀ), Roi is the inﬁltration excess runoﬀ(Horton runoﬀ), Roo is the overﬂow of the uppermost soil layer, and Rob is the groundwater
runoﬀ. The ﬁrst three types are classiﬁed as the surface runoﬀ.
55

The surface runoﬀ Rs calculated by MATSIRO will be:
Rs = Ro − Rob = Ros + Roi + Roo(9.18)

However, when taking snow-fed wetland into account (Nitta et al., 2017), part of the surface
runoﬀ Rs will be stored in the surface storage and part of runoﬀ to rivers will be delayed, whichleads to an increase in land surface wetness and hence evaporation in water-limited regimes.
Then surface runoﬀ Rs will become:Rs = (Ros + Roi + Roo )α(9.19)

here α determines the inﬂow rate into surface tank and is speciﬁed in Wetland section. Pleaserefer to Eq.12.3 for more details.
9.3.1

Estimation of grid average water table depth

The grid average water table is considered to exist in the lowest unsaturated layer, which isexpressed as kW T . The grid average water table depth (z) is estimated by:z = zg(kW T − 1 ) − ψkW T(9.20)

2

where ψkW T is the matric potential in the kW T th soil layer.
The above equation is equivalent to considering the moisture potential on the upper boundary of the unsaturated layer as ψkW T , which denotes that soil moisture distribution reaches
equilibrium state (i.e., the state in which gravity and the capillary force are in equilibrium).
Under unsaturated condition that −ψkW T exceeds soil layer thickness, water table will generate at the lower boundary of soil layer.
When z > zg(kW T + 1 ) , i.e., average water table depth is deeper than the lower boundary of2
kW T th layer, in case kW T is the lowest soil layer, the water table is assumed to not exist; when
kW T is not the lowest soil layer, the layer below (the uppermost layer among the saturated
layers) is assumed to be kW T and water table will generate at zg(kW T + 1 ) .
2
When there is a frozen soil surface in the middle of the soil, estimation of the water table
depth is performed from above the frozen soil surface.
9.3.2

Calculation of groundwater runoﬀ

Because of the quasi-equilibrium assumption, the groundwater runoﬀ is equal to the groundwater
recharge rate in 9.12. In MATSIRO6, Hirabayashi (2004) changed K0 to Ks0 in groundwater
runoﬀ calculation, which denotes a saturation hydraulic conductivity at depth of 2m:
Ks0 = exp(z ” f )K0 Ep(9.21)

where z ” is the depth of 2m, Ep denotes the eﬀect of macropore on groundwater runoﬀ. It isalso worth noting that the value of f1 has changed from 0.6 to 0.33 in MATSIRO6. Therefore,
calculation of groundwater runoﬀ will become:Rob =Ks0 tan βsexp(1 − f z)
f Ls

(9.22)

However, when a frozen soil surface exists under the water table, referring to the case with56

the term exp(−f Z) in 9.2, it is assumed thatRob =Ks0 tan βs[exp(1 − f z) − exp(1 − f zf )]f Ls(9.23)

where zf is the depth of frozen soil surface. Although other relations in TOPMODEL should
also be changed in such a case, they are not changed here for the sake of simpliﬁcation.
When there is an unfrozen layer under the frozen soil surface and a water table exists, thegroundwater runoﬀ from there is added by a similar calculation.
The water content from the groundwater runoﬀ is removed from the kW T th soil layer:

Ro(kW T ) = Rob(9.24)

where Ro(k) denotes the runoﬀ ﬂux from the kW T th soil layer.
9.3.3

Calculation of surface runoﬀ

All of the rainfall that falls on the surface saturated area runs oﬀ as is (saturation excess runoﬀ):
Ros = (P rc∗∗ + P rl∗∗ )Asat(9.25)

The fraction of the surface saturated area Asat is given by 9.16. Here, the correlation between
the rainfall distribution of the subgrid and topography is ignored.
With regard to rainfall that falls on the surface unsaturated area, only the portion thatexceeds the soil inﬁltration capacity runs oﬀ (inﬁltration excess runoﬀ). The soil inﬁltration
capacity is given by the saturation hydraulic conductivity of the uppermost soil layer for simpliﬁcation. The convective precipitation is considered to fall locally, and the fraction of the
precipitation area (Ac ) is assumed to be uniform (0.1 as a standard value). The stratiform
precipitation is also assumed to be uniform.

Roci = max(

P rc∗∗+ P rl∗∗ − Ks(1) , 0)(1 − Asat )Ac(9.26)∗∗
Ronc
i = max(P rl − Ks(1) , 0)(1 − Asat )(9.27)

Roi = Ac Roci + (1 − Ac )Ronc
i(9.28)

where Roci and Ronc
i are Roi in the convective precipitation area and nonconvective precipitation area, respectively; and Ks(1) is the saturation hydraulic conductivity in the uppermost
soil layer.
The overﬂow of the uppermost soil layer, allowing a small amount of ponding wstr (1 mm asa standard value), is assumed to be57

Roo =

max(w(1) − wsat(1) − wstr , 0)ρw ∆zg(1)∆tL(9.29)

This portion is subtracted from the uppermost soil layer later, and therefore should be
remembered as the runoﬀ from the uppermost layer, as follows.

Ro(1) = Ro(1) + Roo

9.4(9.30)

Water ﬂux given to soil

The water ﬂux given to the soil through the runoﬀ process isPr∗∗∗ = P rc∗∗ + P rl∗∗ − Ros − Roi(9.31)

If the wetland scheme is active for the grid, the outﬂow of the wetland should also be addedin this equation (see Eq. 12.4).

10

Soil

The soil temperature, the soil moisture, and the frozen soil are calculated next.

10.1

Calculation of soil heat conduction

10.1.1

Soil heat conduction equations

The prognostic equation for the soil temperature by soil heat conduction is
Cg(k)

∗
τ
Tg(k)
− Tg(k)∆tL= Fg(k+1/2) − Fg(k−1/2)(k = 1, . . . , Kg )(10.1)

with Cg(k) , the soil heat capacity, deﬁned by
Cg(k) = (cg(k) + ρw cpw w(k) )∆zg(k)(10.2)

where cg(k) is the speciﬁc heat of the soil, and is given as a parameter for each soil type; cpw
is the speciﬁc heat of the water; w(k) is the soil moisture (volumetric moisture content); and∆zg(k) is the thickness of the k th soil layer. When including the heat capacity of the soil
moisture in the soil heat capacity in this way, unless the heat transfer accompanying the soil
moisture movement is considered, the energy is not conserved. The heat transfer accompanying
the soil moisture movement is not considered in the MATGND soil submodel at present, and its
introduction is under study. However, it should be noted that unless the heat capacity of suchelements as vapor in the atmosphere, rainfall, etc. is considered, the conservation of energy is
disrupted in certain respects.

58

The heat conduction ﬂux Fg is given by∗∗Fg(1/2) − ∆Fconv− ∆Fc,conv(k = 0)
Tg(k+1) − Tg(k)kg(k+1/2)(k = 1, . . . , Kg − 1)
Fg(k+1/2) =∆zg(k+1/2)0(k = Kg )(10.3)

with kg(k+1/2) , the soil heat conductivity, expressed askg(k+1/2) = kg0(k+1/2) [1 + fkg tanh(w(k) /wkg )](10.4)

where kg0(k+1/2) is the heat conductivity when the soil moisture is 0, and fkg = 6 and wkg = 0.25
are constants.
∆zg(k+1/2) is the thickness between the soil temperature deﬁnition points of the k + 1th layer
and the k + 1th layer (when k = 0, the thickness between the uppermost layer temperaturedeﬁnition point and the soil upper boundary; when k = Kg , the thickness between the lowest
layer temperature deﬁnition point and the soil lower boundary).
In 10.3, the value given to the soil upper boundary condition (Fg(1/2) ) is the value obtainedat the time of solving the ground surface energy balance, with the addition of the energy convergence at the snow lower boundary (including the heat conduction ﬂux at the snow lower
boundary) as well as the allotment to the snow-free portion of the energy correction term due
to phase change of the canopy water. The ﬂux takes an upward (positive) direction, so when
the amount of convergence is added it has a negative sign. The soil lower boundary condition
Fg(Kg +1/2) is assumed to be zero ﬂux.
10.1.2

Solution of heat conduction equations

These equations are solved using the implicit method with regard to the soil temperature fromthe uppermost layer to the lowest layer. That is, for k = 1, . . . , Kg − 1, the heat conduction ﬂux
is expressed as∗τFg(k+1/2)= Fg(k+1/2)+τFg(k+1/2)=∂Fg(k+1/2)∂Fg(k+1/2)∆Tg(k) +∆Tg(k+1)∂Tg(k)∂Tg(k+1)(10.5)

kg(k+1/2)τ(T τ− Tg(k))∆zg(k+1/2) g(k+1)(10.6)∂Fg(k+1/2)
kg(k+1/2)=−∂Tg(k)∆zg(k+1/2)(10.7)∂Fg(k+1/2)
kg(k+1/2)=∂Tg(k+1)∆zg(k+1/2)(10.8)

59

and 10.1 is treated asCg(k)∆Tg(k)∆tL∗∗= Fg(k+1/2)− Fg(k−1/2)∂Fg(k+1/2)∂Fg(k+1/2)∆Tg(k) +∆Tg(k+1)∂Tg(k)∂Tg(k+1)∂Fg(k−1/2)∂Fg(k−1/2)τ− Fg(k−1/2)−∆Tg(k−1) −∆Tg(k)∂Tg(k−1)∂Tg(k−1)τ= Fg(k+1/2)+(10.9)and solved by the LU factorization method as Kg simultaneous equations with respect to∆Tg(k) (k = 1, . . . , Kg ). At this juncture, it should be noted that the equations are solved with
the ﬂuxes at the soil upper boundary and lower boundary ﬁxed as the boundary conditions:
∗
τ
Tg(k)
= Tg(k)+ ∆Tg(k)(10.10)The soil temperature is partially updated by the above equation. By this, as well as throughcorrection of the phase change in the soil moisture mentioned later, the soil temperature iscompletely updated.

10.2

Calculation of soil moisture movement

10.2.1

Soil moisture movement equations

The equation for soil moisture movement (Richards equation) is given by
ρw

τ +1
τ
w(k)
− w(k)=

Fw(k+1/2) − Fw(k−1/2)+ Sw(k)∆zg(k)(k = 1, . . . , Kg )(10.11)


P r∗∗∗ − Et
 (k = 0) (1,1)ψ(k+1) − ψ(k)−1(k = 1, . . . , Kg − 1)K(k+1/2)
Fw(k+1/2) =∆zg(k+1/2)0(k = Kg )(10.12)∆tL

The soil moisture ﬂux Fw is given by

in which K(k+1/2) is the soil hydraulic conductivity that, referring to Clapp and Hornberger(1978), is expressed asK(k+1/2) = Ks(k+1/2) (max(W(k) , W(k+1) ))2b(k)+3 fi(10.13)where Ks(k+1/2) is the saturation hydraulic conductivity and b(k) is the index of the moisturepotential curve, which are given as external parameters for each soil type. W(k) is the degree ofsaturation considered excluding the frozen soil moisture, given byW(k) =

w(k) − wi(k)wsat(k) − wi(k)(10.14)

where wsat(k) is the porosity of the soil, which is also given as a parameter for each soil type.
fi is a parameter that denotes that soil moisture movement is suppressed by the existence of60

frozen soil. Although further study of this point is required, at present it is given byfi = 1 − Wi(k)

1 − Wi(k+1)

(10.15)where Wi(k) = wi(k) /(wsat(k) − wi(k) ).
The soil moisture potential ψ is given as follows from Clapp and Hornberger:−b(k)ψ(k) = ψs(k) W(k)(10.16)

where ψs(k) is given as an external parameter for each soil type.
In 10.11, Sw(k) is a source term which, considering the root uptake and the runoﬀ, is givenbySw(k) = −Froot(k) − Ro(k)(10.17)

In 10.12, the soil upper boundary condition Fw(1/2) is the diﬀerence between the moisture
ﬂux through the runoﬀ process (P ∗∗∗ ) and the evaporation ﬂux from the soil (Et(1,1) ). Separately from this, the sublimation ﬂux portion is subtracted from the frozen soil moisture of theuppermost layer before calculation of the soil moisture movement:

10.2.2ττwi(k)= wi(k)− Et(2,1) ∆tL /(ρ∆zg(1) )(10.18)ττw(k)= w(k)− Et(2,1) ∆tL /(ρ∆zg(1) )(10.19)

Solution of soil moisture movement equations

These equations are solved by using the implicit method for the soil moisture from the uppermost
layer to the lowest layer. For k = 1, . . . , Kg − 1, the soil moisture ﬂux isτ +1τFw(k+1/2)= Fw(k+1/2)+∂Fw(k+1/2)∂Fw(k+1/2)∆w(k) +∆w(k+1)∂w(k)∂w(k+1)τFw(k+1/2)= K(k+1/2)ττ
ψ(k+1)− ψ(k)∆zg(k+1/2)(10.20)!
−1(10.21)ψs(k)∂Fw(k+1/2)
K(k+1/2)−b(k)−1
−b(k)
=−
W
∂w(k)
∆zg(k+1/2)
wsat(k) − wi(k) (k)(10.22)ψs(k+1)∂Fw(k+1/2)
K(k+1/2)−b(k)−1−b(k)=W∂w(k+1)∆zg(k+1/2)
wsat(k+1) − wi(k+1) (k+1)(10.23)and 10.11 is treated as61ρw ∆zg(k)

∆w(k)
∆tL

τ +1
τ +1
= Fw(k+1/2)− Fw(k−1/2)+ Sw(k) ∆zg(k)∂Fw(k+1/2)∂Fw(k+1/2)∆w(k) +∆w(k+1)(10.24)∂w(k)∂w(k+1)∂Fw(k−1/2)∂Fw(k−1/2)τ− Fw(k−1/2)−∆w(k−1) −∆w(k) + Sw(k) ∆zg(k)∂w(k−1)∂w(k−1)

τ= Fw(k+1/2)+

and solved by the LU factorization method as Kg simultaneous equations with respect to∆Tg(k) (k = 1, . . . , Kg ). At this juncture, it should be noted that the equations are solved with
the ﬂuxes at the soil upper boundary and lower boundary ﬁxed as the boundary conditions.
The soil moisture is updated by

τ +1
τ
w(k)
= w(k)+ ∆w(k)(10.25)

As a result of this calculation, if a part appears where the soil moisture become supersaturated, it is adjusted in the vertical direction to eliminate supersaturation. The reason why sucha supersaturated portion is not considered as runoﬀ is that this supersaturation is artiﬁcially
produced because the vertical movement of the soil moisture is solved without saturation data.
First, from the second soil layer downwards, the supersaturated portion of the soil moisture isgiven to the layer below. Next, from the lowest soil layer upwards, the supersaturated portionof the soil moisture is given to the next layer up. With this operation, when the soil moistureis large enough, a saturated layer around the lowest soil layer is formed and the water table ofEq. (278) can be deﬁned.

10.3

Phase change of soil moisture

As a result of calculating the soil heat conductivity, when the temperature in the layer containing
liquid water is lower than Tmelt = 0◦ C, or when the temperature in the layer containing solid
water is higher than Tmelt , the phase change of the soil moisture is calculated. If the amount offreeze (adjustment portion) of the soil moisture in the kth layer is assumed to be ∆wi(k) ,τ +1τ
∗− wi(k)> 0 (frozen):
when Tg(k)< Tmelt and w(k)∆wi(k) = min∗ )Cg(k) (Tmelt − Tg(k)lm ρw ∆zg(k)

!
τ +1
τ
, w(k)− wi(k)(10.26)∗τ
when Tg(k)> Tmelt and wi(k)> 0 (melting):∆wi(k) = max∗ )Cg(k) (Tmelt − Tg(k)lm ρw ∆zg(k)!
τ, −wi(k)The frozen soil moisture and the soil moisture are then updated as follows:

62(10.27)

10.3.1

τ +1
τ
wi(k)
= wi(k)+ ∆wi(k)(10.28)

τ +1
∗
Tg(k)
= Tg(k)+ lm ρw ∆zg(k) ∆wi(k) /Cg(k)(10.29)Ice sheet process

When the land cover type is ice sheet, if the soil temperature exceeds Tmelt , it is returned toTmelt :
τ +1
∗
Tg(k)
= min(Tg(k), Tmelt )(10.30)

The rate of change of the ice sheet Fice is then diagnosed as

Fice = −Et(2,1) −11∗
Cg(k) max(Tg(k)− Tmelt , 0)lm ∆tL

(10.31)Lake

Lake is treated in MATSIRO (lakesf.F, lakeic.F, and lakepo.F), as well as land.
The method for surface ﬂux calculation (section 11.1 and 11.2) is derived from the sea surface
scheme of MIROC-AGCM, while the calculation for the processes below the lake ice (section 11.3
and 11.4) is derived from the ocean model COCO (COCO-OGCM). lakeic.F and lakepo.F werebased on the COCO-OGCM, and the ENTRY statement are used so as to keep the structure of
the original programs. For practical use, note, for example, that the unit of temperature is K
until ﬂux calculation (section 11.1 and 11.2), while it is ◦ C after the ice and inter lake (section
11.3 and 11.4). It is also noted that because the second half part of this section is based on theold version of COCO, hence it is slightly diﬀerent from the MIROC6-AOGCM(Tatebe et al.,2019) and Hasumi (2015).
Dimensions of the lake scheme is deﬁned in include/zkg21c.F. KLMAX is the number of verticallayers set to 5 in MIROC6/MATSIRO6. NLTDIM is the number of tracers, 1:temperature 2:salt.
Since the vertical layers are actually from the ﬁrst number of vertical layer (KLSTR=2) to thelast number of vertical layer (KLEND=KLMAX+1). NLZDIM = KLMAX+KLSTR is also deﬁned for the
management.
Lakes have been calculated as one of the tiles in a grid, as described in section 13.1. Asa standard value, minimum depth of lake of 10 × 102 [cm] is given in matdrv.F, hence any lakes
cannot disappear even in severe dry conditions.

11.1

Calculation of lake surface conditions

In ENTRY LAKEBC (in SUBROUTINE LSFBCS of lakesf.F) lake surface albedo, roughness,and heat ﬂux are calculated. They are calculated supposing ice-free conditions, then modiﬁed
to take into account the eﬀects of ice and snow cover. While the albedo of land snow is aprognostic variable, the lake surface albedo considering with ice and snow above is a diagnostic
variable. The aging eﬀect of the snow is simply depending on its temperature. These methods
are actually same with an old version of COCO-OGCM(Hasumi, 2015). The newest version
of COCO, which is going to be coupled to MIROC7-AOGCM, has been applied a melt pond

63

scheme a snow aging scheme which is basically the same with the treatment in the current land
surface (Komuro, the GCM meeting on 22nd Feb, 2021).
• Output variables of LAKEBC

Variable

Variable in source code

Longname

Unit

αLk
zLk0
G

GRALB
GRZ0
FOGFLXDFGT

Lake surface albedo
Lake surface roughness
Heat ﬂux
Heat diﬀusion coeﬃcient[−][m][W/m2 ]
[W/m2 /K]∂G
∂T

• Input variables of LAKEBC

Variable

Variable in source codes

Longname

Unit

Ts
Tb
Ic
SnLk
RIcLk
ua
va
cosζ

GRTS
GRTB
GRICE
GRSNW
GRICR
GDUA
GDVA
RCOSZ

Surface temperature
Ice base temperaturelake ice amountsnow amountlake ice concentration
Atmospheric 1st layer east-west wind
Atmospheric 1st layer north-south windcos(solar zenith)[K][K][kg/m2 ][kg/m2 ]
[−]
[m/s]
[m/s]
[−]

First, the lake albedo αLk(d,b) is calculated. b = 1, 2, 3 represent the visible, near-infrared,
and infrared wavelength bands, respectively. Also, d = 1, 2 represents direct and scattered light,
respectively. The albedo for the visible bands are calculated in SUBROUTINE LAKEALB,
supposing ice-free conditions. The albedo for near-infrared is set to same as that for the visible
bands. The albedo for infrared is uniformly set to a constant value, as mentioned in section
11.1.1.
When lake ice is present, the albedo is modiﬁed to take into account the ice cover on lakethat can change the albedo.
′
αLk= αLk + (αIcLk − αLk )RIcLk(11.1)

where αIcLk is the lake ice albedo, and RIcLk is the lake ice concentation, respectively. In
addition, snow cover on a lake ice that can change the albedo is also considered. Assuming thatthe snow albedo depends on the surface temperature, we can calculate a function F below.

F (Ts ) =min
Ts − Tmmax − T min
Tmm,(0 ≤ F (Ts ) ≤ 1)

(11.2)

min and T max are the minimum and the minimum
where Ts is the surface temperature, and Tmmtemperature for the albedo change, respectively.
Then, the albedo can be modiﬁed by

64

′′
αLk
= αLk(1,b) + (αSnLk(2,b) − αSnLk(1,b) )F (Ts )(11.3)′′′′
αLk = αLk
+ (αLk− αLk)RSnLk

(11.4)

where the αSnLk(d,b) is the snow albedo covering the lake, and RSnLk is the snow coverage,
respectively.
Second, the lake surface roughness for momentum, heat and vapor are calculated in SUBROUTINE LAKEZ0F, based on Miller et al. (1992), same with COCO-OGCM (Hasumi, 2015),
supposing the ice-free conditions, then modiﬁed to take into account the eﬀects of ice and snow
cover.
When lake ice is present, each roughness is modiﬁed to take into account the lake ice concentration (RIcLk )′= zLk0 + (zLk0,ice − zLk0 )RIcLk
zLk0(11.5)

where zLk0 is surface roughness, and zLk0,ice is surface roughness for lake ice, respectively.
Then, taking into account the snow coverage (RSnLk ), we can express it as′′)RSnLk+ (zLk0,ice − zLk0
zLk0 = zLk0

(11.6)

Third, the heat ﬂux is considered with the temperature diﬀerences between the snow surface
and the ice bottom, because the diﬀerence should be zero in the ice-free conditions.
If the lake ice exists, the heat diﬀusion coeﬃcient is described as ∂G 
∂T

IcLk

=

DIcLk
RIcLk

(11.7)

where DIcLk is the coeﬃcient of lake ice.
If the snow exists, the heat diﬀusion coeﬃcient of snow covered area is ∂G ∂T

SnLk=

DIcLk DSnLk
DIcLk RSnLk + DSnLk RIcLk(11.8)

where DSnLk is the coeﬃcient of snow.
Therefore, the net heat diﬀusion coeﬃcient is ﬁnally ∂G 
∂G  ∂G (1 − RSnLk ) +
RSnLk
=
∂T
∂T IcLk∂T SnLk(11.9)The temperature diﬀerences between the snow surface (TS ) and the ice bottom (TB ) is saved
as heat ﬂux(G).

65

G=
11.1.1∂G(TB − TS )∂T

(11.10)Calculation of lake surface albedo

Processes of SUBROUTINE LAKEALB is the same with SUBROUTINE [SEAALB (in p-sfc.F),
which is the sea surface ﬂux scheme of MIROC6-AGCM(Tatebe et al., 2019). For lake surface
level albedo αLk,L(d) , d = 1, 2 represents direct and scattered light, respectively.
Using the solar zenith angle ζ (cosζ), the albedo for direct light is presented by∗∗

αLk,L(1) = e(C3 A +C2 )A +C1(11.11)

where A∗ = min(max(cosζ, 0.03459), 0.961), and C1 , C2 , C3 is −0.7479, −4.677039, 1.583171
respectively.
On the other hand, the albedo for scattered light is uniformly set to a constant parameter.

αLk,L(2) = 0.06
11.1.2(11.12)

Lake surface roughness

Contents of SUBROUTINE LAKEZ0F is the same with SUBROUTINE SEAZ0F (of pgocn.F),
which is the sea surface scheme of MIROC6-AGCM(Tatebe et al., 2019).
The roughness variation of the lake surface is determined by the friction velocity u⋆ .

u⋆ =qCM0 (ua 2 + va 2 )(11.13)

We perform successive approximation calculation of CM0 , because the variables of the atmosphseric ﬁrst layer (ua , va , Ta , qa ) are required.

zLk0,M = z0,M0 + z0,MR +z0,MR u⋆ 2 z0,MS ν+gu⋆(11.14)zLk0,H = z0,H0 + z0,HR +z0,HR u⋆ 2 z0,HS ν+gu⋆(11.15)zLk0,E = z0,E0 + z0,ER +z0,ER u⋆ 2 z0,ES ν+gu⋆(11.16)

where, ν = 1.5 × 10−5 [m2 /s] is the kinetic viscosity of the atmosphere, z0,M , z0,H and z0,E
are surface roughness for momentum, heat and vapor (0.0, 0.018, and 0.11 as standard values),
and z0,M0 , z0,H0 and z0,E0 are base of them (1.4 × 10−5 , 0.0, and 0.4 as standard values), and66

z0,MR , z0,MR and z0,ER (1.3 × 10−4 , 0.0, and 0.62 as standard values) are rough factor for them,and z0,MS , z0,MS and z0,ES are smooth factor for them, respectively.

11.2

Solution of energy balance at lake surface

In SUBROUTINE [LAKEHB] (of lakesf.F), the energy balance at lake surface is solved.
• Output variables of LAKEHB

Variable

Variable in source code

Longname

Unit

Wf ree/ice
LW ↑F

WFLUXSRFLXLUSFLXBL

Surface water ﬂux
Upward long wave
Flux balance[W/m2 ]
[W/m2 ]
[W/m2 ]• Input variables of LAKEHB

Variable

Variable in source code

Longname

Unit

∂H
∂Ts
∂E
∂Ts
∂G
∂Ts

DTFDS
DQFDS
DGFDS
RFLXSD
RFLXLURFLXLDGRALBL
GRICR

Sensible heat ﬂux coeﬃcient
Latent heat ﬂux coeﬃcient
Surface heat ﬂux coeﬃcient
Downward SW radiation
Upward SW radiation
Downward LW radiation
Lake surface albedo
Lake ice concentration[W/m2 /K]
[W/m2 /K]
[W/m2 /K]
[W/m2 ]
[W/m2 ]
[W/m2 ]
[−]
[−]

SW ↓
SW ↑
LW ↓
αLk
RIcLk

The comments for some variables say “soil”, but this is because the program was adaptedfrom a land surface scheme, and has no particular meaning.
Downward radiative ﬂuxes are not directly dependent on the condition of the lake surface,
and the atmospheric values are simply used in the lake model. Shortwave emission from the lake
surface is negligible, so the upward part of the shortwave radiative ﬂux is accounted for solely
by reﬂection of the incoming downward ﬂux. The upward shortwave radiative ﬂux (SW ↑ ) isrepresented bySW ↑ = −αLk,SW SW ↓(11.17)

where SW ↓ is the downward shortwave radiation ﬂux, and αLk,SW is lake surface albedo forshortwave radiation in the ice-free area, respectively. On the other hand, the upward longwave
radiative ﬂux has both reﬂection of the incoming ﬂux and emission from the lake surface. The
upward shortwave radiative ﬂux is represented by
LW ↑ = −αLk LW ↓ + ϵσTs4(11.18)where αLk be the lake surface albedo for longwave radiation, ϵ be emissivity of the lake67

surface relative to the black body radiation, σ is the Stefan-Boltzmann constant, and Ts issurface temprature, respectively. If lake ice exists, snow or lake ice temperature is consideredby fractions. When radiative equilibrium is assumed, emissivity becomes identical to co-albedo:

ϵ = 1 − αLk(11.19)F ∗ = H + (1 − αLk )σTs4 + αLk LW ↓ − LW ↓ + SW ↑ − SW ↓(11.20)The net surface ﬂux is presented by

The heat ﬂux into the lake surface is presented, with the surface heat ﬂux (G) calculated in
SUBROUTINE [SFCFLX] (in matdrv.F).
G∗ = G − F ∗(11.21)

where G∗ is the net incoming ﬂux (the opposite direction with F ∗ ).
The temperature derivative term is
∂G∗
∂G
∂H
∂R
=
+
+
∂Ts
∂Ts ∂Ts ∂Ts(11.22)

When the lake ice exists, the sublimation ﬂux (ls E) is considered
GIcLk = G∗ − ls E

(11.23)

The temperature derivative term is∂GIcLk
∂G∗
∂E
=
+ ls
∂Ts
∂Ts
∂Ts

(11.24)

Finally, we can update the surface temperature with the lake ice concentration with ∆Ts =IcLk −1
GIcLk ( ∂G∂T)s
Ts = Ts + RIcLk ∆Ts(11.25)

Then, the sensible heat ﬂux (EIcLk ) and latent heat ﬂux (EIcLk ) on the lake ice is updated.

EIcLk = E +

∂E
∆Ts
∂Ts

(11.26)

HIcLk = H +

∂H
∆Ts
∂Ts

(11.27)68

When the lake ice does not existed, otherwise, the evaporation ﬂux (lc E)is added to the net
ﬂux.
Gf reeLk = F ∗ + lc E

(11.28)

Finally each ﬂux is updated.
For the sensible heat ﬂux (H), the temperature change on the lake ice is considered.

H = H + RIcLk HIcLk(11.29)where HIcLk is the sensible heat ﬂux on the lake ice. Then, the heat used for the temperaturechange is saved as:

F = RIcLk HIcLk(11.30)For the upward longwave radiative ﬂux (LW ↑ ), the temperature change on the lake ice isconsidered.
LW ↑ = LW ↑ + 4σ
RIcLk ∆TsTs

(11.31)For the surface heat ﬂux, the lake ice concentration is considered.

G = (1 − RIcLk )Gf reeLk + RIcLk GIcLk(11.32)For the latent heat ﬂux, the lake ice concentration is considered.

E = (1 − RIcLk )E + RIcLk EIcLk(11.33)

Each term above are saved as freshwater ﬂux.

11.3

Wf reeLk = (1 − RIcLk )E(11.34)

WIcLk = RIcLk EIcLk(11.35)

Calculation of lake ice

In this section, the lake ice calculation is described. There are three prognostic variables in thelake ice model described herein: lake ice concentration AI , which is area fraction of a grid coveredby lake ice and takes a value between zero and unity; mean lake ice thickness hI over ice-covered
part of a grid; mean snow depth hS over lake ice. Horizontal ﬂow of ice is not considered in
the lake parts, diﬀerently from the COCO-OGCM. AI , hI and hS are incrementally modiﬁed as69

the prognostic variables. The variables representing the surface state are represented with thesubscript k = 2, since the top of the vertical layers is actually deﬁned as KLSTR=2.
The model also calculates temperature at snow top (lake ice top when there is no snow
cover) TI , which is a diagnostic variable. Density of lake ice (ρI ) and snow (ρS ) are assumed tobe constant. Lake ice is assumed to have nonzero salinity, and its value SI is assumed to be aconstant parameter.
11.3.1

Calculation of heat ﬂux and growth rate

In ENTRY FIHEATL (in SUBROUTINE FIHSTL of lakeic.F), heat ﬂux in lake ice and its
growth rate is calculated.
• Output variables of FIHEATL

Variable

Variable in source code

Longname

Unit

WAO
WAS
WIO

WAO
WAS
WIO

growth rate of the lake ice in ice-free area
Snow growth rate due to heat imbalance
Basal growth rate of lake ice[cm w.e./s][cm w.e./s][cm w.e./s]• Input variables of FIHEATL

70

71

Variable in source code

A
QAI
QIO
SWABS
T(NLTDIM)
TS

Variable

AI
QAI
QIO
SW A
T, S
∆t

Lake ice concentration
Air-ice heat ﬂux multiplied by the factor of lake ice concentration
Vertical heat ﬂux through lake ice and snow
Shortwave radiation absorbed at ice-free lake surface, with the factor of ice-free area multiplied
Lake temperature / Salinity
Time stepLongname[−][−][W/m2 ]
[W/m2 ][◦ C/m2 ]/ [psu]
[s]

Unit

• Internal work variables of FIHEATL

Variable

Variable in source code

Longname

Unit∆T(k=2)
WF Z

TDEV
WFRZ

Freezing point depression
growth rate of the lake ice[◦ C][cm e.w./s]

Temperature at lake ice base is taken to be the lake model’s top level temperature T(k=2) .
In this model, lake ice exists only when and where T(k=2) is at the freezing point Tf , which is adecreasing function of salinity (Tf = −0.0543S is used here, where temperature and salinity aremeasured by ◦ C and psu, respectively). In the heat budget calculation for snow and lake ice, only
latent heat of melting and sublimation is taken into account, and heat content associated withtemperature is neglected. Therefore, temperature inside lake ice and snow are not calculated,
and TI is estimated from surface heat balance.
Nonzero minimum values are prescribed for AI and hI , which are denoted by Aminand hmin,II−61
which are set to 1.0 × 10 and 1.0 × 10 [cm] as default values, respectively. These parametersdeﬁne a minimum volume of lake ice in a grid (VImin ). If a predicted volume AI hI is less thanthe minimum volume (VImin ), AI is reset to zero, and T(k=2) is lowered to compensate thecorresponding latent heat. In this case, the lake model’s top level is kept at a supercooled state.
Such a state continues until the lake is further cooled and the temperature becomes low enoughto produce more lake ice than the minimum volume.
Surface heat ﬂux is separately calculated for each of air-lake and air-ice interfaces in one
grid. The surface temperature of lake ice TI is determined such that

QAI = QIO

(11.36)is satisﬁed, where QIO corresponds to GIcLk and QAI corresponds to GIcLk − WIcLk , respectively. (They are calculated in matdrv.F.) However, when the estimated TI exceeds the melting
point of lake ice Tm (which is set to 0 [◦ C] for convenience), TI is reset to Tm and QAI and QIO
are re-estimated by using it. The heat imbalance between QAI and QIO is consumed to meltsnow (lake ice when there is no snow cover). This heat imbalance changes the growth rate oflake snow (WAS ).

WAS =QAI − QIO
ρO Lf

(11.37)

where ρO is density of water and Lf is the latent heat of melting (the same value is applied
to snow and lake ice), respectively. Otherwise, it is zero when TI < Tm and negative when
TI = Tm .
Although it is assumed that T(k=2) = Tf when lake ice exists, Tk=2 could deviated from Tf
due to a change of salinity or other factors. Such deviation should be adjusted by forming or
melting lake ice. Under a temperature deviation of the top layer of lake,∆T(k=2) = T(k=2) − Tf S(k=2)(11.38)growth rate of the lake ice necessary to compensate it in the single time step (WF Z ) is given
by72

WF Z = −Cpo ∆T(k=2) ∆z1
Lf ∆t(11.39)where Cpo is the heat capacity of lake water and ∆z1 = 100cm is the thickness of the lake
model’s top level (uniformly set to constant in the current lake model.) WF Z is estimated at all
grids, irrespective of lake ice existence, for a technical reason. As described below, this growthrate ﬁrst estimates negative ice volume for ice-free grids, but the same heat ﬂux calculation
procedure as for ice-covered grids ﬁnally results in the correct heat ﬂux to force the lake. Basal
growth rate of lake ice ($ {IO}) is given by

WIO = AI WF Z +

QIO
ρO Lf

(11.40)Lake ice formation also occur in the ice-free area. Air-lake heat ﬂux except for shortwave in
ice-free area is given by
QAO = (1 − AI )[Q − (1 − αLk,SW )SW ↓ ](11.41)

where Q is air-ice heat ﬂux. Shortwave radiation absorbed at ice-free lake surface in ice-freearea (SW A ) is represented by
SW A = (1 − AI )(1 − αLk,SW )SW ↓(11.42)Growth rate of the lake ice in ice-free area (WAO ) is then calculated by

WAO = (1 − AI )WF Z +QAO + I(k=2) SW A
ρO Lf

(11.43)where I(k=2) denotes the fraction of SW A absorbed by the lake model’s top level, which is
calculate in SUBROUTINE [SVTSETL] of lakepo.F, respectively.
Finally, the heat ﬂux for freshwater is

Glake = ∆z1
11.3.2∆T∆t(11.44)

Sublimation and freshwater ﬂux for lake

In ENTRY FWATERL (in SUBROUTINE [FWASTL] of lakeic.F), sublimation (freshwater)
ﬂux, which is practically come from the land ice runoﬀ, is calculated or prescribed over lake ice
cover.
• Output variables (Prognostic variables) of FWATERL73

Variable

VariableLongname

Unit

A′I
h′I
h′S

AX
HIX
HSX

Lake ice fraction
Lake ice thickness
Snow depth

[−]
[cm]
[cm]

• Input variables of FWATERL

Variable

VariableLongname

Unit

F wEv
F wSb
Sof f

WEV
WSB
SOFF

Latent heat ﬂux of evaporation
Latent heat ﬂux of sublimation
Overﬂow snow ﬂux

cm/s
cm/s
cm/s

• Parameters of the lake ice model

Variable
ρS
ρI
RρS
RρI
hmin
I

Variable
rhos
rhoi
rrs
rri
himin

Longname

Unit

Value

Density of snow
Density of lake ice
Ratio of density (ocean/snow)
Ratio of density (ocean/ice)
Minimum thickness of ice[g/cm3 ]0.33
0.9
ρO /ρsρO /ρI
1.0 × 101[g/cm3 ]
[−]
[−]
[cm]

The sublimation ﬂux (F wSb ) is ﬁrst consumed to reduce the depth of lake snow in n-th time
step (hnS ).
ρO F wSb ∆t
ρS AnIh′S = hnS −(11.45)where h′S is the updated depth of lake snow. If h′S becomes less than zero, it is reset to zero.
The melted snow is added to F wSb .

′F wSb= F wSb +ρS AnI (h′S − hnS )ρO ∆t(11.46)′ is updated sublimation ﬂux.
where F wSb′ is not zero, the remain ﬂux is consumed to reduce
When there no remains snow but F wSbnthe thickness of lake ice (hI ).

h′I = hnI −′ ∆t
ρO F wSb
ρI AnI(11.47)

where h′S is the updated thickness of lake ice. If h′I becomes less than hmin
, it is reset toI′ .
zero. The melted ice is added to F wSb

74

′′
′
F wSb
= F wSb−

ρI AnI (hnI − h′I )ρO ∆t(11.48)′′ is updated sublimation ﬂux.
where F wSb′′ is consumed to reduce lake ice concentration (An ).
Nonzero F wSb
I

A′I = AnI −′′ ∆tρO F wSb
ρI hminI(11.49)where A′I is the updated lake ice concentration. If A′I becomes less then 0, it is set to zero.
Even if A′I becomes less than Amin, on the other hand, it is not adjusted here. because such anIadjustment means that the sublimation ﬂux is not used up by eliminating snow and lake ice.
Finally, the evaporation ﬂux F wEv is modiﬁed so as to reduce the lake water.

′
′′
+
F wEv
= F wEv + F wSb(A′I − AnI )hmin
IRρI ∆t(11.50)′′is just as
is updated evaporation ﬂux. When the lake ice does not exist, F wEv
where F wEv′
F wEv= F wEv + F wSb(11.51)

The adjusted evaporation ﬂux is saved as ∆F wEv .
′∆F wEv = F wEv− F wEv(11.52)

When sublimation ﬂux is consumed to the reduce lake ice amount, salt contained in lake ice
has to be added to the remaining lake ice or the underlying water in other to conserve total saltof the ice-lake system. Here, it is added to underlying water, and the way of this adjustment is
described in setcion 11.3. Such an adjustment is not very unreasonable because the ice tends togradually drain high salinity water contained in brine pockets in reality. When A′I is adjusted
to zero, on the other hand, the remaining sublimation ﬂux is consumed to reduce lake water. In
this case, diﬀerence between the latent heat of sublimation and evaporation has to be adjusted,which is also described in section 11.3.
If the ice and/or snow is too thick, they are converted to snow ﬂux. Here, the overﬂow snow
ﬂux Sof f is added to F wS
F wS′ = F wS + Sof f(11.53)where F wS′ is the updated snow ﬂux, and Sof f is the overﬂow snow ﬂux, which is actuallycalculated in SUBROUTINE[MATDRV] (of matdrv.F), then handed to ENTRY FWATER.
11.3.3Updating lake ice fraction• Input variables of PCMCTL
75

Variable

Variable in source code

Longname

Unit

h′I
WAO

HIX
WAO

Thickness of lake ice
Snow growth rate in ice-free area[cm]
[cm w.e./s]

In ENTRY PCMPCTL (in SUBROUTINE CMPSTL of lakeic.F), the lake ice fraction is
updated, using the lake ice thickness (h′I ) and the growth (retreat) rate in ice-free area (WAO ):A′′I = A′I +ρO′ρI hI ϕWAO ∆t(11.54)

where A′′I is the updated lake ice concentration, and ϕ is a factor coeﬃcient, which is set to
4.0 if WAO is greater than 0, otherwise is set to 0.5 as a default, respectively. If A′′I becomes
greater than 1, it is reset to 1, or if A′′I becomes smaller than zero, it is reset to zero.
11.3.4Growth and Melting of lake ice and snow

In ENTRY PTHICKL (in SUBROUTINE OTHKSTL of lakeic.F), the lake ice growth andmelting are calculated.
• Variables of PTHICKL

Variable

Variable in source code

Longname

AI
VI
VS
VIn+1
VSn+1
hI
hS
WAS
WAI
Wres
WIO
F wS
WAO
F wP r
Le
F wSb
FW
∆t
Rof f
–
FS

AX/AZ
AXHIX
AXHSX
AXHIXNAXHSXNHIX HIZ
HSX/HSZ
WAS
WAI
WRES
WIO
SNOW
WAO
PREC
EVAP
SUBI
FT
TS
ROFF
ADJLATFS

Lake ice fraction
Lake ice volume
Lake snow volume
Lake ice volume
Lake snow volume
Lake ice thickness
Snow depth
snow growth rate due to heat imbalance
growth rate of the lake ice due to heat imbalance
Residual heat ﬂux
Basal growth rate of lake ice
Snowfall ﬂux
growth rate of the lake ice in ice-free area
Precipitation ﬂux
Latent heat ﬂux of evaporation
Latent heat ﬂux of sublimation
Freshwater ﬂux
Time step––
Heatﬂux for lake surface

The initial conditions of lake ice volume (VI′ ) and snow volume (VS′ ) in the n-th time step
are presented by76

VI′ = A′′I hnI(11.55)VS′ = A′′I hnS(11.56)

Here after, predicted values of A′′I , VI′ , and VS′ are denoted by A∗I , VI∗ and VS∗ , respectively.
11.3.4.1

Freshwater ﬂux to lake ice and snow

Contribution of snowfall and freshwater ﬂux to the growth of lake ice and snow volume areconsidered.
Changes of snow depth due to snow fall (freshwater) ﬂux (F wS ) (expressed by negativevalues to be consistent with other freshwater ﬂux components) is ﬁrst taken into account.
When the lake ice does not exist, the amount of snow existed before the growth is added tothe snowfall ﬂux (F wS′ ).
F wS′′ = F wS′ +ρS VS′
ρO ∆t(11.57)

where F wS′′ is the updated snowfall ﬂux. In this case, the snow depth and its amount become0.
h∗S = 0(11.58)VS∗∗ = 0(11.59)where h∗S and VS∗∗ are the updated depth of lake snow and lake snow amount, respectively.
When the lake ice exists, the snowfall accumulates over the ice covered region. The depth oflake snow (hnS ) is modiﬁed byh∗S =VS∗∗ ρO F wS′ ∆t
+
A∗I
ρS

(11.60)The snow amount (VS∗ ) is also modiﬁed byVS∗∗ = A∗I h∗S(11.61)The snowfall ﬂux (F wS′ ) is then reduced by VS∗∗ .
F wS′′ = (1 − A∗I )F wS′

77(11.62)

where F wS′′ is the updated snowfall ﬂux. Finally, the snowfall ﬂux (F wS′′ ) is put together
with the precipitation ﬂux.
F wP′ r = F wP r + F wS′′

(11.63)where F wP′ r is the updated precipitation ﬂux.
11.3.4.2

Growth of lake snow

When the lake ice does not exist (A∗I = 0), the snow amount (VS∗∗ ) is converted to ice. In
this case, ﬂux for lake snow growth (WAS ) is used for the basal growth of the lake ice.
∗
WIO
= WIO + WAS(11.64)∗ are the updated basal growth rate of lake ice.
where WIO
When the lake ice exists, a residual heat ﬂux (Wres ) is considered.

Wres =ρO VS∗∗
+ WAS
ρS ∆t

(11.65)

If the residual heat ﬂux is negative (Wres < 0), the snow amount is reduced by

∗
WAS
=−

ρO VS∗∗
ρS ∆t(11.66)∗ is the growth rate of the lake snow. In this case, W
where WAS
res is assumed to reduce thelake ice.

WAI = Wres

(11.67)

where WAI is a growth rate of the lake ice.
The depth of lake snow (h∗S ) is modiﬁed with the accumulation.
h∗∗
S =

∗ ∆t
VS∗∗ + ρO WASρS A∗I(11.68)∗∗where h∗∗S is the updated depth of lake snow. If hS is less than 0, it is set to zero. Whenthe residual heat ﬂux is zero or positive (Wres ≥ 0), the growth rate of the lake ice (WAI ) is
temporally set to 0.
The lake snow amount (VS∗∗ ) is modiﬁed byVS∗∗∗ = A∗I h∗∗
S
where VS∗∗∗ is the updated lake snow amount.

78(11.69)11.3.4.3

Growth of lake ice′ ) for the lake ice growing.
The residual heat ﬂux is updated (Wres

′
Wres
=

ρO VI∗+ WAI
ρI ∆t(11.70)′ <0), W ′When the residual heat ﬂux is negative (Wresres is handed to the lake surface anddeﬁcient ﬂux is used for the lake ice melting.

∗′
WIO= WIO + Wres

∗
WAI
=−

ρO VI∗
ρI ∆t(11.71)(11.72)∗ is the updated basal growth rate of lake ice, and W ∗ is the updated growth rate
where WIOAI
of lake ice, respectively.
The amount of the lake ice (VI∗ ) is modiﬁed byVI∗∗ = VI∗ +∗ +W
ρO (WIO
AO )∆tρI(11.73)where VI∗∗ is the updated lake ice amount. If the lake ice amount becomes equal to or lessthan zero (VI∗∗ ≤ 0),
A∗∗
I = 0,h∗I = hminI(11.74)∗
where A∗∗I and hI are the updated lake ice concentration and its thickness. If the lake ice
amount becomes equal to or less than zero (VI∗∗ > 0),∗A∗∗
I = AI ,

h∗I =

VI∗∗
A∗I

(11.75)

However, if h∗I becomes less than hmin
, they are redeﬁned by
I
A∗∗
I =

VI∗∗
,
hmin
I

h∗I = hminI(11.76)

It is noted that the snow thickness (h∗I ) is not modiﬁed even the A∗∗
I is modiﬁed, so snow on
the disappearing ice is regarded as falling onto the created ice-free lake surface.
min .
A∗∗
I is set to be zero when it is less than AI79

11.3.4.4

Lake ice ﬂows

The lake ice is supposed to ﬂow to cover the lake. A simple method that does not deal withdynamical deformation as in the sea ice model(Hasumi, 2015) is used.
First, a case of the lake ice does not ﬂow because the lake ice amount is small enough is
considered. Ice volume assuming maximum and minimum extent with the minimum thickness) are deﬁned.
of lake ice (hminI
VImax = hmin
Amax
,
I
I

VImin = hmin
Amin
I
I

(11.77)

where Amax
and Aminare maximum and minimum lake ice concentration, which are set toII−6
1 and 1.0 × 10 as default values. When VI∗∗ ≤ (VImax − VImin ), the lake ice does not ﬂow.
Otherwise, the thickness of lake ice (h∗I ) is modiﬁed so as that the lake ice fully covers the lake.

h∗∗
I =

∗
(VImax − VImin ) + A∗∗
I hI
Amax
I

(11.78)where h∗∗I is the updated depth of lake snow. The depth of lake snow and the lake ice
concentration are also updated.
A∗∗ h∗∗ISh∗∗∗= Amax −(V max
S−V min )/h∗∗
I

I

−
An+1
= AmaxI
I

I

VImax −VImin
h∗∗
I

I

(11.79)is the updated depth of lake snow and the lake ice concentration, and An+1
iswhere h∗∗∗SI
the lake ice concentration in the new time step, respectively.
Freshwater and sublimation ﬂuxes aﬀect on salinity (FW (N LT DIM =2) ) are
FW (N LT DIM =2) = F wEv − F wP′ r − Rof f + WS(2) + WI(2)′′ )FS = SI (WI(2) − F wSb(11.80)where WI(2) and WS(2) are the growth rate of the lake ice snow below.
∗
ρS An+1
h∗∗
I −VI
I
ρI ∆t
∗∗
ρ An+1 h∗∗∗S −VS
= S I ρS ∆t

WI(2) =
WS(2)11.3.4.5(11.81)

Sinking snow

If there is a large amount of snow and the lake ice surface is below the lake surface, thesinking snow is converted to the lake ice. The depth of lake snow is updated to take buoyancy
into account.
hn+1= min(h∗∗∗S ,S

ρO − ρI ∗∗
hI )
ρS

80

(11.82)

The ice thickness is also updated.
hn+1
= h∗∗
I +
I

ρS ∗∗∗(h − hn+1
S )ρI S(11.83)The heat ﬂux aﬀects on the lake water temperature (FW (N LT DIM =1) ) is
Lf 
WI(1) + WS(1) − Sn + ∆F wEv
Cp

FW (N LT DIM =1) = −FW (N LT DIM =1) +(11.84)

where WI(1) and WS(1) are the growth rate of the lake ice snow below.
∗∗
ρS An+1
h∗∗
I −VI
I
ρI ∆t
∗∗
ρ An+1 h∗∗∗S −VS
= S I ρS ∆t

WI(1) =
WS(1)

11.4

(11.85)

Physical formulation and processThis sub-section [lakepo.F] introduces the vertical heat and salinity diﬀusion and convection of
the lake water.
11.4.1

Setting the vertical diﬀusion coeﬃcients

Input
Meaning

CharacterIn codeUnit

The default value of vertical diﬀusion coeﬃcient

KV 0

AHVL0

-

OutputMeaning

CharacterIn codeUnit

Vertical diﬀusion coeﬃcient

KV

AHVL

-

In this part, the vertical diﬀusion coeﬃcients are ﬁrst determined for each lake layer, andthe source code is included in ENTRY [VDIFFL] (in SUBROUTINE: [LAKEPO] of lakepo.F).
In the current version of MATSIRO 6, the diﬀusive ﬂux of tracer follows Fick’s law. In Fick’s
ﬁrst law, the vertical diﬀusive tracer ﬂux J of unit time is:
J = KV∂T∂z(11.86)

where T represents temperature (temperature is used as an example here, this equation applies
to salinity as well), KV represents the diﬀusion coeﬃcient, and z represents the vertical distance.
The vertical diﬀusion coeﬃcients for all lake layers are set as the default value, which is KV 0 =1.
Therefore, the vertical diﬀusion coeﬃcient KV 0 of each layer is simply set as:
KV (k) = KV 0

81

(11.87)

It is known that the surface wind and vertical temperature stratiﬁcation will also inﬂuence thevertical diﬀusion coeﬃcients, and they will be considered and added in the future version.
11.4.2Estimate the diﬀusion terms of the tracer equations

Input
Meaning

CharacterIn codeUnit

Vertical diﬀusion coeﬃcient
Water temperature
Water depth

KV
T
h

AHVL
TX(1)
HX

oC

m

OutputMeaning

CharacterIn codeUnit

The vertical component of diﬀusive tracer ﬂux

FD

ADT

o C/s

Internal variables
Meaning

CharacterIn codeUnit

The water depth excluding the ﬁrst layer
Depth of each lake layer

h′
D

HZBOT
DZ

m
cm

This part introduces the calculation of the vertical diﬀusive tracer ﬂux, and the source code
is included in ENTRY [FLXTRCL] (in SUBROUTINE: [LAKEPO] of lakepo.F). Temperature
is shown here as an example, and the tracer equation of salinity is identical.
The thickness of each lake layer is ﬁrst calculated. The water column is divided into 5 layers(k=2, 3, 4, 5, 6) in the current version, and the k=1 and k=7 represent the top surface and thebottom surface of the water column, respectively. The layer of the water column is shown in
Fig. 11-1.
Firstly, the thickness of the ﬁrst layer (k=2) is ﬁxed and represented as D1 , and D1 is set to
100 cm. Therefore, the thickness of the ﬁrst layer can be described as:D(2) = D1(11.88)The thicknesses of the other layers (k=3, 4, 5, 6) are calculated based on the ratio of theremaining thickness of the water column. Using h to represent the total depth of the watercolumn, and the thickness of the remaining water column (h′ ) excluding the top layer is presentedby:h′ = h − D1

(11.89)and thickness (D(k)) of the remaining layer (k=3, 4, 5, 6) can be represented as:
D(k) = S(k)h′(11.90)

The ratio (S(k)) of layer 3 to 6 are 0.1, 0.1. 0.2, and 0.6, respectively. Therefore, the thermalchange of the k-th layer by vertical diﬀusion can be represented as the ﬂux diﬀerence between
82

Solar radiation

Freshwater, heat and salinity flux

Lake surface

K=1 (Surface)
K=2, 1m

Heat and salinity diffusion

Convection when density reversed

K=3, 10%h’

K=4, 10%h

K=5, 20%
K=6, 60%
K=7 (Bottom)

Lake bottom

Figure 11.1: The layer of the water column
upper layer Jk−1,k and lower layer Jk,k+1 :FD (k) = Jk−1,k − Jk,k+1 = KV (k)T (k − 1) − T (k)D(k−1)+D(k)2− KV (k + 1)

T (k) − T (k + 1)
D(k)+D(k+1)
2

The above equations apply to salinity as well.
11.4.3Time integration of the tracer equations

Input
Meaning

CharacterIn codeUnit

The vertical component of diﬀusive tracer ﬂux
Minimum depth of the lake
Heat ﬂux
Freshwater ﬂux
Absorbed shortwave solar radiation
Salt ﬂux
Timestep
Surface-type fraction (lake)

FD
hmin
FT
FW
Sr
FS
∆t
Rlake

ADT
HXMIN
FT (1)
FT(2)
SWABS
FS
TS
LKFRACo C/s

OutputMeaning

CharacterIn codeUnit

Lake water deﬁciency

VD

XHD

cm

83

m
o C/s
cm/s

s
-

(11.91)This part introduces the update of water temperature and salinity of each layer due to thediﬀusion, freshwater ﬂux, and absorption of solar radiation. The source code is included in
ENTRY [SLVTRCL] (in SUBROUTINE: [LAKEPO] of lakepo.F).
Heat diﬀusion is ﬁrst considered. According to Fick’s second law, the temperature (salinity
as well) change of k-th layer follows:∂T∂2T= KV (k) 2∂t∂z(11.92)where z is the vertical distance. The equation is implicitly (backward-in-time) integrated, andobtains:
T n+1 (k) − T n (k)=∆tKV (k) T

n+1 (k−1)−T n+1 (k)D(k−1)+D(k)2− KV (k + 1) T

n+1 (k)−T n+1 (k+1)
D(k)+D(k+1)2

D(k)(11.93)

The Tridiagonal Matrix Algorithm (also known as the Thomas Algorithm) is used to solve theequation, the following equation is constructed:
AA (k)FD (k − 1) + AB (k)FD (k) + AC (k)FD (k + 1) = FD (k)(11.94)

KV (k)AA (k) = − D(k−1)+D(k) ∆t(11.95)

where

2

KV (k + 1)
AC (k) = − D(k)+D(k+1) ∆t(11.96)2

AB (k) = D(k) − AA (k) − AC (k)(11.97)

Then, the new coeﬃcient of the Tridiagonal Matrix Algorithm can be constructed:(
A′c (k) =(A′D (k) =

AC
AB ,
AC
,
AB −AA (k)A′C (k−1)

FD (k)AB (k) ,FD (k)−AA (k)A′D (k−1),
AB −AA (k)A′C (k−1)

k=2
k = 3, 4, 5, 6

k=2

k = 3, 4, 5, 6

(11.98)(11.99)Then, the vertical diﬀusion term can be obtained:FD (k) =A′D (k),′AD (k) − A′C (k)FD (k + 1),

84

k=
6
k = 2, 3, 4, 5

(11.100)

The water temperature of each layer is updated as:
Tk = T (k) + FD (k)∆t(11.101)

After solving diﬀusive changes of tracers, the eﬀect of heat, salinity, and freshwater ﬂuxes atthe lake surface is taken into account. Since the height of the water column changes with the
freshwater ﬂux, tracers at each layer should be re-estimated. It worth noting that a minimum
depth threshold hmin is set, which means the lake depth can not be lower than this value. Whensurface freshwater ﬂux FW (positive upward, i.e., the lake level is lowered when FW > 0) is
imposed, and the change of the water column for each time step can be represented as:
hD = −FW ∆t(11.102)and the lake water deﬁciency (when lake depth h is smaller than the minimum threshold hmin )can be represented as:
VD = max(hmin − h − hD , 0) × Rlake(11.103)

where Rlake is the lake fraction of the grid. Therefore, the change of the lake depth is ﬁnalizedas:hD = max(hD , hmin − h)(11.104)h = h + hD

(11.105)

and the lake depth is updated:

Due to the added freshwater, the depth of each lake layer has been changed, so the watertemperature of each layer is updated as well. The update of depth and water temperature starts
from the bottom layer. When hD > 0, the bottom of the k-th layer in the z coordinate is raised
by:
k
∆zB
= (1 −kX

S(k))hD(11.106)l=2

where S(k) is the vertical depth proportion of the k-th layer, and its top is raised by:∆zTk = (1 −k−1XS(k))hD(11.107)l=2

where the following equation holds:
B∆zkT = ∆zk−1(11.108)

Therefore, the temperature of the k-th (k=3, 4, 5, 6) layer becomes:
T (k) =

T (k)D(k) − T (k)∆zkB + T (k − 1)∆zkT
D(k)

85

(11.109)

The temperature of the ﬁrst layer (k=2) becomes:
T (2) =

T (2)D(2) − T (2)∆z2B + T (2)∆z2T + FT ∆t
D(2)(11.110)where FT is the heat ﬂux of the freshwater added to the lake, and for salinity, the salinity ﬂux
is FS .
When hD < 0, on the other hand, the bottom of the k-th level is lowered by ∆zkB , and its
top is lowered by ∆zkT . In this case, the above tracer equations for the k-th layer (k=3, 4, 5, 6)
and the ﬁrst layer (k=2) become:T (k)D(k) − T (k)∆zkT + T (k + 1)∆zkB
D(k)(11.111)

T (2)D(2) − T (2)∆z2T + T (3)∆z2B + FT ∆t
D(2)(11.112)T (k) =
andT (2) =

Finally, the absorption of solar radiation is considered. The temperature of each layer is updatedas:T (k) = T (k) +Sr × Csr (k)∆tD(k)ρ0 CP 0

(11.113)

where Sr represents the solar radiation, and Csr represents the absorption proportion of eachlayer. ρ0 and Cp0 are density and speciﬁc heat capacity of water, and their values are set as 1
g/cm3 and 3.99 × 107 erg/g/K, respectively. Csr is calculated based on the depth of each layer,
the transitivity Trs (k) at depth Dt (k) can be represented as:−

Trs = Rr × e

Dt (k)z1−+ (1 − Rr ) × e

Dt (k)
z2(11.114)

where Rr , z1 , and z2 are three parameters, and their values are 0.58, 35, and 2300, respectively.
Dt (k) represent the total depth of the bottom of k-th layer:
Dt (k) =

k
X

D(k)

(11.115)

l=2

Therefore, Csr can be represented as:
1 − Trs (k),k=2Csr (k) = Trs (k − 1) − Trs (k), k = 3, 4, 5Pk=6
1 − 5l=2 Csr (k),11.4.4(11.116)

The vertical convection

This part introduced the vertical convection of the water between diﬀerent layers, and the sourcecode is included in ENTRY [OVTURNL] (in SUBROUTINE: [LAKEPO] of lakepo.F).
A classical, still widely used method-convective adjustment, which unstable water column
is artiﬁcially homogenized with conserving heat and salt. MATSIRO 6 employs the convective

86

adjustment for the standard choice, but its algorithm is not the pairwise adjustment. It issummarized as follows:(1) Set the index Ks =2 (the ﬁrst layer), indicating the start level of the convective adjustment,
and set the index Ke =3 (the second layer), indicating the level where instability is currently
judged.
(2) Compare potential density (ρ) of the (Ke -1)-th and Ke -th layers. If instability (ρKe−1 >ρKe ) is found, go to step 3, and if instability is not found, go to step 4.
(3) Mix the water column between the Ks -th layer and Ke -th layer. Increase Ke by 1, and if
the (new) Ke is greater than 6 (the bottom layer), the convection procedure ends. If not,go back to step 2 with the new Ke .
(4) Set Ks =Ke , and increase Ke by 1. Increase Ke by 1, and if the (new) Ke is greater than 6(the bottom layer), the convection procedure ends. If not, go back to step 2 with new Ks
and Ke .
The water density is calculated by:3
ρ = 999.842594 + 6.793952 × 10−2 × T − 9.095290 × 10−3 × T 2 + 1.001685 × 10−4 × T(11.117)
5−1.120083 × 10−6 × T 4 + 6.536332 × 10−9 × T(11.118)

where T represents temperature. When instability happens, the water column between theKs layer and Ke will convect and mix, and the water temperature of these layers will be identicaland are updated as:
PKe
T (Ks , ..., Ke ) =

11.5

k=K

PKes

D(k)T (k)k=Ks D(k)(11.119)

Lake river coupling

The inﬂow rate from the river to the lake, Rin [cm/s], is calculated as follws.
Rin = Hriv /(τriv ∗ 86400)(11.120)

where Hriv [cm] is grid average river water depth, τriv is a time constant, and 86400 is a unit
conversion coeﬃcient. τriv is set to 1 day by default. Hriv is derived the following equation.
Hriv = wriv /ρw ∗ 100(11.121)

where wriv is the river water volume [kg/m2 ], ρw is the density of water, and 100 is thecoeﬃcient for unit conversion.
When water depth H[cm] of a lake becomes deeper than the upper limit Hmax , water outﬂows
into river channel. The outﬂow rate to the river, Rout [cm/s], is diagnosed as follows.
Rout = (H − Hmax )/(τlake ∗ 86400)

87(11.122)

Original MATSIRO
Water into soil
surface (e.g. rainfall,snowmelt)

Infiltration

Outflow

Surface storage

Surface runoff

Yes

Snowmelt?

Surface runoff into
surface storage

No
Surface runoff intoriver

Surface runoff intoriver

Figure 12.1: Flowchart of the wetland scheme in the MATSIRO 6
The Hmax [cm] is, by default, the climatology of lake depth + 10 m. τlake is a time constant,which is also set to 1 day.
Finally, Rout is added to runoﬀ from land and given to the river.

12

Wetland12.1

Outline of wetland scheme

A snow-fed wetland scheme, in which snowmelt can be stored with consideration of sub-grid terrain complexity, is incorporated as a sub-module of TOPMODEL (when # ifdef OPT SW STORAGE
in SUBROUTINE: [MATROF] in matrof.F is active) in MATSIRO 6 to represent the wetlandrelated process in the middle and high latitudes grid with snowmelt (Nitta et al., 2015, 2017,
Fig. 12-1). The wetland scheme has two major eﬀects: 1) the storage of part of the surface
water and delay of runoﬀ to rivers, 2) an increase in land surface wetness thus enhancing theevaporation in water-limited regimes.
With the wetland scheme, when snowmelt occurs, instead of all the generated surface runoﬀ
ﬂows directly to the rivers, only a part of the surface runoﬀ ﬂows into the rivers and the remaining
part of the surface runoﬀ is stored by the added tank (also known as the snow-fed wetland).
Then, the stored water in the wetland is then re-added to the water input of soil combining
with other kinds of water inputs (Fig. 12-1). In the current version, only snow-fed wetlands are
considered, and more types of wetland schemes will be added in the future version.

12.2

Inﬂow and outﬂow of the wetland

The inﬂow of the wetland comes from the fraction of the surface runoﬀ, and its amount isdetermined by the tunable parameter α. The outﬂow from the wetland is calculated using a88

time constant β and the wetland storage S, consequently ﬂowing into the soil surface. Therefore,
the update of the wetland storage S at each time step can be represented as:
S τ +1 − S τS= − + (1 − α)Rs∆t
β(12.1)

where Rs is surface runoﬀ calculated as a total of Ros (saturation excess runoﬀ), Roi (inﬁltration excess runoﬀ), and Roo (overﬂow of the uppermost soil layer); τ is time; and α and β
are parameters related to the inﬂow and outﬂow of the wetland storage, respectively.
β is a spatially dependent time constant, and can be calculated using a function of thestandard deviation of elevation above sea level:β = max(β0 (1 − min(σz (x), σzmax )/σzmax ), ∆t)(12.2)

where β0 is the maximum of the time constant, σz is the standard deviation of elevation above
sea level within each grid at point x, and ∆t is the time step of the model. The parameter σz
is a physical parameter calculated by the topography dataset, with a higher spatial resolution
than the simulation, and β0 , σzmax , and α are tunable parameters. These parameter values
were determined based on sensitivity simulations using an oﬄine land model with perturbed
parameters; 1 month, 200m, and 0.1 were chosen as the most appropriate values for β0 , σzmax ,and α, respectively (Nitta et al., 2015).

12.3

Storage of the surface runoﬀ

The ratio of total surface runoﬀ that ﬂows directly to the rivers is controlled by parameter α.
Therefore, the actual runoﬀ ﬂows into rivers Ro changes to:

Ro = (Ros + Roi + Roo ) × α + Rob(12.3)

where α is the inﬂow parameter (see 12.1); Ros is the saturation excess runoﬀ (Dunne runoﬀ),
Roi is the inﬁltration excess runoﬀ (Horton runoﬀ), and Roo is the overﬂow of the uppermost
soil layer, and all these three kinds of runoﬀ make up the total surface runoﬀ, and Rob is thegroundwater runoﬀ (section 7.3).

12.4

Water input of soil surface

The outﬂow from the wetland storage is re-added to the water input of the soil surface, combining
with the original water input (e.g. precipitation that passes through canopy gaps, water drops
from the canopy, and snowmelt water). Therefore, the updated soil water input W Isoil,total ofeach time step can be represented as:

W Isoil,total = P rc∗∗ + P rl∗∗ − Ros − Roi +S∆tβ(12.4)

where P rc∗∗ represents the convective rainfall, P rl∗∗ represents the non-convective rainfall, S
represents the wetland storage, β represents the outﬂow parameter of the wetland, and ∆t is
the time step.
89

13

Tile scheme

MATSIRO employs a tile treatment of the land surface to represent the subgrid land surface
types, so as to partially mimic the behavior at a higher resolution. The tile scheme is demonstrated in ENTRY:[LNDFLX] and [LNDSTP] (in SUBUROUTINE: [MATSIRO] of matdrv.F),
and the variables and parameters are introduced as follows:• Output variables

Variable

Description

Code

Units

F
Flake
Fi (i = 1, 2)

Fluxes at the land surface
Fluxes at the land surface of lake
Fluxes at the land surface of potential vegetation and cropland

MATFLX--• Parameters

90

91

Description
Fractional weight of lake in grid
Fractional weight of potential vegetation and cropland in gridSub-fractional weight of potential vegetation and cropland on land

PARAMETERflakefi (i = 1, 2)fi” (i = 1, 2)

LKFRACSFFRACSFFRAC1

Code
-

Units

Basically, one land surface grid is divided into three tiles in the control run — lake, potential
vegetation and cropland:
1. There are both snow-covered and snow-free fractions in each tile;
2. The surface heat and water ﬂuxes over lakes have been calculated as one of the tiles in agrid;
3. Both potential vegetation and cropland tiles consist of six soil layers, up to three snow
layers, and a single canopy layer, driving predictions of the temperature and amount ofwater in the canopy, soil, and snow;
4. Potential vegetation is deﬁned according to the vegetation types of the Simple Biosphere
Model 2 (Sellers et al., 1996a) scheme and has 10 categories including land ice. There is
no wetland category for land cover in the original SiB2 vegetation types or soil types.
All the prognostic and diagnostic variables are calculated in each tile, and the ﬂuxes at theland surface F are averaged:
F = Flake flake +

n
X

Fi fi

(13.1)

i=1

fi = fi” (1 − flake )nXfi” = 1

(13.2)(13.3)

i=1

where n is 2, Flake , F1 and F2 denote ﬂuxes at the land surface of lake, potential vegetationand cropland, flake , f1 and f2 denote their corresponding fractional weights (the sum of flake ,f1 and f2 always equals 1), f1” and f2” denote the sub-fractional weights of potential vegetationand cropland on land (note that f1” and f2” have been utilized as boundary condition data in
MATSIRO6), respectively.
By default, tile scheme is applied in land surface type, but it can be used for multiple
purposes.

92

References
Ball, J. T., 1988: An analysis of stomatal conductance. Ph.D. thesis, Stanford University.
Beven, K. J., and M. J. Kirkby, 1979: A physically based, variable contributing area modelof basin hydrology / un modèle à base physique de zone d’appel variable de l’hydrologie dubassin versant. Hydrol. Sci. Bull., 24 (1), 43–69.
Bougamont, M., J. L. Bamber, and W. Greuell, 2005: A surface mass balance model for thegreenland ice sheet: GREENLAND ICE SHEET MASS BALANCE MODEL. J. Geophys.
Res., 110 (F4).
Clapp, R. B., and G. M. Hornberger, 1978: Empirical equations for some soil hydraulic properties. Water Resour. Res., 14 (4), 601–604.
Collatz, G. J., J. T. Ball, C. Grivet, and J. A. Berry, 1991: Physiological and environmental
regulation of stomatal conductance, photosynthesis and transpiration: a model that includes
a laminar boundary layer. Agric. For. Meteorol., 54 (2), 107–136.
Collatz, G. J., J. A. Berry, G. D. Farquhar, and J. Pierce, 1990: The relationship between
the rubisco reaction mechanism and models of photosynthesis. Plant Cell Environ., 13 (3),219–225.
Collatz, G. J., M. Ribas-Carbo, and J. A. Berry, 1992: Coupled Photosynthesis-Stomatal conductance model for leaves of C4 plants. Funct. Plant Biol., 19 (5), 519–538.
Emori, S., 2001: Description of minimal advanced treatments of surface interaction and RunOﬀ(MATSIRO) land surface parameterization. Tech. rep.
Farquhar, G. D., S. von Caemmerer, and J. A. Berry, 1980: A biochemical model of photosynthetic CO2 assimilation in leaves of C 3 species. Planta, 149 (1), 78–90.
Hasumi, H., 2015: CCSR ocean component model (COCO). Atmosphere and Ocean ResearchInstitute, The University of Tokyo.
Hirabayashi, Y., 2004: Global analysis on long term variations of extreme river discharge. Ph.D.
thesis, The University of Tokyo.
K-1 Model Developers, 2004: K-1 coupled GCM (MIROC) description, K-1 tech. rep., 1. Tech.
rep., Center for Climate System Research, the Univ. of Tokyo, Tokyo.
Koirala, S., P. J.-F. Yeh, Y. Hirabayashi, S. Kanae, and T. Oki, 2014: Global-scale land surfacehydrologic modeling with the representation of water table dynamics: GLOBAL LSM WITH
WATER TABLE DYNAMICS. J. Geophys. Res., 119 (1), 75–89.
Kondo, J., and T. Watanabe, 1992: Studies on the bulk transfer coeﬃcients over a vegetatedsurface with a multilayer energy budget model. J. Atmos. Sci., 49 (23), 2183–2199.
Krinner, G., and Coauthors, 2005: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system: DVGM FOR COUPLED CLIMATE STUDIES. Global
Biogeochem. Cycles, 19 (1).
Liston, G. E., 2004: Representing subgrid snow cover heterogeneities in regional and global
models. J. Clim., 17 (6), 1381–1397.

93

Ma, W., and Coauthors, 2021: Applicability of a nationwide ﬂood forecasting system for typhoon
hagibis, 2019. Sci. Rep.
Miller, M. J., A. C. M. Beljaars, and T. N. Palmer, 1992: The sensitivity of the ECMWF modelto the parameterization of evaporation from the tropical oceans. J. Clim., 5 (5), 418–434.
Nitta, T., K. Yoshimura, and A. Abe-Ouchi, 2015: A sensitivity study of a simple wetland
scheme for improvements in the representation of surface hydrology and decrease of surface
air temperature bias. Journal of Japan Society of Civil Engineers, Ser.B1 (Hydraulic Engineering), 71 (4), 955–960.
Nitta, T., K. Yoshimura, and A. Abe-Ouchi, 2017: Impact of arctic wetlands on the climate
system: Model sensitivity simulations with the MIROC5 AGCM and a Snow-Fed wetland
scheme. J. Hydrometeorol., 18 (11), 2923–2936.
Nitta, T., and Coauthors, 2014: Representing variability in subgrid snow cover and snow depth
in a global land model: Oﬄine validation. J. Clim., 27 (9), 3318–3330.
Pokhrel, Y. N., S. Koirala, P. J.-F. Yeh, N. Hanasaki, L. Longuevergne, S. Kanae, and T. Oki,
2015: Incorporation of groundwater pumping in a global land surface model with the representation of human impacts. Water Resour. Res., 51 (1), 78–96.
Rutter, A. J., A. J. Morton, and P. C. Robins, 1975: A predictive model of rainfall interception
in forests. II. generalization of the model and comparison with observations in some coniferous
and hardwood stands. J. Appl. Ecol., 12 (1), 367–380.
Sellers, P. J., and Coauthors, 1996a: The ISLSCP initiative I global datasets: Surface boundary
conditions and atmospheric forcings for land-atmosphere studies. Bull. Am. Meteorol. Soc.,
77 (9), 1987–2005.
Sellers, P. J., and Coauthors, 1996b: A revised land surface parameterization (SiB2) for atmospheric GCMS. part i: Model formulation. J. Clim., 9 (4), 676–705.
Sivapalan, M., K. Beven, and E. F. Wood, 1987: On hydrologic similarity: 2. a scaled model ofstorm runoﬀ production. Water Resour. Res., 23 (12), 2266–2278.
Stieglitz, M., D. Rind, J. Famiglietti, and C. Rosenzweig, 1997: An eﬃcient approach to modeling the topographic control of surface hydrology for regional and global climate modeling.
J. Clim., 10 (1), 118–137.
Takata, K., S. Emori, and T. Watanabe, 2003: Development of the minimal advanced treatmentsof surface interaction and runoﬀ. Glob. Planet. Change, 38 (1), 209–222.
Tatebe, H., and Coauthors, 2019: Description and basic evaluation of simulated mean state,internal variability, and climate sensitivity in MIROC6. Geosci. Model Dev., 12 (7), 2727–
2765.
Watanabe, M., and Coauthors, 2010: Improved climate simulation by MIROC5: Mean states,
variability, and climate sensitivity. J. Clim., 23 (23), 6312–6335.
Watanabe, T., 1994: Bulk parameterization for a vegetated surface and its application to asimulation of nocturnal drainage ﬂow. Bound.-Layer Meteorol., 70 (1), 13–35.
Watanabe, T., and Y. Ohtani, 1995: A simple model of shortwave-radiation transport within
canopy. J. Agric. Meteorol, 51 (1), 57–60.
94

Yamazaki, D., S. Kanae, H. Kim, and T. Oki, 2011: A physically based description of ﬂoodplain inundation dynamics in a global river routing model: FLOODPLAIN INUNDATIONDYNAMICS. Water Resour. Res., 47 (4).
Yang, Z.-L., R. E. Dickinson, A. Robock, and K. Ya Vinnikov, 1997: Validation of the snow
submodel of the Biosphere–Atmosphere transfer scheme with russian snow cover and meteorological observational data. J. Clim., 10 (2), 353–373.
Yoshimura, K., S. Miyazaki, S. Kanae, and T. Oki, 2006: Iso-MATSIRO, a land surface modelthat incorporates stable water isotopes. Glob. Planet. Change, 51 (1), 90–107.

95