Implementation and application of a distributed hydrological model using a component-based approach

Environmental Modelling & Software 80 (2016) 245e258

Contents lists available at ScienceDirect

Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft

Implementation and application of a distributed hydrological model using a component-based approach Yinping Long a, b, Yaonan Zhang a, *, Dawen Yang c, Lihui Luo a a

Cold and Arid Regions Environmental and Engineering Research Institute of Chinese Academy of Sciences, Lanzhou 730000, China University of Chinese Academy of Sciences, Beijing 100049, China c State Key Laboratory of Hydroscience and Engineering, Department of Hydraulic Engineering, Tsinghua University, Beijing 100084, China b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 5 January 2015 Received in revised form 20 November 2015 Accepted 2 March 2016 Available online xxx

The application of a model to different study areas often requires that the model be modified to conform to specific characteristics, but this can be challenging due to the poor readability and reusability of the legacy codes. Component-based programming supported by a modelling framework provides a generic means to develop and modify models. This paper describes the development of a distributed hydrological model using a component-based modelling framework, which is implemented as a set of functional components that are integrated at runtime. The model was applied to runoff simulation in a large scale and data scarce alpine basin, and was further improved by incorporating a simple empirical soil freezing-thawing component. The results show that the componentised model reproduced the daily and monthly flow hydrograph with ‘good’ accuracy. The framework is shown to be flexible enough for model development and model modification. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Component-based programming Modelling framework Distributed hydrological model Runoff simulation

1. Introduction Hydrological models are important in facilitating our understanding of hydrological processes, performing hydrological forecasting, assessing the availability of water resources and making decisions in basin management. A hydrological model is usually developed to solve a specific problem in a specific study area, which limits its application in different regions and different spatial and temporal scales. For example, the Xinanjiang model is applicable in humid and semi-humid areas, but not in arid regions because it merely considers Dunne runoff generation mechanisms, whereas arid regions are dominated by Horton runoff (Hu et al., 2005). Consequently, for specific hydrological modelling purposes, models should be compared, selected and assessed according to the study area characteristics. In most cases (Dechmi et al., 2012; Khakbaz et al., 2011; Ouessar et al., 2009; Reyes et al., 1993), model modification is necessary when including a new process, excluding a process such as reservoir control that does not occur in a natural basin, or replacing one method with another. However, as most

* Corresponding author. Donggang West Road 320, Lanzhou, Gansu 730000, China. E-mail addresses: [email protected] (Y. Long), [email protected] (Y. Zhang). http://dx.doi.org/10.1016/j.envsoft.2016.03.001 1364-8152/© 2016 Elsevier Ltd. All rights reserved.

hydrologists are not professional programmers, model codes have poor readability and reusability, making them difficult to modify. A generic approach for coding hydrological models is needed. Component-based programming offers a nice solution to this problem, as it emphasises the decomposition of a complex system into a set of functional components with manageable complexity (Argent, 2004; Buahin and Horsburgh, 2015; Castronova and Goodall, 2010; Peckham et al., 2013). In a component-based modelling paradigm, the components, which communicate with each other via standard interfaces, can be flexibly and rapidly assembled into a new configuration to solve a specific problem. The ease of incorporating, substituting or deleting components makes it easy to improve the model, compare methods, assess the effects of one process on the whole system and so forth. Component-based programming builds on fundamental object-oriented programming concepts, such as encapsulation, inheritance and polymorphism. A component is implemented as a class, which encapsulates some properties and methods to express a particular physical process or common functionality. The main difference between the two programming paradigms, i.e., component-based and object-oriented programming, is that component-based programming has a modelling framework to link the components together (Peckham et al., 2013). Various frameworks are in development within the water and environmental modelling domains,

246

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

such as the Modular Modelling System (MMS) (Leavesley et al., 1996), the Earth System Modelling Framework (ESMF) (Hill et al., 2004), the Open Modelling Interface (OpenMI) (Moore and Tindall, 2005), the Cold Regions Hydrological Model (CRHM) platform(Pomeroy et al., 2007), the Object Modelling System (OMS) (David et al., 2013), the Community Surface Dynamics Modelling System (CSDMS) (Peckham et al., 2013), etc. These frameworks, which have defined standard interfaces for components and bring suites or libraries of components together, provide an avenue for easier modelling and can reduce the burden of repetitive coding tasks for environmental modellers (Argent et al., 2006). However, in recent years, the environmental modelling framework studies have devoted much more attention to the enhancement of the frameworks themselves, such as analysing the software design requirements of frameworks (Whelan et al., 2014), quantifying the computational overheads introduced by adopting component-based approaches (Buahin and Horsburgh, 2015; Anthony M. Castronova and Goodall, 2013), assessing the framework invasiveness (Lloyd et al., 2011) and extending the framework functionalities (A. M. Castronova and Goodall, 2010; Formetta et al., 2014). In contrast, only a few studies have focused on framework applications in terms of contributing model components to frameworks, building models or comparing different model formulations based on frameworks (Zhou et al., 2014). The following problems may exclude potential framework users: 1) the complexity and enormity of the existing frameworks, which make it difficult to understand the interfaces and structures of frameworks and their method of operation; 2) simple application examples, often involving decomposition and recombination of a simple conceptual model, or standardisation and coupling of two or more models, which fall short of providing complete information on the frameworks and the experience of actual use; and 3) the lack of available framework-based model components. Thus, in addition to improving framework performances, more application examples are required to boost potential users. The objectives of this paper are to present a complete example of implementation and application of a distributed hydrological model based on a componentbased modelling framework and to contribute it to the framework component base. The major challenges of this work lie in understanding the mechanism of hydrological processes and decomposing the hydrologic system into components at the optimum level of granularity. The framework adopted for model implementation and simulation in this study is known as the Heihe river basin Open Modelling Environment (HOME). The Cold and Arid Regions Environmental and Engineering Research Institute of the Chinese Academy of Sciences, with the support of the National Natural Science Foundation of China, has been developing the HOME framework since 2011 to serve as an integrated platform that supports model development, coupling, management, simulation and analysis for ‘Integrated Research on the Eco-hydrological Processes of the Heihe Basin’. HOME uses object-oriented techniques based on JAVA to enable modellers to customise their models by coupling individual components that express different processes. It comprises four parts: a kernel system, a graphic modelling tool, a component base and a database. So far, the kernel system that supports component development, component coupling, model configuration and execution has been realised and can run on any operating system without the other three parts. Our work was based on the kernel system of HOME. We focus on the runoff simulation of the Buha River basin, which is the largest tributary of the largest inland salt-water lake in China, Qinghai Lake, located in the north-eastern Tibetan Plateau. The water level of the lake was reported to have decreased by about 3.7 m overall from 1959 to 2004, but to have increased by nearly

