INTRODUCTION

Tracked and wheeled vehicles used during extensive tillage in Tsumagoi （Gunma Prefecture, Japan） hillslope agricultural fields resulted traffic pan in unsaturated soils. Since densely-packed traffic pan may create conditions that restrict water percolation during the rainfall, hydraulic conductivity contrast between the top and traffic pan layers may lead to the development of the perched saturated zone, and downslope saturated flow ensues above traffic pan. The most baffling observation of Tsumagoi hillslope fields is that gullies and rills are developed mostly by extending their upper end uphill. This upward advance of erosion channels is often caused by undermining, which seems to be due to soil water coming out of the soil. Despite of decades of repeatedly focuses only on hillslope runoff mechanisms （summarized by Anderson and Burt, 1990）, the potential for the Tsumagoi hillslope saturation and erosion seems a perched water table （hereinafter, referred to as PWT） intersected attribute Return Flow Generating Point （hereinafter RFGP） on the slope topography, which is defined as location of the intersection of the PWT and the slope topography attributed to the return of infiltrated rainwater. The unique character of RFGP has been hypothesized （Fig. 1a） analytically proposed such that it would explain erosion phenomena in hillslope saturation processes. The hypothesis offers a fresh look at soil water movement in a uniform hillslope layered with traffic pan in order to the hillslope saturation formation and occurrence of RFGP attributed to erosion initiation on the variably saturated hillslope surface. The objectives of the study are to identify the mechanism of the hillslope saturation formation by observing experimentally and numerically soil water movement in uniform hillslope layered with traffic pan, and to understand the possible starting location of erosion on the variably saturated hillslope surface by observing experimentally and numerically RFGP.

CHARACTERIZATION OF TRAFFIC PAN IN A HILLSLOPE FIELD

Agricultural hillslope field located near the Alpine Crops Research Center, Tsumagoi was chosen for characterizing the presence of traffic pan based on the measurements of bulk density, saturated hydraulic conductivity, soil resistance, and particle size distribution in the soil profiles. Acute increments in bulk density and decrease in saturated hydraulic conductivity at 20 cm depth confirmed the presence of traffic pan. Soil resistance values across the upslope, midslope, and downslope transect of a sub-plot suggested that traffic pan layer was uniformly distributed at 20 cm depth. However, particle size analysis did not reveal any significant changes in fraction distributions versus the soil depth. Preliminary assessment by applying white paint tracer and entrenching tracer distributions on the soil vertical profile suggested that water movement was stopped at traffic pan depth 20 cm.

THEORETICAL CONSIDERATIONS

To derive the two dimensional flow region of interest for Tsumagoi hillsllope, the vertical cross-section of a sloping rectangular region ABCD illustrated in Fig.1 is considered. The saturated hydraulic conductivity of the top （Ks1） is much greater than that of the subjacent pan layer （Ks2）. Boundary AB allows rainfall infiltration. CD is no flux boundary because of the subjacent natural hard brown layers. DA is no flux boundary since it corresponds to the upslope divide of the hillslope. BC is no flux boundary since there is no incised stream channel at the bottom, and the hillslope opposing the hillslope is identical in character so that BC becomes a line of symmetry. It was assumed that runoff removes from this flow region and drains sporadically at the bottom of the hillslope.When the percolating rainwater reaches the interface between the highly permeable top and very low permeable traffic pan layers, the PWT tends to be developed at the downslope above traffic pan layer （Fig.1b）. During the prolonged rainfall, the PWT may then intersect the hillslope surface as RFGP, producing surface seepage as return flow, and also preventing the entry of further rainfall onto the saturated area, which runs off directly as saturation overland flow. With reference to Fig.1b, the intersection of the PWT must coincide with RFGP along the hillslope surface that mass balance at RFGP cross-section （A'D'） under steady state can then be expressed as,

（1）

where r is the steady rainfall flux; Lr is RFGP from the upslope; q is subsurface flux passing through RFGP cross-section; and A is the saturation depth at RFGP section （D1+D2）. The saturated zone is fed by water moving out of the unsaturated soils during the prolonged rainfall. Assuming that soil is isotropic, flux q at any depth of unsaturated soils, governed by Darcy's law, is decomposed into qx and qz ,

