Multiple-tracers-aided surface-subsurface hydrological modeling for detailed characterization of regional catchment water dynamics in Kumamoto area, southern Japan

Abstract

Integrated watershed modeling techniques have been applied in recent years to examine surface and subsurface interactions. Model performance is often evaluated by best fit of the hydrograph, which alone cannot explicitly explain whole catchment dynamics. To overcome this problem, this study incorporated multiple tracers (³H, ⁸⁵Kr, and groundwater temperature) into a physically-based fully distributed modeling framework for characterizing regional-scale hydrological processes in Kumamoto, southern Japan. First, a simulation performed by a hydrometrically calibrated model showed satisfactory performance for river discharge and groundwater level. However, this model showed poor fitting for isotopic composition and temperature due to the structural uncertainty of the model. A new model was established reflecting recent deep bore log data and incorporating tracer data showed acceptable accuracy for hydrographs and tracers. Thus, more reliable estimates of groundwater storage, groundwater age and water flow paths were depicted over the regional catchment. Comparisons between the two models indicate that the model structure of an area with an uncertain lower boundary can be addressed by incorporating multiple tracer data. Tracer-aided models could be applied for a holistic understanding of contaminant transport dynamics besides flow simulation.

Résumé

Les techniques de modélisation intégrée de bassins versants ont été appliquées ces dernières années pour examiner les interactions entre eaux de surface et eaux souterraines. La performance des modèles est souvent évaluée par le meilleur ajustement à l’hydrogramme, ce qui ne peut pas à lui seul, expliquer la dynamique de l’ensemble du bassin. Pour surmonter ce problème, cette étude a intégré des traceurs multiples (³H, ⁸⁵Kr, et température des eaux souterraines) dans un modèle distribué à base physique pour caractériser les processus hydrologiques à l’échelle régionale à Kumamoto, sud du Japon. D’abord, une simulation réalisée par un modèle calibré hydrométriquement a montré des performances satisfaisantes pour les débits de rivières et les niveaux piézométriques. Toutefois, ce modèle montre un mauvais ajustement pour la composition isotopique et la température du fait des incertitudes structurelles du modèle. Un nouveau modèle établi en prenant en compte des données de forages profonds récents et en intégrant des données de traçage a montré une précision acceptable pour les hydrogrammes et les traceurs. Ainsi, des estimations plus fiables du stockage des eaux souterraines, de l’âge des eaux souterraines et des voies d’écoulement de l’eau ont été décrites sur l’ensemble du bassin régional. Des comparaisons entre les deux modèles indiquent que la structure d’un modèle d’un secteur avec des limites inférieures incertaines peut être précisée en introduisant les données de plusieurs traceurs. Les modèles appuyés par des traceurs pourraient être appliqués pour une compréhension holistique de la dynamique du transport de contaminant en plus de la simulation des écoulements.

Resumen

En los últimos años se han aplicado técnicas de modelización integrada de cuencas hidrográficas para examinar las interacciones entre la superficie y el subsuelo. El rendimiento del modelo se evalúa a menudo por el mejor ajuste del hidrograma, que por sí solo no puede explicar explícitamente toda la dinámica de la cuenca. Para superar este problema, este estudio incorporó múltiples trazadores (³H, ⁸⁵Kr, y temperatura de las aguas subterráneas) en un marco de modelado totalmente distribuido de base física para caracterizar los procesos hidrológicos a escala regional en Kumamoto, al sur de Japón. En primer lugar, una simulación realizada por un modelo calibrado hidrométricamente mostró un rendimiento satisfactorio para la descarga del río y el nivel de las aguas subterráneas. Sin embargo, este modelo mostró un pobre ajuste para la composición isotópica y la temperatura debido a la incertidumbre estructural del modelo. Se estableció un nuevo modelo que reflejaba los datos recientes de los registros de perforación profunda y que incorporaba los datos de los trazadores mostraba una precisión aceptable para los hidrogramas y los trazadores. De este modo, se representaron estimaciones más fiables del almacenamiento y la edad de las aguas subterráneas y de las trayectorias del flujo regional de agua en la cuenca. Las comparaciones entre los dos modelos indican que la estructura del modelo de una zona con un límite inferior incierto puede abordarse mediante la incorporación de múltiples datos de trazadores. Los modelos asistidos por trazadores podrían aplicarse para una comprensión holística de la dinámica del transporte de contaminantes, además de la simulación del flujo.

摘要

