Co-evolution of volcanic catchments in Japan

Present day landscapes have evolved over time through interactions between the prevailing climates and geological settings. Understanding the linkage between spatial patterns of landforms, soils, and vegetation in landscapes and their hydrological response is critical to make quantitative predictions in ungaged basins. Catchment co-evolution is a theoretical framework that seeks to formulate hypotheses about the mechanisms and conditions that determine the historical development 5 of catchments and how such evolution affects their hydrological response. In this study, we selected 14 volcanic catchments of different ages (from 0.225 to 82.2 Ma) in Japan. We derived indices of landscape properties (drainage density and slope area relationship) as well as hydrological response (annual water balance, baseflow index, and flow duration curves) and examined their relation with catchment age and climate (through the aridity index). We found a significant correlation between 10 drainage density and baseflow index with age, but not with climate. The intra-annual flow variability was also significantly related to catchments age. Younger catchments tended to have lower peak flows and higher low flows, while older catchments exhibited more flashy runoff. The decrease in baseflow with catchment age is consistent with the existing hypothesis that in volcanic landscapes the major flow pathways change over time from deep groundwater flow to shallow subsurface flow. 15 The drainage density of our catchments decreased with age, contrary to previous findings in a set of similar, but younger volcanic catchments in the Oregon Cascades, in which drainage density increased with age. In that case, older catchments were thought to show more landscape incision due to increasing near-surface lateral flow paths. Our results suggests two competing hypotheses on the evolution of drainage density in mature catchments. One is that as catchments continue to age, the 20 hydrologically active channels retreat because less recharge leads to lower average aquifer levels and less baseflow. The other hypothesis is that the active channels do not undergo much surface dissection after the catchments reach maturity.


Introduction
The hydrological functions that partition incoming water, energy, and carbon at the land surface vary from place to place, reflecting spatial variability in land surface characteristics (landforms, soils, vegetation).Because of this spatial variability, hydrological predictions in ungaged basins are highly uncertain (Wagener et al., 2007;Sivapalan, 2003).Improving these predictions will require a better understanding of the different interactions between landscape features and hydrological response, especially in a changing environment (Mc-Donnell et al., 2007;Wagener et al., 2013).
Present-day land surface characteristics in pristine landscapes have co-evolved over time through various interactions between the prevailing climates and geological settings (bedrock type, tectonic uplift) (Blöschl et al., 2013;Troch et al., 2015).Hydrologists have coined the term catchment coevolution as a theoretical framework that seeks to formulate explanatory hypotheses about spatial patterns of landscape characteristics and the corresponding hydrological response, based on the historical development of catchments (Sivapalan et al., 2012;Troch et al., 2013;Harman and Troch, 2014).The results of historical development are captured in the soil and hillslope catena profiles, the geometry of channel networks and the spatial distribution of plant functional types.This landscape organization determines the hydrological functions that partition incoming water, energy, and carbon.A detailed mechanistic reconstruction of how T. Yoshida and P. A. Troch: Coevolution of volcanic catchments in Japan catchments have evolved to their present state is a daunting task, given the many processes involved and the unknown initial conditions during landscape formation.
Therefore, formulating and testing explanatory hypotheses of catchment coevolution and hydrological response requires a different approach from the classic mechanistic Newtonian approach.Recently, an alternative approach, called Darwinian hydrology, was proposed that seeks to document patterns of variation in populations of hydrological systems, and explain these patterns in terms of the processes and conditions of their historical development (Wagener et al., 2013;Harman and Troch, 2014).Harman and Troch (2014) summarized what the Darwinian approach might look like in watershed science as follows: (1) collect extensive and detailed observations of the pattern to be explained, (2) conceive a hypothesis that accounts for as many of the observations as possible, (3) derive from this hypothesis a subsequent set of circumstances that also must hold if the hypothesis is true, and (4) critically test whether these consequences do hold.In this study, we build upon this framework and explore empirically the links between landscape properties, hydrological response, and time since formation (age) of catchments of similar bedrock geology (volcanic rock deposits of different geologic age).
It may be argued that climate exerts a first-order control on hydrologic partitioning and that physical, chemical, and biological catchment properties only have a second-order impact (e.g., Berghuijs et al., 2014;Sivapalan et al., 2011).However, these second-order controls explain much of the functional patterns that determine flow regimes.Mushiake et al. (1981) observed that bedrock lithology plays a key role in the natural variations of functional patterns in Japanese catchments.They showed that catchments on older sedimentary rocks formed in the Paleozoic and Mesozoic exhibit flashy hydrographs and a low baseflow component compared to younger catchments on volcanic rock.They also found that catchments underlain by younger volcanic rock formed in the Quaternary have a greater groundwater component compared to those on older volcanic rock formed in the Tertiary.Lohse and Dietrich (2005) found that vertical saturated hydraulic conductivity declines as volcanic soils in Hawaii develop with time, due to the formation of secondary clay mineral layers in those soil profiles.This evolution impedes rates of vertical percolation of water and increases shallow lateral subsurface flow.Jefferson et al. (2010) studied landscape and hydrological signatures in basalt catchments in the Oregon Cascades and reported that springs were absent in older catchments, which exhibit well-developed drainage networks and flashy hydrographs.They hypothesized that the differences in hydrological response between young and old catchments was the result of changes in flow pathways, from vertical recharge and groundwater-dominated flow to shallow subsurface lateral flow with sufficient erosive power to create channel networks through landscape incision.
Prior work has suggested that climate controls variations in hydrologic regimes in two ways.One is direct control on the water balance by function of available water and energy (e.g., Budyko, 1974;Milly, 1994;Troch et al., 2009;Sivapalan et al., 2011).The other is indirect control on the hydrologic regime as one of four independent factors (climate, bedrock weatherability, tectonic uplift and time) that control the evolution of landscape properties (landforms, soils, and vegetation;Troch et al., 2015).For example, Gentine et al. (2012) found that catchment-scale rooting depth was strongly correlated with the seasonal distribution of water and energy availability in 431 catchments in the United States.Wang and Wu (2013) found a strong inverse correlation between perennial stream density and climate aridity index in 185 catchments across the United States.Both studies suggested that these catchments evolved by adapting to the long-term climate conditions.
From this overview it is clear that both climate and time influence spatial variations of hydrological response.In this study we investigated the degree of control of climate and time on the differences in hydrological response between catchments that have different internal organizations.More specifically, we asked the following questions: (1) how do landscape features and hydrological response differ with catchment age and climate?(2) What is the dominant control of catchment coevolution in basaltic landscapes?(3) What are the possible mechanisms related to the historical development of volcanic catchments?To address these questions we selected 14 catchments dominated by basaltic rock with different ages and from a range of climatic regions of Japan.The age of the catchments ranges from 0.225 to 82.2 Ma (million years).The aridity index (ratio of annual potential evaporation to annual precipitation) for these catchments ranges from 0.28 to 0.82.For each catchment we derived indices of landscape features (drainage density and slope-area relationship) and hydrological response (annual water balance, baseflow index, and flow-duration curves).We then used regression analysis to explore possible relationships between these features and catchment age and climate.Based on the observed relationships, we proposed a hypothesis that explains as many observations in volcanic catchments, both these and others, as possible.

