A snowmelt runoff forecasting model coupling WRF and DHSVM

This study linked the Weather Research and Forecasting (WRF) modelling system and the Distributed Hydrology Soil Vegetation Model (DHSVM) to forecast snowmelt runoff. The study area was the 800 km 2 Juntanghu watershed of the northern slopes of Tianshan Mountain Range. This paper investigated snowmelt runoff forecasting models suitable for meso-microscale application. In this study, a limited-region 24-h Numeric Weather Forecasting System was formulated using the new generation atmospheric model system WRF with the initial fields and lateral boundaries forced by Chinese T213L31 model. Using the WRF forecasts, the DHSVM hydrological model was used to predict 24 h snowmelt runoff at the outlet of the Juntanghu watershed. Forecasted results showed a good similarity to the observed data, and the average relative error of maximum runoff simulation was less than 15%. The results demonstrate the potential of using a meso-microscale snowmelt runoff forecasting model for forecasting floods. The model provides a longer forecast period compared with traditional models such as those based on rain gauges or statistical forecasting.


Introduction
In some high-altitude mountainous areas of western China, snowmelt water is an important water resource and plays a vital role in management of water resources. Snowmelt water is a primary source for reservoirs and water power sta-Correspondence to: Q. Zhao (dsslab@163.com) tions and plays an important part in controlling the quantity of water in the reservoirs and provision of water used in industry, agriculture and domestic life. Snowmelt water can ease drought in semi-arid and arid areas, but rapid spring snowmelt can cause a flood disaster (Zhao, 2007). Research shows that since the 1980s, the frequency and amount of snowmelt flooding have increased on the northern slopes of Tianshan Mountain Region. The frequency of snowmelt flooding in the 1990s increased 3 times when compared with that in the 1950s, causing serious damage to the national economy, and to arable land and properties in the region (Wu, 2003;Yan, 2003). As population continue to grow, the need for accurate forecasting of flood events is becoming increasingly important.
Traditional flood forecasting models use observed meteorological data. Therefore the forecast period is dependent on the flood routing in a watershed, often predicting floods only several hours in advance. We hope to achieve a longer flood warning forecast period of 1-3 days. High resolution atmospheric models for limited areas offer promisingly accurate regional forecasts of meteorological fields when forced with realistic large-scale conditions. Recent work coupling atmospheric models with hydrological models has shown that forecasting meteorological fields can be used to drive hydrological models to produce hydrographs at selected outlets. The forecast period can thus be extended when compared with traditional methods.
California. Andenson (2002), Lin (2002), and Lu (2006) adopted one-way or two-way atmospheric and hydrological coupling models to successfully forecast rainstorm floods and lengthen the flood prediction time. Kenneth (2001) used the Atmospheric Research Mesoscale Model (MM5) and the Distributed Hydrology Soil Vegetation Model (DHSVM) to simulate a complex rain-on-snow flood event.
This paper focuses on the direct forecasts of 24 h highresolution mesoscale Weather Research and Forecasting (WRF) to drive a distributed hydrological model (DHSVM) to predict the amount of snowmelt runoff. The snowmelt runoff forecasting model was assessed by performing a comparative result analysis between forecasted and observed data.
2 Brief description of the two models

Atmospheric model: WRF
The WRF modelling system is a next-generation mesoscale modelling system (Michalakes et al., 2001;Wang et al., 2004;Skamarock et al., 2005) that serves both operational and research communities. It is designed to be a flexible, state-of-the-art atmospheric simulation system that is portable and efficient on available parallel computing platforms. WRF is suitable for use in a broad range of applications across scales ranging from meters to thousands of kilometres.
The system consists of multiple dynamical cores, preprocessors for producing initial and lateral boundary conditions for simulations, and a three-dimensional variational data assimilation (3DVAR) system. WRF is built using software tools to enable extensibility and efficient computational parallelism. The use of the WRF system has been reported in a variety of areas including storm prediction and research; air-quality modelling; wildfire, hurricane, and tropical storm prediction; and regional climate and weather prediction (Welsh, 2004;Sun, 2003;Zhang, 2004).
The key component of the WRF-model is the Advanced Research WRF (ARW) dynamic solver. The model uses terrain-following, hydrostatic-pressure vertical coordinate with the top of the model being a constant pressure surface. The horizontal grid is the Arakawa-C grid. The time integration scheme in the model uses the third-order Runge-Kutta scheme, and the spatial discretization employs 2nd to 6th order schemes. The model supports both idealised and real-data applications with various lateral boundary condition options. The model also supports one-way, two-way and moving nest options. It runs on single-processor, shared and distributed-memory computers.
There are numerous physics options in the WRF model which are highly modular, transportable, and efficient in the parallel computing environment. There is an advanced data assimilation system developed in tandem with the model it-self. The simulations and real-time forecasting show that WFR model is able to forecast many kind of weather. The WRF model incorporates "online" chemistry; therefore WRF model system has a broad application for not only weather forecasting, but also for air quality forecasting.

