Evaluations of the downward velocity of soil water movement in the unsaturated zone in a groundwater recharge area using δ 18 O tracer: the Kumamoto region, southern Japan

Water and substances from the surface infiltrate the unsaturated zone before reaching ground - water. Yet, little study has been done on the unsaturated zone due to the difficulty of sampling. A lot of studies have been carried out on the top soil down to a depth of one metre and on shal- low aquifers because they are easily accessible for sampling. The unsaturated zone of the Kumamoto region recharge areas is important due to concerns about groundwater pollution from agriculture. The aim of this study was to estimate the downward velocity of soil water movement through the unsaturated zone and the recharge rate using δ18O as a tracer. Five sampling sites were selected and a core was taken from each site. The cores were cut into 0.1 m pieces and soil water was extracted from each to analyze for δD and the δ18O content. Average δD and δ18O compositions of soil water were similar to the isotopic compositions of summer precipita tion. Annual average recharge rate and the downward velocity of soil water in each site were estimated by fitting a vertical δ18O profile pattern to a precipitation δ18O time series as a theo retical water displacement flow model for recharge. An estimated annual average recharge rate in the recharge area ranged from 745 to 1058 mm/yr with the annual average downward velocity of 1.37 to 2.34 m/yr. Based on the estimated downward velocity, the infiltration time for soil water to reach the aquifer was determined as ranging from 9 to 24 years, which corresponds with previous groundwater age estimations presented in an earlier published study on the same area. It was assumed that contaminants will reach aquifers in 9 to 25 years if the effects of dif fusion and microbiological reaction are not taken into account. clarified. DELIN et al. (2017) applied rhodamine and Br tracers to a crude oil contaminated site to detect the influence of oil on the water infiltration rate in the unsaturated zone. These studies show that the plurality of flow processes such as prefer ential flow and matrix flow complicates water infiltration and contaminant behaviour through unsaturated media. Understand ing the behaviour of water in the unsaturated zone is important for making decisions about groundwater usage and formulating mitigation measures against groundwater contamination.


INTRODUCTION
The excessive application of chemical fertilizers and manures to agricultural lands has increased the influx of substances into the ground, which eventually reach and contaminate groundwater (SPALDING et al., 1982;UMEZAWA et al., 2008;ZHANG et al., 2014). In order to identify the type of contamination, determine its extent, trace the source and extract valuable information about the groundwater such as its age, the water is sampled to measure stable isotope ratios, major and minor ion concentrations as well as trace metals. Water and substances reaching aquifers serve as a guide to determining these groundwater characteristics (DEMLIE et al., 2007;CLOUTIER et al., 2008;HOSONO et al., 2014). Usually in groundwater studies, a direct link is created between the characteristics of the groundwater and ions, isotopes and metals on the surface without taking the unsaturated zone into consideration even though substances pass through this zone before reaching groundwater. Although the unsaturated zone has not been as actively studied as groundwater, studies have been carried out over the years that give some insight into the flow processes and solute behaviour through it (BEVEN et al., 1982;GUYMON, 1994;LIN et al., 2006). General unsaturated flow processes have also been studied. Field experiments using Clas a tracer were conducted on a mine waste pile (5 m high) to evaluate the infiltration rate and detect the dominant flow system in the pile (NICHOL et al., 2005). NIMMO et al. (2007) modelled unsaturated zone flow using soil science and hydrologic data. The