Study catchments
The study area is located within the Japanese archipelago, where 28 % of the land is underlain by volcanic rock.The makeup of the volcanic landscape in Japan with respect to time since formation, or geological age, is as follows: 11 % was formed in the Quaternary, 12 % in the Neogene, 1 % in the Paleogene, and 4 % in the Cretaceous (Murata and Kano, 1995).The volcanic landscapes formed in the Quaternary and Neogene are mostly in northeastern and south-western Japan, and those formed in the Paleogene and Cretaceous are predominantly in central and western Japan.The majority of the basement rock of Japan was formed in the Jurassic to Paleogene at the margin of the Asian continent and moved toward the current position during the spreading of the Japan Sea basin in the Neogene (Taira, 2001).The volcanic rocks formed in the Cretaceous, which are prevalent in western Japan, were formed at the margin of the Asian continent and have been exposed to the atmosphere since their formation (Taira, 2001).The volcanic rocks of the Quaternary and Neogene age, on the other hand, were formed after the spreading of the Japan Sea basin and thus have lain in their current position since their formation.
The catchments selected for the study are among the catchments listed in the Japanese Dam Database (National Institute for Land and Infrastructure Management, Japan, 2015).Out of 130 catchments listed in the database, we selected 14 that contain more than 50 % volcanic rock coverage (Table 1).Basalt is the dominant rock type in all catchments, while rhyolite is found in five of the catchments (GOS, YUD, AIM, IKR, and KWM), but at less than 15 % of the volcanic rock area.The catchment areas range from 30.9 to 635.0 km 2 .Daily streamflow observed at each dam was obtained from the Japanese Dam Database.The streamflow observations are inflow to the dams and no major dams are located upstream; therefore, the streamflows examined are not affected by human regulations.
The catchment ages were assigned from the average of the ages of different volcanic rocks, weighted by area of coverage, which was obtained from the Seamless Geological Map of Japan (Geological Survey of Japan, 2012).The map is the product of the merged 1 : 200 000 geologic map sheets covering the entire country of Japan.The age of the catchments ranged from 0.225 to 82.2 Ma.The youngest catchment (SNK) is located in central Japan (Fig. 1).The oldest catchment (HAZ) is covered with basalt formed in the Late Cretaceous (Fig. 2).The other 12 catchments were formed in the Neogene period (2.7-18.5 Ma) and are clustered in northeastern Japan, with the exception of the Shimouke catchment in southern Japan.
To delineate catchments, digital elevation models (DEMs) of 10 m×10 m were downloaded from the website of the Geographic Survey Institute, Japan (2015).The stream networks were derived from digitized geographical information system (GIS) maps from the National Land Numerical Information website (Ministry of Land, Infrastructure and Transportation, Japan, 2015).Tectonic uplift for each catchment was estimated from the map of uplift in the Quaternary with 250 m intervals (National Research Center for Disaster Prevention, Japan, 1973).We imposed the delineated catchments on the tectonic uplift map and determined the uplift from the dominant area of the catchment coverage on the map.