1 m from 2004 to 2009 (Zhang et al., 2014). Modelling the discharge of the Buha River, which accounts for almost half of the lake's inflow, is hence important to assess the water resources of the lake. However, few data are available for this alpine basin, which covers a large area of approximately 15,000 km2, with one hydrometric station at the outlet and seven meteorological stations in and around the area, making it quite a challenge to simulate the hydrological processes. Thus, based on the HOME framework, the GeomorphologyBased Hydrological Model (GBHM) initially developed by Dawen Yang (Yang et al., 2001b) in Fortran, was implemented. This distributed hydrological model (Xu et al., 2008; Xu et al., 2013; Yang et al., 2001a) can be applied to large, complex basins because it adopts the flow interval-hillslope discretisation scheme, is able to represent spatial variability such as topography, land cover and soil type and incorporates vegetation dynamics using the temporal leaf area index (LAI). In particular, this physically-based distributed hydrological model is quite useful in data-poor environments because it implements very few empirical parameters and does not quite rely on the observed data for parameter calibration. Nevertheless, the original code was not well modularised and was difficult to reuse, as the input and configuration parts needed to be modified for different study areas. Additionally, it did not consider the effects of the soil freezing-thawing process, which significantly affects the hydrological processes in the Buha River basin. Consequently, the original GBHM was componentised and a simple empirical soil freezing-thawing component was incorporated. In the following sections, the HOME framework is presented and the reconstruction of the distributed hydrological model based on the HOME is described in detail. The hydrological simulation of the Buha River basin based on the new model and the analysis of the results are then specified. Finally, we summarise the findings, present our conclusions and outline future work. 2. Heihe river basin Open Modelling Environment (HOME) HOME is an integrated framework developed specifically for eco-hydrological model development. Its aims are to 1) enhance the development and evaluation of scientific components; 2) facilitate model development and improvement through the use of well-componentised legacy codes; 3) provide standard interfaces for model importing and coupling for solving complex problems; and 4) provide a wide range of analysis components. To achieve these aims, HOME is implemented in four parts: a graphic modelling tool, a database, a component base and a kernel system (Fig. 1). The graphic modelling tool is mainly used to assist the model configuration, control the model simulation and visualise the modelling outputs. The database and component base provide modelling resources. The physical process components and data processing and analysis tools are managed by the component base. The kernel system, which as previously mentioned is able to run independently, requires at least three parts to work: a configuration file, a model entity and a runtime. The model comprises a set of components, contexts and data pools, which are the key concepts used by HOME. As only the kernel system was used for the model implementation and simulation, only the key elements of the kernel system are introduced in detail in this section. 2.1. Key HOME entities 2.1.1. Component The term ‘component’ in the HOME refers to a modelling entity that implements one of the eco-hydrological processes or a data processing method. Each component extends an abstract class and overrides its three methods, init, run and clear, to represent the

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

247

The context maps the variables of the components to the variables in the data pool. Thus, the incoming values of the components are obtained from the data pool using the getValue method of the reference variables, and the resulting values are assigned to the data pool using the setValue method. To scope the variables, HOME includes the following rules: a context can have only one data pool; components in the parent context cannot access variables in the data pool of a child context, while components in a child context can access the parent's data pool; and a component cannot access variables in other data pools at the same level. 2.2. Configuration file and HOME runtime

Fig. 1. The principle architecture of HOME.

pre-run, runtime and post-run states of the model simulation, as in most component-based systems. The init method can be used for allocating memory for arrays, loading the model parameters and initial conditions from the input files and opening the output files if the resulting values need to be written in real time. Mathematical expressions for each eco-hydrological process are usually coded in the run method. The clear method is usually used to close files opened by the init method. The model developers are responsible for implementing these three methods. The variables in the component that serve as communication media must be defined as ‘reference’ types so that they can pass on values. The HOME implements the basic data types as reference types, and a referencetype variable uses the getValue and setValue methods to get and assign its value. 2.1.2. Context In HOME, the context is a routine in charge of calling the methods of the components and managing their lifetimes. We can consider the context as a container that contains reference-type variables, a set of components or even other contexts, which means that contexts can be nested. A context is essentially a special component that is predefined. The init, run and clear methods of the context are responsible for invoking the init, run and clear methods of the components, respectively. Two types of context, a temporal context and a spatial context, are specially built to control the iterations of the components. The temporal context uses attributes including the start time, end time, time interval and unit of time interval to control the iterations. The spatial context, as the name implies, iterates on each grid of a 2D map that depicts the spatial extent of the study area. In addition, the init and clear methods are invoked once during the whole simulation, while the run method is invoked during each iteration of the context. The model in HOME is also implemented as a context, containing a nested set of contexts, components and variables.

