Solid-phase dissolution is a critical process in subsurface applications, including solution mining of salt caverns, dam foundation stability, and geological carbon sequestration in deep saline aquifers. Under mixed convection, where gravity-driven natural convection and external forced convection coexist, dissolution morphologies are governed by mass transfer at solid-fluid interfaces. However, how buoyancy and flow dynamics influence mass transfer is still not well understood. In this study, by combining micro-particle image velocimetry (micro-PIV) experiments with three-dimensional pore-scale simulations, we investigate the interplay among morphology evolution, velocity fields, and concentration fields in initially flat fracture channels under varying horizontal forced and vertical natural convection conditions. Our results reveal that enhanced forced convection leads to a transition from wedge-shaped to uniform dissolution patterns at the channel scale, while enhanced natural convection induces localized dissolution troughs on the upper surface. These trough-triggered surface roughness, coupled with buoyancy-induced convective rolls, generate eddies that destabilize the concentration boundary layer (CBL), thereby enhancing mass transfer between the troughs and main flow channels. Based on quantitative analysis of solid dissolution, we establish a scaling relationship between the dissolution rate and the mixed convection intensity. This study provides new insights into the complex dynamics of mixed convection-driven dissolution and offers a framework for predicting mineral dissolution rates in large-scale reactive transport modeling.

Understanding particle transport behaviors in porous media during multiphase flow is crucial for a wide range of applications in environmental engineering, energy recovery, and materials science. This paper provides a comprehensive overview of recent advances in experimental techniques, mechanistic studies, and modeling approaches related to particle transport under multiphase flow conditions. We first review recent progress in experimental platforms and imaging techniques that enable visualization and quantification of particle dynamics in porous media. Subsequently, key transport mechanisms, including attachment and detachment, straining, size exclusion, and aggregation, are systematically revisited, along with the influence of various physicochemical conditions. Special attention is given to the coupling between particles and deformable interfaces in, e.g., bubbles, droplets, and thin films. Furthermore, we examine recent developments in numerical modeling approaches for particle transport and multiphase flow. Knowledge gaps in particle transport mechanisms during multiphase flow and challenges associated with multiscale modeling are also outlined. This review can provide a unified framework for understanding and predicting particle transport in porous media under multiphase flow conditions, and guiding future research on interface-mediated particle manipulation and transport control.

Tunnel excavation through large-scale fault zones poses significant risks of water inrush; however, the evolution of unsteady seepage under the coupled effects of progressive excavation and structural heterogeneity (e.g., internal sub-layers) remains insufficiently understood. This study investigates a tunnel section intersecting a major fault zone in the Yangtze-to-Han River Water Diversion Project. A numerical modeling framework integrating the finite element method with the Parabolic Variational Inequality (PVI) approach is developed to simulate transient seepage, explicitly incorporating progressive excavation sequences and seepage control measures according to regional hydrogeological conditions and engineering design specifications. Three-dimensional transient simulations are conducted to characterize groundwater flow evolution and evaluate seepage stability throughout excavation. The results show that excavation alters seepage boundaries and significantly increases hydraulic gradients near the tunnel. As tunneling approaches the fault zone, groundwater disturbance intensifies. The combined application of pre-grouting and advanced drainage effectively reduces the peak tunnel inflow from 26,300 m3/d to 11,397 m3/d (a 56.7% decrease) and decreases the maximum hydraulic gradient from 59.8, initially concentrated within the fault gouge sub-layers and near the tunnel face, to 13 (a 78.2% reduction). Sensitivity analyses further reveal that tunnel inflow increases with the permeability of the grouting zone and with the number of drainage holes, whereas hydraulic gradients decrease correspondingly. These findings highlight the importance of explicitly incorporating internal fault-zone sub-layer structures into transient seepage modeling to accurately identify excavation stages most susceptible to seepage instability. Such consideration enhances hydrogeological risk assessment for deep-buried tunnels crossing large-scale faults with complex internal structures.