Correction of meteorological data
The precipitation in the raw data set accounts for orographic effects based on empirical relationships between precipitation and geographic features, such as elevation and distance from shoreline.It might not represent local orographic effects (e.g., higher rate of increase in precipitation on windward side of mountains), and does not correct for snow undercatch.Thus, we corrected the precipitation data as described in the Appendix.We assessed the validity of the raw and corrected data sets by means of the long-term water balance of the catchments.Given the mean of catchment-normalized annual precipitation for the raw data set P raw (mm yr −1 ) and precipitation-corrected data set P cor (mm yr −1 ), and streamflow Q (mm yr −1 ), we estimated the mean annual evapotranspiration for each data set, E raw and E cor (mm yr −1 ), as follows: To cross-check the acceptability of E raw and E cor , catchment-scale actual evapotranspiration, E a (mm yr −1 ), was directly estimated with the complementary relation method (e.g., Bouchet, 1963;Morton, 1978;Brutsaert and Stricker, 1979).In this study, the complementary relationship method proposed by Otsuki et al. (1984a) was used to estimate E a for each catchment.The method was used because it was rigorously validated with the long-term evapotranspiration in several experimental watersheds across Japan and does not require additional parameters (Otsuki et al., 1984b;Oue et al., 1992;Yoshida et al., 1994).
Thus, we have used three types of meteorological data sets to estimate the long-term water balance: -Raw data set: P raw , E raw ; -Corrected P data set: P cor , E cor ; The estimated data were validated against the Budyko hypothesis (Budyko, 1974).Budyko (1974) postulated an empirical relationship between the long-term evaporation index, E/P (E: evapotranspiration, P: precipitation), and the aridity index, PE/P (PE: potential evaporation) and argued that this relationship can be satisfactorily applied to most catchments.The large bias or discrepancy of estimated evaporation from the prediction of the Budyko curve was interpreted to derive from errors in the data set (Fig. 3).

Geomorphologic signatures
Two measures of landscape characteristics were used: drainage density and slope-area relationship.Drainage density is defined as the total length of streams divided by catchment area (km km −2 ).To obtain the stream length, we used the length of the vectors defining blue-colored streamlines in the digitized version of topographic maps (Geographic Survey Institute, Japan, 2015).The scale of the map is 1 : 25 000.Because our study catchments are all located in humid regions where precipitation is evenly distributed, we assumed that the blue lines of our digital data set represents perennial streams.
To understand the timescale of landscape evolution with climate and geological age, we compared the relationship between local slope and contributing area for all the catchments.The slope-area relationship is the log-log relationship between the local slope and the upslope contributing area for all grid cells within a catchment.The local slope was estimated for each grid cell of a 10 m DEM with the Dinfinity algorithm (Tarboton, 1997).The slope-area relationship can be used to discriminate between different geomorphologic regimes (Tarboton et al., 1992;Tucker and Bras, 1998;Bogaart and Troch, 2006).The slope-area relationship typically shows one or more inflection points, each of which can be interpreted as the threshold points where the dominant geomorphological process switches to other processes (Tucker and Bras, 1998;Bogaart and Troch, 2006).Hillslopes dominated by diffusional transport yield a positive slope-area relationship, i.e., slope increases with area, while hillslopes dominated by landslides or debris flows theoretically shows uniform slope gradient at a value equal to the threshold slope that induces the landslide (Tucker and Bras, 1998).Inflection points are found for a larger contributing area than those at the transition from hillslopes dominated by diffusional transport to the ones dominated by landslide/debris flow.The larger inflection points can be interpreted as the initiation of the hillslopes dominated by fluvial transport, which typically shows a negative slope-area relationship.

Hydrological signatures
Two signatures were used to represent the hydrological response of the study catchments.The first one is the baseflow index, the ratio of long-term baseflow component to total streamflow (e.g., Vogel and Kroll, 1992;Kroll et al., 2004).The higher the baseflow index the higher the groundwater contribution to catchment outflow.We separated the baseflow component from total streamflow as follows: where is the estimated baseflow, and is a low-pass filter parameter (Arnold and Allen, 1999;Eckhardt, 2005).To compare hydrological responses among catchments consistently, the filter parameter was set at 0.925 for all catchments (Sawicz et al., 2011).
To quantify an index of flow variability, the slope of the flow-duration curve was used as the second signature.The flow-duration curve is the cumulative probability distribution of daily streamflow.The slope of the flow-duration curve is calculated between the 33rd and 66th streamflow percentiles, because this range represents a relatively linear part of the flow-duration curve on semi-log scale (Yadav et al., 2007;Zhang et al., 2008).This is calculated from where Q 66 and Q 33 are the streamflows that correspond to 66 % and 33 % exceedance probability, respectively.A high slope value indicates a variable flow regime, while a low slope value means a more damped response.A damped response can arise as a result of a combination of persistent (widespread and year-round) precipitation, dominance of groundwater contribution to streamflow, or a combination of the two.