Evaluations of the downward velocity of soil water movement in the unsaturated zone in a groundwater recharge area using δ 18 O tracer: the Kumamoto region, southern Japan
Azusa Okumura 1 *, Takahiro Hosono 1 , Dennis Boateng 2 and Jun Shimada 1 1 Kumamoto University, Priority Organization for Innovation and Excellence, 2-39-1 Kurokami, Kumamoto 860-8555, Japan; (*corresponding author: osnx.2157@gmail.com) 2 Kumamoto University, Graduate School of Science and Technology, 2-39-1 Kurokami, Kumamoto 860-8555, Japan doi: 10.4154/gc.2018.09 Abstract Water and substances from the surface infiltrate the unsaturated zone before reaching groundwater. Yet, little study has been done on the unsaturated zone due to the difficulty of sampling. A lot of studies have been carried out on the top soil down to a depth of one metre and on shallow aquifers because they are easily accessible for sampling. The unsaturated zone of the Kumamoto region recharge areas is important due to concerns about groundwater pollution from agriculture. The aim of this study was to estimate the downward velocity of soil water movement through the unsaturated zone and the recharge rate using δ18O as a tracer. Five sampling sites were selected and a core was taken from each site. The cores were cut into 0.1 m pieces and soil water was extracted from each to analyze for δD and the δ18O content. Average δD and δ18O compositions of soil water were similar to the isotopic compositions of summer precipitation. Annual average recharge rate and the downward velocity of soil water in each site were estimated by fitting a vertical δ18O profile pattern to a precipitation δ18O time series as a theoretical water displacement flow model for recharge. An estimated annual average recharge rate in the recharge area ranged from 745 to 1058 mm/yr with the annual average downward velocity of 1.37 to 2.34 m/yr. Based on the estimated downward velocity, the infiltration time for soil water to reach the aquifer was determined as ranging from 9 to 24 years, which corresponds with previous groundwater age estimations presented in an earlier published study on the same area. It was assumed that contaminants will reach aquifers in 9 to 25 years if the effects of diffusion and microbiological reaction are not taken into account. study focused on preferential flow and the general transport time was also clarified. DELIN et al. (2017) applied rhodamine and Br tracers to a crude oil contaminated site to detect the influence of oil on the water infiltration rate in the unsaturated zone. These studies show that the plurality of flow processes such as preferential flow and matrix flow complicates water infiltration and contaminant behaviour through unsaturated media. Understanding the behaviour of water in the unsaturated zone is important for making decisions about groundwater usage and formulating mitigation measures against groundwater contamination.
One of the earliest soil water studies using water stable isotopes was performed by ZIMMERMANN et al. (1966). In this study, D 2 O was sprayed on the study area as a tracer. Soil down to 1.2 m depth was sampled 4, 26, 82, 152 and 214 days after spraying for analysis. A time series hydrogen stable isotope profile of the soil water revealed that piston flow, which is the dominant flow mechanism in near saturated formation, was excellent. Several studies have also been carried out with oxygen and hydrogen stable isotopes as tracers in dry and semi-arid areas (BARNERS and ALLISON, 1983;De VRIES et al., 2000). In agriculture, soil at 1 m depth has been actively studied. However, studies on soil at far greater depths have rarely been undertaken (BHARATI et al., 2002).
Many different isotope tracers are used around the world based on the purpose and intent of a study. One of these tracers is tritium, a radioactive isotope of hydrogen. Tritium used to be a good tracer for tracking water movement through the unsatu-rated zone and determining groundwater age because of it abundance in the atmosphere. In the arid climate of Arizona, it was used to determine the mixing degree between shallow and deep groundwater (HARRIS, 2000). In the humid climate of Japan, it was used as a tracer to determine the recharge rate and infiltration time (KAYANE et al., 1980;SHIMADA, 1988). However, a few studies pointed out that tritium cannot be solely relied on because its concentration keeps on dropping year after year in the atmosphere (TSUJIMURA and TANAKA, 1998). Consequently, in later studies, oxygen and hydrogen isotopes were used together instead as tracers to determine recharge rates especially in Japan (YABUSAKI et al., 2011;KUDO et al., 2016). Both recharge rate and unsaturated zone infiltration velocity are important variables that provide insight into contaminant movement through the unsaturated zone.
In this study, the study area as shown in Figure 1, the Kumamoto region, is almost 100% dependent on groundwater for drinking purposes. However, the nitrate concentration has begun to build up in the water, raising concerns about groundwater contamination (HOSONO et al., 2013). The high nitrate buildup occurs at high elevation areas which have lots of agricultural activities and are the recharge zones for Kumamoto groundwater. It is therefore necessary to understand the infiltration mechanism in the unsaturated zone of these areas. Recently, unsaturated zone soil samples (down to 5 m depth) were collected at the slope of the western foot of Mt. Aso and the amount of recharge was estimated using oxygen and hydrogen isotope ratios of the water extracted from the soil (KUDO et al., 2016).
Research on the deep unsaturated zone (up to 15 m in depth) in Kumamoto using isotope tracers is still limited. The aim of this study is to augment the efforts in unsaturated zone studies in the Kumamoto region by trying to clarify the infiltration mechanism through estimation of the recharge amount and downward velocity of soil water using oxygen isotope ratios (water) in the recharge zone.