（2）

（3）

where qx the flux component in the direction of x axis; qz is the component in the direction of z axis; K is the unsaturated hydraulic conductivity and a function of matric head Ψm of the soil water; and α is the slope angle. The K is described by combining the van Genuchten's soil water retention model with the pore-size distribution model of Mualem as VGM model （van Genuchten et al., 1991）. Assuming that the PWT is always almost parallel with the slope bed so that the hydraulic gradient is equal to the hillslope surface at RFGP, lateral subsurface flow passing through RFGP section is given by,

（4） where Ks is the saturated hydraulic conductivity. Hence, RFGP at steady state can be predicted from the Eq. （1）,

（5）

Substituting the Eq. （4） into the Eq. （5） yields REGP from the upslope predictive equation,

（6）

MODEL HILLSLOPE EXPERIMENT

Volcanic ash soils （Andisol） collected from Tsumagoi hillslope plot were used for the experiments. Average values of the saturated hydraulic conductivities and bulk densities for the model hillslope top soils were 0.0225 cm s-1 and 0.3 Mg m-3, respectively. These values for traffic pan were 0.00011 cm s-1 and 0.7 Mg m-3, respectively. Soil water retention behaviors for the soils were measured. The laboratory experimental setup was composed of a rectangular soil packing Plexiglas box, a slope adjustable support, a rainfall simulator, tensiometer assemblages for soil water pressure measurements system, and outflow measurements system. All boundaries of the packing box were impermeable except for the hillslope surface boundary. Field soils sieved to pass 2 mm screen were packed for the top 15 cm and traffic pan 5 cm layers into the box. The 24 experimental runs were performed and the experimental runs were designated as 1, 2, 3, and so forth, up to 24 in the order given. For runs 1-20, rainfall rates 50, 80, 100, 125, and 150 mm h-1 were applied for the 8, 12, 16, and 20・ model hillslopes, respectively. The 15 cm top and 5 cm traffic pan were packed at a 60~70% （mass basis） water content. For runs 21 and 22, 50 mm h-1 was applied for the 8 and 20・ hillslopes with shallow soils, respectively. The 10 cm top and 10 cm traffic pan were packed at a 63~70% water content. For runs 23 and 24, 50 mm h-1 was applied for the 8 and 20・ hillslopes with relatively dry initial conditions, respectively. The 15 cm top and 5 cm traffic pan were packed at a 10~15% water content. Each experimental run was carried out until steady state condition of the outflow prevailed. The outflows were observed at 30 minutes interval. Soil water pressures were observed continuously at 30 minutes interval. Soil water pressures of tensiometric grids were used to assess the PWT and RFGPs. RFGP was visually inspected to identify the location on the hillslope surface at which erosion was considered most likely to initiate.

NUMERICAL SIMULATION

Modified HYDRUS-2D・ （Simunek et al., 1999）, a two-dimensional variably saturated flow model solving Richards' equation with Galerkin finite element method, was used to simulate soil water movement and RFGP against the corresponding model hillslope experiments. The governing water flow equation and associated initial and boundary conditions is given by the following modified form of the Richards' equation,

（7）

where θ is the volumetric water content; Ψ is the matric potential, S is a sink term which was not considered; xi （i=1, 2） are the spatial coordinates （i.e., x1 = x, and x2 = z）; t is time; and KijA are components of a dimensionless anisotropy tensor KA, which was assumed to be isotropic. Three different analytical models Brooks-Corey （BC）, van Genuchten-Mualem （VGM）, and modified VGM for the hydraulic properties were incorporated in modified HYDRUS-2D. The 100 cm ・ 20 cm two-dimensional hillslope domain was discretized into 1 cm ・ 1 cm square finite element grids. As shown in Fig.1b, boundary AB was specified as atmospheric and dynamic boundary condition switching between seepage and atmospheric. The BC, CD, and DA were no flux boundaries. The length of surface saturation （Ls） was unknown a priori for numerical simulation. Soil water returned or rain （r） falling directly on Ls was removed from surface and accounted for in the mass balance. Initial conditions were specified as the soil water pressures of each experimental run prior to rainfall. HYDRUS-2D was calibrated for Ks1, and soil water retention model fitting parameters. To reach a steady state solution, the model was run for a long time until no changes were detected in the solution by monitoring observation nodes. The 24 simulation runs were performed against the corresponding experimental runs, and were also designated as simulation runs 1, 2, 3, and so forth, up to 24 in the similar order of the model hillslope experimental runs were given.