Initial data evaluation
To evaluate data consistency with long-term water partitioning, we plotted evaporation index versus aridity index for all catchments and compared the plotted values against the Budyko curve (Fig. 3).The plotted data points for the raw data set are all well below the Budyko curve.Catchmentscale actual evapotranspiration, E a , estimated from the complementary relation method, on the other hand, showed good agreement with the Budyko curve, suggesting realistic representation of long-term evapotranspiration in each of the catchment.Errors in the estimates of the raw data set could be caused by errors in discharge observations, trans-boundary groundwater flux, and catchment-averaged precipitation.The discharge, or dam inflow, was validated with the dam outflow, and the bias in the total inflow to total outflow was less than 2 % over 20 years.We thus assume that the dam inflow data are reliable.Another potential error is trans-boundary groundwater flux.Quantifying this flux is very difficult; however, by focusing on relatively large catchments (> 30 km 2 ), we argue that the influence of such a flux is negligible.It is more likely that the errors in annual evaporation estimated from the water balance are caused by underestimation of areal precipitation from the gridded data.Indeed, the corrected precipitation data set showed a substantial improvement from the raw data set compared to the Budyko hypothesis although some data still plot below the Budyko curve, likely caused by additional precipitation underestimation due to snow undercatch.
The aridity indices calculated with P cor were used for subsequent analysis (Table 2).All the catchments lie in a energylimited region (i.e., PE/P < 1).Although no correlation between climate aridity indices and catchment ages was found (R = 0.501, p = 0.068), the oldest catchment lies in the driest region, while the youngest catchment lies in a humid region (Fig. 4).Thus, it is not possible to thoroughly investigate the controls of climate and time on catchment coevolution, and the regression of any signatures is substantially affected by the behavior of the youngest and oldest catchments.

Analysis by geological age: controls of time on catchment coevolution
All the landscape and hydrological indices obtained from the analysis are presented in Table 3. Baseflow index, drainage density, and slope of the flow-duration curve were all strongly correlated with catchment age (Figs. 5 and 6). Figure 5 shows a significant correlation between the baseflow index and catchment age (R = −0.756,p = 0.0017), and be- The oldest catchment (HAZ) had a baseflow index of 0.533, while in the youngest catchment (SNK), the baseflow index of 0.835 was the highest of the study catchments, suggesting groundwater-dominated flow regime.The snow-dominated catchments tended to have a larger baseflow index, because we used a low-pass filter to exclude highfrequency signals from daily hydrographs to separate baseflow from total runoff.Two of our catchments (SIM and HAZ) are not snow affected, whereas the others are snow dominated as shown in Figs.A3 and A4.Plot for the catch-ments that receive little snow are negatively deviated from the regression line (Fig. 5).
The decrease in baseflow index with age was supported by the slope of the flow-duration curves (Fig. 6).The slopes of the flow-duration curve were positively correlated with catchment age (R = 0.74, p = 0.002), i.e., as catchments age, the slope of the flow-duration curve increases, indicating flashier runoff.
The flow-duration curves for all the study catchments are presented in the inset of Fig. 7, in which each line represents the data from a certain range in aridity index.It is clear that the aridity index controls the total amount of generated streamflow, i.e., the higher the aridity index is, the less streamflow generated by the catchment.Figure 7 depicts the flow-duration curves for the catchments, whose aridity index are in the range between 0.35 and 0.45.We used this range because it comprises the diversified catchments in terms of age (from 0.225 to 11.7 Ma).The slopes of the flow-duration curves of relatively young catchments, less than 4 Ma, were noticeably less than those of the older catchments.This implies that flow paths in older catchments with a similar aridity index have changed from deep groundwater to surface or near-surface flow paths.
It should be noted that there is no statistically significant log-linear fit between age and drainage density (p = 0.28) or age and slope of the flow-duration curve (p = 0.30) without the youngest catchment.We will later discuss the uniqueness of the youngest catchments and how it influences the interpretation of our results.

Analysis by climate: climate controls on catchment coevolution
The drainage density and baseflow index shows no significant correlation with aridity index (Fig. 8), although there are decreasing trends in both.
To highlight the seasonality of the intra-annual distribution of precipitation in each catchment, we calculated the seasonal index (SI) (Walsh and Lawler, 1981).
where P is the mean annual precipitation and X m is the mean precipitation of month m.The maximum value of SI is 1.83 if all the precipitation fell in a single month, while the minimum is 0 where the precipitation is perfectly evenly distributed throughout the year.The calculated SI of our catchments ranged from 0.182 to 0.506, which means they are classified as "very equable", "equable but with a definite wetter season" or "rather seasonal with a short drier season" (Walsh and Lawler, 1981).We therefore assumed the season- ality of precipitation of our catchments does not influence the results.

