Integrated land-surface hydrological and biogeochemical processes in simulating water, energy and carbon fluxes over two different ecosystems

Journal Pre-proofs Research papers Integrated land-surface hydrological and biogeochemical processes in simulating water, energy and carbon fluxes over two different ecosystems Sidong Zeng, Jun Xia, Xiangdong Chen, Lei Zou, Hong Du, Dunxian She PII: DOI: Reference:

S0022-1694(19)31125-4 https://doi.org/10.1016/j.jhydrol.2019.124390 HYDROL 124390

To appear in:

Journal of Hydrology

Received Date: Revised Date: Accepted Date:

3 July 2019 2 October 2019 19 November 2019

Please cite this article as: Zeng, S., Xia, J., Chen, X., Zou, L., Du, H., She, D., Integrated land-surface hydrological and biogeochemical processes in simulating water, energy and carbon fluxes over two different ecosystems, Journal of Hydrology (2019), doi: https://doi.org/10.1016/j.jhydrol.2019.124390

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

1

Integrated land-surface hydrological and biogeochemical processes in simulating water,

2

energy and carbon fluxes over two different ecosystems

3

Sidong Zenga, c, Jun Xiaa,b,d,*, Xiangdong Chene, Lei Zoud, Hong Duf, Dunxian Sheb

4

aChongqing

5

Chongqing, 400714, China

6

bState

7

University, Wuhan 430072, China;

8

cChangjiang

9

dKey

Institute of Green and Intelligent Technology, Chinese Academy of Sciences,

Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan

Institute of Survey, Planning, Design and Research, Wuhan 430010, China

Laboratory of Water Cycle and Related land Surface Processes, CAS, Beijing 100101

10

China;

11

e China

12

f College

13

Wuhan 430074, China

14

Abstract: Understanding the mechanism of the interactions between vegetation dynamics and

15

the water cycle is rather important for determining global and regional water and carbon

16

budgets. In this paper, a physically-based model integrating land-surface hydrological and

17

biogeochemical processes by coupling a hydrological model and a biogeochemical model is

18

developed to simulate the water, energy and carbon fluxes. The model is validated against

19

observed biometric, eddy-covariance flux, soil moisture and temperature data over two

20

different ecosystems. Results show that the model could simulate the vegetation physiological

Water Exchange, Beijing, 100053, China of Resources and Environmental Science, South-Central University for Nationalities,

* Corresponding author: Address: State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072. E-mail address: [email protected] 1

21

and hydrological processes satisfactorily including net radiation, latent heat, gross primary

22

production (GPP), net ecosystem exchange and soil moisture and temperature. Sensitivity

23

analysis illustrates that evapotranspiration, GPP, net ecosystem production are quite sensitive

24

to plant physiological controls such as the maximum electron transport rate, quantum

25

efficiency of electron transport and runoff parameter. Moreover, results show that a close

26

relationship between photosynthesis and transpiration. Water use efficiency shows a U-bend

27

curve at different time steps and it also indicates a higher value in forest ecosystem than that

28

of grass ecosystem. In the forest ecosystem, evapotranspiration is higher and the surface

29

runoff is lower. Meanwhile it is almost carbon sink during the whole year, while in the grass

30

ecosystem, it shows a carbon source during May to September. This study could provide an

31

effective model for the simulation of water-carbon cycles.

32

Keywords: Hydrological; Biogeochemical; DTVGM; CASACNP; Coupling

33

1. Introduction

34

The bidirectional effects of vegetation dynamics and hydrologic balance are well known

35

and have received much attention in both hydrological and ecological research (Arora, 2002;

36

Wagener et al., 2010; Shen et al., 2013). For example, regarding the interactions between

37

vegetation and soil moisture via transpiration, water extraction by plant roots affects soil

38

moisture transport and then runoff generation, and the influences are various from different

39

species. Moreover, the feedback of soil moisture through its stress effects on transpiration,

40

photosynthesis, and respiration affect plant growth and the carbon cycle. The leaf area index

41

(LAI) also affects solar radiation transfer and rainfall interception by the canopy, subsequently

42

affects the evapotranspiration, runoff, etc. In addition, nutrient limitation on carbon

2

43

assimilation affects vegetation growth and the water balance. Several studies have shown that

44

under climate change condition with increasing carbon dioxide (CO2) levels, plants tend to

45

reduce stomatal conductance, suppress transpiration and enhance carbon uptake, which

46

implies increases in plant biomass (Pritchard et al. 1999), water use efficiency and soil

47

moisture (Volk et al., 2000). Moreover, the increases in plant biomass may not be sustainable

48

in the long-term due to nutrient limitation (McMurtrie and Comins, 1996). These findings

49

suggest that it is necessary to consider the interactions between carbon/nutrient dynamics and

50

hydrological processes to better understand vegetation growth and hydrologic balance.

51

Physically based hydrological models derived deductively from fundamental physical

52

laws (Beven, 2002) have successfully been used to explain complex hydrological processes,

53

such as the rainfall-runoff relationship, exchanges between surface and subsurface water, and

54

flood routing. However, these physically based models often have not considered the effects

55

of biogeochemical changes in vegetation such as LAI dynamics, root growth and stomatal

56

conductance on hydrological balance. Biogeochemical models have used to characterize the

57

carbon, nitrogen, and phosphorus cycles and simulate carbon, nitrogen, phosphorus transfer

58

between different pools. However, the hydrological processes in most biogeochemical models

59

are simplified, and the effects of the water cycle on biogeochemical processes are not well

60

understood.

61

From a modeling perspective, the coupling of physically based hydrological models and

62

biogeochemical models that characterize vegetation dynamics and related carbon/nutrition

63

cycles is required to explore the interactions between vegetation and hydrology (Park et al.,

64

2011), and it is currently receiving much attention (Manzoni and Porporato, 2011;

3

65

Riveros-Iregui et al., 2012; Wu et al, 2016; Zhao et al., 2018, 2019). For example, Park et al.

66

(2011) developed a coupled hydrological and biogeochemical model and applied it to a larch

67

forest in eastern Siberia. The authors found that soil water is the determinant that influences

68

CO2 fluxes. Shen et al. (2013) coupled the Process-based Adaptive Watershed Simulator and

69

Community Land Model and studied the complex effects of key controls of nitrogen,

70

groundwater and soil water content on evapotranspiration, net primary productivity and other

71

important water/carbon flux variables in a humid continental climate watershed in the Great

72

Lakes region of North America. They found that nitrogen significantly controls transpiration

73

and then influences other hydrologic fluxes. However, hydrological models still have rarely

74

considered biogeochemical processes.

75

The Distributed Time Variant Gain Model (DTVGM) developed by Xia et al. (2005) is

76

an extension of a nonlinear approach derived from Volterra function (Xia et al., 1997; Xia,

77

2002) coupled hydrological processes based on geographic information system and remote

78

sensing platform. The DTVGM has been applied to many basins in China and has achieved

79

good performance (Wang et al., 2009; Li et al., 2010; Zhan et al., 2013). The model includes

80

several components of hydro-information analysis and modeling, such as data processing,

81

snow melt, the evapotranspiration and runoff generation on each grid and flow routing

82

between adjacent concentration belts. The model has smaller parameter sets for runoff

83

generation, and the parameters can be estimated in terms of system identification approach to

84

reduce uncertainty (Xia et al., 2005).

85

CASACNP (Carnegie-Ames-Stanford Approach Carbon Nitrogen Phosphorous) is a

86

process-based biogeochemical model that was developed by Wang et al. (2007) and Houlton

4

87

et al. (2008) from the Carnegie-Ames-Stanford Approach (CASA) model (Randerson et al.,

88

1997). This model has been well calibrated and used in several previous studies to study

89

global and regional carbon cycles (Randerson et al., 2002; Wang et al., 2010). The model can

90

simulate vegetation dynamics, including those of leaves, wood and root, and the carbon cycle

91

in litter and soil organic matter. The model also considers the coupling of the biogeochemical

92

cycles of carbon, nitrogen and phosphorus. Based on the coupling, the model can be used to

93

examine the nutrition limit on vegetation productivity and growth, which is not considered in

94

most other terrestrial biogeochemical models (Wang et al., 2007). However, CASACNP

95

rarely considers hydrological processes or energy transfer. For example, soil moisture and soil

96

temperature should be prepared as inputs to drive the model (Shi et al., 2011), that is the

97

bidirectional effects of water and carbon are not well simulated. Besides, the accuracy of the

98

simulation for photosynthesis affected by physiological characteristics and environmental

99

factors is significant to carbon assimilation. In the CASA model, the net primary production

100

is estimated simplified by the absorbed photosynthetically active radiation which ignores the

101

complex physiological and environmental effects, which may bring certain error in the

102

simulation of carbon assimilation. From the above descriptions, neither of the hydrological or

103

biogeochemical models could meet the requirement of studying interactions between

104

vegetation and hydrological balance. For these reasons, we integrated CASACNP with carbon