近年来, 综合流域模拟技术已应用于研究地表水和地下相互作用。通常通过水位曲线的最佳拟合来评估模型的性能, 仅靠水文曲线无法单独解释整个集水区的动力学。为了克服这个问题, 本研究将多个示踪剂(³H, ⁸⁵Kr和地下水温度)纳入了基于物理基础的完全分布式建模框架, 以表征日本南部Kumamoto市的区域尺度的水文过程。首先, 通过水文校准模型进行的模拟能较好地模拟河流的排泄和地下水位。但是, 由于该模型的结构不确定性, 该模型显示出对同位素组成和温度的拟合性较差。建立了一个新模型, 该模型反映了最近的深孔岩心记录数据, 并结合示踪剂数据, 结果表明水位曲线和示踪剂的准确性可以接受。因此, 在该区域汇水区估算了更可靠的地下水储量, 地下水年龄和水流路径。两个模型的比较表明, 具有较低下边界的区域的模型结构可以通过合并多个示踪数据来解决。示踪剂辅助模型除了可以用于水流模拟外, 还可以用于全面了解污染物的运移机制。

Resumo

Técnicas integradas de modelagem de bacias hidrográficas têm sido aplicadas nos últimos anos para examinar as interações entre superfície e subsuperfície. O desempenho do modelo é frequentemente avaliado pelo melhor ajuste do hidrograma, que sozinho não pode explicar explicitamente toda a dinâmica da bacia. Para superar esse problema, este estudo integra diferentes traçadores (³H, ⁸⁵Kr e temperatura da água subterrânea) em uma estrutura de modelagem totalmente distribuída com base física para caracterizar processos hidrológicos em escala regional em Kumamoto, sul do Japão. Em primeiro lugar, uma simulação realizada por um modelo calibrado hidrometricamente mostrou um desempenho satisfatório para a descarga do rio e nível do lençol freático. No entanto, este modelo apresentou ajuste pobre para composição isotópica e temperatura, devido à incerteza estrutural do modelo. Um novo modelo foi estabelecido refletindo dados recentes de registro de um poço profundo e incorporando dados de traçadores, mostrando uma precisão aceitável para hidrogramas e traçadores. Assim, estimativas mais confiáveis de armazenamento de água subterrânea, idade da água subterrânea e caminhos de fluxo de água foram representadas sobre a bacia hidrográfica regional. As comparações entre os dois modelos indicam que a estrutura do modelo de uma área com um limite inferior incerto pode ser abordada incorporando vários dados de traçadores. Modelos auxiliados por rastreadores podem ser aplicados para uma compreensão holística da dinâmica de transporte de contaminantes, além da simulação de fluxo.

Appendices

Appendix 1: governing equations of GETFLOWS

The governing equations of GETFLOWS simulator include fluid flow, transport of dissolved materials (e.g., concentrations and isotope ratios), and heat. Fluid flow is modeled as coupled surface and subsurface flows, and the associated conjunctive behavior of the dissolved materials and heat can be traced simultaneously. The individual governing equations are based on the mass and energy conservation laws. Table 1 summarizes the major parameters shown in the governing equations.

Table 1 The parameters used in the governing equations for operation of the GETFLOWS

Coupled surface and subsurface fluid flows

The governing equations are derived from the mass conservation laws for each fluid phase (surface and subsurface). Manning’s law was adapted for surface water flow, and diffusive wave approximation of Saint Venant equations was applied to estimate the surface water velocity in 2D space. Subsurface fluid flows (i.e., groundwater) are characterized by generalized Darcy’s law in 3D space. The law of mass conservation for both surface and subsurface fluid flows can be expressed as follows:

$$ -\nabla \bullet {M}_p+{\rho}_p{q}_p=\frac{\partial \left(\phi {\rho}_p{S}_p\right)}{\partial t},\left(p=\mathrm{water},\mathrm{air}\right) $$

(2)

where, ∇ and ρ_p denote differential operator and density (kg/m³) of fluid, respectively. Mass flux of fluid (kg/m²/s) and volumetric flux of source and sink (m³/m³/s) are represented by M_p and ϕ, and q_p and t represent effective porosity (m³/m³) and time (s), respectively. The saturation of fluid (m³/m³) is embodied by S_p, and the subscript p is used to denote the air (a) or water (w) phase of the fluid. The effective porosity value is 1.0 in Eq. (1) for the surface fluid.

Fluid mass flux (M_p) per unit area can be given by adapting Manning’s law for surface flow and Darcy’s law for subsurface as follows:

Surface fluid:

• Water

$$ {M}_p=-\frac{\rho_p{R}^{\frac{2}{3}}}{n_l}\sqrt{\left|\frac{\partial h}{\partial l}+\frac{\partial z}{\partial l}\right|}\operatorname{sgn}\left(\frac{\partial h}{\partial l}+\frac{\partial z}{\partial l}\right),\left(p=\mathrm{water}\ l=x,y\right) $$