STUDY AREA
The Kumamoto region is about 1041 km 2 in area and bordered by the Chikushi-mountains to the North, the Mashiki-mountains to the South and the slope of the western foot of Mt. Aso to the East (Figure 1). There are three main rivers in the region; the Kikuchi, Shira and Midori Rivers which all flow into the Ariake Sea. The hydrogeological basement of the Kumamoto region consists of Mesozoic metamorphic and igneous rocks. The basement rocks are overlain by Quaternary pyroclastic flow deposits with high permeability (MIYOSHI et al., 2009). These deposits also cover the three major plateaus in the region which are the Kikuchi and Takayubaru plateaus to the east and the Ueki plateau to the North (Figure 1). They have altitudes of approximately 100 m. The unsaturated zone of these plateaus has a thickness of about 30 m. New volcanic ash layers are widely distributed on the top soil. The soil profiles of cores taken from the various sites studied are shown in Figure 2. Soil types include sandy clay, clay sand, volcanic ash sand, volcanic ash clay, gravelly clay, nonwelded tuff and gravel. There are variations in soil type among the cores.
Groundwater in the Kumamoto region is recharged on the slope of Mt. Aso and around the plateaus, and flows from the northeastern part of the region to the southwest (TANIGUCHI et al., 2003). The aquifer formed by Quaternary pyroclastic flow deposits are divided into two main types by an aquiclude layer of lacustrine sediments; unconfined and confined aquifers (SHI- A R T I C L E I N P R E S S MADA, 2012). The land-cover is 31% farmland, 29% forest and 11% paddy. The plateaus, which are groundwater recharge areas are mainly used as farmland. Average annual precipitation and temperature of Kumamoto from 1985 to 2015 are 1986 mm and 16.9 ˚C, respectively.

SAMPLING AND METHODS
In this study, 5 sampling sites were selected on the plateaus in the region ( Figure 1). Each sampling site shows a different soil type, different fertilizer application method as well as treatment, different crops cultivated and different drilling times as shown in Table 1 and Figure 2. Core samples were taken from the unsaturated zone of each site using an excavator without adding water (ECO-3V, YBM, Japan).
The core was cut into 0.1 m samples and vacuum packed onsite to prevent evaporation of soil moisture before being taken to the laboratory for analysis. The 0.1 m fractions were split into two and each placed in a 100 cc sampling can. From one, soil water was centrifuged at pF 4.2 using centrifuge (MODEL SS-2200, SAKUMA, Japan) which spins for 2 hours. Although different for each depth, 2 to 20 ml of soil water was extracted. For the other sample, a three-phase meter that operates by Charles Boyle's law was used to measure the total volume of the liquid and solid phases after weighing (HORITA, 1985). The gas phase content (air within the pores) was obtained by subtracting the liquid and solid phase volumes from the total volume of the can.
Precipitation data used in this study was gathered at the Kumamoto University from Jan, 2005 to Jul, 2016. The data was acquired by sampling rainfall two to three times every month in a 2000 ml plastic bottle equipped with a funnel. A ping-pong ball was placed in the funnel to prevent evaporation. Samples were analyzed for hydrogen and oxygen isotope ratios. The results were converted to monthly data by evaluating the weighted average for each month.
Hydrogen and oxygen isotope ratios were measured by mass spectrometer (Delta-V, Thermo Fisher Scientific, USA). In the mass spectrometer, the gas achieved equilibrium with the oxygen and hydrogen isotopes of the water samples that were introduced.   Table 1. Details of sampling points. For the three-phase distribution, samples at 0.9 -1.0 m, 5.0 -5.1 m and 9.9 -10.0 m are displayed. Details of all the water content data of the cores are shown in Table A1.

A R T I C L E I N P R E S S
The hydrogen and oxygen isotopes are expressed as relative differences (δ values) from Standard Mean Ocean Water in per mil (‰). The accuracy of the measured value was within ± 0.50‰ for δD and within ± 0.05‰ for δ 18 O.

