The improved Noah land surface model based on storage capacity curve and Muskingum method and application in GRAPES model

Noah‐LSM, applied in GRAPES model, cannot express the changes of runoff‐yield area effectively and describe hydrological cycle completely. It is improved as follows: (1) Add saturation excess runoff generation module, which considers the change in runoff‐yielding area based on a storage capacity curve for the better grid runoff simulation; (2) introduce routing module based on Muskingum method, to consider the redistribution of soil moisture in two‐dimensional level. The improved model application results indicated that the meteorological elements in the surface layer simulated are influenced by land surface water cycle process, such as soil moisture, 2 m temperature, is closer to the observed.


Introduction
Land surface can remember weather events through variations in its heat and water storages. In turn, land heat and water storage can affect atmospheric predictability through their effects on surface energy and water fluxes (Roesch et al., 2001;Jiang et al., 2009). As the water exchange in the different forms and the different scales, the relationship between atmosphere and hydrology is mutually related to each other and mutually restricts each other. Hydrological cycling between the atmosphere, vegetation, and soil plays an important role in land surface processes. The changes of elements in hydrological cycle are feedback to the atmosphere, through the latent heat flux and sensible heat flux. Therefore, hydrological process in land surface model is a critical link between land and atmosphere. It is crucial to understand interactions between hydrological cycle and atmospheric system, which includes controls of atmospheric conditions to regional hydrologic system and feedbacks of hydrological cycle to atmospheric system. The simulation of water balance is important in atmospheric system which is directly influence the evapotranspiration, surface runoff, interflow, and underground runoff in hydrological processes (Cao and Liu, 2005). However, in most land surface model, the implemented hydrological process simulation is too simple to predict soil moisture accurately, even the water cycle. Several important problems exist in most land surface models. First, heterogeneity of the runoff-yielding area is not fully considered. For example, infiltration capacity is usually assumed to be homogeneous within a grid cell, but it is likely highly heterogeneous within a cell. Second, land surface models usually lack sophisticated runoff concentration. It is true that physical mechanisms and theories of surface hydrologic processes have been well studied for a long time; hydrologists are able to model these processes adequately -even though the systems and processes are spatially heterogeneous -if many important topographic, soil, vegetation, and climatic parameters are known. In China, the Xin'anjiang model is widely applied in flood forecasting and runoff simulation (Zhao, 1983;Li et al., 2006). It is a well-tested hydrological model. Therefore, it can be used to improve the modelling of hydrological processes in land surface models in terms of providing better simulation runoff. As a result, it should be able to largely improve the simulation of land surface water cycle in land surface models.
Noah-LSM is a well-known land surface model, which couples the diurnally dependent Penman potential evaporation approach of Mahrt and Ek (1984), the multilayer soil model of Mahrt and Pan (1984), and the primitive canopy model of Pan and Mahrt (1987). It has been extended by Chen et al. (1996) and Chen and Jimy (2001) to include the modestly complex canopy resistance approach of Noilhan and Planton (1989) and Jacquemin and Noilhan (1990). During the 1990s, National Centers for Environmental Prediction (NCEP) greatly expanded its land surface modelling collaborations via several components of the Global Energy and Water Cycle Experiment (GEWEX), most notably, the GEWEX Continental-Scale International Project (GCIP) and the Project for Intercomparison of Land-surface Parameterization Schemes (PILPS). These collaborations included the Office of Hydrological Development (OHD) of the National Weather Service, National Environmental Satellite Data and Information Service (NESDIS), NASA, National Center for Atmospheric Research (NCAR), the US Air Force, and OSU and other university partners. As an outgrowth of these collaborations and their broad scope of LSM testing in both uncoupled and coupled mode over a wide range of space scales and timescales, NCEP substantially enhanced the OSU LSM, now renamed the Noah-LSM in recognition of the broad partnership above (Ek et al., 2003).
The Global/Regional Assimilation and Prediction System (GRAPES, Chen and Shen, 2006;Chen et al., 2008;Xue and Chen, 2008) is an important operational numerical weather prediction system in the China Meteorological Administration (CMA). This model adopts a structure of standardized and module-based software in accordance with the strict requirements of software engineering.A preliminary study (Wu et al., 2005) shows that the application of GRAPES model meets the requirement for sustainable development of the numerical prediction system of China. Nowadays, GRAPES model has been expanded to applications in various fields, such as GRAPES_Meso for mesoscale weather prediction, GRAPES_TCM for typhoon prediction, GRAPES_SDM for sandstorm forecast, and GRAPES_SWIFT for short-time weather forecasting (Zhao and Li, 2006;Zhu et al., 2007). Further development of GRAPES is undergoing. So far, however, the operational GRAPES model in the CMA has not yet been able to directly predict runoff and flood events. Therefore, it is necessary to improve the Noah-LSM model that accommodates the whole process of the closed land surface water cycle, and is able to predict the runoff or discharge. But there is obviously insufficient that is the Noah-LSM land surface model for hydrological processes especially runoff description. It cannot simulate the complete description of the hydrological cycle, on the basis of the land surface hydrological using parameterization scheme is a simple water balance model SWB (Simple Water Balance, Entekhabi et al., 1999).
The aims of this article, therefore, are to (1) improve Noah-LSM with introducing saturation excess runoff generation to descript the inhomogeneity of the yield runoff area and routing module, which can represent better land surface water cycle processes and (2) apply the newly improved Noah-LSM to describe hydrological process of GRAPES_MESO model and compare with the forecast results of the near-surface meteorological elements of the original model.