Geomorphological signatures of landscape evolution
The rate of landscape evolution is a function not only of time but also of climate, geology, and tectonics that drive the evolution process (Troch et al., 2015).To understand the timescale and governing processes of landscape evolution, an estimation of tectonic uplift is essential, especially in dynamic landscapes such as those of Japan.Figure 9 shows a strong correlation between the maximum slope captured in the mean slope-area relationship for the study catchments and the rate of tectonic uplift in the Quaternary period in Japan (National Research Center for Disaster Prevention, Japan, 1973).Numerical experiments with landscape evolution models also suggested that higher tectonic uplift rate yields steeper slopes (e.g., Yetemen, 2015).We thus assumed that the effect of the difference in tectonic uplift rate can be eliminated by comparing the catchments with similar maximum slope.
We then performed slope-area analysis (Tarboton, 1997;Jefferson et al., 2010) to infer the evolution pathways based on the current evidence.To quantify the activity level of climate on the landscape evolution, we used effective precipitation P eff = P cor − E a , because Jefferson et al. (2014) showed that mean annual precipitation has a strong link with erosion rate in basaltic islands and is the major driving factor for landscape evolution.
The mean slope-area relationship for three catchments with the smallest range of uplift rates, less than 250 m in the Quaternary, are presented in Fig. 10: SNK (0.225 Ma), ASE (3.52 Ma), and GOS (6.39 Ma).The mean slope-area relationship for the oldest catchment HAZ (82.2 Ma) is also shown for comparison.Note that the maximum slopes of the mean slope-area relationship were similar in SNK, ASE, and GOS as were the effective precipitation values; the mean slope-area relationship for HAZ differed in both respects.The smaller inflection points occurred at the contributing area value of approximately 300 m 2 in all four catchments in Fig. 10.The mean slope-area relationship in the younger two catchments (SNK and ASE) was relatively straight from the inflection points down to the outlet of the catchments, whereas a larger inflection point, which can be interpreted as the point of initiation of a fluvial channel, was found in the GOS catchment at a contributing area value of 4.4 × 10 4 m 2 (green dashed line in Fig. 10), suggesting the development of fluvial channels had occurred between 3 and 6 Ma since formation.Although the activity level of climate on the landscape evolution, i.e., effective precipitation, of HAZ is the lowest among the catchments in Fig. 10, the mean slopearea relationship for HAZ exhibited the gentlest slope among them where the contributing area is larger than 10 4 m 2 .Figure 11 depicts the mean slope-area relationship of catchments with uplift rate of 750 m in the Quaternary: TMG (2.78 Ma), NAR (8.43 Ma), ISB (11.8 Ma), YUD (18.5 Ma) and HAZ (82.2 Ma).The maximum slopes were similar in all of these catchments, while the effective precipitation value of HAZ was substantially less than those of the other catchments.The smaller inflection points were found at a contributing area value of approximately 300 m 2 , which is similar to the catchments with the lowest uplift rate (Fig. 10), whereas clear larger inflection points were found in catchments except HAZ (Fig. 11).The larger inflection point in the youngest catchment in Fig. 11, TMG, at a contributing area of 2.1 × 10 6 m 2 (red dashed line) suggested that fluvial erosion had emerged earlier (less than 3 Ma of age) than in the catchments with the lowest uplift rate.With the only exception of ISB, the inflection points chronologically shifted from a lower position of the catchments (larger contributing areas) towards upper areas (smaller contributing areas); the inflection points emerged at 4.5 × 10 5 m 2 for NAR (dashed light blue line) and at 4.2 × 10 4 m 2 for YUD (dashed deep blue line), respectively.Note that the contributing areas of the inflection points do not necessarily coincide with those of the averaged contributing areas for the channel heads of digitized stream networks: these were 5.2×10 6 m 2 for TMG, 9.2×10 5 m 2 for NAR, 5.3×10 5 m 2 for ISB, 7.7×10 5 m 2 for YUD, and 7.0 × 10 5 m 2 for HAZ.
The shift of the inflection points towards an upper position of the catchments implies the extension of channels dominated by fluvial erosion.Although there was no noticeable inflection point in the relationship for HAZ, the lowest slope in HAZ suggested that the catchment was substantially eroded.The inflection point for ISB was observed at a contributing area of 2.9 × 10 6 m 2 , which was the largest value for the inflection point in all study catchments.The exceptional catchment, ISB, may invalidate our interpretation that the channels dominated by fluvial erosion extended chronologically.
Last, we compared the slope-area relationships for the catchments with tectonic uplift of 1000 m in the Quaternary: KWM (3.69 Ma), IKR (10.7 Ma), and AIM (11.7 Ma) (Fig. 12).Again, the slope-area relationship for the oldest catchment, HAZ, was added to the figure for reference.Clear inflection points were observed at contributing area values of 2.8 × 10 5 m 2 for KWM, 1.5 × 10 5 m 2 for IKR, and 3.7 × 10 5 m 2 for AIM.
Interpretation of slope-area relationship is not straightforward when tectonic uplift rates are incorporated.The inflection points for the catchments with the highest uplift rate do not corroborate the findings in the catchments with modest uplift rates.One possible explanation for this is that the steep slopes in the high uplift catchments caused frequent landslides or debris flows, and subsequent rate of landscape evo-