105

assimilation modules modified to DTVGM in our research to enhance the ability and

106

adaptability in simulating water, energy and carbon fluxes.

107

In this study, the coupled DTVGM-CASACNP model is introduced to study water,

108

energy and carbon fluxes over two different ecosystems. The objectives of this study are (1)

5

109

to describe the physical processes of the model, (2) to evaluate the performance of the

110

coupled model in simulating water, energy and carbon fluxes, soil moisture and temperature

111

and (3) to study the coupling characteristic of vegetation and hydrological processes.

112

2. Model description

113

The coupling of the hydrological model (i.e., DTVGM) and biogeochemical model (i.e.,

114

CASACNP) is described below, and the schematic of the coupled model is shown in Fig.1.

115

The coupled model consists of three main modules: the land surface hydrological module

116

included in the DTVGM, the carbon assimilation module based on Farquhar et al. (1980), and

117

biogeochemical cycle module based on CASACNP.

118

-------------------------------------------------------------------------

119

Place Fig.1 here

120

-------------------------------------------------------------------------

121

2.1 Land surface hydrological processes

122

Water budget The coupled model solves the water budget for both the canopy and soil,

123

including rainfall interception, evaporation of intercepted water from the canopy,

124

transpiration, soil evaporation, surface runoff, subsurface runoff, vertical movement of soil

125

water, groundwater recharge and base flow.

126

Interception by the canopy plays an important role in redistributing water components

127

including evapotranspiration, infiltration and surface runoff, and even the deposition of

128

solutes such as dissolved organic carbon and organic nitrogen (Lohse et al., 2009). The

129

intercepted water (Wc) can be obtained by determining difference in gross precipitation (P)

130

arriving at the top of the canopy and the throughfall (Pd), stemflow (Pa) and evaporation (Ec)

6

131

from the intercepted water following Rutter et al. (1972) as follows.

132

dWc  P  Pd  Pa  Ec dt

133

Evapotranspiration (ET) is the link between the water budget and energy balance. In this

134

study, it consists of three components: evaporation from canopy intercepted water (Ec),

135

transpiration (Tr) and soil evaporation (Es). Ec is considered to be the potential evaporate from

136

water surface with the limitation of water storage in the wet leaf area. Tr, which is the main

137

force for water uptake by plant roots, is assumed to be the accumulation of water extraction

138

from each soil layer where plant roots appear. It is influenced by the degree of stomatal

139

opening and the atmosphere surrounding the canopy, including the absorbed net irradiation

140

and aerodynamics factors. Soil water content and soil temperature also affect transpiration by

141

controlling stomatal conductance. Es is controlled by the net radiation absorbed by the soil

142

surface, aerodynamics factors and soil surface resistance determined by soil moisture and soil

143

texture. The major equations for estimating ET are follows(Mo et al., 2005):

144

145

146

Tr 

Rnc   C p D0 / rac    1  rc / rac 

 Ec 

 Es 

1  f w 

Rnc   C p D0 / rac

   

fw

  Rns  G    C p D0 / ras    1  rs / ras 

(1)

(2)

(3)

(4)

147

where Rnc and Rns are the net radiation absorbed by the canopy and soil surface (Wm-2),

148

respectively. λ is the latent heat of vaporization of water (Jkg-1),  is the slope of the curve of

149

saturated vapor pressure to temperature (PaK-1),  is the psychrometric constant (hPaK-1),  is

150

air density (kgm-3), Cp is the specific heat at a constant air pressure (Jkg-1K-1), D0 is the vapor

151

pressure deficit at the source height (hPa), fw is the fraction of wet leaf area. rc, rs, rac, and ras 7

152

are the canopy resistance, soil resistance, leaf boundary aerodynamic resistance and

153

aerodynamic resistance from the soil surface to the source height, respectively (sm-1).

154 155

Surface runoff (Rs) is calculated by a nonlinear method which is a special formulation of the second order Volterra nonlinear functional series (Xia et al., 2002), as follows:

156

Rs   g1  g 2 API  Pe

157

API   exp   / K e  Pe  t    / K e d

t

0

(5) (6)

158

where g1 and g2 are parameters, API is the Antecedent Precipitation Index which indicates

159

underlying surface wetness, Pe is the effective rainfall arriving at the ground surface (mm), is

160

a time integrating variable, Ke is a parameter indicating the rate of soil moisture recession, and

161

t is the time (s).

162

Many previous studies have shown that soil moisture affects both the water and carbon

163

cycles including photosynthesis, plant and soil respiration, evapotranspiration, etc. (Reynolds

164

et al, 2004; Davidson and Janssens, 2006). The vertical transport of soil water limits the water

165

available for water extraction by plant roots in different soil layers. In this paper, soil moisture

166

transport is simulated by using the Richards’ equation as follows:

167

     K  K R  t z  z  z

(7)

168

where  is the volumetric soil water content (m3m-3), K is the soil hydraulic conductivity

169

(cm/s),  is the matric potential (cm), R is a source term that includes the uptake water by

170

roots and subsurface runoff, and z is the soil depth (m). The relationship between K,  and 

171

follows van Genuchten (1980):

8

s  r    r m  1   n           s

172

(  0) (8)

(  0)

m 2   r 1 K ( )  K s Se0.5 1  1  Se1/ m   , m  1  , n  1, Se    n s  r

173

(9)

174

where s and r are the saturated and residual volumetric water contents, respectively, Ks is

175

the saturated soil hydraulic conductivity (cms-1), and , m and n are parameters. The

176

Richards’ equation can be solved numerically with the given boundary conditions and initial

177

soil moisture at sub-daily timestep. The implicit difference equation of Equation (7) is as

178

following:

179

zi

i j 1,m 1  i j t

K

j 1 i 1/2



j 1, m 1 i 1

 i j 1,m 1 

zi 1  zi

K

j 1 i 1/2



j 1.m 1 i

 i j11,m 1 

zi  zi 1

(10)

1, m 1, m  K i j 1/2  K i j 1/2  Ri zi

180 181 182

In this study, the implicit difference equation is solved according to Celia et al.（1990） with the following boundary conditions.

K

 K z

 Pe  Es  RS , z  0,t  0

 z

0

(11)

z  z L ,t  0

183

Energy balance Energy balance is considered for both the canopy and soil surface in this

184

model. In this study, a “root-finding” numerically iterative algorithm is used to calculate the

185

canopy and soil surface temperature at sub-daily timestep. Energy balance equations are

186

solved iteratively until the energy balance error is smaller than a given permissible error. The

187

energy balance equations are as follows, and considering the heat storage in the canopy and

188

the soil surface layer are in the left terms in equation (10) (Sellers et al., 1986):

9

Tv  Rnv  H v  LEv t T Cgs gs  Rns  H s  LEs  G t Cv

189

(12)

190

where Cv and Cgs are the bulk heat capacity per unit area of the canopy and upper soil layer,

191

respectively (Jm-2K-1); Tv and Tgs are the canopy and soil surface temperature, respectively

192

(K), Rn, H and LE are the net radiation, sensible and latent heat fluxes, respectively, and the

193

subscripts v and s indicate the canopy and soil surface respectively (Wm-2). G is the soil heat

194

flux (Wm-2) and obtained by the following equation.

195

T T   Tgs  Ts ,1 G  0.5  k1  k2 s ,1 s ,2  z1 z2  

196

where k is the soil thermal conductivity (Wm-1K-1), Ts is the soil temperature. The subscript 1

197

and 2 indicates the first and second layers under the ground surface, respectively.

(13)

198

Soil temperature is also an important factor that affects root growth, photosynthesis,

199

respiration and transpiration. In this coupled model, soil heat transfer is simulated by

200

numerically solving the following heat flow equation with the given boundary conditions and

201

initial soil temperature.

Cs

202

Ts   Ts   k  t z  z 

(14)

203

where Cs is the soil volumetric heat capacity (Jm-3K-1), and k is the soil thermal conductivity

204

(Wm-1K-1). Cs is calculated following Flerchinger and Saxton (1989) and k is calculated

205

following Jonhansen (1975) which considering the soil water content.

206

2.2 Carbon assimilation

207

Photosynthesis The photosynthetic rate can be obtained by both empirical and physical

208

models. In this coupled model, the photosynthetic rate is calculated for C3 and C4 species 10

209

based on Farquhar et al. (1980) and Collatz et al. (1991, 1992), respectively. The

210

photosynthetic rate is calculated for both sunlit and shaded leaves at the canopy scale

211

following De Pury and Farquhar (1997). For both kinds of leaves, the net photosynthetic (A)

212

rates are limited by Rubisco-limited photosynthesis (Av), electron transport-limited

213

photosynthesis (Ae), export-limited carboxylation (for C3 species) or PEP-carboxylase-limited

214

carboxylation (for C4 species) (As), and dark respiration (Rd) as follows:

215

A  min  Av , Ae , As   Rd

216

The units of all the terms in the above equation are mol m-2s-1. In the equation, Av, Ae and As