A configuration file with a specific XML scheme is designed for modellers to customise their models. Based on this file, modellers are able to decide which contexts and components will be included, define the mapping relationships between the component variables and the data pool variables and assign initial values to the variables (see Table 1). The runtime is HOME's scheduler, which handles the model simulation from reading the configuration file to loading, initialising, running and finalising the model. At the start of the simulation, the configuration file is parsed by the runtime. According to the information provided by the configuration file, instances of contexts and components are generated by the Java reflection mechanism. Then, the variables and these instances are respectively added to the data pool and the component list of the contexts that they belong to. After completing the model loading, the init, run and clear methods are invoked by the runtime. The components then run in sequence in accordance with the control structure of the contexts and the order of the components in the XML file. 3. Model description and implementation 3.1. Model introduction The geomorphology-based hydrological model (GBHM) is a distributed hydrological simulator. It operates on hourly time steps at the basin scale and uses the hillslope-flow interval discretisation scheme to represent the catchment. In the GBHM, the drainage network is extracted from the Digital Elevation Model (DEM) and the basin is delineated by a set of subcatchments coded according to the Pfafstetter numbering system (D Yang et al., 2000). In each subcatchment, the water flow pathway is simplified as a single main channel with the maximum flow distance of the subcatchment (Yang et al., 2001b). Along the flow direction, each subcatchment is subdivided into a set of flow intervals, the area between two flow distance contours. The main channel is correspondingly subdivided into stream segments. Each grid in a flow interval is abstracted as a rectangular inclined plane, called the hillslope element, which is symmetrically distributed on each side of the stream segment with a slope of b and a length of l, calculated P P as l ¼ a=2 L, where a is the area of the grid and L is the total length of the stream network in the grid (see Fig. 2) (Yang et al., 2004). The hillslope element is the basin's basic modelling unit, in which hydrological processes, including canopy interception, snowmelt, evapotranspiration, infiltration, layer interflow, subsurface flow, the exchange between soil water and ground water, and sheet flow, are calculated (Yang, 1998). 1) Interception

2.1.3. Data pool A common storage space, referred to as the data pool, is used to pass values between the components. The data pool, which stores the values of the reference variables, is implemented in the context.

Actual interception depends on precipitation and the amount and deficit of the canopy water storage. The deficit of the canopy water storage SCd at time t is the difference between interception

248

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

Table 1 Key elements and their expressions in the XML configuration file. Element

Format

Roles

model context component variable



Assigning the name of the model Assigning the class name and object name of the context Assigning the class name and object name of the component Defining the data type, size and initial value of the component variable and mapping the relationship between the data pool and component variables, where the ‘attribute’ element indicates the name of the data pool variable, the ‘context’ element indicates the data pool, and the ‘name’ element indicates the name of the component variable

Fig. 2. Hillslope-flow interval discretisation scheme of the GBHM.

capacity SC0 and canopy water storage SC, given by

SCd ðtÞ ¼ SC0 ðtÞ  SC ðtÞ

 Ecanopy ðtÞ ¼ (1)

The interception capacity SC0 is a function of the leaf area indext(LAI), calculated as

SC0 ðtÞ ¼ 0:2LAIðtÞ

(2)

Kv Kc Ep SC ðtÞ=Dt

M ¼ Sf ðTi  Tb Þ

(3)

where M is the depth of the melting water (mm); Sf is the snow melting factor; Ti is the mean air temperature ( C); and Tb is the base temperature when the snow starts to melt ( C).

Etr ðtÞ ¼ Kv Kc Ep f1 ðzÞf2 ðqÞ

LAIðtÞ LAImax

(5)

where Etr is the transpiration rate (mm/hr); f1(z) is the root distribution function, treated as the inverted triangular distribution with the base located at the soil surface; z is the average depth of the soil layer; f2(q) is the soil moisture function, treated as the linear decrease from the field capacity to the wilting point; q is the soil moisture content; and LAImax is the maximum leaf area index of the year. The evaporation rate from the surface water storage in time interval [t, tþDt] is calculated as

3) Evapotranspiration Actual evapotranspiration is assumed to take place only during the daytime 12 h. It is calculated as the sum of the evaporation from the canopy water storage, transpiration from the root zone, evaporation from the surface water storage and evaporation from the soil surface. The actual evaporation rate from the canopy water storage is given by

(4)

where Ecanopy is the evaporation rate from the canopy (mm/hr); Kv is the fraction of the vegetation cover; Kc is the crop coefficient, defined as the ratio of the crop evapotranspiration under standard conditions to the reference evapotranspiration (Allen et al., 1998); Ep is the potential evapotranspiration rate (mm/hr); and Dt is the time interval (hr). The vegetation transpiration rate from the root zone in time interval [t, tþDt] is calculated as

2) Snowmelt Snowmelt is estimated using the temperature index approach, expressed as

SC ðtÞ  Kv Kc Ep Dt SC ðtÞ < Kv Kc Ep Dt

 Esurface ðtÞ ¼

ð1  Kv ÞEp Ss ðtÞ=Dt

Ss ðtÞ  ð1  Kv ÞEp Dt Ss ðtÞ < ð1  Kv ÞEp Dt

(6)

where Esurface is the evaporation rate from the surface water storage (mm/hr) and Ss(t) is the surface water storage (mm) at time t. The evaporation rate from the soil surface in time interval [t, tþDt] is given by

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

6) Exchange between saturated zone and river

h i Es ðtÞ ¼ ð1  Kv ÞEp  Esurface ðtÞ f2 ðqÞ

(7)

where Es is the evaporation rate from the soil surface (mm/hr). 4) Unsaturated zone water flow The unsaturated zone is vertically divided into a maximum of 10 layers on the hillslope, and its depth changes when the level of the groundwater table increases. Vertical water movement in the unsaturated zone is described using the one dimensional Richards equation, given as

vqðz; tÞ vqv ¼ þ sðz; tÞ vt vz

(8)

where s(z, t) is the soil moisture source and sink, usually a negative value capturing the evaporation from the soil surface and transpiration from the root zone. qv refers to the vertical soil moisture fluxes, calculated according to Darcy's Law:

qv ¼ Kðq; zÞ

  vJðqÞ 1 vz

(9)

where K(q, z) is the soil hydraulic conductivity and J(q) is the soil suction. The relationship between soil suctionJ(q) and water content q is expressed by the Van-Genuchten model:

8  > > > < Se ¼

1 1 þ ðaJÞn

m

> > q  qr > : Se ¼ qs  qr

(10)