RESULTS

Although results of the 24 experimental and simulation runs have been obtained, it is appropriate to provide a demonstrative description of the flow field conditions for the experimental run 1 of the 8・ model hillslope under 50 mm h-1. Fig.2 will be referred to in this qualitative description.

Soil Water Movement

Evolution of Matric Potentials and the PWT

Matric potential profiles of the experimental and simulation runs 1 at steady state are shown in Fig.2, indicating that the PWT above traffic pan rises in turn to an increasing depth of saturation downslope and eventually intersects the surface as RFGP. Steady state condition is attained at t = 750 and t = 420 minutes for experimental and simulation runs, respectively. Nearly the entire flow region is saturated and the PWT intersects the surface. Even at steady state, matric potential patterns indicate the infiltration zone at the upper end of the upslope. Traffic pan eventually undergoes positive pressures, being caused by the weight of the column of the PWT above this layer and slow percolation into pan layer during the prolonged rainfall.

Flux Distributions

Corresponding to Fig.2 （a, b）, flux distributions of the representing experimental and simulation runs 1 at steady state are illustrated in Fig.2 （c, d）. Simulated profiles show relatively higher fluxes, attributing to the optimized relatively higher Ks1. The infiltrating water is split into two regions with one region at the unsaturated upslope indicating nearly vertical flow, and the other a saturated region at the downslope above traffic pan indicating lateral subsurface flow at the saturated zone. Lateral unsaturated subsurface flow is observed at the midslope right above the PWT. Flux distributions indicate that saturated conditions extend up only so far as it necessary to accommodate fluxes from the upslope. Maximum flux is observed at RFGP section.

Observed and Simulated RFGPs

The Observed and simulated RFGPs for the 8・ model hillslope runs are demonstratively shown in Fig.3. RFGPs were assessed from the matric potentials at hillslope surface along the distance from upslope, and zero ordinates in Fig.3 indicate the occurrence of RFGP. Simulated RFGPs show agreement with the observed. Observed and simulated RFGPs of the 8 and 20・ hillslope with shallow top soils under 50 mm h-1 are shown in Fig.4 for the experimental and simulation runs 21 and 22. Observed and simulated RFGPs of the 8 and 20・ hillslope with relatively dry initial conditions under 50 mm h-1 are shown in Fig.5 for the experimental and simulation runs 23 and 24.

Observed and Simulated Outflows

Surface runoff was not occurred since saturation as the PWT was always moved upward from the below, and time lag of the onset of outflows for the experimental and simulation runs assured saturation overland flow. The discharge hydrographs, demonstratively shown in Fig.6 for the 8・ model hillslope experimental and simulation runs. Time lags to onset of indicate that runoff does not occur for the model hillslope show saturation overland flows, suggesting return flow and saturation-excess overland flow that may accelerate erosion. Visual inspections were performed to identify the breakout location of the PWT on the slope surface at which erosion was likely to initiate. For a demonstrative purpose, Fig.7 depicts that RFGP is visualized on the slope surface at 82 cm distance from the upslope for the experimental run 1.

ANALYSIS AND DISCUSSION

The Observed, simulated and analytically-estimated RFGP

It is more likely that all the terms in the Eq. （6） have the same proportional effect on RFGP （Lr） at steady state. Since Ks, A, and α are known parameters for model hillslope, Lr could also be predicted. Subsequently, a qualitative comparison among the observed, simulated and the predicted （Lrc） is demonstratively made for the 8・ model hillslope runs under the rainfalls of 50, 80, 100, 125, and 150 mm h-1 （Fig.8）. The analytically-estimated RFGPs are within the agreed variability of the observed and simulated RFGPs at steady.

Factors affecting RFGP

Rainfall Rates and Slope Angles on RFGP