217

are calculated as follows:

218

219

220

ci    (C3 ) Vc max c  K 1  o / K Av   i c i o V (C4 )  c max  J ci    Ae   4 ci  2  Q  par

(C3 )

(15)

(16)

(17)

(C4 )

(C3 )  0.5Vc max As   3  4 10 Vc max ci (C４)

(18)

221

where Vcmax is the photosynthetic Rubisco capacity (mol m-2s-1), ci and oi are the intercellular

222

CO2 and oxygen (O2) concentrations, respectively (mol mol-1). Kc and Ko are the

223

Michaelis-Menten constants of Rubisco for CO2 and O2, respectively (mol mol-1).  is the

224

CO2 compensation point of photosynthesis (mol mol-1). J is the electron transport rate (mol

225

m-2s-1).  is intrinsic quantum efficiency (mol CO2 mol-1), and Qpar is the absorbed

226

photosynthetically active radiation (mol m-2s-1).

227

Stomatal conductance Stomatal behavior is a controlling factor for both transpiration

228

and photosynthesis. Stomatal conductance is influenced by both plant physiology and the 11

229

surrounding weather conditions. Following Leuning (1995) and Wang and Leuning (1998),

230

the relationship between the photosynthesis rate and stomatal conductance is as follows:

g s  g0 

231

a1 f sm A  cs   1  Ds / D0 

(19)

232

A  g c  ca  ci 

(20)

233

g c 1  1.56 g s1  1.37 gb1  g a1

(21)

234

where gs, g0, gb, and ga are the bulk stomatal, residual, boundary and aerodynamic

235

conductance for water vapor, respectively (mol m-2s-1), gc is the total conductance for CO2

236

from the intercellular level to the reference source height above the canopy (mol m-2s-1). a1

237

and D0 (Pa) are empirical parameters, cs and ca are the CO2 concentrations at the leaf surface

238

and source height above the canopy, respectively (mol mol-1). Ds is the vapor pressure

239

deficit (Pa), and fsm is the effect of soil water content on stomatal conduction, which is

240

simulated using an empirical function (Wang and Leuning, 1998). In this study, the stomata

241

model and photosynthesis model are coupled to simulate photosynthesis and stomatal

242

conductance by solving equations (14) to (19) iteratively until convergence. For the accuracy

243

of the simulation, the stomata model and photosynthesis model should be run at sub-daily

244

timestep.

245

2.3 Biogeochemical cycles

246

Carbon cycle In the CASACNP model, there are nine carbon pools consist of three plant

247

pools: leaf, root and wood, three litter pools: structural, metabolic and coarse woody and three

248

soil organic pools: microbial biomass, slow and passive pools. The carbon is transferred from

249

plant pools to litter pools and from litter pools to soil pools. Meanwhile carbon is also

250

turnover among different pools in litter and soil pools individually. The major equations of 12

251

carbon cycle are as followings (Wang et al., 2010):

dCi  ac ,i Fc  i Ci , i  leaf,wood or root dt

252

253

254

dC j dt

(22)

  b j ,i i Ci  mn  j C j j  metabolic, structural or coarse woody (23) i

dCk   ck , j mn  j C j   d k ,kk kk Ckk  k Ck , kk  k , dt j kk k  microbial, slow or passive,  ac ,i  1,  b j ,i  1 i

(24)

j

255

Where C is pool size (g C m-2) and  is turnover rate (d-1) with the subscript i, j and k

256

denoting for plant, litter and soil, respectively. Fc is the net primary productivity (NPP) and

257

ac,i is the fraction of NPP allocated to three plant pools. bj,i is the fraction of carbon from plant

258

pool i to litter pool j and ck,j is the fraction of carbon from litter pool j to soil pool k. dk,kk is the

259

fraction of carbon from soil pool kk to soil pool k. mn is the limitation of nitrogen on litter

260

decomposition. Leaf area index (LAI) is simulated by the function of carbon storage in leaf:

261 262

LAI  SLA0  Cleaf

(25)

where SLA is the leaf specific area (m2 gC-1).

263

Phenology Phenology dynamics are an important component when relating carbon

264

modeling and ecosystem function (Parmesan and Yohe, 2003). Phenology dynamics,

265

including the time of budburst and blossom of different species are measured and simulated

266

and have been widely used to study responses to climate change at different site. However,

267

individual phenology dynamics can only work well at the local scale and for a few kinds of

268

vegetation and do not satisfy the differences for different ecosystems (Zhang et al., 2006).

269

According to Moulin et al. (1997) and Zhang et al. (2006), using vegetation index

270

measurements can be an alternative approach to identify vegetation phenology. With this

13

271

approach, phenology is divided into four phases by greenness, which is quantified using a

272

vegetation index: (1) greenup: the date of onset of greenness increase; (2) maturity: the date

273

when greenness reaches the maximum degree (3) senescence: the date when greenness begins

274

to decrease; (4) dormancy, the date when greenness is lowest. In this study, the phenology

275

characteristics are derived from the remote sensing observations by Zhang et al. (2006). In

276

different phenology phases, the allocations of NPP to plant pools differ. In phase 1, the

277

fraction of NPP allocated to leaf is set to 0.8, root is 0.1, and wood is 0 and 0.2 respectively

278

for woody and non-woody vegetation. In phase 2, the fractions of allocation are constant but

279

vary between different vegetation types, and the values follow Fung et al. (2005). During

280

phases 3 and 4, the fraction of allocation to leaf is set to 0 and the allocation coefficient in

281

phase 2 is allocated to wood and root.

282

3. Sites and model initialization

283

3.1 Sites dataset

284

This research was performed in the Loblolly Pine (LP) site and Open Field (OF) site in

285

the Blackwood Division of the Duke Forest near Durham, North Carolina (35.98N,

286

79.09W). This area has a moderate climate with warm and humid summers and mild winters.

287

The mean annual precipitation is approximately 1117.6 mm and is well distributed throughout

288

the year. The mean annual temperature ranges from 21.7 C to 9.0 C in the growing and

289

nongrowing seasons, respectively. The soils are low-fertility acidic Hapludalf in the Enon

290

Series, with clay loam in the upper 0.3 m and clay below down to the bedrock at 0.7 m.

291

In the LP site, the surrounding dominant stand is a uniform-age loblolly pine forest that

292

was planted in 1983. The understory is diverse and includes more than 26 species dominated

14

293

by sweetgum, red maple, winged elm and flowering dogwood, which utilize the C3

294

photosynthesis pathway. The growing season occurs from March to mid-October and the leaf

295

area index ranges from 2.0 to 5.0 in all phases. In the OF site, the dominant stand is the C3

296

grass Festuca arundinacea shreb, and few plants use the C4 photosynthesis pathways. The leaf

297

area index ranges from 1.0 to 3.0 in all growth phases. It should be noted that the field is

298

mowed once or twice each summer since 1992.

299

Information of the two sites is shown in Table 1.

300

-------------------------------------------------------------------------

301

Place Table 1 here

302

-------------------------------------------------------------------------

303

The data used in this study were downloaded from the AmeriFlux Site and Data

304

Exploration System at the site http://ameriflux.ornl.gov/. The LP site is equipped with seven

305

20 m aluminum walkup towers. The flux and radiation instruments are set at 20 m from the

306

ground surface and other instruments such as temperature probe are set at 20 m and relative

307

humidity probe is set at 7 m. The OF site is equipped with 4.8 m walkup towers and the eddy

308

covariance flux observation system is set at 2.8m. The meteorological data including

309

precipitation, atmospheric pressure, air temperature, relative humidity, wind speed, and

310

incoming solar radiation at 30 min time step are used to drive the coupled model to simulate

311

the water, energy and carbon fluxes. Other meteorological data, such as the net radiation,

312

sensible heat, latent heat, and evapotranspiration and biological data, such as gross primary

313

productivity are used to validate the model.

314

The “goodness” of energy balance closure (i.e. net radiation equals the sum of the sensible

15

315

heat, latent heat and ground heat flux) is an indicator of observed data quality. Because the

316

ground heat flux was not available, therefore, the observed net radiation is compared to the

317

sum of the sensible heat and latent heat, as shown in Fig.2. The correlation coefficients are

318

0.76 and 0.91 for the LP site and OF site, respectively. The average relative error of energy

319

ARE 6.2Wm-2 and 15.2Wm-2 for the two sites.

320

-------------------------------------------------------------------------

321

Place Fig.2 here

322

-------------------------------------------------------------------------

323

3.2 Model initialization

324

In order to get more accurate water, energy and carbon fluxes at higher time resolution,

325

the simulation is conducted at a time step of 30 minutes using half-hourly measured

326

meteorological data including air pressure, temperature, vapor pressure, wind speed,

327

precipitation, and downward short and long wave radiation to drive the coupled model. The

328

CO2 concentration is dynamically set during the simulation periods using the measured

329

ambient CO2 concentration above the vegetation. The average ambient CO2 concentration is