where qr is the residual soil moisture; qs is the saturated soil moisture; a and n are constants determined by experiments or field measurements; and m ¼ 1/n. The soil hydraulic conductivity is calculated as 1=2

Kðq; zÞ ¼ Ks Se

h   i2 1=m m 1  1  Se

(11)

where Ks is the saturated hydraulic conductivity. The Richards equation is solved using an implicit finite difference scheme. The surface runoff, which can be interpreted as Dunne saturated runoff and Horton runoff, is also generated by solving the Richards equation. 5) Overland flow After filling the surface depressions, the surface runoff flows through the hillslope into the stream. The overland flow along the hillslope is described using the one dimensional kinematic wave model:

8 vh vqs > > < vt þ vx ¼ i > > : qs ¼ 1 S1=2 h5=3 ns 0

249

(12)

where qs is the discharge per unit width (m3/(s$m)); h is the water depth after subtracting surface water detention (m); t is time (s); x is the length of the hillslope element; i is the effective precipitation (m); S0 is the gradient of the hillslope element; and ns is Manning's roughness parameter.

The mass balance and Darcy's law are used to estimate the exchange rate between the saturated zone and river, given as

8 vS ðtÞ 1000 > > > G ¼ rechðtÞ  LðtÞ  qG ðtÞ < vt A > H  H h þ h > 1 2 1 2 > : qG ðtÞ ¼ KG l=2 2

(13)

where qG(t) is the discharge rate to the river per unit width (m3/ (hr$m)); vSG(t)/vt is the change of water storage in the unconfined aquifer (mm/hr); rech(t) is the recharge rate from the upper unsaturated zone (mm/hr), also calculated using Darcy's law; L(t) is the leakage to the deep aquifer (mm/hr); A is the area of the hillslope element per unit width (m2/m); KG is the hydraulic conductivity of the unconfined aquifer (m/hr); l is the length of the hillslope element (m); H1 and H2 are the groundwater tables at the last and current time steps (m), respectively; and h1 and h2 are the water levels of the river at the last and current time steps (m), respectively. 7) River flow In the hillslope-flow interval system, the sheet flow and groundwater flow from all of the hillslope elements in a flow interval are assumed to drain directly into the main river. The river flow is essentially described by the one-dimensional kinematic wave equations used for overland flow:

8 vA vQ > > þ ¼q > < vt vx 1=2 > S > > : Q ¼ 0 2=3 A5=3 nr p

(14)

where q is the lateral inflow (m3/(s$m)), calculated as q ¼ qs þ qG; x is the distance along the longitudinal axis of the river (m); A is the cross-section area of the river (m2); Q is the discharge at x (m3/s); S0 is the slope of the river bed; nr is Manning's roughness parameter; and p is the wetting perimeter (m).

3.2. Components of the model Based on the HOME framework, a modular version of the geomorphology based hydrological model was constructed. The components of the distributed hydrological model can be divided into the core unit components and the auxiliary components. Nine unit components are first developed to perform the calculations for canopy interception, snow melt, actual evaporation (including evaporation from canopy storage, transpiration from crops and evaporation from the soil surface), the recharge rate from the unsaturated zone to groundwater, infiltration, layer interflow, subsurface flow (including groundwater and river exchange and the flow from the top saturated zone above the groundwater), sheet flow and channel flow. A unit component includes neither spatial nor temporal iterations to enhance its reusability; i.e. where and when the behaviour of the component is not known. In particular, a simple empirical soil freezing-thawing component, which accounts for the hydrological effects of soil freezing and thawing by calibrating the soil hydraulic conductivity according to the daily mean air temperature, was implemented for applications in high-elevation regions. The calibrated hydraulic conductivity is calculated as

250

8 0 > >   > < 10  T k ¼ k0 $exp  f 15 þ T > > > : k0

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

T  15 15 < T < 10

(15)

T  10

Where k is the modified hydraulic conductivity(mm/hr), k0 is the initial hydraulic conductivity(mm/hr), T is the air temperature( C) and f is a constant parameter. This component was only applied from mid-November to the end of May, according to the daily flow hydrograph. Auxiliary components used for data input and output and data processing were redefined, including the base maps (e.g., elevation, land-use type, slope and slope length) importing components, climate data reading component, reading and setting components for parameters (i.e., soil, vegetation and other model parameters), reading and setting components for initial conditions (e.g., soil moisture, river water level and groundwater depth), outputting components for modelling results (e.g., discharge, evapotranspiration and soil moisture), and other data processing components for spatial interpolation and statistical tasks. In addtion, to facilitate the spatial and temporal display of the modelling results, such as a spatiotemporal variation of the evapotranspiration and soil moisture in the basin, most of the results were exported using the ESRI ASCII raster format at a daily time step. All of the GBHM components were compliant with the HOME's component standard. The components implemented the init, run and cleanup methods. Their exchanged items were defined as pool data types, and the getValue and setValue methods were used to input and output values. To avoid overuse of the getValue and setValue methods of the ‘reference’ data types, the algorithms in the HOME-based component were implemented first through native language data types. The getValue methods of the input variables were then called at the front part of the execution code, and the setValue methods of the output variables were called at the end part. For example, the temperature index-based snowmelt component was implemented as shown in Table 2. Annotations of the component and variables are not mandatory, as they simply play their parts in the graphic modelling tool of HOME. Thus far, HOME has not provided any functionality for model code implementation. Hence, the Netbeans integrated development environment was used for component code development, which also supported the development of HOME. To be identified by HOME's kernel system, components of the model have to be manually packaged into a JAR file and added to the HOME library. 3.3. Model configuration A simplified XML configuration file with component variables omitted defines the structure of the distributed model (see Table 3), which is also depicted in Fig. 3. The first part is an input context, where the components for reading the base maps, parameters and initial conditions are defined. These components run once during the whole simulation and hence all of the algorithms are realised in the init method. The second part is a daily temporal context, in which the time-dependent input components that perform daily climate data readings, spatial interpolations and 4-day LAI readings are defined, a soil freezing-thawing component for calibrating the overall hydraulic conductivities of the basin is added, an hourly temporal context for simulating the core of the model is nested and daily output components are incorporated. The hourly temporal context contains a temporal downscaling component for calculating hourly values from daily weather data, as the model simulates on an hourly time scale; a spatial context for iterating each