RESULTS
Plots of δD and δ 18 O for all precipitation and soil water samples are shown in Figure 3. The elevation effect on water isotopes was evaluated. The maximum elevation difference between the cores was calculated to be 70 m. This elevation effect on the hydrogen isotope in the study area is reported as -0.027‰ per metre (KAGABU et al., 2017). Hence, the maximum elevation effect of the sampling points is about 1.9‰. However, it can be seen from were divided into two groups depending on season; summer (April to September) and winter (October to March). The slopes of the regression lines of the summer and winter data are 8.5 (r 2 = 0.87, p < 0.0001, n = 70) and 8.0 (r 2 = 0.84, p < 0.0001, n = 65) respectively, which are almost equal to 8, the slope of the World Meteoric Water Line (CRAIG, 1961). However, the intercept of the regression line for summer and winter data are 12.8 and 19.2, respectively. The intercept for all soil water was 8.8 (r 2 = 0.83, p < 0.01, n = 566) which was rela-  A R T I C L E I N P R E S S tively close to the intercept for summer ( Figure 3). Moreover, the averages of δD and δ 18 O of soil water samples were -41.3‰ and -6.7‰, respectively which are similar to the isotopic composition for summer precipitation (δD = -45‰, δ 18 O = -7‰) compared with that of winter (δ 18 O = -5‰、δD = -32‰). Evaporation occurring in the soil is generally considered to be non-equilibrium evaporation, in which case the slope of a regression line in the soil water tends to be smaller ranging from 2 to 5 (BARNES & TURNER, 1998;CLARK & FRITZ, 1997). In this study, the slope for all soil water sample plots was 8.1 (r 2 = 0.83, p < 0.01, n = 566), which implies evaporation in the top soil was not considerable. Figure 4. shows the vertical profiles of δ 18 O in soil water. These profiles display several relatively high isotopic composition peaks, e.g., at depths of 1.2 m and 3.5 m for the S1 core. It is interpreted that the vertical change of δ 18 O is reflective of the seasonal isotopic composition of the recharging precipitation which is high in summer and low in winter. These peaks hold the seasonal information of precipitation.
The plots in Figure 4. show that it takes one year for new peaks to appear on the profile. It can be inferred from the S1 core profile that the infiltration rate in the top 3.5 m would be 2.2 m/ yr approximately. The next peak after 3.5 m is at 10 m depth. Therefore, the depth interval from peak to peak is observed to be increasing. Moreover, there is a gradual change in peak compared with the peak at 3.5 m depth. The profiles of the other cores show similar characteristic peaks. It is assumed that as soil water infiltrates deeper, the more dispersion may occur when the peaks are averaged. Therefore, there is a limit to understanding the infiltration rate from the peak position. Table 1. shows the liquid and solid phases at depths of 1 m, 5 m and 10 m for each core. In all cores, the solid phase varies between 15% and 27% at 1 m depth, while porosity ranges between 73 and 88%. According to SHIMOZU (1986), Mt. Aso volcanic ash soil porosity ranges from approximately 70 to 90%. Hence, the sampling points have the characteristics of volcanic ash soil. The S1 liquid at 10 m depth is small compared with other cores. On the other hand, in C3, the liquid phase occupies 96% of the pore volume at a depth of 10 m. In particular, C3 has a high liquid phase at all depths. Therefore, C3 soil is considered to have a relatively high water holding capacity. It is important to determine the recharge rate and downward velocity of soil water by incorporating a water content profile.

DISCUSSION
In order to estimate the average downward velocity of soil water from 0-15 m depth using stable isotope ratios, the isotope profile of soil water was compared with the model isotope profile of the monthly infiltration amount obtained by subtracting evapotranspiration from precipitation. In building the model profile based on precipitation, the time-series data of precipitation was converted into depth data. To do this, the Displacement Flow Model (DFM) proposed by Anderson and SEVEL (1974) was applied. DFM was used earlier by SHIMADA (1988) and KUDO et al. (2016) to understand the behaviour of water in the unsaturated zone. The model assumes a piston flow. Here, a piston flow is also assumed in an attempt to estimate the recharge rates and the downward velocity.
The methods used to calculate DFM are described below. Equation (1) which is the water balance equation was used to estimate Infiltration (I) into the ground.
Where P, R and E are amounts of precipitation, runoff and evapotranspiration, respectively. In estimating E, Equation (2) was used. This is a modified form of Thornthwaite's equation (THORNTHWAITE, 1948) used to handle wide temperature ranges. This is based on the principle that the amounts of precipitation and evapotranspiration are directly proportional to the atmospheric pressure (TAKAHASHI, 1979).
Where T is monthly temperature, P is monthly precipitation.
Here, R in equation (1) was replaced with αP by setting the direct runoff coefficient α, assuming that a constant rate of P directly flows out. α was varied from 0 to 0.3 at an interval of 0.02 with reference to Kudo et al. (2016), and the highest coefficient of determination between the actual profile and the modeled profile was selected. The fraction of monthly precipitation used for each depth was determined to successively represent infiltration from the surface which corresponds to the soil water content. Based on the fraction of the monthly precipitation for each depth, weighted averages of the oxygen isotope ratios in precipitation were estimated to create a model isotope vertical profile (TSU-JIMURA et al., 1994).
DFM was used to create a depth profile by applying the time variation of the tracer concentration. The downward velocity of soil water and recharge rate were estimated based on the depthtime information obtained when the profile was created using DFM. Monthly temperature and precipitation data for Kumamoto City collected by the Japan Meteorological Agency were used in the DFM. Time series data used in the DFM calculations are shown in Figure 5 and include precipitation, evapotranspiration, temperature and δ 18 O in precipitation.
The profile of δ 18 O in soil water and that of the modeled profile by DFM are shown in Figure 6. The coefficient of determination (r 2 ) between the measured δ 18 O and modeled δ 18 O by DFM was lower than 0.5 for all cores (black line vs gray line in Figure  6). This indicates that direct results by DFM could not trace the behaviour of measured δ 18 O profile successfully. This may be due to the fact that input precipitation δ 18 O data are based on a weighted average of several precipitation events and did not cover precisely all actual precipitation δ 18 O records. Moreover, the concept of the DFM model might be the other reason which does not account for dispersion phenomena. To solve this problem, variance due to these factors was considered by calculating moving averages of DFM values and the condition that resulted in the highest correlation between measured δ 18 O and modeled δ 18 O was employed (KUDO et al., 2016). Consequently, we found the best condition resulting in the highest correlation between measured δ 18 O and modeled δ 18 O (moving average DFM), as shown by the black and dotted lines in Figure 6. The coefficient of determination for all cores was (r 2 > 0.8, p <0.02), suggesting that our DFM approach satisfactorily reflects water behaviour at the studied sites.
The estimated result for each site is shown in Table 2. In S1, it was observed that the downward velocity was the highest of all the sites. This is because the soil has a low liquid phase around 10 m depth. Comparing C3 and S1 which have the same soil type, S1 has a higher recharge rate and velocity. This may be due to the A R T I C L E I N P R E S S difference in soil conditions. For C3 and S2, it can be seen that the downward velocity of C3 is relatively lower, even though the recharge rate is higher. The most likely reason is that, C3 soil has the highest water content (Table 1). Therefore, it is evident that the difference in the average recharge rate and downward velocity for each core is due to the difference in soil conditions such as soil type and water content at each core depth (Table 1 and Figure 2).