330

about 375 ppmv.

331

The main parameterization of the model consists of soil physical parameters, runoff

332

parameters, eco-physiological and biogeochemical parameters. The key parameters for the

333

simulation are listed in Table 2. The soil physical parameters are based on the United States

334

Department of Agriculture (USDA) classification system. The soil moisture and temperature

335

are calculated for 50 layers from the ground surface to a depth of 2 m. Most of the

336

eco-physiological and biogeochemical parameters are according to the previous studies. In

16

337

this study, we mainly calibrate the runoff, photosynthetic and stomatal conductance

338

parameters manually.

339

To evaluate the model performance, the total root mean square error (RMSE), mean bias

340

error (MBE), and coefficient of determination (r2) are used. The lower of the values of MBE

341

and RMSE, and the higher of the values of r2 near 1.0, the better the model performs.

342

-------------------------------------------------------------------------

343

Place Table 2 here

344

-------------------------------------------------------------------------

345

4. Results

346

4.1 Diurnal cycle of energy fluxes

347

Ten-day time series of the diurnal cycle of simulated and observed net radiation (Rn),

348

latent heat (LE) and sensible heat (H) in the two ecosystems are shown in Fig.3. Results show

349

that the energy is mostly used in the form of evaporative flux and the LE is much greater than

350

H for the two sites. The coupled model simulates the magnitude and pattern of diurnal cycle

351

for Rn, LE and H relatively well for both the two sites, and extremely well for Rn. The

352

simulated Rn captures the diurnal variation quite well agree with the observation, while Rn is

353

overestimated at noon in some days with the maximum error of 61.82 Wm-2. The LE is also

354

simulated quite well except for some errors at noon and overestimation after 6:00 pm. The H

355

is not modeled as satisfactorily as LE. In some of the days, the H is overestimated or

356

underestimated at noon and the maximum bias is around 50%. The simulation error of H

357

may be related to the estimation of canopy and ground surface temperature for the reason that

358

H is determined by the difference of air temperature and canopy or ground surface

17

359

temperature. Besides, the energy imbalance of flux measurements may also contribute to

360

some bias in the simulation. Mo et al. (2012) indicate that the possible reason is that the

361

sensible heat is usually fluctuates irregularly over canopy and less efficient to capture its

362

variation. The simulated Rn is not imbalance with the combination of simulated LE and H. The

363

possible reason is that the residual energy is partitioned into soil heat flux which increases the

364

soil temperature and the overestimated Rn due to low albedo.

365

-------------------------------------------------------------------------

366

Place Fig.3 here

367

-------------------------------------------------------------------------

368

The half-hourly values of observed energy fluxes are compared with the simulated

369

values by the coupled model in the LP and OF sites shown in Fig.4. In the LP site, results show

370

that net radiation is predicted quite satisfactorily. Comparing to the observed Rn, the MBE is

371

10.15 Wm-2. The r2 between the observed and simulated series is quite high with the value of

372

0.99 and the RMSE value is 56.08 Wm-2. For the heat components, the MBE is -0.34 Wm-2 for

373

the simulated mean LE with the RMSE value 40.18 Wm-2, while the r2 value is a little lower

374

than that of Rn with the value about 0.93. The sensible heat H is modeled not as well as the

375

LE. The MBE is about 1.78 Wm-2 and the r2 and RMSE values are 0.90 and 45.47 Wm-2,

376

respectively.

377

In the OF site, the coupled model performs quite well in simulating Rn and LE, while H is

378

not simulated as well as the other two fluxes. The MBE of Rn is 8.86 Wm-2 comparing to the

379

observed value. The r2 value is 0.99 and the RMSE value is 34.52 Wm-2. For the heat

380

components, the simulated LE is also generally consistent with the observation with a low

18

381

MBE value of -13.85 Wm-2, the r2 value is 0.92 and the RMSE value is 40.26 Wm-2.

382

Comparing to the observation, the MBE of the simulated H is higher with a value of 30.35

383

Wm-2, while the r2 value is 0.94 and the RMSE value is 43.43 Wm-2.

384

-------------------------------------------------------------------------

385

Place Fig.4 here

386

-------------------------------------------------------------------------

387

4.2 Diurnal cycle of net ecosystem exchange

388

The coupled model simulates the carbon fluxes of photosynthesis, plants including leaf,

389

wood and root respiration and litter and soil respiration. For the lack of the detail observation,

390

The half-hourly values of observed gross primary production (GPP) and net ecosystem

391

exchange (NEE) are compared with the simulated values by the coupled model in the LP and

392

OF sites shown in Fig.5.

393

In the LP site, the MBE value between simulated and observed NEE is -0.43 mol C

394

m-2s-1, the r2 and RMSE are 0.76 and 5.54 mol C m-2s-1, respectively. The MBE for the

395

simulated GPP is 2.90 mol C m-2s-1, and the r2 and RMSE are 0.83 and 6.25 mol C m-2s-1,

396

respectively. In the OF site, the MBE is higher than that of LP site, and the value is -1.47

397

mol C m-2s-1, the r2 and RMSE are 0.88 and 2.68 mol C m-2s-1, respectively. For GPP, the

398

MBE is lower than that of LP site with the value of 1.74 mol C m-2s-1, the r2 value is 0.89

399

and the RMSE value is 3.57 mol C m-2s-1.

400

As shown in Fig.5, the agreement between the simulation and observation of NEE and

401

GPP is generally satisfactory. From the results of NEE, we can see it is underestimated in the

402

night which is due to the bias of the simulation of respiration, which was reported in numbers

19

403

of studies (Baldocchi, 1997, Baldocchi, 2003). For GPP, it is overestimated in the daytime,

404

especially at noon in the OF site. This may be due to the overestimation of observed solar

405

radiation for photosynthesis.

406

-------------------------------------------------------------------------

407

Place Fig.5 here

408

-------------------------------------------------------------------------

409

4.3 Diurnal cycle of evapotranspiration

410

Evapotranspiration (ET) including transpiration, soil evaporation and canopy interception

411

evaporation are simulated and the comparisons of the results to observed evapotranspiration

412

are shown in Fig.6. In the LP site, the MBE is 0.01mm, the r2 and RMSE between predicted

413

and observed ET series are 0.95 and 0.02 mm, respectively. In the OF site, the MBE value is

414

-0.01mm. The r2 and RMSE comparing the simuated and observed ET series are 0.92 and

415

0.03mm, respectively. The results indicate that the model could capture the ET diurnal

416

variation quite well.

417

Fig.6 also shows that the coupled model overestimates ET at noon in some days, which

418

maybe relating to a higher absorption simulation of solar radiation or lower canopy resistance

419

for water flux.

420

-------------------------------------------------------------------------

421

Place Fig.6 here

422

-------------------------------------------------------------------------

423

-------------------------------------------------------------------------

424

Place Table 3 here

20

425

-------------------------------------------------------------------------

426

4.4 Seasonal cycle of energy fluxes

427

Fig.7 shows the seasonal cycle of daily averaged energy fluxes in the LP site. The

428

agreement between predicted and observed Rn is quite good with a small RMSE value of

429

17.01 Wm-2, MBE value of 2.36 Wm-2 and a high r2 value of 0.97. The simulated Rn is lightly

430

overestimated with a 2.31% bias. The seasonal variation in Rn throughout the whole year

431

indicates that the variations of observed and simulated values are quite consistent. The

432

performances of LE and H are not as good as Rn , while LE is simulated better than H. The

433

MBE between the simulation and observation of LE and H are -1.75 Wm-2 and 11.52 Wm-2,

434

respectively. The RMSE values for the two energy components are 15.02Wm-2 and

435

25.07Wm-2 respectively and the r2 values are 0.94 and 0.70 respectively. Fig.7 also shows that

436

LE is generally overestimated in nongrowing season, while it underestimated in the growing

437

season. As less energy is allocated to LE in the simulation, H is overestimated particularly

438

from August to October.

439

-------------------------------------------------------------------------

440

Place Fig.7 here

441

-------------------------------------------------------------------------

442

Fig.8 shows the seasonal cycle of daily averaged energy fluxes in the OF site. The

443

seasonal simulated Rn is also quite consistent with the observation with a small RMSE value

444

of 11.80 Wm-2, a high r2 of 0.98 and MBE value of 5.30 Wm-2. The simulation of LE and H is

445

also relatively well during the period. Comparing to the observed LE, the MBE is -5.51 Wm-2,

446

the r2 is 0.84 and RMSE value is 19.18 Wm-2. For the sensible heat H, the MBE is 15.77

21

447

Wm-2, the r2 is 0.75 and RMSE value is 21.65 Wm-2. Fig.8 also shows LE is generally

448

underestimated in non-growing season but overestimated during August to October, while H

449

is generally overestimated during the simulation period.

450

-------------------------------------------------------------------------

451

Place Fig.8 here

452

-------------------------------------------------------------------------

453