Relation between catchment age and hydrologic response
This study was predicated on the assumption that there is a connection between geological age, climatic conditions, and hydrological response.
Our results show that there are significant negative correlations between geologic age, baseflow index, and drainage density.Variability in the baseflow index and flow-duration curves suggest that younger catchments tend to have a larger groundwater flow component, while older catchments exhibit more flashy runoff.The decrease in baseflow in older catchments suggests that the major flow pathways have changed over time from deep groundwater aquifers to shallow subsurface flows.
The observed decrease of baseflow index with catchment age is consistent with the relationship observed in the Oregon Cascades in the western United States (Fig. 13, Jefferson et al., 2010).To compare the baseflow index values of our study with those presented by Jefferson et al. (2010), we separated the baseflow from total streamflow by using the Hysep procedure (Sloto and Crouse, 1996) presented in Jefferson et al. (2010).The combined data set of our catchments and with those of Jefferson et al. (2010) showed a steeper decline of the baseflow index with catchment age and more significant correlation between them (R = −0.879,p < 0.001) than those for the catchments in Japan alone.Jefferson et al. (2010) concluded that this decrease is the result of decreasing groundwater recharge due to a development of clay layers in the near subsurface.The emergent negative trend of baseflow index for the integrated data sets suggests that the clay layers have continued to develop over longer periods than studied by Jefferson et al. (2010).
The consistency in the rate of decline of baseflow index with time since formation suggests the similarity in catchment evolution rate in both studies.Troch et al. (2015) proposed the concept of hydrological age in catchment coevolution; i.e., "the age of a catchment is not only a function of time since formation, but also a function of how 'active' the climate, geology, and tectonics are to drive coevolution".The aridity index for the catchments in the Oregon Cascades ranged from 0.5 to 1.0 (Zomer et al., 2008), and annual precipitation from 1800 to 2700 mm (Jefferson et al., 2010).These conditions are in the same range of our study (aridity index: 0.281-0.826;annual precipitation: 1513-3305 mm), which allows the two studies to be compared.Studies addressing the tectonic uplift in the western Cascades suggests that the portions of the landscape are uplifted during the early Pliocene (5.3-3.6 Ma; Priest, 1990); and the cross section of the western Cascades suggested that the rate of tectonic uplift is on the order of 700 m in the period (Conrey et al., 2002), which corresponds to 350-500 m in the Quaternary.The rates of tectonic uplift of our catchments, ranging from 250 to 1000 m, suggests some of our study catchments are more dynamic than the catchments in the Oregon Cascades.However, the significant correlation between the catchment age and baseflow index for the combined data set implies that the differences in tectonic uplift rate do not exert significant influence on the rate of catchment coevolution.
No correlation was found between baseflow index and aridity index, not even in the United States (Fig. 14).These results indicate that the partitioning of total runoff into baseflow and quick flow was not dominantly controlled by climate, at least in the narrow range of aridity index of our catchments, suggesting that the catchment age is a strong control on hydrologic response, which is in agreement with the findings of Jefferson et al. (2010).Significant correlations between catchment age and hydrological responses allow us to enhance the predictive capabilities in ungaged environments (e.g., Sivapalan et al., 2011;Troch et al., 2013;Harman and Troch, 2014).