Hydrological model: DHSVM
The DHSVM is a physically based, distributed hydrological model developed for use in complex terrain (Wigmosta et al., 1994). The model accounts explicitly for the spatial distribution of land-surface process, and can be applied over a range of scales, from a small plot to large watershed at sub-daily to daily timescales.
The DHSVM model includes a two-layer canopy model for evapotranspiration, an energy balance model for snow accumulation and melting, a two-layer rooting zone model and a saturated subsurface flow model. Digital elevation data are used to model topographic controls on incoming shortwave radiation, precipitation, air temperature, downslope surface water and soil moisture movement. At each time step the model provides a simultaneous solution to the energy and water balance equations for every grid cell in the watershed. Individual grid cells are hydrologically linked through a quasi-three-dimensional saturated subsurface transport scheme. The effects of topography on flow routing are obtained through the direct use of Digital Elevation Model (DEM) data. Each grid cell can exchange water with eight adjacent neighbours. Local hydraulic gradients are approximated by ground surface slope (kinematic approximation). Thus a given grid cell will receive water from upslope neighbours and discharge to the downslope.
Water confluence processes in DHSVM involve three parts: surface slope flow, road/route flow and soil moisture flow. Surface slope flow is the downslope surface runoff flow of remaining water following the processes of evapotranspiration, vertical infiltration and evaporation. However, infiltration will continue along the flow route during the surface slope flow process. Road/route flow occurs when the land surface category is road or route. During this type of flow evaporation may occur, however downward infiltration does not occur. Soil moisture flow is the lateral flow of water in the soil layer, in which horizontal diffusion in soil is accounted for.
There is a perfect snow accumulation and melt algorithm in the DHSVM model. DHSVM models the processes associated with snowpack morphology as described by Storck andLettenmaier (1999, 2000) and Storck (2000) using a two-layer ground snowpack representation of snow accumulation and melt. The snowpack model utilizes separate energy and mass balance components to represent the various physical processes affecting the snowpack. It also accounts for energy exchanges taking place between the atmosphere, overstory canopy, and main snowpack. The energy balance components of the model address snowmelt, refreezing, and  changes in snowpack heat content, while the mass-balance equations address the snow accumulation and ablation processes, transformations in the snow water equivalent, and snowpack water yield (Wigmosta, 2001).

Juntanghu watershed: the study area
The Juntanghu River is located on the northern slope of Tianshan Mountain Range, Xinjiang China ( Fig. 1). It is a small river, originating from the Tenniscar Glacier, starting with the Terssi. According to the statistical analysis from Geographic Information System (GIS), the headstream elevation is approximately 3400 m, and the main section is between 1000 m and 1500 m. Multiple streams converge at Mazal, located in the middle reach of Juntanghu River. The river then flows into the Red-Mountain Reservoir at the outlet of mountain area before entering the plains. The catchment area is approximately 800 km 2 , the catchment length is 45 km. The average elevation is approximately 1500 m, the slope of the upriver is 62.5‰, and the slope of downstream is 5.26%. The average annual runoff of this basin is approximately 3.89×10 8 m 3 . The watershed has some obvious hydrological characteristics of an arid area river, and can be divided into a runoff forming region and a runoff dissipation region, the boundary located at the outlet of the mountain area. One reason for choosing this basin as our study area is that it is relatively small with a close hydrological circumscription, and that snowmelt flood damage in the watershed is serious.

A limited-region 24-h Numeric Weather Forecasting System
Currently every Chinese meteorological station and meteorological service system can access the fourth-generation medium-term global numerical weather prediction system T213L31 forecast of the Chinese National Weather Service. In this paper, the T213L31 provided at 00:00:00 was used for the initial field and lateral boundaries of WRF v2.2. For this study the forecast period was 144 h: at 3 hourly intervals from 0 to 72 h and 12 hourly intervals from 72 to 144 h The Numeric Weather Forecasting System was run for 24-h meteorological forecasting everyday.