4.5 Seasonal cycle of net ecosystem exchange

454

Fig.9(a) and Fig.9(b) presents the comparison of daily observed and simulated GPP and

455

NEE in the LP site. The MBE is -0.03mol C m-2s-1 comparing to the observed GPP values.

456

The r2 and RMSE are 0.86 and 1.68 mol C m-2s-1, respectively. The MBE value is -1.56mol

457

C m-2s-1 between simulated and observed NEE. The r2 and RMSE are 0.62 and 1.81 mol C

458

m-2s-1, respectively.

459

The comparison of daily observed and simulated GPP and NEE in the OF site are shown

460

in Fig.9(c) and Fig.9(d). The MBE of simulated and observed GPP is -0.31mol C m-2s-1. The

461

r2 is relatively lower than that of LP site with the value of 0.83, while the RMSE value is 1.23

462

mol C m-2s-1. Compare of the simulated and observed mean NEE, the MBE is -0.15mol C

463

m-2s-1, the r2 value is 0.64 and the RMSE value is 1.27 mol C m-2s-1, respectively.

464

From the results, the model could capture the transition between the carbon source and

465

carbon sink quite well with a good agreement with the observation except for some large

466

errors in certain month. In January and March, NEE at the LP site is underestimated which

467

may relate to the underestimation of soil temperature and soil moisture that stressed NEE,

468

which has been well known in forest ecosystem (Ohta et al., 2008), while in October,

22

469

November and December, the NEE is overestimated is also related to the overestimation of

470

soil moisture that reduce the stress on NEE. In addition, uncertainties in the plant physiology

471

and soil carbon parameters would also bring errors in the simulated NEE. Generally, NEE in

472

the OF site is underestimated in nongrowing season and overestimated in the growing season.

473

-------------------------------------------------------------------------

474

Place Fig.9 here

475

-------------------------------------------------------------------------

476

4.6 Seasonal cycle of evapotranspiration

477

The variations of daily simulated and observed evapotranspiration in the two sites are

478

shown in Fig.10. In the LP site, the MBE between the simulated and observed ET is 0.11mm

479

and the r2 and RMSE values are 0.94 and 0.48 mm, respectively. The comparison between

480

observed and simulated ET series show overestimation in non-growing season and

481

underestimation in growing season. In the OF site, the MBE for ET is -0.02mm, the r2 and

482

RMSE values are 0.90 and 0.51 mm, respectively. Different from the LP site, the seasonal

483

changes show the ET is underestimated in nongrowing season and overestimated in growing

484

season. The bias of the simulation of ET may be attributed to the estimation stomatal

485

conductance, LAI and soil resistance.

486

-------------------------------------------------------------------------

487

Place Fig.10 here

488

-------------------------------------------------------------------------

489

4.7 Soil moisture and soil temperature

490

The simulation of soil moisture and soil temperature is important for determining the

23

491

water and carbon budgets such as evapotranspiration and soil respiration. The simulated soil

492

moisture and temperature at different soil depths in the two sites are shown in Fig.11. From

493

the results, the most sensitive soil layers are in the upper 0.3 m. The soil moisture and

494

temperature changes exhibit a time lag as the soil depth increases.

495

The seasonal changes of soil moisture 30cm and soil temperature Ts,30cm in 30cm depth at

496

daily time step for the two sites are shown in Fig. 12. The results show that the simulation for

497

soil moisture and temperature are quite effective for the two sites. In the LP site, the MBE

498

values are about -0.01 m3 m-3 and -0.39 C, the RMSE values are also small at 0.03 m3 m-3 and

499

1.91 C, and the r2 values are 0.90 and 0.99, respectively. In the OF site, the MBE values are

500

about -0.03 m3 m-3 and -0.47 C, the r2 values are 0.92 and 0.99, the RMSE values are also

501

small at 0.08 m3 m-3 and 2.02C, respectively.

502

In the two sites, 30cm is underestimated during January, February and March, which is

503

related to the bias of simulation of evapotranspiration and runoff in these months from the

504

perspective of water balance. In addition, the heterogeneity in the soil structure and physical

505

properties also result in complex distributions of soil water (Cuenca et al., 1997) which

506

increase the uncertainties of soil moisture simulation. There is a satisfactory agreement

507

between the observed and simulated Ts,30cm from March to October. In other months, Ts,30cm is

508

underestimated especially in January and February. The maximum error is -6.1C, which

509

result in less soil respiration and may be the reason for the underestimation of NEE. Overall,

510

the model reasonably captures the variations in soil moisture and temperature.

511

-------------------------------------------------------------------------

512

Place Fig.11 here

24

513

-------------------------------------------------------------------------

514

-------------------------------------------------------------------------

515

Place Fig.12 here

516

-------------------------------------------------------------------------

517

-------------------------------------------------------------------------

518

Place Table 4 here

519

-------------------------------------------------------------------------

520

5. Discussion

521

5.1 Sensitivities to parameters

522

In general, the coupled model could simulate the water, energy and carbon fluxes well over

523

the period in this study. However, previous study (Baldocchi and Wilson, 2001; Kothavala et

524

al., 2005) argues that despite close agreement between simulated and measured data, one

525

cannot claim complete success of the model for the reason that there may be numerous

526

sources of errors including errors in the model parameters, the forcing variable, inadequate

527

representation of physical processes and the flux underestimation problem. Therefore, the

528

coupled model needs to be investigated to further analyze the applicability. In this section, the

529

model sensitivities to the parameters due to parameterization are analyzed.

530

The mono-factor analysis is used to analyze the sensitivities of ET, GPP and net

531

ecosystem production (NEP, NEP=-NEE) to certain parameters, i.e. the analyzed parameter is

532

selected in a reasonable range from the lower value to upper value while other parameters are

533

fixed. The analyzed parameters are shown in Table 1. The sensitivities of ET, GPP and NEP

534

to the main runoff, photosynthetic and stomatal conductance parameters including g1, g2,

25

535

Vcmax, Jmax,  and a1 are shown in Fig.13. The results shown that while g1 varies from 0 and 1,

536

the total ET and GPP changes about 1.3% and 4.6%, while the total NEP shows a big change.

537

The runoff parameter g2 varies from 0 to 1, the total ET, GPP and NEP all show small

538

changes lower than 5%. The photosynthetic parameter Vcmax varies between 50 and 150, all

539

the three elements show small changes, however, the photosynthetic parameter Jmax varies

540

from 50 to 100 could result in large changes for ET about 17.4% and GPP about 37.3%,

541

especially for NEP nearly six times higher. The stomatal conductance parameter  changes

542

from 0.05 to 0.2 could also result in large changes approximately 131%, 287% and 126% for

543

ET, GPP and NEP, respectively. The parameter a1 ranges between 15 to 18 could result in

544

6.8%, -2.0% and 6.9% changes for ET, GPP and NEP, respectively. The results above show

545

that the water and carbon fluxes are quite sensitive to plant physiological controls such as the

546

maximum electron transport rate, quantum efficiency of electron transport and runoff

547

parameter.

548

-------------------------------------------------------------------------

549

Place Fig.13 here

550

-------------------------------------------------------------------------

551

5.2

Anaysis of water use efficiency (WUE)

552

Water use efficiency (WUE), the ratio of water loss to carbon gain, is a key characteristic

553

of ecosystem function that is central to the water, energy and carbon cycles (Beer et al. 2009).

554

WUE is defined as the ratio of photosynthesis and transpiration rates. The relationships

555

between photosynthesis and transpiration are investigated firstly. Fig.14 shows the

556

relationship between photosynthesis and transpiration in the LP and OF sites. The results

26

557

show a close relationship between the two elements with high correlation coefficients of

558

0.9742 and 0.9259, respectively. The photosynthesis increases and reaches a convergence as

559

the increase of transpiration at both sites.

560

-------------------------------------------------------------------------

561

Place Fig.14 here

562

-------------------------------------------------------------------------

563

Fig.15 shows the WUE in diurnal and seasonal changes for the two sites. The diurnal variation

564

of the WUE at half-hourly time step shows a U-bend curve in Fig.15(a). In the night time,

565

WUE is zero without carbon assimilation. It reaches a peak value after the sunrise and then

566

reduces to a valley value at noon; after that it increases and reaches another peak value which

567

is mostly lower than the peak value in the morning. This is due to the reason that changing

568

rate of transpiration is faster than that of photosynthesis. Fig.15(b) shows that seasonal

569

variation of the daily WUE with the values mostly between 5~10 mmol C/mol H2O. The results

570

show that the WUE reduces firstly and then increases around July, indicating a high value in

571

winter and a low value in summer, showing an opposite changing trends of temperature,

572

evapotranspiration and net radiation. Comparing of the two sites in Fig.16, generally, the

573

WUE is higher in the Loblolly Pine site than that of the Open Field site, which indicates a

574

higher water use efficiency of forest ecosystem than grass ecosystem.

575

-------------------------------------------------------------------------

576

Place Fig.15 here