Fig.9 shows that RFGPs decrease because of high rainfall rates, indicating more slope surface saturation. However, incremental variation for different rainfalls does not emphasize any certain extent, suggesting the significant nonlinearities in the occurrence of RFGP with respect to the rainfall rates. As shown in Fig.9, increase in slope angles increases RFGPs. Since each of the experimental and simulation runs approached downslope subsurface flux ensues at steady state, RFGP could be represented demonstratively as function of the rainfall rates （r） based on the experimental or simulation runs 1, 2, 3, 4, and 5 by using the Eq. （5）. Subsequently, prediction for RFGPs （Lrc） could be extended theoretically to a range of the rainfall rates. Fig.10 shows r as a function of Lrc derived from the experimental runs 1, 2, 3, 4, and 5.

Hillslope with Shallow Soils and Relatively Dry Initial Conditions

Fig.4 shows that shallow soils provide less profile storage, causing the PWT advanced upward and decreased RFGP from the upslope, i.e., more hillslope saturation. Such condition occurred when shallow traffic pan was developed within the plow top soils. Hillslope with relatively unsaturated initial conditions （Fig.5） provided prolonged experimental and simulation runs, and thereby decreasing RFGPs, i.e., increase hillslope saturation. However, simulation showed the threshold of drier initial conditions for the flow domain, indicating that steady state solution for simulating RFGP was more suited for the initial conditions close to saturation.

Key State Variables in modified HYDRUS-2D Configuration

Sensitivity showed that soil water retention models VGM, modified VGM or Brooks-Corey would not influence significantly on RFGPs. However, Ks1 was very sensitive to simulate RFGPs. Optimized Ks1 （0.04 cm s-1） was slightly higher than those typically seen in field （0.014~0.038 cm s-1）, but were, however, physically realistic for Tsumagoi hillslope with high infiltration capacity.

Implications on Hillslope Length

In all the experimental and simulation runs, RFGPs were observed such that Lr ・ L. However, referred to Fig.1, slope length L （100 cm in this study） does not experience RFGP unless the PWT intersects the surface. Although a long hillslope increases the available storage profiles （since 20 m long hillslope simulation shows even hillslope saturation formation due to occurrence of RFGP）, RFGP defined in this study would be observed on the long uniform hillslope wherever the PWT intersects the hillslope surface during the prolonged rainfall episodes.

Soil Water Movement and Erosion

Each of the model hillslope experimental and simulation runs, conceptualized for Tsumagoi hillslope, extends the norm of the generalized soil water movement in uniform hillslope layered with traffic pan （Fig.11）. Build-up of the saturated zone is triggered with the hydraulic conductivity contrast between the top and traffic pan layers. Traffic pan causes the perching of soil water as the PWT, leading to saturated-unsaturated subsurface flow. Fig.7, RFGP may explain hillslope field erosion phenomena on the acceptance of two processes: （i） outcoming water erodes the soil right at RFGP, and （ii） water comes out of the soil as return flow and contributes to surface runoff. The first interpretation is undoubtedly outward seepage of water that, by itself, may be an erosive mechanism. The second interpretation is that saturation overland flow is initiated at RFGP. So RFGP concept extends advances for soil conservation and erosion.

COLCLUSIONS

From Experiment

1.Traffic pan causes hillslope saturation because of the PWT fed by infiltration into unsaturated zone, and by laterally inflowing unsaturated-saturated subsurface flow during the prolonged rainfalls, resulting RFGP.

2.RFGP extends the potentiality to identify the location at which erosion （e.g., rills or ephemeral gullies） are most likely to initiate on hillslope.

From Numerical Simulation

1.Modified HYDRUS-2D shows the applicability of reproducing soil water movement phenomena in hillslope layered with traffic pan and RFGP as well.

From Experiment and Numerical Simulation

1.High rainfall rates decrease RFGPs, i.e., increase hillslope saturation

2.Uniform hillslope with smaller slope angles, shallow top soils, and relatively dry initial conditions tend to decrease RFGPs, i.e., more hillslope saturation

From Experiment, Numerical, and Analytical Considerations

1.Observed and simulated RFGPs shows agreement with the approached analytically-estimated steady state RFGPs, generalizing RFGP occurrence

2.For a long uniform hillslope layered with traffic pan, it is envisaged that RFGP defined in this study must occur wherever the PWT intersects hillslope topography during the prolonged rainfall episodes