Hydrological process on Noah-LSM
In Noah-LSM, the Simple Water Balance (SWB) model is used to calculate the surface runoff (R). The SWB model (Schaake et al., 1996) is a two-reservoir hydrological model. A thin upper layer consists of the vegetation canopy and the soil surface. A lower layer includes both the root zone of the vegetation and the groundwater system. The surface runoff, Rs, is difference between precipitation (P d ) and maximum infiltration (I max ): I max is formulated as where D x is the soil water diffusivity; Θ s is the maximum soil moisture (porosity); Θ i is the soil moisture content of the ith layer; ΔZ i is the thickness of the ith soil layer; kdt = 3.0 and K ref = 2 × 10 − 6 m s −1 ; i is the conversion of the current model time step t (in terms of seconds) into daily values.
The underground runoff is as computed in Noah-LSM model as: Others (5) where Q max is the maximum underground runoff, S max is the critical value of soil water moisture deficit, D b is the soil water shortage in the lower soil layer, the maximum value is the lower soil water content.

Problems in hydrological process modelling of Noah-LSM
In the rainfall-runoff process, runoff yielding area is an area that contributes to the runoff flow in the basin. It changes during the precipitation process. According to Equations (1) and (4), only vertical movement of water is considered in the runoff module in Noah-LSM model. It cannot effectively simulate spatial changes in runoff-producing area during the rainfall (Abbott et al., 1986). Therefore, runoff module of NOAH-LSM should be improved in the runoff generation and routing.

Improvement on hydrological process modelling of Noah-LSM
The concept of runoff formation on repletion of storage holds true so far as one single point, in the sense of a very small area, of a basin is concerned. For an entire watershed, things are more complicated. The moisture deficit often varies from place to place. This non-uniform distribution effects the runoff production of a whole basin significantly. To solve the problem, the tension water storage capacity curve (Zhao, 1983;Zhao and Liu, 1995), which represents the concept of runoff formation on repletion of Storage, is suggested to improve the runoff module in Noah-LSM.
where f /F represents the proportion of the pervious area of the basin whose tension water capacity is less than or equal to the value of the ordinate; W ' is point tension water capacity, varies from zero to a maximum W MM according to the relationship (mm); W MM is the maximum tension water capacity (mm); B is an exponent of the tension water capacity distribution curve. Based on the relationship between storage capacity curve and exchange between rainfall and runoff, the improved runoff generation in Noah-LSM is computed as follows: When P − E + A < W MM ,only part of the grid cell is producing runoff: When P − E + A ≥ W MM ,the whole grid cell is producing runoff: where P is areal mean rainfall in the grid (mm); E is areal mean evapotranspiration in the grid; W 0 is the initial areal mean tension water storage in the grid (mm); W M is the areal mean tension water capacity in a grid cell (mm), W M and W MM are related through the parameter B; R is the runoff generation in the grid cell (mm); ] .
In the tension water storage capacity curve and runoff generation model, P is the input, E and W 0 can be calculated with according to measured pan evaporation and soil moisture, the relation between W M and W MM is easy to be shown by W MM = W M (1 + B), and B can be calibrated by basin area and underlying surface factors (Zhao and Liu, 1995). Therefore, the parameter estimate of W M is very important and necessary. W M corresponds to available soil water defined as the difference between field capacity and wilting point: where fc is field capacity; wp is wilting point; L a is thickness of aeration zone (mm). The values of fc and wp were determined from soil texture following Williams et al. (1998) and Anderson et al. (2002). To derive the spatial distribution of W M , thicknesses of aeration zone and humus soil layer should be given beforehand. Soil thickness, which is an important factor in hydrologic processes, can vary spatially because of complex interactions of many different factors; practical measurement of soil thickness is laborious and time consuming (Tesfa et al., 2009). Thus, we developed a simple method to predict roughly the spatial pattern of soil thickness from topographic, soil texture, and land cover attributes.  Assuming that topography is the dominant source of heterogeneity in thickness of aeration zone L a , a simple scheme based on topographic index (Beven and Kirkby, 1979) was developed: where TI is topographic index, a and b are two coefficients.
The TI reflects the tendency of flow to accumulate in the watershed and the tendency for gravitational forces to move that flow downslope (Quinn et al., 1991;Yao et al., 2012). A location with high TI value usually has a high water table, i.e. a thin aeration zone. It is reasonable to expect relatively low moisture capacity W M in such a location. It is thus valid to assume that hydrologic unit with the minimum or maximum value of TI correspond to the unit with the minimum or maximum value of W M (Zhao and Liu, 1995;Yao et al., 2012).
In terms of the soil textures, a and b are calculated by (Yao et al., 2012): where W M max and W M min could be easily pre-set by referring to the literature (Zhao, 1992;Zhao and Liu, 1995) or the experience with implementation of the Xinanjiang model.