577

-------------------------------------------------------------------------

578

-------------------------------------------------------------------------

27

579

Place Fig.16 here

580

-------------------------------------------------------------------------

581

5.3 Comparisons of water, energy and carbon fluxes over the two sites

582

Comparisons of water, energy and carbon fluxes over the LP site and OF site are shown

583

in Fig.17. The hydrological components in the LP site and OF site including transpiration

584

(Tr), evaporation from canopy intercepted water (Ec), soil evaporation (Es) and surface runoff

585

(Rs) are shown in Fig.17(a) and Fig.17(b), respectively. Total surface runoff and soil

586

evaporation in the OF site are much higher than those in the LP site, while transpiration and

587

the total ET in the LP site is much higher than that in the OF site.

588

Fig.17(c) and Fig.17(d) are the energy components including sensible (H) and latent heat

589

(LE) fluxes, soil heat flux (G) in the LP site and OF site, respectively. The LE in the LP site is

590

higher than that in the OF site, which is caused by the higher evapotranspiration absorbing

591

more heat in the LP site. Correspondingly, the H is higher in the OF site from the perspective

592

of energy balance. Th soil heat fluxes have small differences in the two sites are. The monthly

593

variances of soil heat fluxes show that during September to January, the soil heat flux are

594

negative which indicating a output heat to the air.

595

Fig.17(e) and Fig.17(f) are the carbon fluxes including ecosystem respiration (RE) and

596

net ecosystem production (NEP, NEP=-NEE, GPP=NEP+RE) in the LP site and OF site,

597

respectively. Positive value of NEP indicates carbon sink and the negative value indicates

598

carbon source. From the results, we can find that the total GPP and RE are higher in the LP

599

site than those in the OF site. The LP site ecosystem are almost carbon sink during the whole

600

year, while in the OF site, it shows a carbon source during May to September.

28

601

-------------------------------------------------------------------------

602

Place Fig.17 here

603

-------------------------------------------------------------------------

604

5.4 Comparisons between coupled model with single models

605

In order to validate the importance of the coupling of hydrological processes and

606

biogeochemical processes, comparisons of the simulation of water, energy and carbon fluxes

607

by coupled model and single model are conducted to analyze the improvement of the coupled

608

model. The Loblolly Pine site is taken as the example in this section. Results of the single

609

model simulation (SIM-S) and the coupled model simulation (SIM-C) are compared and

610

analyzed. The comparisons of simulated evapotranspiration (ET), soil moisture and GPP,

611

which are three main elements of water carbon cycles, are conducted by the SIM-S and

612

SIM-C.

613

Fig.18 shows the comparisons of simulated ET, soil moisture and GPP by SIM-S and

614

SIM-C via observed values. The results show that the simulated ET and GPP are

615

overestimated about 35.7% and 29.6%, while soil moisture is underestimated about 10.3% by

616

SIM-S. And the relative error reduce to 5.6%, 0.6% and 2.7% for ET, GPP and soil moisture

617

by SIM-C. The r2 between simulated and observed values of ET, GPP and soil moisture by

618

SIM-S are 0.83, 0.75 and 0.82, respectively. And the r2 values increase to 0.94, 0.86 and 0.90

619

by SIM-C. The RMSE between simulated and observed values of ET, GPP and soil moisture

620

by SIM-S are 1.2mm, 2.9 mol C m-2s-1 and 0.06 m3 m-3, respectively. And the RMSE values

621

reduce to 0.48mm, 1.68 mol C m-2s-1 and 0.03 m3 m-3 by SIM-C. The MBE between

622

simulated and observed values of total ET, GPP and soil moisture by SIM-S are 0.72mm, 1.7

29

623

mol C m-2s-1 and 0.04 m3 m-3, respectively. And the MBE values reduce to 0.11mm, -0.03

624

mol C m-2s-1 and -0.01m3 m-3 by SIM-C. From the results above, we find that the coupled

625

model could improve the performances in simulating water and carbon fluxes obviously than

626

the single models. The main reason is that the coupled model could characterize the

627

relationship of transpiration and photosynthesis. The single model could not capture the

628

effects of stomatal behavior on evapotranspiration and evapotranspiration is mainly estimated

629

via climate conditions and soil moisture without considering the vegetation effects.

630

Meanwhile, the effects of land surface hydrological processes on carbon assimilation could

631

also not be considered in the single model, which would also bring systematic error in

632

simulating the carbon flux.

633

-------------------------------------------------------------------------

634

Place Fig.18 here

635

-------------------------------------------------------------------------

636

6. Conclusions

637

The hydrological cycle and vegetation dynamics/carbon cycle interact strongly with

638

each other. We developed a physical process-based model by coupling a hydrological model

639

(DTVGM) and a biogeochemical model (CASACNP) to simulate the water, energy and

640

carbon fluxes and analyze their relationships over two different ecosystems. The coupled

641

model DTVGM-CASACNP is tested well at both sites. The results show that the coupled

642

model could simulate net radiation, latent and sensible heat fluxes, evapotranspiration, GPP,

643

NEE, soil moisture and soil temperature quite well at both the diurnal and seasonal scales by

644

comparing the measured and simulated water, energy and carbon fluxes. Sensitivity analysis

30

645

illustrates that evapotranspiration, GPP, net ecosystem production are quite sensitive to plant

646

physiological controls such as the maximum electron transport rate, quantum efficiency of

647

electron transport and runoff parameter. The coupling characteristics of water and carbon

648

cycles are also investigated by analyzing the water use efficiency (WUE) at different time

649

scales and the water, energy and carbon fluxes are compared over the two ecosystems.

650

Since the complexity of the hydrological and ecological processes, the coupled model

651

still needs to be further improved by including the phenology and dynamic root growth at

652

different soil depths etc. Parameterization is still need to be further investigated for the

653

applications of the coupled model. In addition, the coupled model is only tested at two sites,

654

the results obtained in this study may not represent sites with different environment

655

conditions. At other sites and regional verification is still need to be conducted over different

656

ecosystems in further studies.

657

Acknowledgements

658

This study was supported by National Natural Science Foundation of China (No. 51809008),

659

the Fundamental Research Funds for the Central Universities, South-Centtral University for

660

Nationalities(CZY18042), Hubei Provincial Natural Science Foundation of China

661

（ 2018CFB123), Geology and Mineral Resources Survey Project: Ecological Configuration

662

and Global Strategy of China Water Resources (DD20190652). We acknowledge the Duke

663

Forest-Loblolly Pine and Duke Forest-Open Field sites of AmeriFlux for the data to test the

664

coupled model.

665

References

666

Arora, V., 2002. Modeling vegetation as a dynamic component in soil-vegetation-atmosphere

31

667

transfer schemes and hydrological models. Reviews of Geophysics 40 (2), 3-1-3-26.

668

Baldocchi, D.D., 1997. Measuring and modeling carbon dioxide and water vapor exchange

669

over a temperate broadleaved forest during the 1995 summer drought. Plant Cell and

670

Environment 20, 1108-1122.

671

Baldocchi, D.D., 2003. Assessing the eddy covariance technique for evaluating carbon

672

dioxide exchange rates of ecosystems: past, present and future. Global Change Biology. 9,

673

479-492.

674

Beer, C., Ciais, P., Reichstein, M., Baldocchi, D., Law, B.E., Papale, D., Soussana, J.-F.,

675

Ammann, C., Buchmann, N., Frank, D., Gianelle, D., Janssens, I.A., Knohl, A., Köstner,

676

B., Moors, E., Roupsard, O., Verbeeck, H., Vesala, T., Williams, C.A., Wohlfahrt, G.,

677

2009. Temporal and among-site variability of inherent water use efficiency at the

678

ecosystem level. Global Biogeochemical Cycles 23 (2), DOI: 10.1029/2008GB003233.

679

Beven, K., 2002. Towards an alternative blueprint for a physically based digitally simulated

680 681 682

hydrologic response modelling system. Hydrological Processes, 16 (2), 189-206. Celia M A, Bouloutas E T, Zarba R L, et al. 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research, 26(7): 1483-1496

683

Collatz, G.J., Ball, J.T., Grivet, C., Berry, J.A., 1991. Physiological and environmental

684

regulation of stomatal conductance, photosynthesis and transpirationda model that

685

includes a laminar boundary-layer. Agricultural and Forest Meteorology 54 (2-4),

686

107-136.

687

Collatz, G.J., Ribas-Carbo, M., Berry, J.A., 1992. Coupled photosynthesis-stomatal

688

conductance model for leaves of C4 plants. Australian Journal of Plant Physiology 19,

32

689 690 691 692 693 694 695 696 697

519-538. Cuenca, R.H., Stangel, D.E., Kelly, S.F., 1997. Soil water balance in a boreal forest. Journal of Geophysical Research, 102 (D24), 29355-29365. Davidson, E.A., Janssens, I.A., 2006. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. Nature 440, 165-173. De Pury, D.G.G., Farquhar, G.D., 1997. Simple scaling of photosynthesis from leaves to canopies without the errors of big-leaf models. Plant, Cell and Environment 20, 537-557. Farquhar, G.D., von Caemmerer, S., Berry, J.A., 1980. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149, 78-90.