A R T I C L E I N P R E S S
Although there are regional differences, the annual average recharge rate in the recharge area of the plateau was estimated to range from 745 to 1058 mm/yr, and the annual average downward velocity from 1.37 to 2.34 m/yr. Furthermore, the average recharge rate and average downward velocity calculated from the estimation results in Table 2 for S1, C2 and C3 in the northeastern recharge area (Figure 1) are 893 mm/yr and 1.7 m/yr, respectively. These somehow match the recharge amount of 1137 mm/ yr and infiltration rate of 2.3 m/yr of the Kumamoto region estimated using tritium (SHIMADA & UENO, 2016). However, our estimated value was smaller than that of SHIMADA & UENO (2016). The difference may be due to regional heterogeneities.
In estimating the infiltration time for water to reach the aquifers, the unsaturated zone thickness was divided by the downward velocity as shown in Table 2. The unsaturated zone thickness was obtained by subtracting the groundwater level elevation in the lowest water level period at the sampling points from the ground surface elevation. The transportation time from the ground surface to the aquifer was estimated to range from 9 to 24 years. KAGABU et al. (2017) estimated the ages of groundwater in the recharge zone at Kikuchi plateau using 85 Kr age tracer to range from 15 to 25 years (n = 3). Although different sample points were selected in this study, the water transportation time from the surface to the aquifer estimated by DFM is equivalent to the recharge area's groundwater ages calculated by KAGABU et al., (2017). Therefore, it can be reliably concluded that soil water reaches the aquifer over a period of about 9 to 25 years in the main recharge area of the Kumamoto region.

CONCLUSION
The recharge rate and downward velocity in the unsaturated zone of the recharge area of the Kumamoto region were estimated for 5 soil cores using oxygen isotope ratios in soil water as a tracer. Consequently, the annual average recharge rate was estimated to range from 745 to 1058 mm/yr, and the annual average downward velocity from 1.37 to 2.34 m/yr. Differences in the estimated recharge rate and downward movement velocity are due to the variation in soil properties of each site. Water transportation time from the surface to the aquifers was estimated to be 9 to 24 years, corresponding well to groundwater ages in the same recharge area. With the assumption that nitrogen loaded on the surface will be transported with soil water, nitrate should leach to the aquifer in 9 to 24 years once fertilizers and manures are applied. Such information is primarily important for the management and preservation of groundwater resources. However, nitrate transportation and behaviour in the unsaturated zone should be affected by diffusion, microbial activities, etc. The results of this study therefore may not be able to fully account for the transportation of solutes through the unsaturated zone to the aquifers. To better understand nitrate infiltration characteristics, it is necessary to promote future studies on nitrate concentrations and isotope ratios to aid further understanding of soil water infiltration mechanisms.     Three phase distribution (%) S1 core depth (m) S1 core depth (m) water amount (mm)