The construction of a large-scale hydropower project inevitably leads to disruption of hydrogeological conditions at the site, primarily through extensive underground excavations and reservoir impoundment. Understanding the resultant groundwater dynamics is therefore essential to safety and environmental assessment of the project. This study investigated the transient seepage behavior around an underground powerhouse currently under construction. The sources of leakage into each tunnel during construction were analyzed via hydrochemical analysis. The permeability of rock formations was determined by packer tests and calibrated against available groundwater level and tunnel inflow observations. The groundwater dynamics during the entire process of construction and reservoir filling were assessed with transient flow modeling, with an emphasis on the initial recharge pattern, phreatic surface evolution, and leakage contribution from each source. The results show that the groundwater in the cavern area is primarily recharged from a regional-scale fault F34 at the rear margin and secondarily from two perennial gullies. The earliest-excavated bypass access tunnel (BAT) has acted as a permanent sub-horizontal drain to the groundwater in the cavern area, with its discharge peaking at 79.25 L/s during excavation, stabilizing at 47.98 L/s during subsequent construction, and increasing to 53.65 L/s after reservoir filling. The subsequent construction of the underground caverns and drainage system leads to a groundwater drawdown of about 120 m at maximum and 55 m on average, with an influence radius of about 250 m. Groundwater level recovery is predicted to occur only locally near the reservoir, with the depletion cone in the cavern area being maintained by the drainage system. The results deepen our understanding of groundwater flow dynamics at the site, which assists in seepage control design for the hydropower project with improved performance evaluation accuracy and reliability.

Abstract Coupled mineral dissolution‐precipitation (CMDP) in fractured rocks alters pore structure and flow pathway in many reactive subsurface systems. Although previous studies have shown that simultaneous dissolution of primary minerals and precipitation of secondary phases can significantly modify the macroscopic porosity‐permeability relationship, the pore‐scale mechanism controlling precipitate distribution, texture, and their feedback on macroscopic properties remains poorly understood. Here, we combine real‐rock microfluidic experiments with high‐resolution imaging to investigate CMDP in limestone fractures exposed to injection of acidic sulfate solutions under different concentrations and flow rates. High‐speed microscopy tracks the dynamic evolution of the channels during dissolution, and post‐experiment scanning electron microscopy resolves the texture and spatial patterns of gypsum precipitation. We identify three CMDP regimes: a porous precipitation regime with localized supersaturation; a compact dissolution‐precipitation regime with inlet‐focused dissolution and outlet‐focused precipitation; and a uniform precipitation regime characterized by continuous, dense precipitate layers formed under pervasive supersaturation. A dimensionless timescale ratio of precipitation to advective transport is introduced to delineate regime boundaries. We quantify the passivation of calcite dissolution by secondary gypsum, and evaluate the evolution of hydraulic transmissivity under the combined effects of dissolution, precipitation, and gas bubble dynamics. Notably, we observe a counterintuitive hydraulic response where hydraulic transmissivity can stagnate or even decrease despite an increase in mean fracture aperture, due to low‐transmissivity bottlenecks formed by localized precipitation. These findings elucidate pore‐scale mechanisms controlling hydraulic property evolution and provide a regime‐based framework for parameterizing dissolution rates and transmissivity in reactive transport modeling.

Spontaneous imbibition at the nanoscale plays a pivotal role in both natural phenomena and industrial applications. Despite numerous theoretical models have been proposed, their validity remains to be further validated, especially in the cases with strong solid-liquid interactions. To address this knowledge gap, we employ molecular dynamics simulations with atomic resolution to characterize spontaneous water imbibition in calcite nano-channels, motivated by their prevalence and intense solid-liquid interactions. We find that the three-phase contact line on the calcite surface moves forward in a rolling mode. We demonstrate that this rolling mode is primarily attributed to the layered structures and the restricted mobility of water molecules on the calcite surface due to the strong water-rock interactions. We show that the classical theoretical models fail to fully capture nanoscale imbibition dynamics. We thus propose an improved Lucas-Washburn model to incorporate the dynamic contact angle, the advancement mode of the three-phase contact line, and changes in water properties. The simulation results show that the proposed model can predict the imbibition process in the calcite nano-channels with sufficient accuracy. The findings of this work not only enrich the understanding of nanoscale spontaneous imbibition dynamics, but also offer theoretical guidance for the related engineering applications.