Numerical experimental plan
Basic parameters of simulated area: The central longitude and latitude was 86.5 • E and 44.0 • N respectively. The horizontal resolution was 1 km and the grid numbers in North-South direction and East-West direction were 130 and 121, respectively. There were 18 vertical layers. The total simulated time length was 24 h with a time step of 3 s. A forecast of meteorological fields was produced every hour.
Terrestrial data: Data included terrain elevation, land-use/vegetation, landwater mask, soil type, vegetation fraction and deep soil temperature obtained from USA AVHRR data. Soil class was based on United States Department of Agriculture texture. Terrain elevation was Global 30s DEM data. Vegetation category was US Geological Survey standard.
Physical process options: There are many physical process options in WRF for every parameterisation scheme. In this study, the schemes were selected as follows: the cumulus parameterisation was New Kain-Fritsch scheme; the microphysics scheme was WRF Single-Moment 3-class (WSM3); a rapid and accurate radiative transfer model (RRTM) long wave and Dudhia scheme were adopted for long-wave radiation and short-wave radiation; the planetary boundary layer scheme was Yonsei University (YSU); the 5-layer thermal diffusion surface physics scheme was selected.

Data-processing of meteorological fields
Temperature, humidity, wind speed, incident short and long wave radiation, and surface pressure at 2 m are required by DHSVM. Humidity and wind speed are not available directly from the WRF grid. However, the wind speed of every sigma lever in the simulation data of WRF model is filed. In this study, the wind speed at the lowest sigma level was used instead of the wind speed at 2 m. Humidity can be calculated by simulated water vapour mixing ratio and simulated surface air pressure. The formula is as follows: Where es is the saturation vapour pressure (pa) (Murray, 1967), t is temperature (K) at 2 m; RH is humidity (%), qv is water vapour mixing ratio (g/g) at 2 m, prs is surface pressure (pa).

Hydrological model initialisation
The DHSVM parameters can be broadly divided into two major steps. The first involves the assembly of surface characteristics data. This included digital elevation data, soil characteristics, vegetation, snow data and stream network information. The model requires attributes derivable from surface characteristics data for each pixel. This step was facilitated by use of GIS, with appropriate overlays for each of the attributes. The second step was to assemble the model forcing data, which consists of time series of meteorological variables and spatial overlays used to distribute these forcing data. This information was provided by the WRF model.

DEM
Elevation data taken directly from the DEM is used by DHSVM. Other topography attributes (e.g., surface slope and drainage patterns) were also derived from the DEM. DEM for the catchment was obtained from a 1:50 000 contour map at a spatial resolution of 30 m. The DEM was used to delineate the catchments. This procedure, which was implemented using an algorithm described by Jensen and Domingue (1998), is coded in most GIS programmes, including Arc/INFO routing flow-direction. Additional processing was performed to preserve general flow characteristics.

Soils
The DHSVM soil data were based on three types of information: soil type, soil physical parameter (e.g., lateral conductivity, exponential decrease, maximum infiltration, porosity, bubbling pressure, field capacity, wilting point, bulk density) and soil depth. The soil type data was obtained directly from the Chinese 1:1 000 000 soil type classification map, which was interpolated at a spatial resolution of 30 m in Arc/INFO (Fig. 2). Soil physical parameters were defined firstly according to the FAO global 17-catergory soil physical parameters data and the book "Soil in Xinjiang (Agricultural Bureau of Uygur Autonomous Region of Xinjiang, 1996)". Table 1 lists some of the soil parameters. Soil depth data were calculated by defining the maximum soil depth (1.3 m) and the minimum soil depth (0.25 m) according to observations based on DEM and a program provided by Washington University.

Vegetation
Five vegetation classes (grassland, farmland, water, bare land and evergreen needle leaf) were derived from Enhanced Thematic Mapper (ETM) classified satellite imagery (Fig.3). Data were processed to be similar in form to the data sets used by Kirschbaum (1997) and Matheussen, et al. (2000). In addition to land classification type, the DHSVM vegetation parameters (e.g., height, maximum resistance, minimum resistance, moisture threshold, vapour pressure deficit, monthly Leaf Area Index (LAI), monthly albedo) were defined according to field observation and USGS 25-category vegetation physical parameters. Table 2 lists some of the vegetation parameters.

