We present our new multiscale pairwise-force smoothed particle hydrodynamics (PF-SPH) model for the characterization of flow in fractured porous media. The fully coupled multiscale PF-SPH model is able to simulate flow dynamics in a porous and permeable matrix and in adjacent fractures. Porous medium flow is governed by the volume-effective Richards equation, while the flow in fractures is governed by the Navier–Stokes equation. Flow from a fracture to the porous matrix is modeled by an efficient particle removal algorithm and a virtual water redistribution formulation to enforce mass and momentum conservation. The model is validated by (1) comparison to a finite element model (FEM) COMSOL for Richards-based flow dynamics in a partially saturated medium and (2) laboratory experiments to cover more complex cases of free-surface flow dynamics and imbibition into the porous matrix. For the laboratory experiments, Seeberger sandstone is used because of its well-known homogeneous pore space properties. The saturated hydraulic conductivity of the permeable matrix is estimated from a pore size and grain size distribution analysis. The developed PF-SPH model shows good correlation with the COMSOL model and all types of laboratory experiments. We employ the proposed model to study preferential flow dynamics for different infiltration rates. Here, flow in fracture is associated with the term “preferential flow”, providing rapid water transmission, while flow within the adjacent porous matrix enables only slow and diffuse water transmission. Depending on the infiltration rate and water inlet location, two cases can be distinguished: (1) immediate preferential/fracture flow or (2) delayed preferential flow. In the latter case, water accumulates at the surface first (ponding), then the fracture rapidly transmits water to the bottom system outlet. For the immediate fracture flow response, ponding only occurs once the fracture is fully saturated with water. In all cases, preferential flow is much more rapid than diffuse flow even under saturated porous medium conditions. Furthermore, infiltration dynamics in rough fractures adjacent to an impermeable or permeable matrix for different infiltration rates are studied as well. The simulation results show a significant lag in arrival times for small infiltration rates when a permeable porous matrix is employed, rather than an impermeable one. For higher infiltration rates, water rapidly flows through the fracture to the system outlet without any significant delay in arrival times even in the presence of the permeable matrix. The analysis of the amount of water stored in permeable fracture walls and in a fracture void space shows that for small infiltration rates, most of the injected water is retarded within the porous matrix. Flow velocity is higher for large infiltration rates, such that most of the water flows rapidly to the bottom of the fracture with very little influence of matrix imbibition process

The physical mechanisms that govern preferential flow dynamics in unsaturated fractured rock formations are complex and not well understood. Fracture intersections are critical relay points along preferential flow paths and control the partitioning behavior, leading to temporal delay and intermittent flow. In this work, a three-dimensional Pairwise-Force Smoothed Particle Hydrodynamics (PF-SPH) model is being applied in order to simulate gravity- driven droplet flow at synthetic fracture intersections. SPH, as a mesh-less Lagrangian method, is particularly suitable for modeling deformable interfaces, such as three-phase contact dynamics of droplets. The static and dynamic contact angle can be recognized as the most important parameter of gravity-driven free-surface flow. In SPH, surface tension and adhesion naturally emerges from the implemented pairwise fluid-fluid (s_f f ) and solid- fluid (s_sf ) interaction force. The model was calibrated to a contact angle of 65 ◦ , which corresponds to the wetting properties of water on Poly(methyl methacrylate). The accuracy of the SPH simulations were validated against an analytical solution of Poiseuille flow between two parallel plates and against laboratory experiments. Using the SPH model, the complex flow mode transitions from droplet to rivulet flow of an experimental study were repro- duced. Additionally, laboratory dimensionless scaling experiments of water droplets were successfully replicated in SPH. Finally, SPH simulations were used to investigate the partitioning dynamics of single droplets into syn- thetic horizontal fractures with various apertures (∆d_f = 0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 mm) and offsets (∆d_of f = 1.5, 1.0, 0.5, 0, 1.0, 2.0, 3.0 mm). The perfect conditions of ideally smooth surfaces and the SPH inherent advantage of particle tracking allow the recognition of small scale partitioning mechanisms and its importance for bulk flow behavior. The aim of this study is to derive an analytical correlation and interpretation of partitioning dynamics, droplet height and aperture