The capillary pressure curve provides fundamental insights into the dynamics of fluid-fluid displacement and phase distributions. Capillary scaling is crucial for extrapolating capillary pressure-saturation data from laboratory tests to field applications. However, the classic scaling method fails to capture the effect of wettability as the pore surface approaches neutral wetting. Here, inspired by the role of pore-filling events in controlling fluid-fluid displacement, we perform a theoretical analysis of the burst events occurring during drainage processes. We find that the median threshold capillary pressure, which corresponds to the occurrence of burst events for the median pore throat, is closely correlated with the capillary pressure curve across various contact angles. Using this concept, we propose a new scaling method for capillary pressure curves under various wetting conditions. We conduct microfluidic experiments and pore-network modeling across different contact angles, porosities, and disorders to evaluate the new scaling methods, indicating that the new scaling method performs better than the Leverett J-function as the contact angle approaches 90°. We further perform geometry analysis on the critical radius of curvature for burst events in an ideal tetrahedral arrangement and extend the new scaling method to 3D (three-dimensional) porous media. Model evaluation shows that the 3D version of the scaling method also performs well but requires fewer parameters compared to the Leverett J-function. Our work enhances the prediction and interpretation of experimental data for capillary pressure curves under various wet conditions, and more importantly, establishes a methodology that relates Darcy-scale flow behavior to pore-scale fluid displacements.

Immiscible displacement in rough-walled fractures governs multiphase transport in many subsurface energy and environmental systems, yet the coupled influence of fracture roughness and wettability under favorable viscosity-contrast conditions remains unclear. Here, we investigate how roughness, wettability, and structural anisotropy jointly control immiscible displacement under a favorable viscosity ratio (M = 10) using a three-dimensional multicomponent Shan-Chen lattice Boltzmann model. An explicit wetting boundary condition with bias improvement is implemented to accurately impose contact angles on rough surfaces while suppressing spurious currents over a wide wettability range (15(degrees) <= theta <= 175(degrees)). Simulations in self-affine rough fractures reveal a non-monotonic dependence of breakthrough saturation on wettability, with a universal maximum near theta approximate to 80(degrees). This behavior is resulted from the competition between corner-driven wetting film transport and cooperative meniscus filling. Increasing fracture roughness progressively weakens wettability-controlled mechanisms and promotes roughness-induced channelization, leading to reduced displacement efficiency under strongly non-wetting conditions. Three-dimensional interface analyses further show that roughness regulates wetting film continuity and pinch-off, whereas structural anisotropy constrains transverse connectivity and accelerates preferential flow development. These results reconcile previously reported discrepancies regarding wettability effects in fractured media and indicate that immiscible displacement gradually transitions to a geometry-dominated regime once wetting film transport is suppressed (i.e., theta > 90 degrees). Further quantitative analysis shows that when fracture roughness reaches lambda >= 0.25 and the structural anisotropy satisfies a >= 4, the displacement structure evolves into a channelized invasion pattern controlled by fracture geometry.

Rocks release subtle geochemical warning signals before breaking. These signals, coming from naturally occurring nuclides (e.g., radon, helium, argon, and thoron), have often been reported before earthquakes, volcanic eruptions, landslides, and rock and ice avalanches. However, despite their high sensitivity to deformation, their detectability, as well as myriad promising observations over half a century, nuclide signals are still far from being applied to geohazard prediction or widely used for monitoring. Here, we first develop a decomposition and interpretation method for nuclide signals. By analyzing nuclide signal time series observed from a month-long laboratory rock failure experiment and year-long slope deformation in a field setting, we identify a universal paradigm unit of nuclide signal evolution. We find that this paradigm unit is characterized by two core characteristics: a transient pulse and equilibrium fluctuation which are intrinsically correlated to rupture area and crack aperture, respectively. Through analytical derivation and pore-scale simulations, we establish the constitutive equations that link these characteristic nuclide signals to key rupture structural parameters. Rooted in these constitutive relations, we further develop a diagnostic theory of rock rupture via nuclide signals. We apply the model to track rock failures at the laboratory and field scale. The proposed nuclide signal decomposition and rupturing model enable the unification of discrete signal units emitted by individual microrupturing events, with the integrated signal evolution observed during macroscopic failure. This integration may serve as a foundation for both the mesoscopic assessment of rock damage and the early warning of geohazards induced by rock ruptures.

Two-phase flow driven by coupled gravitational and viscous instabilities is widely encountered in many natural systems and engineering applications. The interfacial dynamics involved in such processes are highly complex and remain insufficiently explored experimentally in three dimensions (3D). In this study, a 3D visualization platform based on planar laser-induced fluorescence is developed to investigate the effects of injection rate and pore structure on immiscible displacement behavior in a non-wetting phase displacing a wetting phase under coupled gravitational-viscous instability conditions. We identify four distinct unstable displacement morphologies, namely, discrete buoyant droplet, slug, straight continuous fingering, and dendritic fingering cluster. We establish an experimental phase diagram of the displacement pattern in the space of capillary number and the Bond number. We also introduce a set of geometric descriptors to quantify morphological evolution and phase connectivity. This work reveals how the competition and coupling among capillary, viscous, and buoyancy forces govern fluid fragmentation, coalescence, and fingering development in 3D porous media. The findings provide mechanistic insights into multiphase flow instabilities and have direct implications for a number of subsurface applications.