hillslope element of the study area; a routing component for accumulating flow through river networks; and a statistical component for modelling the results. In the spatial context, hydrological processes including interception, snowmelt, evapotranspiration, percolation, infiltration, subsurface flow and sheet flow are simulated on each hillslope element. All of the time- and space-dependent algorithms are implemented in the run method. The execution sequence of the components and contexts is indicated by the arrows in Fig. 3. 4. Hydrological modelling 4.1. Study area and data availability The origin of the Buha River (Fig. 4) is in Shulenan Mountain. The river is 286 km long and is the largest tributary of Qinghai Lake, flowing into the lake from the northwest to the southeast, with annual mean runoff of 799.2 million m3. The precipitation period is between May and September, whereas runoff is concentrated from June to October, lagging a month behind precipitation. There is hardly any precipitation in winter. The mean air temperature is above 0  C from May to September, with the maximum in July and August. The basin's elevation varies from 3188 to 5279 m above sea level, with an average of 3979 m. The major types of land cover are alpine meadow (55.73% of the study area), bare soil, everglade, alpine desert and tundra. The primary soil types include felty, thin dark felty, castanozems, peaty bog, dark frigid calcic, calcic litho, humid felty and bog soils. The distributed hydrological model requires data on the DEM, land use, soil types, LAI and climate. The DEM was obtained from the SRTM database (http://srtm.csi.cgiar.org/) with a spatial resolution of 90 m, from which the slope, aspect and length of hillslope were derived. The land use data were from the 1:100,000 land use map of China (Liu et al., 2003), downloaded from the Cold and Arid Regions Science Data Centre in Lanzhou (http://westdc.westgis.ac. cn/), and reclassified as nine land cover types, as shown in Table 4. The soil map was taken from the 1:1,000,000 soil database of China (Shi et al., 2002), whereas the soil properties data were taken from the China dataset of soil hydraulic parameters (Dai et al., 2013). The LAI was taken from the MCD15A3 product (https://lpdaa c.usgs.gov/products/modis_products_table/mcd15a3), i.e. the level-4 combined (Terra and Aqua) MODIS global LAI and FPAR product, which is composited every 4 days at a spatial resolution of 1 km. The climate data, including daily precipitation, air temperature, relative humidity, wind speed and sunshine duration, were downloaded from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn/). The daily observed stream flow in the Buha River estuary, which was used as the calibrated data, was provided by Hydrology and Water Resources Survey Bureau of Qinghai Province. To ensure relatively high computing efficiency while retaining adequate geographical information from the data sources and maintaining sufficient accuracy in extracting the stream network from the DEM, the gridded system on which the model was based was set to a resolution of 900 m. Thus, all datasets were processed to a 900-m resolution using Albers projection with the GCS_Krasovsky_1940 geographic coordinate system. 4.2. Model parameters There are four parameter categories in the model: 1) soil parameters, including saturated and residual soil moisture, saturated hydraulic conductivity and parameters of the Van-Genuchten model; 2) vegetation and underlying surface parameters, such as the LAI, fraction of vegetation coverage, crop coefficient, root depth,

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258 Table 2 Implementation of the snowmelt component based on the HOME framework.

251

252 Table 3 Simplified XML configuration file for GBHM.

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

253

Fig. 3. Framework of the distributed hydrological model in the HOME.

maximum surface water detention and Manning's roughness of the hillslope; 3) river parameters, including Manning's roughness and shape parameters, such as length, width, depth and bed slope of the river; and 4) empirical parameters such as the snow melting factor, snowmelt base temperature and the storage coefficient and hydraulic conductivity of the groundwater. Section 4.1 details how the soil parameters and LAI were accessed. In this study, the root depth, maximum surface water detention and Manning's roughness of the hillslope were estimated according to land use types, the shape parameters of the river were estimated from measured data and DEM, Manning's roughness of the stream segments was set between 0.01 and 0.014, and the empirical parameters were estimated via model calibration. In particular, actual evapotranspiration is a major term of water

budgets and must be estimated as accurately as possible. In this model, using the fraction of green vegetation cover (Kv), a hillslope element is split into a transpirative surface, on which evaporation from the canopy water storage and transpiration from the root zone are calculated, and an evaporative surface, on which evaporation from the surface water storage and evaporation from the soil surface are estimated. The crop coefficient (Kc) is combined with potential evapotranspiration for the prediction of canopy evaporation and transpiration. The crop coefficient and fraction of vegetation cover are hence key parameters for estimating actual evapotranspiration. In this study, the MODIS Normalized Difference Vegetation Index (NDVI) for the period 2003e2010, taken from the MOD13Q1 product (https://lpdaac.usgs.gov/products/modis_products_table/

254

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

Fig. 4. Locations of the study area and the hydrological and meteorological stations.

Table 4 Crop coefficient for each type of land cover. Land cover

Kc

Land cover

Kc

Water Urban and built-up Bare soil Forest Croplands

1.05 0.2 0.2 1.0 1.1

Grass Shrub Wetlands Snow and ice

1.05 1.0 1.2 0

mod13q1) every 16 days at a spatial resolution of 250 m, was used to estimate the Kv per grid cell with a simple linear mixing equation (Gutman and Ignatov, 1997):

Kv ¼ ðNDVI  NDVI0 Þ=ðNDVI∞  NDVI0 Þ

(16)

where NDVI0 and NDVI∞ are the minimum and maximum NDVI values. These values were assumed to be independent of the vegetation/soil type and were calculated for each grid cell as the lower and upper 5% of the NDVI histogram for the period of 2003e2010. The 250-m NDVI images were initially projected to Albers projection with 900-m resolution and then used to derive the Kv values, which were updated every 16 days during the model simulation. Estimation of the seasonal distribution of crop coefficient is a difficult task, especially when monitoring large areas dominated by natural landscapes, as the conventional crop coefficient is specific