Characterization of karst systems, especially the assessment of structure and geometry of conduits along with forecast of state-variables, are essential for groundwater quality/quantity management and implementation/rehabilitation of large-scale engineering projects in karst regions. These objectives can be fully met by utilizing process-based discrete-continuum models, such as MODFLOW-2005 CFPv2, as employed here. However, such tools should be used with the caveat of the potential non-uniqueness of results. This research focuses on the joint-inversion of discharge, water temperature, and solute concentration signatures of Freiheit Spring in Minnesota, USA, in response to a spatiotemporally small-scale hydraulic and transport experiment. Adopting the multi-model concept to address conceptual uncertainty, seven distinctive model variants were considered. Spring hydro-chemo-thermo-graphs for all variants were simultaneously simulated, employing joint-inversion by PEST. Subsequently, calibrated models were compared in terms of calibration performance, parameter uncertainties and reasonableness, as well as forecast capability. Overall, results reveal the reliability of the discrete-continuum flow and transport modeling, even at a spatiotemporally small-scale, on the order of meters and seconds. All conceptualized variants suggest almost identical conduit tracer passage sizes which are close to the flood-pulse method estimates. In addition, the significance of immobile conduit-associated-drainable storages in karst hydrodynamic modeling, which is uniquely provided in our model code, was highlighted. Moreover, it was demonstrated that the spring thermograph and hydrograph carry more information about the aquifer characteristics than the chemograph. However, this last result can be site-specific and depends on the scale of the experiment and the conceptualized variants of the respective hydrological state.

Unsaturated fractured aquifer systems offer a domain for complex gravity-driven flow dynamics leading to the development of preferential flow along fracture networks that often strongly contributes to rapid mass fluxes. This behaviour is difficult to recover by volume-effective modeling approaches (e.g. Richards equation) due to the non-linear nature of free-surface flows and mass partitioning processes at unsaturated fracture intersections. The application of well-controlled laboratory experiments enables to isolate single aspects of the mass redistribution process that ultimately affects travel time distributions across scales. We use custom-made acrylic cubes (20 cm x 20 cm x 20 cm) in analogue percolation experiments to create simple fracture networks with single or multiple horizontal fractures. A high precision multichannel dispenser produces gravity-driven free surface flow (droplets; rivulets) at flow rates ranging from 1 ml/min to 5 ml/min. Hereby, total inflow rates are kept constant while the fluid is injected via 15 (droplet flow) or 3 inlets (rivulet flow) to reduce the impact of erratic flow dynamics. Normalized fracture inflow rates (Q_f/Q_0) are calculated and compared for aperture widths d_f of 1 mm and 2.5 mm. A higher efficiency in filling an unsaturated fracture by rivulet flow observed in former studies can be confirmed. The onset of a capillary driven Washburn-type flow is determined and recovered by an analytical solution. In order to upscale the dynamics and enable the prediction of mass partitioning for arbitrary-sized fracture cascades a Gaussian transfer function is derived that reproduces the repetitive filling of fractures, where rivulet flow is the prevailing regime. Results show good agreement with experimental data for all tested aperture widths.

Infiltration processes in fractured-porous media remain a crucial, yet not very well understood component of recharge and vulnerability assessment. Under partially-saturated conditions flows in fractures, percolating fracture networks and fault zones contribute to the fastest spectrum of infiltration velocities via preferential pathways. Specifically, the partitioning dynamics at fracture intersections determine the magnitude of flow fragmentation into vertical and horizontal components and hence the bulk flow velocity and dispersion of fracture networks. In this work we derive an analytical solution for the partitioning processes based on smoothed particle hydrodynamics simulations and laboratory studies. The developed transfer function allows to efficiently simulate flow through arbitrary long wide aperture fracture networks with simple cubic structure via linear response theory and convolution of a given input signal. We derive a non-dimensional bulk flow velocity ($\widetilde{v}$) and dispersion coefficient ($\widetilde{D}$) to characterize the system in terms of dimensionless horizontal and vertical time scales $\tau_m$ and $\tau_0$. The dispersion coefficient is shown to strongly depend on the horizontal time scale and converges towards a constant value of $0.08$ within reasonable ranges for the fluid and geometrical parameters, while the non-dimensional velocity exhibits a characteristic $\widetilde{v} \sim \tau_m^{-1/2}$ scaling. Given that hydraulic information is often only available at limited places within (fractured-porous) aquifer system, such as boreholes or springs, our study intends to provide a rudimentary analytical concept to potentially reconstruct internal fracture network geometries from external boundary information, e.g., the dispersive properties of discharge (groundwater level fluctuations).