(3)

• Air

$$ {M}_p=-\frac{\rho_pK{k}_{r,p}}{\mu_p}\nabla \left({P}_p+{\rho}_p gz\right),\left(p=\mathrm{air}\right) $$

(4)

Subsurface fluid:

$$ {M}_p=-\frac{\rho_pK{k}_{r,p}}{\mu_p}\nabla \left({P}_p+{\rho}_p gz\right),\left(p=\mathrm{water},\mathrm{air}\right) $$

(5)

where, n_l indicates Manning’s roughness coefficient (m^-1/3 s). The R, h and l represent hydraulic radius (m), surface water depth (m), and surface water flow distance (m), respectively. The z, g, μ_p and P_p denote elevation from datum surface (m), gravitational acceleration (m/s²), viscosity coefficient (Pa s), and pressure (Pa), respectively. K and k_{r, p} are intrinsic permeability (m²) and relative permeability, respectively.

The surface-water depth h was computed from the water saturation and the height of the grid block in the surface environment. The difference between air and water pressure is associated with capillary pressure P_c (Pa):

$$ {P}_{\mathrm{a}}={P}_{\mathrm{w}}+{P}_{\mathrm{c}} $$

(6)

The relative permeability and capillary pressure are a function of the water saturation in the two-phase flow system. Additionally, the fluid density and viscosity coefficient are a function of fluid pressure and temperature. The total of fluid saturations is treated as 1.0 in the two-phase flow system as follows:

$$ {S}_{\mathrm{w}}+{S}_{\mathrm{a}}=1.0 $$

(7)

Transport of dissolved materials in water

This study coupled the water flow with mass transport of groundwater age tracers (³H and ⁸⁵Kr) as dissolved materials in water. The mass balance of dissolved material i in water can be expressed as the following advection-diffusion/dispersion equation considering radioactive decay as:

$$ -\nabla \bullet \left({M}_{\mathrm{w}}{C}_{\mathrm{w},i}\right)+\nabla \bullet {D}_{\mathrm{w},i}\nabla \left({\rho}_{\mathrm{w}}{C}_{\mathrm{w},i}\right)+{\rho}_{\mathrm{w}}{q}_{\mathrm{w}}{C}_{\mathrm{w},i}-\phi {\lambda}_i{C}_{\mathrm{w},i}=\frac{\partial \left(\phi {\rho}_{\mathrm{w}}{S}_{\mathrm{w}}{C}_{\mathrm{w},i}\right)}{\partial t} $$

(8)

where, C_{w, i} is the mass fraction of dissolved material i in water (kg/kg) (defined as the mass of dissolved material i per unit mass of water in each grid). D_{w, i} and λ_i are the hydrodynamic dispersion coefficient (m²/s) and the decay constant (1/s) of dissolved material i, respectively. The hydrodynamic dispersion coefficient D_{w, i} is computed from the molecular diffusion coefficient, tortuosity factor, dispersion length, and effective water velocity in porous media. The λ_iis derived from the half-life of each radioactive material (³H and ⁸⁵Kr). Adsorption and desorption are not considered for these isotopes. The molecular diffusion coefficient is considered as a constant value in this study.

Heat transport

The governing equation for heat transport modeling consists of energy conservation equations for fluid and solid phases in addition to mass balance equations for air and water phases under nonisothermal conditions. The equation can be expressed as:

$$ -\nabla \bullet \left({M}_{\mathrm{w}}{H}_{\mathrm{w}}\right)-\nabla \bullet \left({M}_a{H}_{\mathrm{a}}\right)+\nabla \bullet {\kappa}_{\mathrm{f}}\nabla T+\nabla \bullet {\kappa}_{\mathrm{s}}\nabla T+{\rho}_{\mathrm{w}}{q}_{\mathrm{w}}{H}_w+{\rho}_a{q}_a{H}_a+{E}_{\mathrm{s}}=\frac{\partial \left(\phi {\rho}_{\mathrm{w}}{S}_{\mathrm{w}}{U}_{\mathrm{w}}+\phi {\rho}_{\mathrm{a}}{S}_{\mathrm{a}}{U}_a+\left(1-\phi \right){\rho}_{\mathrm{s}}{U}_{\mathrm{s}}\right)}{\partial t} $$

(9)