to the crop planted in each field and its calendar (Maselli et al., 2014). Thus, the crop coefficient was assumed to be constant for each type of land cover. The Kc values were referenced from the Food and Agricultural Organization of the United Nations Irrigation and Drainage Paper No. 56 values for full-grown crops (Allen et al., 1998), except for urban and built-up land, bare soils and snow and ice. The Kc for bare soils and urban and built-up land was referenced from Maselli et al. (2014). For snow and ice, it was assumed that there was no transpiration and the Kc was reduced to 0. The Kc values are summarised in Table 4. 4.3. Model calibration and validation A threshold area of 2000 ha was adopted for watershed delineation, and the Buhakou hydrometric station was set as the outlet of the whole basin to enable comparison of the simulated and observed stream flow. Consequently, 161 subcatchments were identified. Stream segments in each subcatchment were abstracted according to flow distance with a contour interval of 2000 m. To obtain the initial state variables, a warm up run was carried out by cycling the climate forcing over the period 2005e2009, until the basin's hydrological state of equilibrium was achieved. An additional spin up run was carried out from July to December 2002, as the MCD15A3 product was not available before that time. The discharge observed at the Buhahekou hydrometric station was used in the calibration and validation of the model. The period 2003e2004 was used as the calibration period, and period

Table 5 Descriptions and values of the calibrated parameters in GBHM. Parameters

Description

Method

Value

Sf Tb Ks KG GWcs

Snow melting factor Base temperature at which snowmelt starts ( C) Saturated hydraulic conductivity (mm/hr) Hydraulic conductivity of the groundwater (mm/hr) Groundwater storage coefficient

Basin Basin Hillslope element Hillslope element Basin

0.1 5.0 70% 2.5%Ksat 0.1

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

255

Fig. 5. Simulated and observed discharge at Buhahekou hydrometric station during the whole simulation period.

Table 6 Model performance for the calibration and validation period based on daily and monthly discharge at Buhahekou hydrometric station. Period

Ens-daily

R2-daily

RE-daily

Ens-monthly

R2-monthly

RE-monthly

Calibration: 2003e2004 Validation: 2005e2010

0.75 0.56

0.79 0.70

8.15% 5.67%

0.95 0.75

0.96 0.80

8.15% 5.67%

Table 7 Mean annual water balance of the whole study area during the whole simulation period. Precipitation (mm)

Evapotranspiration (mm)

Runoff (mm)

Change in water storage (mm)

378.0

302.8

65.7

9.5

2005e2010 was used as the validation period. Three metrics were used to evaluate the stream flow modelling performance (Legates and McCabe, 1999): the Nash-Sutcliffe efficiency (Ens), the coefficient of determination (R2) and the relative error (RE). Although HOME provides automated calibration components, the model was manually calibrated because running the model was time consuming and only the empirical parameters and some sensitive soil parameters required calibration. Table 5 lists the calibrated parameters and their values. In particular, the saturated hydraulic conductivity of each soil layer was calibrated to 70% of the initial value from the soil dataset. The hydraulic conductivity of the groundwater was set to 2.5% of the saturated hydraulic conductivity of the last soil layer. The preceding two parameters were calibrated based on the hillslope element, and the other three were set to constants for the whole basin. Fig. 5 compares the simulated and observed discharge at Buhahekou hydrometric station. The model shows good performance in simulating the daily and monthly discharge. Table 6 lists the values of the three metrics. For monthly discharge, the Ens values were greater than 0.7, the R2 values were greater than 0.8

and the absolute values of RE were less than 10% for both the calibration and validation periods. For daily discharge, the RE values showed no change. The Ens and R2 values decreased slightly but still within an acceptable range. Fig. 5 also shows that the daily discharge was highly overestimated during some summer days. The overestimation may be due to the extremely poor density of the rain gauges and errors carried over from the spatial interpolation processes, whereas the Buha River water supply is mainly from rainfall. Uncertainties in the physical parameters of the soil, such as the depth of the soil profiles, saturated water content and saturated hydraulic conductivity, which are very sensitive parameters in the model, could also explain this overestimation. Nevertheless, the modelling results were within acceptable limits, indicating that the model had the capacity to simulate runoff generation appropriately. 4.4. Water balance characteristics The annual water balance of the basin was analysed using the simulated results from 2003 to 2010 (with monthly Ens ¼ 0.77,

256

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

Fig. 6. Spatial distributions of the mean annual values from 2003 to 2010: (a) evapotranspiration and (b) runoff.

Table 8 Comparison of different criteria produced by model simulations with and without the soil freezing-thawing component over the total period.

Ens-daily R2-daily Ens-monthly R2-monthly RE Evapotranspiration (mm) Total runoff (mm) Soil moisture (mm)

Without soil freezing-thawing

With soil freezing-thawing

0.63 0.64 0.73 0.74 7.24% 316.7 58.9 865.2

0.58 0.71 0.77 0.82 2.93% 302.8 65.7 878.9

R2 ¼ 0.82 and RE ¼ 2.93%). Table 7 lists mean annual values of the water balance terms, from which it could be observed that evapotranspiration was the largest output item, accounting for around 80% of precipitation. The runoff coefficient was 0.17, indicating that the basin was dominated by a semi-arid climate. The water storage changes in the basin, including changes in the soil water, groundwater, snowpack and surface water, was a very small fraction of the precipitation. As such, the changes were always neglected in the mean annual water balance analysis. Thus, the mean annual water balance equation was typically simplified so that the runoff was the difference between precipitation and actual evapotranspiration. Hence, it was crucial to accurately estimate precipitation and evapotranspiration for the runoff simulation. The model simulated actual evapotranspiration in three parts as mentioned before: evaporation from canopy storage, and transpiration from crop and soil evaporation, which were estimated, respectively, as 4% (12.9 mm), 39% (117.4 mm) and 57% (172.5 mm) of the total annual evapotranspiration. These figures suggest that evapotranspiration was mainly determined by soil evaporation and crop transpiration. Fig. 6 shows high spatial variability in the water balance components. The mean annual actual evapotranspiration ranged from 131.7 to 419.2 mm, increasing from upstream to downstream (see Fig. 6(a)). The annual runoff generation ranged from 0 to 246.2 mm, which has the opposite spatial tendency to that of evapotranspiration, i.e. decreasing from upstream to downstream (see Fig. 6(b)). The major source areas of runoff were located in the northeast of the basin. This spatial variability was a good reflection of the vegetation distribution and topography, impling that the model