Adding routing module to hydrological process modelling of Noah-LSM
The runoff flow produced by the rainfall could flow with the terrain. If there is no routing module, runoff redistribution in the horizontal direction cannot be accounted for; it will cause the runoff to accumulate in the grid cell without flowing to the downstream areas, and an incomplete simulation of water cycle. In land surface model, the runoff directly affects the sensible heat flux and the latent heat flux, which have important feedbacks to the meteorological model. By adding the routing module, it consider the movement of the runoff in the horizontal direction, and making the simulation closer the reality and potentially more accurate. The Muskingum routing method (Bates and De Roo, 2000) is adopted into Noah-LSM in this study. As shown in Figure 1, the flow in the grid cells a, b, and c can all flow into the grid d; the outflow are Q a , Q b , and Q c , respectively. Q a ', Q b ', and Q c ' are obtained by the Muskingum flow routing method, where where x e and k e are the two parameters of the Muskingum flow routing method.
As a result, the outflow Q d of the grid d is the sum of Q a ', Q b ', Q c ' and the runoff Q d0 produced at the grid d: The depth of flow and discharge in reach can be derived from empirical relationships recommended by the US Reclamation Service (Chow, 1959), hence, the Muskingum k e and x e parameters can be estimated with simplifying overland routing for wide rectangular channel routing, channel routing for parabolic channel routing (Tewolde and Smithers, 2006;Bao, 2009;Li et al., 2008).
The Muskingum k e is estimated from Equation (17) as follows (Chow, 1959;Fread, 1993): where ΔL is reach length (m), V w is celerity (m s − 1 ). For a wide rectangular channel, the celerity may be estimated by V w = 5 3 V av . For a parabolic channel, the celerity may be estimated by V w = 11 9 V av , where V av is the average velocity, which can be calculated from the Manning equation (Chow, 1959;Tewolde and Smithers, 2006;Bao, 2009).
The Muskingum x e is estimated by (Fread, 1993;Tewolde and Smithers, 2006): where Q 0 is reference flow, which can be estimated by Q 0 = I b + 0.5(I P − I b ). I b is minimum discharge and I P is peak discharge. S is slope, P is wetted perimeter (m),which can be calculted by P = c √ Q 0 . c is coefficient (between 4.71 and 4.81) (Bao, 2009).

Experiments and analysis
In order to compare the performance of the improved Noah-LSM model with the original one, an experiment is conducted for simulating the precipitation from 1 August 08 BJT (Beijing Time), corresponded for 1 August 00 UTC, to 1 October 08 BJT. The prediction time is 48 h for GRAPES_MESO model within the domain (15 ∘ -64.5 ∘ N, 70 ∘ -145.3 ∘ E), on 08 BJT daily. The experiment are conducted with the initial fields and lateral boundaries provided by NCEP forecast fields with 1 ∘ × 1 ∘ resolution, while the GRAPES_MESO model has a 0.15 ∘ × 0.15 ∘ horizontal resolution. Because the improvement is on the hydrological processes of Noah-LSM, soil moisture, especially the soil moisture of the first level (in 10 cm below the ground), is directly affected, and 2 m temperature is a  after the improvement. Therefore, these variables are selected for comparison. More importantly, all these variables were also measured. It is worth to note that the distribution of the soil moisture stations is different from that of the 2 m temperature stations.

