Notes and Correspondence
Investigating Runoff Sensitivity in the Land-Surface Model MATSIRO to Reduce Low Runoff Bias
2021 Volume 99 Issue 3 Pages 685-695
Details
2021 Volume 99 Issue 3 Pages 685-695
The Minimal Advanced Treatments of Surface Interaction and RunOff (MATSIRO), which has been used as a land-surface scheme in the global climate model, the Model for Interdisciplinary Research on Climate (MIROC), calculates Dunne runoff and base runoff using the TOPography-based MODEL (TOPMODEL). In past experiments that used MATSIRO, the runoff and its response to precipitation were too low compared to observation. We conjectured that those biases could be attributed to the water table's excessive depth. Its depth was diagnosed based on grid-mean soil moisture using a saturation threshold that was originally set to almost equal the porosity. In this study, sensitivity experiments, in which the threshold was decreased to 75, 50, 25, and less than 13 % of the porosity, were conducted, and the subsequent effects on river flow were investigated in the Chao Phraya River basin, Thailand, as a case study. As a result, both Dunne and base runoff increased along with the response of river flow to precipitation. The simulated river flow matched observations most closely with a threshold of 50 % saturation. In addition, soil moisture and the Bowen ratio changed significantly with the runoff changes induced by the threshold changes. These results suggested the importance of the relationship between grid-mean soil moisture and groundwater level for the TOPMODEL. Preliminary global experiments indicate that runoff sensitivity might be dependent on climate zone.
The importance of the role of terrestrial hydrological processes in climate change has been widely recognized since the 1980s (Manabe and Wetherald 1987), and the effects of vegetation canopy have been incorporated (Sellers et al. 1996). The Minimal Advanced Treatments of Surface Interaction and RunOff (MATSIRO; Takata et al. 2003; Nitta et al. 2014) model represents land-surface processes in a global climate model, the Model for Interdisciplinary Research on Climate (MIROC; e.g., Tatebe et al. 2019). MIROC is a core model used for climate change projections in the Intergovernmental Panel on Climate Change (IPCC) assessment reports, and MATSIRO has been used in various comparisons of land-surface hydrological models (e.g., Veldkamp et al. 2017).
In land-surface models used for climate studies, such as MATSIRO, spatiotemporal distributions of surface temperature, evapotranspiration, runoff, soil moisture and temperature, snow cover, and other variables are calculated from energy and water balances at the ground surface. In the past, MATSIRO significantly underestimated runoff compared to observations. A parameter for base runoff was calibrated (Hirabayashi et al. 2005), and the algorithm for diagnosing the depth of the groundwater table was improved (Yoshimura et al. 2006) at a global scale. Yet, MATSIRO still belongs to a group of land-surface models with low runoff (e.g., Fig. 2a of Haddeland et al. 2011).
In this study, we investigated the behavior of MATSIRO, focusing specifically on the runoff generation mechanism (i.e., the relationship between the water table depth and soil moisture). We conducted offline experiments using meteorological data for the Chao Phraya River basin in Thailand as a case study to see whether daily runoff values can be close to observation.
We performed a numerical experiment for the Chao Phraya River basin in Thailand, for which high-resolution input and verification data are available. The average annual temperature in Thailand from 1981 to 2010 was 26.9°C, and the annual precipitation was 1,588 mm, with 1,265 mm during the rainy season (May to October) and 322 mm during the dry season (November to April; Thai Meteorological Department 2015). The Chao Phraya River basin has an area of approximately 160,000 km2 and extends from central Thailand to the northern part of the country. The northern part of the basin is mountainous (with an altitude of about 1,500 m), and the central and southern parts are made up of vast plains (Fig. 1). There are two major dams: the Bhumibol dam on the Ping River in the western upstream area and the Sirikit dam on the Nan River in the eastern upstream area. These dams are responsible for controlling floods in the rainy season and supplying water for irrigation in the dry season (Shintani et al. 1994). The Bhumibol dam was selected for analysis in this study, and it has a catchment area of 26,214 km2 and a total water storage capacity of approximately 13,462 MCM (GRanD v1.1, https://globaldamwatch.org/grand/).
Map of the Chao Phraya River basin. Colors indicate altitude (m). The catchment of the Bhumibol dam is hatched with “1”. White dots indicate the locations of the Bhumibol dam, Kog Ma D site, and Tak.
MATSIRO was developed as a model of land-surface hydrological processes for a global climate model. It calculates surface temperature and heat and moisture fluxes using energy and water balances at the ground surface. The following processes are taken into consideration: the effects of the vegetation canopy on radiation, interception of precipitation, turbulent fluxes, and transpiration (modeled with a single layer), the accumulation and ablation of snow (with a variable number of model layers between one and three), and the moisture and temperature conduction in soil. Soil is divided into six layers with depths of 0–5, 5–25, 25–100 cm, 1–2, 2–4, and 4–14 m from the surface to the bottom. See Takata et al. (2003) and Nitta et al. (2014) for details.
The four types of runoff considered in our calculations are overland runoff when the soil moisture in the uppermost layer exceeds the porosity, infiltration excess runoff (known as Horton runoff) when the precipitation intensity exceeds the infiltration capacity at the ground surface, saturation excess runoff (known as Dunne outflow) which is generated from the runoff contributing area or saturated zone, and base runoff (Fig. 2). The sum of these four types of runoff is taken as the total runoff at each grid point. Dunne and base runoff are calculated based on the TOPography-based MODEL (TOPMODEL; Beven and Kirkby 1979; Sivapalan et al. 1987; Stieglitz et al. 1997) by simplifying the subgrid-scale topography in a grid box. Dunne runoff is generated from the runoff contributing area, taken as the saturated zone, calculated using the TOPMODEL. On the other hand, overland runoff and infiltration excess runoff (Horton runoff) are considered to be generated outside the contributing area.
Schematic figure of runoff in MATSIRO when the water table depth (zWTD) is (a) shallow and (b) deep. As is contributing area, Qdunne is Dunne runoff, Qbase is base runoff, Qo is overland runoff, Qi is infiltration excess runoff, and P is precipitation intensity at the ground surface.
Dunne runoff, Qdunne, is calculated by multiplying precipitation by the contributing area (As), and As is calculated from the water table depth (zWTD) on the basis of the TOPMODEL using Eq. (1):
 |
 |
The base runoff (Qbase) is also calculated based on the TOPMODEL using Eq. (3):
 |
zWTD is diagnosed using two steps. In the first step, the grid-mean matric potential (unit: m) in the k-th soil layer (numbered from the surface) is examined to assess whether the water table can exist in it. If the water table is expected to be in the k-th layer, in the second step, saturation judgment is conducted with the grid-mean volumetric soil moisture content in the k-th layer (θk, unit: m3 m−3) using Eqs. (4a) and (4b):
 |
 |
The total runoff, which is the sum of all runoff components calculated at each grid point, flows down along the river channel network scheme, the Total Runoff Integrating Pathways (Oki and Sud 1998). In this study, the spatial resolution of the channel network is 5 arc min (modified from Kotsuki and Tanaka 2013), and a constant flow velocity (0.3 m s−1) is assumed.
2.3 Experiments and analysesWe used the surface atmospheric data from the Integrated study on Hydro-Meteorological Prediction and Adaptation to Climate Change in Thailand (IMPAC-T) Forcing Dataset (IFD-K10) of Kotsuki et al. (2014). This dataset includes precipitation, pressure, and wind speed at a temporal resolution of 1 h, temperature and downward longwave radiation at a temporal resolution of 3 h, and humidity and downward shortwave radiation at a daily temporal resolution. Precipitation, which is one of the most important climate variables for runoff simulation, was derived from the gauge-based daily data in combination with the satellite-based mean monthly diurnal variation. The precipitation was thus rather representative as actual precipitation on a daily or longer time scale because it was based on gauged data (Kotsuki and Tanaka 2013), though its diurnal changes could include uncertainty. The time step of the model was set to 60 min. When the time of the input data and the time of the model simulation do not match, time interpolation is performed using the before-and-after input data, except that the daily short-wave radiation was divided into an hourly value accounting for the solar zenith angle. The study area was bounded by the coordinates 97–102°E and 13–20°N, with a spatial resolution of 5 arcmin (approximately 9 km) for the period 1981–2004. The initial condition was taken by conducting a separate spin-up experiment with recursive forcing of 1981, and it was used for all sensitivity experiments. In our analyses, we used results for 1993 and 2002 as case studies of lowand high-precipitation years, respectively. Vegetation distribution and monthly averaged leaf area index in 2000 (Takata and Hanasaki 2020) were used.
In the control experiment (denoted EP.0), εGW was set to the standard value of 0.001, and in our sensitivity experiments EP.1, EP.2, EP.3, and EP.4, this value was changed to 0.1, 0.2, 0.3, and 0.4, respectively (Table 1). Considering that θs ranged from 0.404 to 0.465 depending on soil type, a soil layer was judged to be “saturated” when θk was more than approximately 0.3, 0.2, 0.1, and 0.004–0.06 in EP.1, EP.2, EP.3, and EP.4, respectively (Eq. 4). In other words, it was judged to be “saturated” when θk was approximately 75, 50, 25, and only 1–13 % of the porosity in the sensitivity experiments. Consequently, zWTD is expected to become shallow if the vertical profile of soil moisture is the same because zWTD is diagnosed as being within the deepest “unsaturated” layer (see Section 2.2). The saturation area (As) is then expected to increase (Eq. 2), and Qdunne (Eq. 1) and Qbase (Eq. 3) are also expected to increase (Eq. 3). Note that the soil moisture calculation remained unchanged in this study. Only the threshold (εGW) for saturation judgment in diagnosing zWTD was changed in the examinations of runoff sensitivity.
We compared the observed river flow at the Bhumibol dam inlet to the calculated river flow at the model grid point nearest to the Bhumibol dam and examined the sensitivity to εGW. In addition, the sensitivity of runoff, zWTD, and soil moisture to εGW were investigated in the catchment area of the Bhumibol dam. Modeled river flow was compared directly to observations because there is no large reservoir or irrigation water intake upstream of the dam, and the river flow scheme in MATSIRO does not consider such human operations.
In addition, the sensible heat flux, the latent heat flux, and the Bowen ratio (i.e., the ratio of sensible heat fluxes to latent heat fluxes) were compared with the observation to examine the effects of the sensitivity experiments on surface energy and water balances. The observation data from Kim et al. (2014) recorded at Tak (16°56.390′N, 99°25.793′E; located about 50 km downstream of the Bhumibol dam) using the eddy correlation method were provided. The data were compared with the calculated results at the grid point nearest to Tak.
Figure 3 shows the time series of daily rainfall (used as model input values) averaged over the catchment of the Bhumibol dam and the time series of daily river flow at the inlet of the dam. The results from a low-precipitation year (1993, 724 mm yr−1) and a high-precipitation year (2002, 1221 mm yr−1) are shown, and the mean precipitation during the study period (1981–2004) was 968 mm yr−1. For both years, the control experiment EP.0 showed that river flow was significantly underestimated compared to observations, and there was almost no response of river flow to precipitation. The Nash–Sutcliffe efficiency (NSE), one of the most widely used validation indices for hydrology models, was calculated to compare the modeled daily river flow with the observation using Eq. (5):
 |
Daily precipitation (blue bars, mm day−1) averaged over the catchment area of the Bhumibol dam and hydrograph daily river flow (lines, m3 s−1) at the inlet of the dam in (a) 2002 (rainy year) and (b) 1993 (dry year). Observed river flow is indicated by a dotted black line, and calculated river flows are shown by solid lines (color legend shown above the panels).
Figure 4 shows the time series in 2002 of the modeled zWTD, Qdunne, and Qbase averaged over the Bhumibol dam catchment area for the control and sensitivity experiments. The results in the high-precipitation year are shown because seasonal changes are generally marked in high-precipitation years. Note that the arrival of peaks for Qdunne and Qbase was earlier than that of river flow at the inlet of the dam because of the time for passing through the river channels. In EP.0 and EP.1, zWTD (Fig. 4a) showed a similar seasonal change throughout the year, and it deepened gradually from January to May, reaching a maximum depth of about 6 m, subsequently rose, and became shallowest at about 5 m in December. zWTD in EP.2 was much shallower, and it deepened from January to May, became deepest at about 3 m, subsequently rose, and became shallowest at about 0.5 m in early September, followed by gradual deepening to December. Some smoothed increases occurred in response to precipitation events in November at intra-monthly time scales. In EP.3 and EP.4, zWTD was even shallower and did not deepen from January to May, remaining at a depth of 1.8 m in EP.3 and 1 m in EP.4. Quick rises in response to precipitation (at a time scale of approximately 5 days) were more marked than seasonal changes.
(a) Daily water table depth (m) averaged over the catchment area of the Bhumibol dam for 2002 (rainy year). (b) As in (a) but for Dunne flow (mm day−1). (c) As in (a) but for base flow (mm day−1). See the color legend above the panel.
The groundwater level measured from July to December in 2006 at the Kog Ma D site in the Bhumibol dam catchment area (Shiraki et al. 2017) rose from July to early September, and the shallowest depth was approximately 1.5 m. It deepened toward December, reaching about 2 m deep. Although it is not possible to make a direct comparison owing to the difference between the observation and experimental periods and the difference in the spatial scales, the seasonal change of zWTD in EP.2 qualitatively agreed with the observation, and it showed a gradual rise until early September, followed by gradual deepening until December.
Qdunne and Qbase both increased with increasing εGW. Qdunne (Fig. 4b) markedly increased between EP.1 and EP.2 and between EP.2 and EP.3, which suggests that the sensitivity of Qdunne was great near εGW = 0.2 (as in EP.2). Qbase (Fig. 4c) hardly changed between EP.0 and EP.1 but increased significantly in EP.2, EP.3, and EP.4. In EP.2, Qbase increased gradually from the beginning of the rainy season (from May), peaking in September, with responses to precipitation until November. In EP.3, Qbase increased in May, subsequently decreased in June and July, and then increased again with a peak in September with marked responses to precipitation during the rainy season. By contrast, in EP.4, Qbase increased in May and remained high throughout the rainy season.
Our results show that as the saturation threshold εGW increased, zWTD became shallower, Qdunne and Qbase became larger, and river flow increased. The increase in Qdunne was led by the increase in the saturated zone As, which induced decreases in infiltration of precipitation into the soil, reducing soil moisture. Figure 5 shows the vertical profiles of monthly soil moisture averaged over the catchment of the Bhumibol dam in 2002 for the experiments EP.0, EP.2, and EP.4. In EP.0, As was small, and infiltration into soil was large. Thus, soil moisture was high, and its seasonal amplitude at the surface was large, reaching to the fifth layer (2–4 m) (Fig. 5a). Because εGW was 0.001 for EP.0, the bottom layer (4–14 m) was the deepest “unsaturated” layer, and the water table averaged over the catchment was diagnosed as varying between 4.8 m and 5.9 m deep in 2002 (Fig. 4a). In EP.2 (εGW = 0.2), As became larger, and infiltration into the soil decreased. Consequently, soil moisture and its seasonal amplitude tended to be lower and smaller than EP.0 (Fig. 5b). The deepest “unsaturated” layer varied from the third to the fifth layers (0.25–4 m), and the water table averaged over the catchment was diagnosed as varying between 0.5 m and 3.1 m deep in 2002 (Fig. 4a). In EP.4, As was much larger, infiltration was much smaller, soil moisture and its amplitude were much smaller (Fig. 5c), and the water table averaged over the catchment was much shallower (varying between 0.1 m and 1.1 m in 2002; Fig. 4a).
Vertical profiles of monthly mean volumetric soil moisture content (m3 m−3) (lines with marks, legends of the colored marks for each month are shown right to the panel) averaged over the catchment area of the Bhumibol dam for 2002 in (a) EP.0, (b) EP.2, and (c) EP.4.
Such changes in soil moisture influence not only runoff (i.e., water balance at the surface) but also energy balance at the surface. In examining the effects of the sensitivity experiments on the surface energy and water balance, the sensible heat fluxes, the latent heat fluxes, and the Bowen ratios were examined along with comparison to the observation at Tak. Figure 6 shows the monthly values in 2003 and 2004 when the periods of observation and numerical experiments overlapped. Both fluxes and the Bowen ratios showed seasonal changes similar to the observations. The latent heat fluxes were overall smaller than the observations, though, because of the difference in spatial scale and meteorological conditions between the experiments and the observations. The values of both fluxes and the Bowen ratios for EP.0 and EP.1 were almost the same, and their seasonal amplitudes were smaller than the observations. The seasonal amplitudes of both fluxes became large for EP.2, and the amplitude of the Bowen ratio was as large as the observations. In EP.3 and EP.4, the seasonal amplitudes of the Bowen ratio were larger than the observations because of the high sensible heat fluxes and the low latent heat fluxes particularly in the dry seasons. These results imply that the change in εGW has significant effects not only on runoff but also on the energy and water balance at the ground surface.
Monthly mean (a) sensible heat fluxes (W m−2), (b) latent heat fluxes (W m−2), and (c) Bowen ratios measured (black line with dots) and modeled (all other lines; see legend) at Tak in 2003–2004.
MATSIRO is a land-surface hydrology scheme used in the global climate model MIROC, a major model used by the IPCC to make global warming projections since the fourth assessment report (Randall et al. 2007). Although the capability of MIROC is rated highly because of, for example, its ability to reproduce an El Niño event (Watanabe et al. 2010), runoff in MATSIRO has been underestimated in various offline experiments at global and regional scales (e.g., Haddeland et al. 2011). The method used in this study's sensitivity experiments is extremely simple but has not been tried since the model was developed, and it markedly changed runoff in the study area. We speculate that this is because of the conceptual gaps between the assumptions of the TOPMODEL and the application in MATSIRO.
Specifically, the TOPMODEL calculates Qdunne from the contributing area according to topographical structure and Qbase from groundwater with quasi-stationary recharge on slopes in small watersheds. In MATSIRO, these concepts are applied to a grid point whose size ranges from a regional scale (e.g., 10 km) to a global scale (e.g., 100 km). Thus, treatment of subgrid heterogeneity in topography and water table depth is essential to apply the TOPMODEL concept to a relatively large grid cell such as MATSIRO. The threshold change in this study is believed to exemplify such treatment.
The seasonal changes in river flow, which increased from the latter half of August to October, agreed with the observation in EP.2 (Fig. 3). However, the peaks in response to the intra-seasonal peaks of precipitation in September and November 2002 were approximately 10 days later than those observed. Moreover, the observed small peaks in May in response to precipitation were not represented in either year. Considering that the delays of the observed peaks from precipitation were very small, investigations into the river channel scheme and the uncertainty of observation would also be needed, but other runoff components and parameters should be further examined to sophisticate representation of the runoff response to precipitation.
Subgrid topography parameters and the vertical attenuation parameter for hydraulic conductivity can also affect runoff. In addition, the assumption of the uniform topographic index would be a tentative treatment. Lastly, possible introduction of other types of runoff should be considered, including subsurface stormflow, macropores, and other preferential flow paths (Brutsaert 2005). A more realistic form of the function should be used in future land-surface models. It should also be noted that runoff sensitivity also depends on the model's temporal and spatial resolutions.
4.2 Regionality of runoff sensitivityRunoff depends on climate, geology, topography, soil characteristics, vegetation, and land use, and the predominant runoff process varies in time and space depending on the precipitation intensity and ground-surface conditions preceding precipitation (Sivapalan et al. 1987). Runoff sensitivity may therefore vary regionally.
Next, global experiments were conducted offline with global surface meteorological data (ELSE; Kim et al. 2009). For experiment G-EP.0, εgw = 0.001, and for experiment G-EP.2, εgw = 0.2. The horizontal resolution was 1° × 1° (about 100 km), and the study period was 1979–2004.
The global mean annual runoff (averaged from 1995 to 2004) was 199 mm yr−1 for G-EP.0 and 327 mm yr−1 for G-EP.2. Thus, changing εgw from 0.001 to 0.2 increased the annual runoff by 64 %. Figure 7 shows the global distribution of the rate of change in annual runoff: (QG-EP.2 − QG-EP.0)/(0.5 × QG-EP.2 + 0.5 × QG-EP.0). The runoff sensitivity to εgw was high in regions with a savanna climate (Africa, India, Indochina Peninsula, northern Australia, and South America) and regions with a subtropical dry winter climate (east China, central North America, and Western Europe). In cold regions (north of 50°N), the runoff sensitivity to εgw was moderate. By contrast, runoff was not sensitive in regions with a semiarid climate (the Middle East, western China, and western North America) or regions with a rainforest (Indonesia). In semiarid regions, changing εgw did not change zWTD significantly, whereas in regions with a rainforest, where it is sufficiently humid, changes in zWTD did not noticeably change the total runoff.
Difference in annual runoff between G-EP.0 (εGW = 0.001) and G-EP.2 (εGW = 0.2) calculated as (QG-EP.2 - QG-EP.0)/(0.5 × QG-EP.2 + 0.5 × QG-EP.0). Mean values from 1995–2004 are shown. The regions where the difference between QG-EP.0 and QG-EP.2 has been smaller than 1 mm yr−1 are uncolored.
These results show that runoff sensitivity to εgw is not globally uniform. Therefore, it is important to fully survey modeled river flow by comparing it to observations to understand the relationship between soil moisture, zWTD, runoff components, and river flow.
Aiming to reduce the low biases in runoff with MATSIRO by focusing on the runoff generation mechanism, sensitivity experiments were performed to raise the water table by changing the saturation judgment threshold εGW in diagnosing zWTD (i.e., without changing the soil moisture). Specifically, the threshold was set to approximately 100 % (EP.0), 75 % (EP.1), 50 % (EP.2), 25 % (EP.3), and less than 13 % (EP.4) saturation. Consequently, zWTD became shallow in accordance with decreasing the judgment level of “unsaturated” soil moisture, Qdunne and Qbase increased, river flow increased, and the response of river flow to precipitation increased in the catchment area of the Bhumibol dam in the Chao Phraya River basin. The NSE for the observations was highest in EP.2 in both the high-precipitation and low-precipitation years. Soil moisture and the Bowen ratio changed significantly in accordance with the changes in εGW. The results of global experiments imply that runoff sensitivity to εGW may vary by region.
This research was supported by the Science and Technology Research Partnership for Sustainable Development (SATREPS) program of the Japan Science and Technology Agency and the Japan International Cooperation Agency, and by the Integrated Research Program for Advancing Climate Models (TOUGOU), Grant No. JPMXD0717935457, from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. The authors thank Dr. Seita Emori for providing the draft of Fig. 2.