had the capacity to simulate the effects of topography and land cover on the hydrological processes. 4.5. Effects of frozen soil on runoff simulation To investigate the effects of the soil freezing-thawing processes on hydrological processes, an additional simulation with the soil freezing-thawing component eliminated was conducted and compared to the former simulation. Table 8 provides an overview of the 8-year (2003e2010) statistical results from the two simulations. The runoff simulation was clearly improved when the effects of frozen soil were considered. The model with the soil freezing-thawing component generated more runoff, evaporated less and was able to store more moisture. Fig. 7 shows the details of the 8-year monthly mean results. In winter, the total runoff and evapotranspiration decreased and more water was conserved in the soil profile compared with the scenario without soil freezing-thawing, as the soil evaporation, infiltration, soil water movement and exchange rate between the unsaturated soil layer and unconfined aquifer were restricted due to the reduced hydraulic conductivity. In summer, the total runoff simulated in the experiment with soil freezing-thawing was much higher due to the higher soil moisture level, but with almost no difference in evapotranspiration. The runoff difference between the two experiments became smaller in later periods as the soil moisture levels evened out. Overall, the incorporation of the soil freezing-thawing component altered the magnitude and temporal scheme of the runoff.

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

257

2) The HOME-based distributed hydrological model performed well in the monthly runoff simulation despite a lack of adequate data, thus demonstrating the robustness of the distributed modelling scheme and the helpfulness of the remote sensing based parameterisation scheme. 3) Including the empirical soil freezing-thawing component improved the performance of the runoff simulation, which implies that the soil freezing-thawing process has a significant effect on runoff simulation, especially during the cold season. This comparison also demonstrated the flexibility of the component-based programming scheme for model modification. Overall, we have produced an efficient way of developing models and understanding hydrological processes in high altitude areas with scarce data availability. However, this component-based approach involves a significant execution time overhead compared with the initial Fortran-coded GBHM. The time spent by running the Fortran-based GBHM over period 2003e2010 was 8638 s (about 2.4 h). In contrast, the HOME-based GBHM spent 18,472 s (nearly 5.1 h), more than twice the time of the former. Sources of simulation inefficiency may include 1) the HOME's data exchange scheme, which may have a considerable effect on the simulation time; 2) the fine-grained components of the model, which may increase run-time data exchange; 3) the computational distinction between JAVA and FORTRAN; or 4) the greater number of writing tasks than seen in the original Fortran version. Further research is urgently needed to identify the specific causes of the computational overhead and improve simulation efficiency. Acknowledgements This work was supported by the National Natural Science Foundation of China under grant No. 91125005, J1210003/J0109 and 41301508 and the Research cloud of alpine joint observation of Academy of Chinese Academy of Sciences (grant No. XXH12503-0507).We give special thanks to the development team of the Heihe river basin Open Modelling Environment (HOME). Fig. 7. The average monthly basin quantities produced by the two simulations.

References From the perspective of the yearly water balance, the main difference resulted from greater evapotranspiration in the cold season when frozen soil was not considered.

5. Summary, discussion and conclusions Component-based programming provides a generic approach for constructing hydrological models, which reduces the amount of repetitive coding and speeds up tasks such as improving the model and comparing and assessing the components. This paper outlines HOME, a new component-based modelling framework that was used as the basis for the implementation of a distributed hydrological model. A runoff simulation of a large-scale alpine basin with scarce data was carried out to test the componentised hydrological model. The effects of the soil freezing-thawing process on runoff generation were assessed through the inclusion and exclusion of a simple empirical soil freezing-thawing component. The major conclusions drawn from this study are as follows. 1) The HOME-based component development and model construction are proven to be feasible, which will boost the confidence of potential users.

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements. Irrigation and drainage paper NO. 56, 300 pp., FAO, Rome, Italy. Argent, R.M., 2004. An overview of model integration for environmental application - components, frameworks and semantics. Environ. Modell. Softw. 19 (3), 219e234. Argent, R.M., Voinov, A., Maxwell, T., Cuddy, S.M., Rahman, J.M., Seaton, S., Vertessy, R.A., Braddock, R.D., 2006. Comparing modelling frameworks e a workshop approach. Environ. Modell. Softw. 21 (7), 895e910. Buahin, C.A., Horsburgh, J.S., 2015. Evaluating the simulation times and mass balance errors of component-based models: an application of OpenMI 2.0 to an urban stormwater system. Environ. Modell. Softw. 72, 92e109. Castronova, A.M., Goodall, J.L., 2010. A generic approach for developing processlevel hydrologic modeling components. Environ. Modell. Softw. 25 (7), 819e825. Castronova, A.M., Goodall, J.L., 2013. Simulating watersheds using loosely integrated model components: evaluation of computational scaling using OpenMI. Environ. Modell. Softw. 39, 304e313. Dai, Y., Wei, S., Duan, Q., Liu, B., Fu, S., Niu, G., 2013. Development of a China dataset of soil hydraulic parameters using pedotransfer functions for land surface modeling. J. Hydrometeorol. 14 (3), 869e887. David, O., Ascough, J.C., Lloyd, W., Green, T.R., Rojas, K.W., Leavesley, G.H., Ahuja, L.R., 2013. A software engineering perspective on environmental modeling framework design: the object modeling system. Environ. Modell. Softw. 39, 201e213. Dechmi, F., Burguete, J., Skhiri, A., 2012. SWAT application in intensive irrigation systems: model modification, calibration and validation. J. Hydrol. 470, 227e238. Formetta, G., Antonello, A., Franceschi, S., David, O., Rigon, R., 2014. Hydrological modelling with components: a GIS-based open-source framework. Environ. Modell. Softw. 55, 190e200.