Fig.1. Definition illustrations of （a） RFGP hypothesis, and （b） the flow domain used for model hillslope experiments and numerical simulations.

Fig.2 Illustrations of matric potential profiles in cm and flux distributions in cm h-1 for the 8・ model hillslope and simulations at steady state applying 50 mm h-1 rainfall: （a） matric potentials at experimental run 1, t = 750 minutes; and （b） matric potentials at simulation run 1 applying 50 mm h-1, t = 420 minutes; zero contour levels represents the positions of the PWT and thus RFGP, （c） flux distributions for the experimental run 1, t = 750 minutes; （d） flux distributions for the simulation run 1, t = 420 minutes.

Fig.3 Comparisons of matric potential at hillslope surface along the distance from the upslope at steady state between the observed and simulated profiles for the 8・ model hillslope runs 1, 2, 3, 4, and 5 at steady state. Top and pan layers are 15 cm and 5 cm, respectively, and soil was packed at 60~70% （mass basis） water content. Ordinate zero potential represents the occurrence of RFGP at surface from the upslope.

Fig.4 Comparisons of matric potential at hillslope surface along the distance from the upslope between the observed and simulated profiles for the 8 and 20・ model hillslope with shallow soils steady state, i.e., observed and simulated runs 21 and 22 applying 50 mm h-1. Top and pan layers are 10 and 10 cm, respectively, and soil was packed at 63~70% water content. Zero potential represents RFGP at surface from the upslope.

Fig.5 Comparisons of matric potential at hillslope surface along the distance from the upslope between the observed and simulated profiles for the 8 and 20・ model hillslope with relatively dry initial conditions at steady state, i.e., observed and simulated runs 23 and 24 applying 50 mm h-1. Top and pan layers are 15 cm and 5 cm, respectively, and soil was packed at 10~15% water content. Zero potential represents the occurrence of RFGP at surface from the upslope.

Fig.6 Discharge hydrograph as saturation overland flow （return flow coupled with the direct rainfall onto the saturated areas） for the 8・ model hillslope experimental and simulation runs 1, 2, 3, 4 and 5 under the rainfalls of 50, 80, 100, 125, and 150 mm h-1, respectively. Time lags to onset rainfall are inevitable for outflow initiation.

Fig.7 Visual inspection of RFGP on the downslope surface for the 8・ model hillslope experimental run 1 under 50 mm h-1 rainfall at steady state. Picture shows the top view of the downslope. RFGP is measured at 82 cm （indicated by arrowhead） from the upslope.. Erosion is most likely to start at 82 cm and is followed up to the outflow end. Erosion mechanism at RFGP is shown in the schematic.

Fig.8 The Observed and simulated RFGP versus the predicted RFGP for the 8・ model hillslope runs 1, 2, 3, 4 and 5 under the rainfalls of 50, 80, 100, 125, and 150 mm h-1, respectively. The 1 to 1 relationship is represented by the solid line.

Fig.9 The Observed and simulated RFGPs for the 8, 12, 16 and 20・ model hillslope runs 1-20 under the rainfalls of 50, 80, 100, 125, and 150 mm h-1 for each slope angle, respectively.

Fig.10 The analytically-estimated RFGP as the function of the different steady rainfalls derived （using Eq. 5） from the 8・ model hillslope experimental runs 1, 2, 3, 4 and 5 under the rainfalls of 50, 80, 100, 125, and 150 mm h-1, respectively, and extended to a range of the rainfall rates.

Fig.11 Schematic illustrations of soil water movement in uniform hillslope layered with traffic pan, which causes erosion initiated most likely at RFGP on the hillslope surface.

傾斜地における畑作は、機械走行の累積効果により耕盤層が形成され、しばしば雨水浸透を阻害したり、作物根の伸張に支障をもたらしたりする。特に、耕盤層を有する斜面における浸透阻害は、雨水の地表面流去水の発生を促し、土壌侵食の引き金になることが少なくない。然るに、従来の土壌水分移動と地表面流去水との関連を扱った研究は、水文学分野の流出解析などマクロスコピックのものが多く、地表面近傍に人為的に形成された耕盤層の影響については、ほとんど研究例が見当たらない。