Snow information
Snow information is very important for stimulating snowmelt runoff, and for the DHSVM model is requires this as an initial snow state file. Snow information included snow cover and spatial distribution of snow water equivalent. In this  Satellite reflectance in MODIS bands 4 and 6 were used to calculate the NDSI for the snow cover map. A pixel is mapped as snow if the NDSI value is >0.4 and the reflectance in MODIS band 2 is >11% (Barton, 2001;Klein, 2003;Salomonson, 2004).
The spatial distribution of snow water equivalent was obtained from the National Centre for Atmospheric Research (NCAR) Final analysis data.

Stream network
The DHSVM stream network was based on three types of information: a mapping table which located a portion of stream reaches within its appropriate grid cell and described the depth, width, and aspect of the channel cut into the soil; a reach table describing the length, slope, and class of a reach connected with the next reach downstream; and a class file with routing characteristics of width, depth, and roughness for each stream class. These files were derived from the DEM using an algorithm described by Wigmosta and Perkins (2001).
In essence, the DEM topology defines the stream locations, while the extent of the network is specified by the model user via a given support area (minimum area below which a stream channel is assumed to exist). For Juntanghu catchment, the contributing area is 324 000 m 2 , which was in part based on field observations. Stream order was used to perform an initial classification of reaches. This produced a manageable number of types to which channel characteristics could be indexed. Manual adjustments based on limited field observation were conducted where necessary (Fig.4). Class characteristics were defined according to field observations (where available) and the relative descriptive size of the classes. The channel depth, width, and roughness were classified according to the reach classifications using GISWA algorithms (Wigmosta and Perkins, 1997). The average slope of each reach was calculated using the DEM.

Analyses of forecasting meteorological field results
The time period for the case study was from 01:00:00 GMT on 29 February 2008 to 00:00:00 GMT, 6 March 2008. This study used the WRF with the initial fields and lateral boundaries provided by the T213L31 model to realise the 24 h-numerical weather forecast from 29 February 2008 to 6 March 2008. Figure 5 shows the comparative map of observed meteorological data and corresponding grid forecasts.
From Fig. 5, we can see that: (1) There are certain deviations between forecasted and observed temperature at the highest and lowest points. The average error is 1.2 K; (2) The forecasted wind speed is higher than the observed data. This is because WRF only supplies wind speed in every sigma layer, so the wind speed in the lowest sigma layer was adopted. Because the wind speed is a relatively weak influence on the snowmelt runoff and the error is only 1.54 m/s, the forecasted data is acceptable; (3) The relative humidity forecasted error is larger when the humidity fluctuates significantly. When the relative humidity is stable, the forecast is more accurate. The overall average error is 6%. (4) The forecasted results of solar radiation are generally good. Forecasted data is significantly higher than observed data at midday when the water vapour is higher and cloud activity is greater. It is difficult to consider the effect of clouds due to WRF low spatial resolution.
Overall the forecasting errors are relatively small, proving that limited regional numerical weather forecasting precision can meet the requirements of accurate snowmelt runoff forecasting.

Improvement of DHSVM parameters
Hydrology, vegetation and soil parameter schemes have been successfully developed for simulation in North America. A total of 33 parameters were calculated and adjusted in terms of basin climatic and natural conditions. To apply DHSVM model system to snowmelt runoff modelling, the parameter scheme must be improved and reviewed. In this study, all 33 parameters were recalculated and reset by using up-todate hydrometeorology theory and methodology, focussing on critical parameters such as soil porosity, vertical saturated hydraulic conductivity, maximum infiltration rate and coefficient of roughness for each layer of soil type.
During the spring melt season, there is little Evergreen Needle leaf coverage within the study area. As deciduous foliage is not yet present or is covered by snow, LAI and height of vegetation (except Evergreen Needle leaf) were classifieds as bare land. There are two important soil parameters: Maximum Infiltration rate and Manning's n which need to be adjusted in snowmelt runoff modelling.