258

Y. Long et al. / Environmental Modelling & Software 80 (2016) 245e258

Gutman, G., Ignatov, A., 1997. Satellite-derived green vegetation fraction for the use in numerical weather prediction models. Adv. Space Res Ser. 19 (3), 477e480. Hill, C., DeLuca, C., Suarez, Balaji, M., da Silva, A., 2004. The architecture of the earth system modeling framework. Comput. Sci. Eng. 6 (1), 18e28. Hu, C.H., Guo, S.L., Xiong, L.H., Peng, D.Z., 2005. A modified Xinanjiang model and its application in northern China. Nord. Hydrol. 36 (2), 175e192. Khakbaz, B., Imam, B., Sorooshian, S., Koren, V.I., Cui, Z.T., Smith, M.B., Restrepo, P., 2011. Modification of the national weather Service distributed hydrologic model for subsurface water exchanges between grids. Water Resour. Res. 47. Leavesley, G.H., Markstrom, S.L., Brewer, M.S., Viger, R.J., 1996. The modular modeling system (MMS) - the physical process modeling component of a database-centered decision support system for water and power management. Water Air Soil Poll. 90 (1e2), 303e311. Legates, D.R., McCabe, G.J., 1999. Evaluating the use of “goodness-of-fit” measures in hydrologic and hydroclimatic model validation. Water Resour. Res. 35 (1), 233e241. Liu, J., Liu, M., Zhuang, D., Zhang, Z., Deng, X., 2003. Study on spatial pattern of landuse change in China during 1995e2000. Sci. China Ser. D Earth Sci. 46 (4), 373e384. Lloyd, W., David, O., Ascough, J., Rojas, K.W., Carlson, J.R., Leavesley, G., Krause, P., Green, T.R., Ahuja, L., 2011. Environmental modeling framework invasiveness: analysis and implications. Environ. Modell. Softw. 26 (10), 1240e1250. Maselli, F., Papale, D., Chiesi, M., Matteucci, G., Angeli, L., Raschi, A., Seufert, G., 2014. Operational monitoring of daily evapotranspiration by the combination of MODIS NDVI and ground meteorological data: application and evaluation in Central Italy. Remote Sens. Environ. 152 (0), 279e290. Moore, R.V., Tindall, C.I., 2005. An overview of the open modelling interface and environment (the OpenMI). Environ. Sci. Policy 8 (3), 279e286. Ouessar, M., Bruggeman, A., Abdelli, F., Mohtar, R.H., Gabriels, D., Cornelis, W.M., 2009. Modelling water-harvesting systems in the arid south of Tunisia using SWAT. Hydrol. Earth Syst. S. C. 13 (10), 2003e2021. Peckham, S.D., Hutton, E.W.H., Norris, B., 2013. A component-based approach to integrated modeling in the geosciences: the design of CSDMS. Comput. GeosciUk 53, 3e12. Pomeroy, J.W., Gray, D.M., Brown, T., Hedstrom, N.R., Quinton, W.L., Granger, R.J., Carey, S.K., 2007. The cold regions hydrological model: a platform for basing process representation and model structure on physical evidence. Hydrol.

Process 21 (19), 2650e2667. Reyes, M.R., Bengtson, R.L., Fouss, J.L., Rogers, J.S., 1993. Gleams Hydrology submodel modified for shallow-water table conditions. T Asae 36 (6), 1771e1778. Shi, X., Yu, D., Pan, X., Sun, W., Gong, Z., Warner, E., Petersen, G., 2002. A Framework for the 1: 1,000,000 Soil Database of China. paper presented at the 17th world congress of soil science, Bangkok. Whelan, G., Kim, K., Pelton, M.A., Castleton, K.J., Laniak, G.F., Wolfe, K., Parmar, R., Babendreier, J., Galvin, M., 2014. Design of a component-based integrated environmental modeling framework. Environ. Modell. Softw. 55, 1e24. Xu, J.J., Yang, D.W., Yi, Y.H., Lei, Z.D., Chen, J., Yang, W.J., 2008. Spatial and temporal variation of runoff in the Yangtze River basin during the past 40 years. Quatern Int. 186, 32e42. Xu, X.Y., Yang, H.B., Yang, D.W., Ma, H., 2013. Assessing the impacts of climate variability and human activities on annual runoff in the Luan River basin, China. Hydrol. Res. 44 (5), 940e952. Yang, D., 1998. Distributed hydrologic model using hillslope discretization based on catchment area function: Development and Applications. University of Tokyo, Tokyo, Japan. Yang, D., Koike, T., Tanizawa, H., 2004. Application of a distributed hydrological model and weather radar observations for flood management in the upper Tone River of Japan. Hydrol. Process 18 (16), 3119e3132. Yang, D., Musiake, K., Kanae, S., Oki, T., 2000. Use of the Pfafstetter Basin Numbering System in Hydrological Modeling. paper presented at 2000 Annual Conference, Japan Society of Hydrology and Water Resources. Yang, D.W., Herath, S., Musiake, K., 2001a. Spatial resolution sensitivity of catchment geomorphologic properties and the effect on hydrological simulation. Hydrol. Process 15 (11), 2085e2099. Yang, D.W., Herath, S., Oki, T., Musiake, K., 2001b. Application of distributed hydrological model in the Asian monsoon tropic region with a perspective of coupling with atmospheric models. J. Meteorol. Soc. Jpn. 79 (1B), 373e385. Zhang, G., Xie, H., Yao, T., Li, H., Duan, S., 2014. Quantitative water resources assessment of Qinghai Lake basin using Snowmelt Runoff Model (SRM). J. Hydrol. 519 (Part A(0)), 976e987. Zhou, J., Pomeroy, J.W., Zhang, W., Cheng, G., Wang, G., Chen, C., 2014. Simulating cold regions hydrological processes using a modular model in the west of China. J. Hydrol. 509, 13e24.