Reading the landscape
In addition to studying the pathways in terms of individual catchment forming factors independently, Troch et al. (2015) proposed that the landscapes could be read to infer their evolutional pathways from the current geomorphological evidence.
In the current set of catchments, we found a strong negative correlation between the drainage density and catchment age.This finding contradicts that of Jefferson et al. (2010), who found a positive correlation between drainage density and catchment age.The disparity may arise from the differences in the range of catchment ages: 0.225-82.2Ma for our catchments, versus 0.017-15.5 Ma for the Oregon Cas- cades catchments of Jefferson et al. (2010).From the global data set of visible dissection on volcanic islands, Jefferson et al. (2014) observed that the visible dissection begins after about 0.5 Ma, that the dissection rate is much faster in humid landscapes, and that there is a transition from weak dissection to substantial dissection between 0.5 and 2 Ma, even in arid environments.Conversely, Suzuki (1969) observed that the dissection rate in volcanic landscapes in Japan declined markedly after 0.1 Ma after initiation.These data suggest that the increase in drainage density in the Oregon Cascades with age arose as a result of intense surface dissection, whereas the catchments of our study had possibly reached maturity in regard to dissection.We thus hypothesize that the decline in drainage density with catchment age is caused by a mechanism other than the one described in Jefferson et al. (2010).Wang and Wu (2013) examined 185 catchments in the United States with aridity index ranging from 0.26 to 5.50.They found a systematic decrease in perennial stream density with increasing aridity index (Fig. 15).They hypothesized that the aridity index is the first-order indicator of the partitioning of precipi-  3).The correlation efficient between them was 0.74 without the study catchments (Wang and Wu, 2013) and 0.72 for the entire data set.tation into baseflow, and found a strong correlation between perennial stream density and the baseflow coefficient, Q b /P , which follows the complementary Budyko curve (Fig. 16).We assume that our digital stream data set represents perennial streams because all the study catchments are located in a humid climate; hence, the drainage density can be compared with the perennial stream density of Wang and Wu (2013).However, the perennial stream density-aridity index curve derived by Wang and Wu (2013) does not explain the decreasing trend in drainage density of our catchments (Fig. 15, inset) even though the plotted data of our catchments cluster around the fitted curve (Fig. 15).Based on the general trend that catchments with lower baseflow coefficients tend to have lower drainage densities, we hypothesize the decline in drainage density with age therefore suggests that the proportion of precipitation that recharges groundwater aquifers declines with age.
However, plotting the drainage density of our study together with those observed in the Oregon Cascades leads to a different interpretation of the data (Fig. 17).First, the youngest catchment in our data set, SNK, has a much higher drainage density (1.08 km km −2 ) than those of similar age in the Oregon Cascades.Second, if we remove SNK, no significant correlation between drainage density and catchment age can be found for the study catchments, although they still exhibit a negative trend (R = −0.32,p = 0.28).From lines (b) and (c) in Fig. 17, we can formulate an alternative hypothesis: the volcanic catchments reach maturity with respect to surface dissection at approximately 2 Ma of development, and the drainage density does not undergo significant change after that.
The mapping standards of the United States and Japan are both not available; hence, the difference in the mapping standards might be a reason for the anomalous drainage density of SNK.However, it is not convincing given the drainage densities of the catchments older than 2 Ma are similar.We thus assume that the difference in the mapping standards are negligible for comparison although there still remains uncertainty.
To test the alternative hypotheses, we investigated why the SNK catchment has such a high drainage density.The digital elevation model of SNK reveals distinct landscape differences between the upper and lower part of the catchment; i.e., less dissected plateau in the upper and highly dissected valleys in the lower catchments (Fig. 18).The longitudinal profile of the longest stream in SNK shows a point of convexity, or a knick point, defined as a steep reach between two low-gradient reaches (Fig. 19).The knick point is interpreted as a hallmark of the transient state of a river system caused by changes in the lithology across which the river flows, changes in the rate of tectonic uplift or spatial variability in climate forces (Anderson and Anderson, 2010).A detailed investigation of the cause of the knick point is beyond the scope of this paper; however, the basalt in the upper catchment is younger than that in the lower catchment.By dividing the upper and lower catchment along the 1700 m contour line, which is where the stream gradient is steepest (Fig. 19), we obtain the weighted average of the basalt age as 72 ka for the upper and 380 ka for the lower catchment areas and the drainage densities for those areas as 0.51 and 1.35 km km −2 , respectively.The plotted drainage density for the upper part of SNK is comparable with those observed in the Oregon Cascades; however, the high drainage density for the lower part of the catchment does not fit with either relationship (squares associated with SNK in Fig. 17).One possible explanation for the high drainage density in the lower part of the catchment is the acidity of stream water.The pH of the stream water in this catchment is as low as 2-3 due to its high sulfuric acid content (Ossaka et al., 1980).Acid precipitation significantly promotes weathering and soil formation by the input of rainfall and H + (Huang et al., 2013).Acidic stream water is neutralized by the chemical weathering of primary silicate minerals in humid forested watersheds (Johnson et al., 1981).The lower part of the SNK catchment is known for numerous hot springs of high acidity, which may have dissected at a much faster rate than neutral water.
If the extremely high drainage density in the lower part of SNK has been caused by weathering due to acidic streamflow, it supports the latter hypothesis that the drainage density reaches maturity and does not change over time, but it does not reject the former on that the drainage density declines with age.
From the above analysis, the drainage densities reached a steady state or slightly decreasing trend 2 Ma after catchment formation.However, this contradicts the slope-area relationship, which indicates that the channels dominated by fluvial erosion emerged between 2 and 6 Ma since their formation, depending on the tectonic uplift rate, and extended towards the upper portion of the catchments until they become substantially eroded gentle hillslopes.It should also be noted that the contributing areas of the inflection points in mean slope-area relationship are generally smaller than those of the channel heads defined by the digitized stream lines.In other words, the signatures of fluvial erosion emerge in the upper portion of the catchments.The discrepancy be- tween the inflection points and the channel heads implies that some portion of the fluvial channels are ephemeral rather than perennial.The extension of ephemeral streams in mature catchments supports the argument that the major flow pathways have changed over time from deep groundwater to shallow subsurface flow.