Mineral precipitation is ubiquitous in subsurface environments, influencing processes such as karst evolution and carbon mineralization. While precipitation dynamics have been extensively studied in porous media, the pattern formation and transitions in fractured media remain underexplored. In this study, we develop a novel experimental system to investigate precipitation dynamics in microfluidic fractures by integrating charge-coupled device camera, micro-PIV, and confocal microscopy. Flow-through experiments on CaCO3 precipitation are performed by co-injecting Na2CO3 and CaCl2 under controlled flow rates Pe and saturation indices SI. Real-time imaging of precipitation dynamics and velocity fields revealed two distinct patterns. At low Pe and low SI, mineral precipitation exhibits precipitation band acting as a barrier that inhibits mixing of the reactants and results in minimal permeability reduction. In contrast, at high Pe and high SI, the system transitions to precipitation clusters, characterized by widespread particle distribution that significantly reduces permeability. We demonstrate that this pattern shift is governed by the balance between fluid shear forces and repulsive forces between CaCO3 particles. When repulsive forces dominate, particles cannot aggregate, leading to band formation, whereas shear-induced aggregation promotes cluster growth. Theoretical analysis is developed to interpret the regime transition, and a phase diagram mapping precipitation regime as a function of Pe and SI shows well agreement with experimental results. These findings provide critical insights into fracture mineral precipitation, which are crucial for predicting injectivity and permeability in CO2 mineralization processes.

Weak imbibition in permeable media is crucial for CO2 geological sequestration, oil and gas extraction, and groundwater contamination and remediation. Previous studies predicted that weak imbibition under favorable viscosity ratios produces stable displacement, but this assumption may fail in fractured porous media. Here, we combine high-resolution microfluidic experiments with theoretical analysis to investigate water-air imbibition with four aperture-throat ratios at seven capillary numbers (Ca). We identify three weak imbibition patterns: capillary-dominated stable pattern at low Ca, viscous-dominated stable pattern at high Ca, and a previously unrecognized unstable pattern at intermediate Ca. The region of unstable pattern, in terms of Ca, expands with the aperture-throat ratio. Pore-scale imaging reveals that the unstable pattern arises from the formation and breaching of capillary barriers at fracture intersections, coupled with burst events in the porous matrix and preferential flow in the fractures. Based on this mechanism, we perform scaling analysis and derive the critical capillary numbers that capture the transitions from stable to unstable imbibition patterns. The predicted phase diagram aligns well with our microfluidic experiments. We further find that the interactions between fractures and matrix are co-current across all Ca. For unstable patterns, the capillary barriers locally enhance the co-current imbibition while viscous forces suppress it, causing a reduction in dimensionless imbibition rate within the matrix. Our work provides a basic understanding of the role of capillary barriers in controlling imbibition patterns. It offers new predictive capability for fluid-fluid displacement in fractured porous systems, with implications for subsurface resource recovery and contaminant mitigation.

Fluid-rock interactions involving chemical dissolution, mechanical erosion, and multiphase flow are central to a wide range of geological and engineering processes, yet they remain poorly understood due to the lack of integrated in situ observation tools. Existing methods often compromise between spatial resolution and temporal dynamics. Here, we develop a real-rock microfluidic platform that enables simultaneous visualization and quantification of erosion dynamics in multiphase reactive systems. The platform integrates fluorescence microscopy, micro-particle image velocimetry, and ion chromatography to monitor the coupled evolution of solid-liquid-gas interfaces and flow velocity fields at micrometer-scale resolution. Microfluidic chips fabricated directly from limestone preserve natural mineral heterogeneity, and the platform enables direct observation of rock surface evolution and multiphase flow behavior. This facilitates decoupled analysis of chemical dissolution and mechanical erosion-two processes often difficult to isolate in traditional systems. Using this system, we investigate erosion during acid-rock interactions and identify a transition between two regimes-transport-limited and reaction-limited-controlled by CO2 bubble mobility. In the transport-limited regime, immobile bubbles confine flow to thin films, enhancing dissolution and particle detachment. In the reaction-limited regime, surface-adhered bubbles shield reactive areas and reduce shear stress, suppressing erosion. We derive scaling laws that distinguish chemical and mechanical erosion rates and validate a theoretical model for the critical Péclet number marking the regime transition. This study advances understanding of erosion under multiphase flow and introduces a versatile experimental framework for probing pore-scale reactive transport. The platform can be extended to other rock types and fluids, offering a powerful tool for studying geochemical, physical, and biological processes in complex subsurface environments.