Soil moisture
Comparing to the original Noah-LSM, not only the vertical variation of water but also its horizontal redistributionhas been changed in the improved Noah-LSM. The overland flow affects the horizontal distribution of the soil moisture, especially the soil moisture of the first level. Land surface water cycle caused the soil moisture change at first. The observed soil moisture records at 10 cm of Nanchong station and Linfen station were chosen for the comparison. Table 1 shows the latitude and longitude of the stations. The observed soil moisture at Nanchong station starts to rise at 08 BJT 1 August as shown in Figure 2(a). It becomes saturated from 3 August to 30 September (Figure 2(b)) in Linfen station. Comparing to the observation, the GRAPES_MESO results are generally low, but the improved results are much closer to the observations. This showing that improving the new runoff module and adding the routing module help improve the simulation of soil moisture indeed. In land surface model, change in soil moisture directly affects the simulation of soil temperature. Furtherly, this effect can provide feedback to the atmospheric model, to impact the simulation of the temperature in 2 m, wind speed in 10 m, and other variables Peters-Lidard et al., 1998;Chen and Jimy, 2001).

Temperature in 2 m
The change of the soil moisture can affected the 2 m temperature, which is another key variable for GRAPES_MESO model in our continuous experiments. In this study, the 2-m observation stations include Shouxian, Shangqiu, and Macheng stations ( Table 2). It is clear that the 2-m temperature series in GRAPES model are different from those in the revised model in Figure 3. The root mean square error of temperature in 2 m simulated by the improved one is lower than the GRAPES_MESO model (Table 3). This indicates that the improved hydrological process has also help to improve the simulation of temperature in 2 m.

Wind speed in 10 m
The soil moisture can also affect the simulated on the wind speed in 10 m. Observations at three stations were selected to compare with the simulations (Table 1). In Figure 4, the 10-m wind speed modelled by the improved model is slightly different to the results of the original model. According to the root mean square error on wind speed in 10 m, the improved one is slightly better than the original one (Table 4). The better results in 10-m wind speed simulation highlight the important consequence of the improvement on the hydrological processed conducted in this study.

Precipitation
Threat score (TS) as the point-based precipitation verification method is commonly used in weather forecast operation and model research (Huang, 2001). In this article, 2513 stations are chosen for precipitation verification. The improved and original models have close TS scores and bias errors as shown in Table 5. We further selected two stations for precipitation verification (in Table 6). It is clear that the simulated precipitation by the improved model is different from this by the original model ( Figure 5), which proves that the changes in the hydrological process has a direct consequence on the simulation of the rainfall regime.

Discharge
In this study, due to the limited stream flow observations, the observed flow at Wangjiaba station from 13 to 16 August area chosen for verification in this experiment. The contributing area of Wangjiaba stationis about 30 672 km 2 and is the upper stream of Huai River, where channel bed slope and flow velocity are generally large. The first key flood control gate of the catchment is located at Wangjiaba station. Behind this gate is the Mengwa flood retention zone, with a design capacity of 750 million m 3 and a design maximum discharge of 1626 m 3 s −1 . The area, during drier periods, usually serves as farmland of approximately 12 000 ha for a local population of about 157 800. The water level at Wangjiaba station is a key flooding indicator for the entire catchment. It is therefore important to obtain reliable forecasting of discharge at Wangjiaba station. Seen from Figure 6, the accuracy of runoff flow simulation has been greatly improved. This is shown that the improved hydrological process in GRAPES_MESO model has the ability of the simulation on the runoff flow.

Conclusions
The Noah-LSM in GRAPES_MESO model does not consider the heterogeneity distribution within grid cell. For example, infiltration capacity is usually assumed to be homogeneous within a grid cell, but it is likely highly heterogeneous within a cell. Land surface models usually lack sophisticated runoff concentration, because hydrological module is too simple. In this study, we improved the hydrological process modelling in the Noah-LSM model. First, we incorporated saturation-excess runoff generation module, which considers the change of runoff-yielding area based on a storage capacity curve. Second, we introduced a routing module based on the Muskingum method into the model to consider the redistribution of soil moisture in two-dimensional level and to improve the simulation accuracy of the surface runoff. To evaluate the performance of the improved model, we applied both original and improved models to conduct simulations from August to September 2009 and further investigated the feedback of changes in the land surface water modelling to meteorological simulation. The results show that the revised land surface water modelling has affected the results of simulated meteorological variables such as soil moisture, air temperature, and wind speed and generally led to better results. This suggests that the model improvement conducted in this study is valid and effective. It also highlights the important connection between surface hydrology and meteorology.