Hypothesis on mechanisms of coevolution in volcanic catchments
Table 4 summarizes all the hydrological and geomorphological signatures of the study catchments.From this summary, we propose a hypothesis that explains as many observations of catchment evolution in volcanic landscapes as possible.
In the early stages of development in volcanic landscapes, most of the infiltrated water flows vertically due to the high permeability of the young basalt, which typically contains many cracks and fissures (Lohse and Dietrich, 2005).In catchments younger than 2 Ma, subsurface flow is dominant to the extent that signatures of fluvial erosion in the slopearea relationship are weak (Fig. 20a).As time progresses, subsurface impermeable layers develop due to chemical weathering and mineral precipitation, and the major flow pathways change from vertical flow to shallow subsurface flow The slope-area relationship suggests that from 2 to 6 Ma of age, fluvial channels emerge in the lower part of catchments (Fig. 20b).As chemical weathering further continues (6-20 Ma), recharge to deep aquifers decreases, which would result in a shorter transit time of water and flashier hydrological responses (Fig. 20c).Increased pore pressure in steep slopes causes landslides.Mass wasting from landslides gradually decreases the local elevation until a threshold slope is reached below which landslides cannot be triggered.The lower gradient of the channel probably cannot support a high perennial stream density, leaving most of the stream network as ephemeral.The channels dominated by fluvial erosion extended towards the upper portion of the catchments until they become substantially eroded gentle hillslopes (Fig. 20d).As chemical weathering further continues, recharge to deep aquifers decreases, and further disconnection of the channel network from aquifers might result in a decrease in the extent of the perennial stream network.
It should be noted that the interconnected nature of landscape evolution and hydrologic coevolution make it difficult to determine causality.The limitation of our data set is also a confounding influence, as the oldest catchment lies in the driest region and the youngest catchment lies in a humid region.Separating the effects of climate and base rock geology is therefore difficult if we are restricted to empirical analysis.One way to test the hypothesis in such an environment T. Yoshida and P. A. Troch: Coevolution of volcanic catchments in Japan is to use process-based hydrological models in an attempt to derive responses that are not observable with direct measurement.For instance, Troch et al. (2013) decoupled the landscape and climate properties using a hydrological model developed by Carrillo et al. (2011) and illustrated how catchment characteristics varied with climate gradients.

Conclusions
We examined 14 volcanic catchments of different ages, ranging from 0.225 to 82.2 Ma.We derived indices describing landscape features (drainage density and slope-area relationship) and hydrological responses (baseflow index, annual water balance, and slope of the flow-duration curve).The age of the volcanic rock was significantly related to the partitioning of deep groundwater recharge and shallow subsurface flow.Younger catchments tend to be associated with a larger component of groundwater flow, older ones with more flashy hydrographs.The decrease in baseflow with age suggests that the major flow pathways have changed from deep groundwater aquifers to shallow subsurface flows.The drainage densities of the catchments also decreased with catchment age.While aridity index does not explain the trend in drainage density given a quite narrow range of aridity index of our data set, catchment age has a more significant impact on the changes in drainage density; alternatively, drainage density has stayed constant in mature catchments.Our results suggests two hypotheses to be tested in future studies on the evolution of drainage density in matured catchments.One is that as catchments further evolve, hydrologically active channels retreat as less recharge leads to lower average aquifer levels and less baseflow; the other is that it does not significantly change after catchments reached maturity in terms of surface dissection.

Figure 3 .
Figure 3. Raw, corrected P , and estimated E data sets plotted with Budyko curve.

Figure 4 .Figure 5 .
Figure 4. Relation between aridity index and age for the study catchments.

Figure 8 .
Figure 8. Climate controls on baseflow index and drainage density.The relationships with aridity index were not significant for both drainage density (R = −0.451,p = 0.105) and baseflow index (R = −0.372,p = 0.190).

Figure 9 .
Figure 9. Tectonic uplift in the Quaternary period and maximum slope of mean slope-area relationship.

Figure 14 .Figure 15 .
Figure 14.Relation between baseflow index and aridity index in Japan and MOPEX (Model Parameter Estimation Experiment) catchments in the US aridity index does not explain variability in baseflow index.

Figure 16 .
Figure 16.Perennial stream density and baseflow coefficient (Q b /P ).The baseflow, Q b , in their study was estimated by using the same equation and parameter value as we used in Eq. (3).The correlation efficient between them was 0.74 without the study catchments(Wang and Wu, 2013) and 0.72 for the entire data set.

Figure 17 .
Figure 17.Drainage density with catchment age.Line (a) is the regression line for the study catchments (R = −0.603,p = 0.0222).Line (b) is the mean drainage density for the catchments older than 2 Ma in Japan and the Oregon Cascades (0.873 km km −2 ).Line (c) is the regressed line for younger than 2 Ma in the Oregon Cascades.The drainage densities and basalt ages of the upper and lower part of SNK are plotted as squares (the upper part: drainage density is 0.51 km km −2 and basalt age is 72 ka; the lower part: 1.35 km km −2 and 380 ka).

Figure 20 .
Figure 20.Conceptual model describing coevolution in volcanic catchments.In the early stages of evolution, deep groundwater flow in the bedrock (deep blue arrows) is the dominant flow pathways.As time progresses, subsurface impermeable layers develop (deepening of color in soil layer), and the major flow pathways change from vertical flow to shallow subsurface flow (light blue arrows).

Table 1 .
List of study catchments

Table 4 .
Summary of hydrological and geomorphological analysis.