698

Flerchinger, G.N. and K.E. Saxton. 1989. Simultaneous heat and water model of a freezing

699

snow-residue-soil system I. Theory and Development. Transactions of the ASAE,

700

32(2):565-571

701

Fung, I.Y., Doney, S.C., Lindsay, K., John, J., 2005. Evolution of carbon sinks in a changing

702

Climate, Proceedings of National Academy Sciences of USA. 102 (32), 11201-11206.

703

Houlton, B.Z., Wang, Y.P., Vitousek, P.M., Field, C.B., 2008. A unifying framework for

704

dinitrogen fixation in the terrestrial biosphere. Nature, 454 (7202), 327-330.

705

Jonhansen O. 1975. Thermal Conductivity of Soils. Ph D. thesis, University of Trondheim.

706

Leuning, R., 1995. A critical appraisal of a combined stomatal-photosynthesis model for C3

707

plants. Plant, Cell and Environment 18, 339-355.

708

Li, L., Xia, J., Xu, C.Y., Singh, V.P., 2010. Evaluation of the subjective factors of the GLUE

709

method and comparison with the formal Bayesian method in uncertainty assessment of

710

hydrological models. Journal of Hydrology, 390, 210-221.

33

711

Lohse, K.A., Brooks, P.D., Mclntosh, J.C., Meixner, T., Huxman, T.E., 2009. Interactions

712

Between Biogeochemistry and Hydrologic Systems. Environment and Resources 34,

713

65-96.

714 715 716 717

Manzoni, S., Porporato A., 2009. Soil carbon and nitrogen mineralization: Theory and models across scales. Soil Biology and Biochemistry, 41 (7), 1355-1379. McMurtrie, R.E., Comins, H.N., 1996. The temporal response of forest ecosystems to doubled atmospheric CO2 concentration. Global Change Biology, 2 (1), 49-57.

718

Mo, X., Liu, S., Lin Z., Xu, Y., Xiang, Y., McVicar T.R., 2005. Prediction of crop yield,

719

water consumption and water use efficiency with a SVAT-crop growth model using

720

remotely sensed data on the North China Plain. Ecological Modelling, 183, 301-322.

721

Mo, X., Liu, S., Lin, Z., 2012. Evaluation of an ecosystem model for a wheat–maize double

722

cropping system over the North China Plain. Environmental Modelling & Software 32,

723

61-73.

724

Moulin, S., Kergoat, L., Viovy, N., Dedieu, G., 1997. Global-scale assessment of vegetation

725

phenology using NOAA/AVHRR satellite measurements, Journal of Climate 10 (6),

726

1154-1170.

727

Ohta, T., Maximov, T.C., Dolman, A.J., Nakai, T., van der Molen, M.K., Kononov, A.V.,

728

Maximov, A.P., Hiyama, T., Iijima, Y., Moors, E.J., Tanaka, H., Toba, T., Yabuki, H.,

729

2008. Interannual variation of water balance and summer evapotranspiration in an eastern

730

Siberian larch forest over a 7-year period (1998-2006). Agricultural and Forest

731

Meteorology 148 (12), 1941-1953.

732

Park, H., Iijima, Y., Yabuki, H., Ohta, T., Walsh, J., Kodama, Y., Ohata, T., 2011. The

34

733

application of a coupled hydrological and biogeochemical model (CHANGE) for

734

modeling of energy, water, and CO2 exchanges over a larch forest in eastern Siberia.

735

Journal of Geophysical Research, 116, D15102, doi:10.1029/2010JD015386.

736 737 738 739

Parmesan, C., Yohe, G., 2003. A globally coherent fingerprint of climate change impacts across natural systems, Nature 421, 37-42. Pritchard, S.G., Rogers, H.H., Prior S.A., Peterson, C.M., 1999. Elevated CO2 and plant structure: A review. Global Change Biology 5 (7), 807-837.

740

Randerson, J.T., Thompson, M.V., Conway, T.J., Fung, I.Y., Field, C.B., 1997. The

741

contribution of terrestrial sources and sinks to trends in the seasonal cycle of atmospheric

742

carbon dioxide. Global Biogeochemical Cycles, 11, 535-560.

743

Randerson, J.T., Enting, I.G., Schuur, E.A.G., Caldeira, K., Fung, I.Y., 2002. Seasonal and

744

latitudinal variability of troposphere 14CO2: post bomb contributions from fossil fuels,

745

oceans, the stratosphere, and the terrestrial biosphere. Global Biogeochemical Cycles, 16

746

(4), 59-1-59-19.

747

Reynolds, J.F., Kemp, P.R., Ogle, K., Fernández R.J., 2004. Modifying the 'pulse-reserve'

748

paradigm for deserts of North America: precipitation pulses, soil water, and plant

749

responses. Oecologia 141 (2), 194-210.

750

Riveros-Iregui, D.A., McGlynn, B.L., Emanuel, R.E., Epstein, H.E., 2012. Complex terrain

751

leads to bidirectional responses of soil respiration to inter-annual water availability.

752

Global Change Biology, 18 (2), 749-756.

753

Rutter, A.J., Kershaw, K.A., Robins, P.C., Morton, A.J., 1972. A predictive model of rainfall

754

interception in forests, 1. Derivation of the model from observations in a plantation of

35

755

Corsican pine. Agricultural Meteorology, 9, 367-384.

756

Sellers, P.J., Mintz, Y., Sud, Y.C., Dalcher, A., 1986. A simple biosphere model (SiB) for use

757

within general circulation models. Journal of the Atmospheric Sciences 43 (6), 505-531.

758

Shen, C., Niu, J., Phanikumar M.S., 2013. Evaluating controls on coupled hydrologic and

759

vegetation dynamics in a humid continental climate watershed using a subsurface-land

760

surface processes model. Water Resources Research 49, 1-21, doi:10.1002/wrcr.20189.

761

Shi, X., Mao, J., Wang, Y., Dai, Y., Tang, X., 2011. Coupling a terrestrial biogeochemical

762

model to the common land model. Advances in Atmospheric Sciences, 28 (5), 1129-1142.

763

Thornton, P.E., Rosenbloom, N.A., 2005. Ecosystem model spin-up: Estimating steady state

764

conditions in a coupled terrestrial carbon and nitrogen cycle model. Ecological Modelling

765

189, 25-48.

766

van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic

767

conductivity of unsaturated soils. Soil Science Society of America Journal 44 (5),

768

892-898.

769 770

Volk, M., Niklaus, P.A., Korner, C., 2000. Soil moisture effects determine CO2 responses of grassland species. Oecologia 125 (3), 380-388.

771

Wagener, T., Sivapalan, M., Troch, P.A., McGlynn, B.L., Harman, C.J., Gupta, H.V., Kumar

772

P., Rao, P.S.C., Basu, N.B., Wilson, J.S., 2010. The future of hydrology: An evolving

773

science

774

doi:10.1029/2009WR008906.

for

a

changing

world.

Water

Resources

Research

46,

W05301,

775

Wang, G., Xia, J., Chen, J., 2009. Quantification of effects of climate variations and human

776

activities on runoff by a monthly water balance model: A case study of the Chaobai River

36

777

basin

778

doi:10.1029/2007WR006768.

in

northern

China.

Water

Resources

Research,

45,

W00A11,

779

Wang, Y., Leuning, R., 1998. A two-leaf model for canopy conductance, photosynthesis and

780

partitioning of available energy. I: Model description and comparison with a multi-layered

781

model. Agricultural and Forest Meteorology 91, 89-111.

782

Wang, Y.P., Houlton, B.Z., Field, C.B., 2007. A model of biogeochemical cycles of carbon,

783

nitrogen and phosphorus including symbiotic nitrogen fixation and phosphatase

784

production.

785

doi:10.1029/2006GB002797.

786 787

Global

Biogeochemical

Cycles,

21

(1),

GB1018,

Wang, Y.P., Law, R.M., Pak, B., 2010. A global model of carbon, nitrogen and phosphorus cycles for the terrestrial biosphere. Biogeosciences, 7, 2261-2282.

788

Wu, Y., Liu, S., Qiu, L., Sun, Y., 2016. SWAT-DayCent coupler: an integration tool for

789

simultaneous hydro-biogeochemical modeling using SWAT and DayCent. Environmental

790

Modelling & Software, 86, 81–90.

791 792

Xia, J., 2002. A system approach to real-time hydrologic forecast in watersheds. Water International, 27 (1), 87-97.

793

Xia, J., O’Connor, K.M., Kachroo, R.K., Liang, G.C., 1997. A non-linear perturbation model

794

considering catchment wetness and its application in fiver flow forecasting. Journal of

795

Hydrology, 200 (1-4), 164-178.

796