本研究の目的は、傾斜畑作地帯の耕盤層の実態を反映した室内実験と数値シミュレーションにより、耕盤を有する斜面における土壌水分移動メカニズムを明らかにすること、および、地表面近傍に形成される水分飽和層と復帰流発生に関する知見を土壌侵食の予測と防止に適用すること、とした。

本論文は8章で構成されている。第1章は序論であり、研究の背景として傾斜畑作地帯に生じる地表面流去水とそれに伴う土壌面侵食の問題があることを述べた。

第2章は、既往の研究レビューである。ここでは、傾斜地における水文流出現象のモデリングや実験的研究に多数の事例が存在すること、しかしながら、復帰流（Return Flow）に関する知見は必ずしも多くないこと、傾斜畑地に形成される耕盤が雨水浸透と地表流去に及ぼす影響の研究は、ほとんど行われていないことを述べた。

第3章は、調査地として選んだ群馬県嬬恋村のキャベツ畑地帯における土壌特性の調査結果を述べ、この地帯は全体が傾斜地畑であること、降雨期や凍結土壌融解期において、地表面流去水が増加し、土壌侵食の発生が顕著になること、深さ約20cmの位置に機械走行に起因する耕盤層が形成されていて、この層の飽和透水係数が著しく低く、浸透阻害を起こしている可能性が高いこと、などを指摘した。

第4章は、耕盤を有する斜面における復帰流発生位置（Return Flow Generation Point, RFGP）概念を提起し、その定義を明確にし、この理論に基づく復帰流発生位置（RFGP）予測法を述べた。

第5章は、耕盤を有する斜面を再現するモデル土槽を用いた実験について述べた。モデル土槽サイズは長さ100cm、深さ22cm、奥行き5cmとし、表層土厚さを15cm、耕盤土層厚を7cmとして、群馬県嬬恋村のキャベツ畑黒ボク土壌を充填した。モデル土槽の勾配は、8度、12度、16度、20度の4種類、人工降雨強度は50, 80, 100, 125, 150 mm/hの5種類とし、表土厚さを変えた実験も含めて24種類の条件を設定して反復実験を行った。土槽則壁面からは等間隔に15点のマトリックポテンシャル測定センサーを埋め込み、その測定値から、ジオスタティスティックス（地理的統計学）で用いられるクリッジング手法によって則壁面全体のマトリックポテンシャル分布を推定した。なお、実験に供試した土壌の保水性、透水性は、別途測定した。この実験から、降雨浸透水がモデル斜面の地表に復帰する現象（復帰流発生）が確かめられ、その発生位置（RFGP）を測定することができ、更に、RFGPと土壌侵食発生位置との一致現象についても確かめることができた。

第6章は、第5章で行われた実験と同一の初期条件、境界条件のもとで、2次元の非定常水分移動数値計算を行った。プログラムとしては、現在急速に世界中に普及中の有限要素法プログラム、HYDRUD-2Dを適用したが、このプログラムは2次元斜面の非定常流解析用には完成していないため、HYDRUS-2Dの開発者であるSimunek博士と共同でプログラム改良を行い、利用可能なプログラム開発を行った。24種の実験に対応する24種の数値計算を行った結果、本研究で初めて提起されたRFGP概念の妥当性が確認され、さらにRFGPと土壌侵食発生メカニズムとの関連性が明確にされた。数値計算では、実験的に実現が困難である長大斜面のRFGPについても検討し、フィールドスケールでのRFGPと土壌侵食との関連性について予測することができた。

第7章では、RFGP概念が、理論的にも、実験的にも、数値計算においても妥当であることが検証されたことを再確認し、この新しい概念がフィールドの問題により多く適用されることの重要性を主張すると共に、フィールドスケールでの理論の妥当性の検証が必要であることも示唆し、第8章においてこれらの結論をまとめた。

以上要するに、本論文は理論、室内実験、数値解析を通じて、斜面畑に一度浸透した降雨が地表面に戻る復帰流について、その発生位置を予測する概念である、復帰流発生位置（Return Flow Generation Point, RFGP）概念を提唱し、その妥当性を実証すると共に、この概念が土壌侵食発生位置の予測に直結することの重要性も証明したものであり、学術応用上寄与するところが大きい。よって、審査委員一同は本論文が博士（農学）の学位論文として価値あるものと認めた。