where, H_w and H_aare the enthalpies of water and air phases (J/kg), respectively, κ_f and κ_s are the saturation weighted average thermal conductivities of fluid and solid phases (W/m/K), respectively, T is the average temperature for both fluid and solid phases (K), E_s is the volumetric heat flux of the sink and source in the solid phase (J/m³/s), ρ_s is the density of solid phase (kg/m³), and U_w, U_a and U_R are the internal energies of water, air, and solid phases (J/kg), respectively. The heat transport parameters such as heat capacity and thermal conductivity are treated as fixed values in this simulation.

Appendix 2: verification and validation

In the early stage of the development of the GETFLOWS simulator, Tosaka et al. (2000) verified surface flow simulation by comparing with experimental results, and its application was proved in some local fields (Itoh et al. 2000). Studies (Tosaka et al. 1996, 2000) also compared results of some analytical solutions for surface–subsurface flow coupling simulation, e.g., relative permeability and capillary rise of water, using GETFLOWS (Tosaka et al. 1996, 2000). More recently, Mori et al. (2015) have shown inter-comparisons with another simulator (i.e., InHM, VanderKwaak 1999) for surface flow and sediment discharges. The detailed verification and validation (V & V) of GETFLOWS are available through the website of the Geosphere Environmental Technology Corp. (2020) or direct email. However, some basic information on V & V procedures and model discretization are provided also in this paper’s Appendices and ESM, since these manuals are not completely accessible through their website.

In the development of the GETFLOWS systems, a lot of verification was undertaken for relevant cases, mainly for the subsurface domain, by comparing results of analytical solutions, inter-comparisons with other simulators, and using results of experimental study. Detailed descriptions of the analytical solutions are provided in Sections S3–S8 of the ESM. Briefly, an analytical solution for the conceptual model is shown in Fig. S14 of the ESM and numerical simulation for 1D saturated horizontal groundwater flows (Fig. S15 of the ESM) and 2D vertical flows (Fig. S16 of the ESM) are described in Section S3 of the ESM. Details of the numerical model and analytical solutions are provided in Tables S7–S10 of the ESM, and their results are presented in Tables S11 and S12 of the ESM. Similarly, the simulator was also verified using experimental results of a falling head test to monitor water level changes over a given time (Tosaka 2007; Section S4, Figs. S17–S19 and Tables S13–S15 of the ESM). Furthermore, verifications for the pumping test of a confined aquifer (Figs. S20–S22, Tables S16–S18 of the ESM) and capillary pressure of unsaturated zone (Figs. S23 and S24, Tables S19–S21 of the ESM) were provided in sections S5 and S6. These comparisons showed excellent agreement between analytical and numerical solutions and ensured that there is no bug in the GETFLOWS numerical code for water flow simulation.

Comparison of simulation results between GETFLOWS and TOUGH2 (Pruess et al. 1999), another multiphase flow simulator, is documented in Section S7 of the ESM, for simulating groundwater flows in porous media with heterogenous permeability (Fig. S25 of the ESM). This inter-comparison example was taken from the international HYDROCOIN project for groundwater flow simulation evaluation (Thunvik 1987). A 2D heterogeneous porous medium placed between impermeable layers was considered for the simulation (Fig. S26 of the ESM). Initially, the system was kept in an unsaturated condition; then, water was injected to allow fluids flow through the media replacing the air (Fig. S26 of the ESM). Further details on model domains and properties are given in Tables S22–S24 of the ESM. Comparisons of transient simulation results between GETFLOWS and TOUGH2 (Fig. S27, Fig. S28 of the ESM) showed almost identical results over the time steps. Overall, the GETFLOWS simulations showed an excellent agreement with TOUGH2 for transient simulation (Fig. S27, Fig. S28 and Table S25 of the ESM).

Several previous studies have already demonstrated the capability of GETFLOWS for simultaneous simulations of water flows, isotope tracers and sediment transport (e.g., Mori et al. 2015; Kitamura et al. 2016; Sakuma et al. 2017, 2018; Hosono et al. 2019; Tawara et al. 2020). However, number of studies simulating for heat transport processes along with groundwater flow systems are still small on certain field sites using GETFLOWS. In principal, performance of GETFLOWS simulation is verified with the theoretical solutions (Domenico and Palciauskas 1973) for simultaneous transport of groundwater and heat. All the model parameters (Table S26) relevant to the Domenico and Palciauskas 1973) problem (Fig. S28 of the ESM) and numerical simulation details are provided in Section S8 of the ESM. All numerical simulation results extracted both for vertical and horizontal plains, and the comparisons between analytical and numerical solutions, showed identical results (Fig. S30 of the ESM). Hence, GETFLOWS is capable of simulating heat transport processes.