Xia, J., Wang, G., Tan, G., Ye, A., Huang, G.H., 2005. Development of distributed

797

time-variant gain model for nonlinear hydrological systems. Science in China Ser. D

798

Earth Sciences, 48 (6), 713-723.

37

799

Zhan, C.S., Song, X.M., Xia, J., Tong, C., 2013. An efficient integrated approach for global

800

sensitivity analysis of hydrological model parameters. Environmental Modelling &

801

Software, 41, 39-52.

802

Zhang, X., Friedl, M.A., Schaaf, C.B., 2006. Global vegetation phenology from Moderate

803

Resolution Imaging Spectroradiometer (MODIS): Evaluation of global patterns and

804

comparison with in situ measurements. Journal of Geophysical Research 111, G04017,

805

doi:10.1029/2006JG00021.

806

Zhao, F., Wu, Y., Qiu, L., Sivakumar, B., Zhang, F., Sun, Y., et al., 2018. Spatiotemporal

807

features of the hydro-biogeochemical cycles in a typical loess gully watershed. Ecological

808

Indicators, 91, 542–554.

809

Zhao, F., Wu, Y., Bellie, S., Long, A., Qiu, L., Chen, J., et al., 2019. Climatic and hydrologic

810

controls on net primary production in a semiarid loess watershed. Journal of Hydrology,

811

568, 803–815.

812

Abstract: Understanding the mechanism of the interactions between vegetation dynamics and the

813

water cycle is rather important for determining global and regional water and carbon budgets. In

814

this paper, a physically-based model integrating land-surface hydrological and biogeochemical

815

processes by coupling a hydrological model and a biogeochemical model is developed to simulate

816

the water, energy and carbon fluxes. The model is validated against observed biometric,

817

eddy-covariance flux, soil moisture and temperature data over two different ecosystems. Results

818

show that the model could simulate the vegetation physiological and hydrological processes

819

satisfactorily including net radiation, latent heat, gross primary production (GPP), net ecosystem

820

exchange

and

soil

moisture

and

temperature.

38

Sensitivity

analysis

illustrates

that

821

evapotranspiration, GPP, net ecosystem production are quite sensitive to plant physiological

822

controls such as the maximum electron transport rate, quantum efficiency of electron transport and

823

runoff parameter. Moreover, results show that a close relationship between photosynthesis and

824

transpiration. Water use efficiency shows a U-bend curve at different time steps and it also

825

indicates a higher value in forest ecosystem than that of grass ecosystem. In the forest ecosystem,

826

evapotranspiration is higher and the surface runoff is lower. Meanwhile it is almost carbon sink

827

during the whole year, while in the grass ecosystem, it shows a carbon source during May to

828

September. This study could provide an effective model for the simulation of water-carbon cycles.

829 830 831

Declaration of interests

832 833 834

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

835 836 837

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

838 839 840 841 842 843

Figure captions

844

Fig.1 Schematic of the coupled DTVGM-CASACNP model

39

845

Fig.2 Comparisons of observed net radiation and the sum of latent heat and sensible heat in the LP

846

(a) and OF (b) sites

847

Fig. 3 Time series of simulated net radiation, Rn (black line), latent heat LE (dash line), sensible

848

heat H (dot line) and observed Rn (diamonds), LE (squares) and H (circles) in the LP and OF sites

849

Fig. 4 Simulated vs. observed net radiation , latent heat, sensible heat in the LP and OF sites

850

Fig. 5 Time series of simulated and observed GPP and NEE in the LP and OF sites

851

Fig. 6 Time series of simulated ET (black line) and observed ET (circles) in the LP and OF sites

852

Fig. 7 Time series of simulated and observed daily net radiation, Rn (a), latent heat LE (b),

853

sensible heat H (c) in the LP site

854

Fig. 8 Time series of simulated and observed daily net radiation, Rn (a), latent heat LE (b),

855

sensible heat H (c) in the OF site

856

Fig. 9 Time series of simulated and observed daily GPP and NEE in the LP and OF sites

857

Fig. 10 Time series of simulated and observed daily ET in the LP and OF sites

858

Fig. 11 Simulated soil moisture and temperature at different soil depth in the LP and OF sites

859

Fig. 12 Time series of simulated and observed soil moisture and temperature in the LP and OF

860

sites

861

Fig. 13 Sensitivities of total accumulated ET, GPP and NEP to model key parameters

862

Fig.14 Relationships between transpiration (Tr) and photosynthesis (An) in the LP (a )and OF (b)

863

Fig.15 Diurnal and seasonal changes of WUE in the LP and OF sites

864

Fig.16 Comparison of the WUE in the LP and OF sites

865

Fig.17 Comparisons of water, energy and carbon fluxes in the LP and OF sites

866

Fig.18 Comparisons of ET, Soil moisture and GPP by SIM-S and SIM-C

40

867 868

Highlights:

869



developed

870 871



Water, energy and carbon fluxes are simulated satisfactorily by the coupled model

872 873

A coupled hydrological and biogeochemical model (DTVGM-CASACNP) is



Water and carbon fluxes are sensitive to maximum electron transport rate, quantum efficiency of electron transport and runoff parameter

874 875



Water use efficiency of forest ecosystem is higher than that of grass ecosystem

876



Coupled model could improve the performances of the simulation of water and

877

carbon fluxes than single models

878 879

Table 1 Descriptions of the two types of ecosystem sites Description

Loblolly Pine site

Open Field site

Latitude, Longitude

35.9782, -79.0942

35.9712, -79.0934

Elevation (m)

163

168

Vegetation IGBP

ENF (Evergreen Needleleaf Forests)

GRA (Grasslands)

Dominant stand

loblolly pine forest

Festuca arundinacea shreb

photosynthesis

C3

C3, a few C4

Canopy height

19.0

0.1～1.0

Leaf Area Index

2.0～5.0

1.0～3.0

880

41

881 882

Table 2 The key parameters in the coupled DTVGM-CASACNP model Description

LP site

OF site

runoff parameter, g1

0.2

0.5

runoff parameter, g2

0.95

0.5

Photosynthetic Rubisco capacity, Vcmax (mol m-2s-1)

81

65

maximum electron transport rate, Jmax (mol m-2s-1)

75

50

quantum efficiency of electron transport,  (mol electrons mol-1)

0.20

0.20

Parameter related to the intercellular CO2 concentration, a1

18

16

Parameter for stomatal sensitivity to vapor pressure, D0 (Pa)

1250

1250

Plant respiration, Rm (gC gN-1 day-1)

0.066

0.04

the fraction of NPP allocated to leaf, ac,leaf

0.42

0.30

the fraction of NPP allocated to wood, ac,wood

0.33

0.0

the fraction of NPP allocated to root, ac,root

0.25

0.70

turnover rate of leaf, 1/ leaf (year)

2.0

1.0

turnover rate of wood, 1/ wood (year)

70

1.0

turnover rate of root, 1/ root (year)

18

3.0

saturated volumetric water content (s)

0.45, 0.46

0.54, 0.46

residual volumetric water content (r)

0.07, 0.10

0.07, 0.10

Saturated hydraulic conductivity (mm/day) (Ks)

9.60, 6.15

9.60, 6.15

maximum canopy interception parameter (Kcmax)

1.6

0.0

soil resistance parameter a

3.5

3.5

42

soil resistance parameter b

2.3

2.3

soil resistance parameter c

400

1000

883 884 885

Table 3 The performances of the model

886

in simulating diurnal cycle of water, energy and

carbon fluxes over the two ecosystems LP site

OF site

r2

RMSE

MBE

r2

RMSE

MBE

Rn

0.99

56.08

10.15

0.99

34.52

8.86

H

0.90

45.47

1.78

0.94

43.43

30.35

LE

0.93

40.18

-0.34

0.92

40.26

-13.85

GPP

0.83

6.25

2.90

0.89

3.57

1.74

NEE

0.76

5.54

-0.43

0.88

2.68

-1.47

ET

0.95

0.02

0.01

0.92

0.03

-0.01

887 888 889

Table 4 The performances of the model

890

in simulating seasonal cycle of water, energy and

carbon fluxes over the two ecosystems LP site

OF site

r2

RMSE

MBE

r2

RMSE

MBE

Rn

0.97

17.01

2.36

0.98

11.80

5.30

H

0.70

25.07

11.52

0.75

21.65

15.77

43

LE

0.94

15.02

-1.75

0.84

19.18

-5.51

GPP

0.86

1.68

-0.03

0.83

1.23

-0.31

NEE

0.62

1.81

-1.56

0.64

1.27

-0.15

ET

0.94

0.48

0.11

0.90

0.51

-0.02

Soil moisture

0.90

0.03

-0.01

0.92

0.08

-0.03

Soil temperature

0.99

1.91

-0.39

0.99

2.02

-0.47

891 892

893

44

894

895

45

896

897

898 46

899

900

47

901

902

48

903

904

49

905

906

50

907

908

51

909

910

52