The coupled mineral dissolution-precipitation process plays a critical role in natural and engineered systems with fluid-rock interaction including geological CO2 sequestration and contamination treatment of groundwater. A representative scenario is the reaction between limestone and acidic, sulfate-bearing solutions, where the dissolution of carbonate minerals drives gypsum precipitation. Previous studies have investigated the impact of this coupled process on hydraulic parameters through core-flooding experiments. However, the interplay between calcite dissolution (CaCO3) and secondary precipitation (CaSO4) as well as its influence on the reaction rate of CaCO3 remains unclear. Here, we combine flow visualization experiments in limestone-based microfluidics with post-experiment SEM-EDS characterization to probe how sulfate concentration affects both precipitation morphology and CaCO3 dissolution rates under the same pH condition (pH=1). Two distinct dissolution-precipitation coupling regimes have been identified. At low sulfate concentrations, the dissolution-dominated regime causes continuous dissolution of limestone throughout experiments because the porous gypsum coating permits continuous acid infiltration into the rock surfaces. Conversely, high sulfate concentrations induce a precipitation-dominated regime, reducing dissolution rates by an order of magnitude as a dense, thick coating completely shields the reactive surface area. We find a critical sulfate concentration threshold that triggers regime transition, wherein the released calcium concentration equilibrates with the supplied sulfate in solution, facilitating the formation of a dense coating layer. Quantified results indicate that the average dissolution rate decreases linearly with the saturation index by one order of magnitude as sulfate concentration increases from 0 to 0.2mol/L. These findings provide insights into coupled geochemical processes and offer practical guidance for predicting reaction rates in reactive transport modeling.

Natural fractures are commonly filled with materials such as sediments and mineral cements. Under hydrodynamic conditions, these infillings may be scoured, eroded, or removed, leading to alterations in pore structure and bulk permeability. However, the mechanisms driving these changes, particularly the interactions between hydrodynamic conditions, particle migration, and permeability evolution, remain insufficiently understood. In this study, we conducted a series of visual hydrodynamic erosion experiments on fully-filled, rough-walled fractures with varying apertures and roughness characteristics. Using well-calibrated image monitoring and processing techniques, we tracked the erosion process in real time and quantified the resulting eroded flow channels. The results identify five distinct stages across the entire erosion process: particle incipient motion, erosion initiation along with channel penetration, erosion acceleration, deceleration, and depletion. The compiled phase diagrams indicate that the Reynolds number plays a decisive role in erosion dynamics, with fracture aperture serving as the primary geometric control while roughness having a comparatively weaker impact under the identical hydrodynamic condition. We further developed two phenomenological models to predict the variations of erosion ratio and bulk permeability throughout the erosion process. These models capture the effects of Reynolds number, aperture, and roughness on the initiation, growth, and stabilization of erosion and permeability changes. These findings offer a deeper understanding of how hydrodynamic forces drive erosion in complex fracture systems and provide valuable insights into various fields concerned with the coupled issues of seepage and erosion.

Abstract Natural or engineered microparticles are often encountered in subsurface multiphase flow systems. This introduces a complex flow scenario with a variety of applications. However, the influence of particles on multiphase flow dynamics and the underlying mechanism remain elusive. Here, we investigate particle transport behavior and fluid phase distribution within 3D porous media through direct visualization utilizing laser scanning confocal microscopy. The mechanisms and factors governing particle aggregation during two‐phase flow are elucidated. We identify a previously overlooked pore‐scale phenomenon: fragmentation of wetting liquid induced by spontaneous particle aggregation in localized regions. This process dramatically increases the number of wetting clusters while reducing the fluid connectivity, resulting in a change in the relative permeability of fluids. These findings reveal the dynamic coupling mechanism between pore‐scale particle aggregation and fluid flow and transport properties at larger scales, providing critical insights for predicting and controlling particle mediated geophysical flow processes.