(1) Maximum infiltration rate:
Seasonal ground frost is widespread in the catchment during spring melt season. The spatial distribution of frozen soil and snow cover at the start of the spring melt season plays an important role in the generation of spring runoff. Many field studies on snowmelt infiltration into frozen soils have been reported in the literature (Kane and Stein, 1983;Burn, 1991;Gray, Toth and Zhao, 2001;Cherkauer, and Lettenmaier, 2003;Niu and Yang, 2006;Zhang and Sun, 2007;Ye, 2009).
Hydrologically, frozen soil suppresses infiltration and encourages surface runoff. In this paper, we empirically hypothesise that as the seasonal frozen soil is distributed  under snow cover regions, the maximum infiltration rate of frozen soil is 1.0 e −6 .
(2) Manning's n: When compared to the time of energy input at snow surface, the delay of snowmelt runoff is due to the water holding capacity of the snowpack and the horizontal travel time of melt water along the ground. This results in a delay in the peak time of daily runoff. In this paper the soil parameter, Manning's n (coefficient of roughness), was adjusted so that the simulated daily flood-peak time matched the observation data.
There is no hydrological and meteorological station in the study area. We have observed the snowmelt process for 3 years (2006, 2007, and 2008). We have observed the daily flood-peak time during the spring melt season in 2006 and 2007. However for the purpose of this study there was insufficient time series runoff data available. Several parameters (e.g., coefficient of roughness, stream network parameters) were adjusted based on observations from 2006 and 2007. With the new model parameter schemes, the forecasted snowmelt runoff agreed with the record database. Modelling efficiency was better than that with original parameter schemes (Fig. 7).

Analyses of snowmelt runoff results
The DHSVM model was forced by forecasted meteorological fields at a spatial resolution of 1 km. However, the DHSVM model is initialized at a spatial resolution of 30 m. Therefore there is an interpolated programme embedded within the DHSVM model, which is based on the DEM data (30 ) used by WRF model, the high-resolution DEM data (30 m) and temperature gradient. Figure 6 shows the spatial change of snow water equivalent. Figure 7 shows the comparative map of 24 h-forecasted discharge model with the new model parameter schemes, forecasted discharge with the old model parameter schemes and observed discharge from the outlet of Juntanghu basin from 29 February 2008 to 6 March 2008. The model efficiency coefficient and the relative error of maximum value are used to evaluate the efficiency of the snowmelt runoff forecasting model. Table 3 shows the test of runoff forecasted with the new model parameter schemes and observed results.
The model efficiency coefficient: Where Q obs is the observed discharge (m 3 /s), Q obs is the average observed discharge (m 3 /s), Q fore is forecasted discharge (m 3 /s).
The relative error of maximum value: Where Q obs.m , and Q fore.m is the maximum observed and forecasted discharge (m 3 /s) respectively.
From Fig. 7 and Table 3, the following results can be observed: (1) The average efficiency coefficient is 0.8, which shows the forecasted data is in strong agreement with the observed data. The same trends are present in hydrological processes; and (2) The maximum relative error of maximum data, which is very important for flood warning, is 13.2%. This means the snowmelt runoff forecasting model is able to meet the needs of snowmelt flood forecasting and flood warning.

Conclusion
Based on the latest development of atmospheric science and hydrology, using the features of snowmelt flooding on the northern slopes of Tianshan Mountains, this study has built a snowmelt runoff forecasting model by coupling WRF and DHSVM. The forecasted results of this model have been verified. This was achieved as follows: The limited-region 24-hour Numeric Weather Forecasting System was established by using the new generation atmospheric model system WRF2.2 with the initial fields and lateral boundaries provided by the T213L31. Overall, weather predictions were in accordance with observational data. The atmospheric and hydrologic models were coupled and a 24-h snowmelt runoff forecasting model was run using forecasted meteorological fields to force the DHSVM model. With the new parameter schemes taking seasonal ground frost and snow cover into consideration, the simulated data showed strong agreement with the observed data. The average absolute relative error of the maximum runoff in simulation is below 15%. The model has successfully achieved practical snowmelt runoff forecasting.
This study provides safeguards for flood early warning systems, flood prevention and disaster reduction, and water resource management through the forecasting of the mesomicroscale snowmelt runoff forecasting model. Coupling the atmospheric and hydrologic models can offer useful reference for hydrological forecasting and water resources management in areas where there is no observed data or incomplete data.
The results demonstrate the potential of using mesomicroscale snowmelt runoff forecasting model for flood forecasting. The model can provide a longer forecast period compared to traditional models such as those based on rain gauges, or statistical forecasting.