Reactive flow through geological fractures under stress is fundamental to subsurface hydrology. In this flow-stress-dissolution system, free-surface dissolution within voids increases permeability, whereas pressure dissolution at contacting asperities combined with mechanical deformation seals the fracture. The interplay between these mechanisms controlling fracture sealing or opening remains to be explored. Here, we integrate numerical simulations with theoretical analysis to examine the coupled mechanical-chemical processes. We develop a flow-stress-dissolution modeling approach with the incorporation of pressure dissolution, free-surface dissolution and mechanical deformation for fractures in homogeneous mineral systems. Comparison with previously experiments demonstrates the ability of this approach and highlights the significant role of pressure dissolution for fracture sealing. Through extensive simulations with a wide range of flow rates, reaction rates and normal stresses, we find that stress induced-mechanical compaction enhances the positive feedback mechanisms between free-surface dissolution and solute transport, promoting the development of dissolution channels. Our analysis of permeability evolution under various stresses and dissolution patterns reveals that fracture sealing can occur in uniform and compact patterns depending on stress, whereas fracture opening spontaneously occurs in wormholes. Through theoretical analysis of the interplay between pressure dissolution and free-surface dissolution, we propose a theoretical model for the occurrence of fracture sealing. A generalized phase diagram is then established that can effectively describe fracture sealing under varying flow rates, reaction rates, and normal stress conditions. Our work provides a foundation for assessing fracture leakage risks in subsurface engineering and holds practical implications for geological carbon and hydrogen storage.

Particle transport and clogging in porous media is a critical process for a number of engineering applications, including membrane filtration, groundwater remediation, and hydraulic fracturing. However, the fundamental mechanisms behind the particle transport and clogging behaviors remain to be fully understood. In this study, we combine microfluidic experiments and numerical simulations to investigate the particle transport and clogging phenomenon at the pore scale. We perform a large set of visualized experiments by systematically considering the impact of critical parameters such as particle volume fraction, particle size, and flow rate. Four particle clogging regimes, including non-clogging, depth filtration, caking-depth filtration, and caking, are identified. The influence of critical parameters on clogging area and clogging probability is quantified, revealing that larger particle sizes and higher particle volume fractions significantly increase both clogging area and probability. Additionally, results indicate that particle clogging alters its preferential flow paths within porous media. Pore scale direct numerical simulations reveal that the permeability reduction strongly depends on the particle clogging regime. In the caking and caking-depth filtration regimes, clogging can cause 2-3 orders of magnitude reduction in permeability. This work enhances our understanding of particle clogging in porous media and may provide insights on controlling particle clogging in porous media in practical applications.

Dams are often constructed on thick overburden deposits when a complete removal of the overburden becomes cost-prohibitive. This may induce an intense coupling between seepage and deformation in the deposits owing to their high permeability and poor mechanical properties. In this study, a coupled transient flow and nonlinear elastic deformation model is used to evaluate the design of the cutoff wall and foundation treatment for a concrete sluice dam of 42 m height located on 133-m-thick overburden. It is found that the depth of the cutoff wall has a significant effect on seepage and deformation control in the dam foundation, and a closed cutoff wall of 120 m depth performs much better than suspended ones with smaller depths in terms of the control of flow rate, uplift pressure, and seepage-induced uplift deformation. The backfill treatment of the top two overburden layers is crucial for controlling both flow rate and settlement, and the bored piles help reduce local settlement at the powerhouse foundation. Due to the relatively smaller height of the dam, the foundation deformation is less sensitive to permeability models. The results provide an important guidance for optimization design of the dam foundation treatments.

Seepage control is a critical issue for construction of underground powerhouse caverns, especially when the caverns are located near large-scale, water-conductive faults that provide flow channels into the cavern area. Grouting and draining have been considered as the most effective measures for regulating leakage into the caverns, which usually calls for an optimization design to achieve a balance between the performance and cost. This study proposes to optimize both the location of an underground cavern system and the design parameters (depth and/or spacing) of seepage control system located near a large-scale fault by numerical simulations. A comprehensive site characterization is performed to quantify the permeability of the surrounding rocks, and the excavation-induced permeability variation is characterized with a strain-dependent model. The groundwater flow is described by a steady-state flow model with unilateral boundary condition for rigorous simulation of drains. It is found that both the location of caverns and the layout of grout curtains and drains can be effectively optimized in terms of the metrics including the pore water pressure distribution, the discharge into the cavern area, and the stability of fault against seepage erosion. This work underscores the importance of optimization design in regulating the groundwater flow around underground caverns located near permeable faults.