Porosity Fields Research Articles

An alternative finite element formulation to predict ductile fracture in highly deformable materials

Abstract An alternative finite element formulation to predict ductile damage and fracture in highly deformable materials is presented. For this purpose, a finite-strain elastoplastic model based on the Gurson-Tvergaard-Needleman (GTN) formulation is employed, in which the level of damage is described by the void volume fraction (or porosity). The model accounts for large strains, associative plasticity and isotropic hardening, as well as void nucleation, coalescence and material failure. To avoid severe damage localization, a nonlocal enrichment is adopted, resulting in a mixed finite element whose degrees of freedom are the current positions and nonlocal porosity at the nodes. In this work, 2D triangular elements of linear-order and plane-stress conditions are used. Two systems of equations have to be solved: the global variables system, involving the degrees of freedom; and the internal variables system, including the damage and plastic variables. To this end, a new numerical strategy has been developed, in which the change in material stiffness due to the evolution of internal variables is embedded in the consistent tangent operator regarding the global system. The performance of the proposed formulation is assessed by three numerical examples involving large elastoplastic strains and ductile fracture. Results confirm that the present formulation is capable of reproducing fracture initiation and evolution, as well as necking instability. Convergence analysis is also performed in order to evaluate the effect of mesh refinement on the mechanical response. In addition, it is demonstrated that the nonlocal parameter alleviates damage localization, providing smoother porosity fields.

A model of mechanical loading of the lungs including gravity and a balancing heterogeneous pleural pressure.

Recent years have seen the development of multiple in silico lung models, notably with the aim of improving patient care for pulmonary diseases. These models vary in complexity and typically only consider the implementation of pleural pressure, a depression that keeps the lungs inflated. Gravity, often considered negligible compared to pleural pressure, has been largely overlooked, also due to the complexity of formulating physiological boundary conditions to counterbalance it. However, gravity is known to affect pulmonary functions, such as ventilation. In this study, we incorporated gravity into a recent lung poromechanical model. To do so, in addition to the gravitational body force, we proposed novel boundary conditions consisting in a heterogeneous pleural pressure field constrained to counterbalance gravity to reach global equilibrium of applied forces. We assessed the impact of gravity on the global and local behavior of the model, including the pressure-volume response and porosity field. Our findings reveal that gravity, despite being small, influences lung response. Specifically, the inclusion of gravity in our model led to the emergence of heterogeneities in deformation and stress distribution, compatible with in vivo imaging data. This could provide valuable insights for predicting the progression of certain pulmonary diseases by correlating areas subjected to higher deformation and stresses with disease evolution patterns.

A Novel Approach to Continuum FE Multiphysics Simulation of Batteries via Heterogeneous Homogenization

Advancements in modern technology necessitate batteries that can be safely charged and discharged at high rates with minimal aging. A promising strategy to achieve these ends is through optimization of porous microstructures of existing electrode materials to promote better electrical and ionic transport without compromising energy density, as has been shown experimentally by introducing laser ablation channels to the electrodes. The present work proposes a microstructure up-scaling methodology that efficiently accounts for material heterogeneity by assigning a porosity field to a continuum scale finite element mesh, effectively homogenizing materials while maintaining the heterogeneous porosity distribution inherent to the original microstructure. Constitutive relationships relating solid electric conductivity and pore tortuosity to local porosity are directly computed from fully resolved anode, cathode, and separator microstructure 3D models and assigned in the continuum as material properties of a finite element mesh. This allows for utilization of the computationally efficient, homogenized Newman model while maintaining dependency on the heterogeneous local porosity. Then, coupled electrochemical-thermal-mechanical-pore pressure finite element simulations are performed on the battery cell model with this heterogeneous porosity field, and associated local variations on transport properties are prescribed as material properties by the precomputed constitutive relationships. Performance metrics and lithium plating potential are compared to the results of simulations of the fully homogenized cell with the same constitutive relationships and to the base case of a fully homogenous cell with typical modelling assumptions applied to porosity dependent constitutive relationships. These comparisons are made both as a means of validating the methodology and to demonstrate the vast differences that arise when simplifying assumptions are used instead of the constitutive relationships of the actual material microstructure. Finally, the methodology is used in its intended optimization application to compare the performance of a typical battery microstructure to a similar porous microstructure with micro holes added via laser ablation. The results indicate that differences in the microstructure of materials with identical average porosities and similar constitutive relationships can have significant differences in both cell performance and aging that would typically be impossible to detect using fully homogeneous continuum modeling. Figure 1

An Improved Porosity Calculation Algorithm for Particle Flow Code

The widely used discrete-element particle flow software PFC’s (PFC 7.0 and previous versions) algorithm for calculating porosity is not sufficiently accurate. Because of this, when the particles are densely packed, the solution to the equation produces an algorithm exception for odd calculations of porosity, which results in the inability to calculate the results. This paper, based on a Darcy seepage model of fluid flow through a granular bed, analyzed the shortcomings of the two porosity calculation methods of PFC and the function analysis method. Combining this analysis with the theory of computer graphics, a new and efficient porosity calculation algorithm was proposed. The result showed that the new proposed porosity calculation algorithm calculated a more accurate and reasonable porosity field and made the iterative solution of the CFD equation more stable. This method makes porosity-related models of PFC more accurate. The algorithm can be not only used to calculate porosity, but also applied to other fields.

An improved weighted topology optimization lattice Boltzmann model for porous structures of advection–diffusion chemical reaction systems

An improved topology optimization lattice Boltzmann model for porous structures is proposed for advection–diffusion chemical reaction systems. Nonlinear porous medium drag models related to the mesoscale porosities and velocities are introduced for the momentum and mass transport coupled with the topology structure evolution. The evolution of the temporal-spatially varying porosity field is governed by a mesoscale porosity transfer equation. At each evolution step, the porous structure is sharpened by the porosity adjustment of a transient deviation of the objective function relative to porosity, and is slightly smoothed by the porosity diffusion simultaneously. Adjustable weights for chemical reaction, mass diffusion, kinetic energy, and viscous shears are applied for calculation of the transient deviation fields, which control the evolution speeds of the topology structures due to local chemical reaction, mass diffusion, flow convection, and interface drag, respectively. The effects of governing parameters and objective weights on structure patterns are illustrated. Three topology structure patterns with natural shapes are observed: (i) trunk shaped convective flow induced channels, (ii) fractal-branch shaped diffusion enhanced sub-channels, and (iii) leaf shaped chemical reaction fins. The trunk–branch–leaf topology structures result from the dynamic balance of the mesoscale convection, diffusion, and reaction.

A parameter-free and monolithic approach for multiscale simulations of flow, transport, and chemical reactions in porous media

Traditional approaches for transport in porous media rely on empirical hydrodynamic/transport parameters to quantify the interfacial exchange terms between fluid and solid, which leads to considerable uncertainties and also to a narrow application range of each method for realistic investigations of porous media. A parameter-free approach for flow, transport, and chemical reactions in porous media is presented in this paper. This approach adopts a novel, single-domain set of equations that are rigorously derived based on the elementary conservation laws describing the physical, continuous domain. The impact of the solid on nearby fluid being directly quantified in terms of the quantity to be solved, this avoids the introduction of empirical parameters. Various transport problems including external, conjugate, and reactive processes are properly formulated in this context, with the coupling between fluid and porous structures at different scales naturally solved in a monolithic manner. All these equations remain simple in terms of their numerical implementation, with solid structures represented by the porosity field on a background mesh. The application of the present method in both pore-scale and representative elementary volume (REV-scale) simulations are checked by performing a series of test cases, including interfacial momentum/heat/mass transport, conjugate heat/mass transfer, chemical reactions, and two practical applications containing reacting interfaces. An overall first-order convergence rate has been observed in both spatial and temporal accuracy tests.

Homogenized color-gradient lattice Boltzmann model for immiscible two-phase flow in multiscale porous media

In this paper, a homogenized multiphase lattice Boltzmann (LB) model is established for parallelly simulating immiscible two-phase flow in both solid-free regions (pore scale) and porous areas (continuum scale). It combines the color-gradient multiphase model with the Darcy–Brinkman–Stokes method by adding a term that includes surface force and drag force of porous matrix to multiple-relaxation-time LB equation in moment space. Moreover, an improved algorithm is proposed to characterize and implement the apparent wettability in the locally homogenized porosity field. Validations and test cases are given to demonstrate the accuracy and robustness of this new model, as well as its applicability for trans-scale fluid simulation of transport and sorption behavior from porous (Darcy flow) area to free (Stokes flow) area. For practicality, the two-phase seepage flow in a composite rock structure with multiscale pores is simulated by this new model, and the effects of viscosity ratio and wettability on the displacement process are discussed.

Upscaling reactive transport models from pore-scale to continuum-scale using deep learning method

Reactive transport modeling of subsurface environments plays an important role in addressing critical problems of geochemical processes, such as dissolution and precipitation of minerals. Current transport models for porous media span various scales, ranging from pore-scale to continuum-scale. In this study, we established an upscaling method connecting pore-scale and continuum-scale models by employing a deep learning methodology of Convolutional Neural Networks (CNNs). We applied Darcy-Brinkmann-Stokes (DBS) method to simulate the fluid flow and reactive transport in pore-scale models, which would act as constituents of a continuum-scale model. The datasets of spatial pore distribution of subvolume samples were used as the input for the deep learning model, and the continuum (Darcy)-scale parameters such as permeability, effective surface area, and effective diffusion coefficient were figured out as outputs (i.e., labels). By applying the trained models of the subvolumes in the entire sample volume, we generated the initial field of porosity, permeability, effective diffusion coefficient, and effective surface area for continuum-scale simulation of a mineral dissolution problem. We took an acid dissolution case as an example to utilize the outcomes of trained deep learning models as input data in the continuum-scale simulation. This work presents a comprehensive upscaling workflow, as bridging the findings of microscale simulations to the continuum-scale simulations of a reactive transport problem.

Modelling a gas injection experiment incorporating embedded fractures and heterogeneous material properties

This study focuses on the modelling of a gas injection experiment to assess the effects of incorporating heterogeneous material properties. The numerical model considers a two-phase flow coupled hydro-mechanical problem, and includes embedded fractures that open with deformation, thereby enhancing permeability. The approach used is integrated in the CODE_BRIGHT software, which allows for the consideration of geomaterials with a spatially correlated heterogeneous field of porosity that follows a normal distribution. This spatial correlation can be either isotropic or anisotropic. A key aspect of this approach is that material properties such as intrinsic permeability, diffusivity or cohesion are defined as a function of porosity. Consequently, these properties also exhibit heterogeneity with spatial correlation and, eventually, anisotropy. The results derived from the numerical model align well with in-situ measurements. The study also includes sensitivity analyses to the variation of critical variables. The calibration of the model has been validated through a similar experiment. The findings indicate that the consideration of heterogeneous material properties can have a significant influence on gas injection problems, particularly when a hydraulic fracture is formed.

Probabilistic Forecast of Multiphase Transport Under Viscous and Buoyancy Forces in Heterogeneous Porous Media

AbstractWe develop a probabilistic approach to map parametric uncertainty to output state uncertainty in first‐order hyperbolic conservation laws. We analyze this problem for nonlinear immiscible two‐phase transport in heterogeneous porous media in the presence of a stochastic velocity field. The uncertainty in the velocity field can arise from incomplete descriptions of either porosity field, injection flux, or both. This uncertainty leads to spatiotemporal uncertainty in the saturation field. Given information about spatial/temporal statistics of spatially correlated heterogeneity, we leverage the method of distributions to derive deterministic equations that govern the evolution of pointwise cumulative distribution functions (CDFs) of saturation for a vertical reservoir, while handling the manipulation of multiple shocks arising due to buoyancy forces. Unlike the Buckley‐Leverett equation, the equation describing the fine‐grained CDF is linear in space and time. Ensemble averaging of the fine‐grained CDF results in the CDF of saturation. Thus, we give routes to circumventing the computational cost of Monte Carlo simulations (MCS), while obtaining a pointwise description of the saturation field. We conduct a set of numerical experiments for one‐dimensional transport, and compare the obtained saturation CDFs, against those obtained using MCS as our reference solution, and the statistical moment equation method. This comparison demonstrates that the CDF equations remain accurate over a wide range of statistical properties, that is, standard deviation and correlation length of the underlying random fields, whereas the corresponding low‐order statistical moment equations significantly deviate from the MCS results, except for very small values of standard deviation and correlation length.

A novel semi‐resolved CFD‐DEM method with two‐grid mapping: Methodology and verification

AbstractThe semi‐resolved Computational Fluid Dynamics coupled with the Discrete Element Method (CFD‐DEM) method has emerged as approach to modeling particle‐fluid interactions in granular materials with high particle size ratios. However, challenges arise from conflicting requirements regarding the CFD grid size, which must adequately resolve fluid flow in the pore space while maintaining a physically meaningful porosity field. This study addresses these challenges by introducing a two‐grid mapping approach. Initially, the porosity field associated with fine particles is estimated using a coarse CFD grid, which is then mapped to a dynamically refined grid. To ensure conservation of total solid volume, a volume compensation procedure is implemented. The proposed method has been rigorously verified using benchmark cases, showing its high computational efficiency and accurate handling of complex porosity calculations near the surface of coarse particles. Moreover, the previously unreported impact of the empirical drag correlation on fluid‐particle force calculations for both coarse and fine particles has been revealed.

Augmenting Deep Residual Surrogates with Fourier Neural Operators for Rapid Two-Phase Flow and Transport Simulations

Summary Accurate numerical modeling of multiphase flow and transport mechanisms is essential to study varied, complex physical phenomena including flow in subsurface oil and gas reservoirs and subsurface aquifers subject to CO2 sequestration. State-of-the-art complete physics-based solvers suffer from many computational challenges. High-fidelity data-driven surrogate models that solve the governing partial differential equations (PDEs) have the potential to optimize the time to solution and increase confidence in critical business and engineering decisions through better quantification of solution statistics. We leverage the recently proposed Fourier neural operators (FNOs) with quasilinear time complexity to capture the spectral information from feature maps to solve the coupled porous flow and transport PDEs. Embedding Fourier layers within the residual blocks results in a highly effective structure that, while achieving competitive accuracy, also enables efficient training of deeper networks with a dramatically reduced number of trainable parameters. The resulting novel deep-learning (DL) architecture is coined as FResNet++. FResNet++ uses squeeze and excitation blocks, atrous spatial pyramid pooling (ASPP), and attention blocks to increase its sensitivity to the relevant features and capture multiscale information, and it is specifically tuned to operate optimally to learn from and predict numerically simulated flow (pressure and saturation) fields. We demonstrate the ability of FResNet++ to generalize over multiple high-dimensional input parameter spaces that describe subsurface permeability and porosity heterogeneity. The resulting DL architecture accurately captures the complex interplay between viscous forces and highly heterogeneous permeability and porosity fields. We investigate two-phase flow in porous media, which is the archetypal problem for reservoir simulation giving rise to a system of nonlinearly coupled PDEs with highly heterogeneous coefficients. We show in blind tests that FResNet++ predicts saturation fields more accurately compared to ResU-Net and original FNO with fully connected linear layers. We additionally investigate the effects of using alternative loss functions and an alternative way of utilizing FResNet++ to increase its effectiveness. For the first time in the literature, we show that the spatiotemporal evolution of pressure and saturation fields can be jointly predicted with good accuracy using a single FResNet++ network over long time horizons in response to previously unseen permeability and porosity fields. After a moderate training investment on graphics processing units (GPUs), FResNet++ yields a speedup of at least four orders of magnitude compared to a conventional numerical PDE solver and operates with notably fewer trainable parameters compared to the original FNO. Our numerical experiments validate that FNOs can be utilized in various convolutional neural network-based architectures and can effectively substitute for repetitive physics-based forward simulations for scenario testing.

Nonlinear Ultrasonic Imaging for Porosity Evaluation.

The influence of porosity on the mechanical behaviour of composite laminates represents a complex problem that involves many variables. Therefore, the evaluation of the type and volume content of porosity in a composite specimen is important for quality control and for predicting material behaviour during service. A suitable way to evaluate the porosity content in composites is by using nonlinear ultrasonics because of their sensitivity to small cracks. The main objective of this research work is to present an imaging method for the porosity field in composites. Two nonlinear ultrasound techniques are proposed using backscattered signals acquired by a phased array system. The first method was based on the amplitude of the half-harmonic frequency components generated by microbubble reflections, while the second one involved the frequency derivative of the attenuation coefficient, which is proportional to the porosity content in the specimen. Two composite samples with induced porosity were considered in the experimental tests, and the results showed the high accuracy of both methods with respect to a classic C-scan baseline. The attenuation coefficient results showed high accuracy in defining bubble shapes in comparison with the half-harmonic technique when surface effects were neglected.

Solving Geophysical Inversion Problems with Intractable Likelihoods: Linearized Gaussian Approximations Versus the Correlated Pseudo-marginal Method

A geophysical Bayesian inversion problem may target the posterior distribution of geological or hydrogeological parameters given geophysical data. To account for the scatter in the petrophysical relationship linking the target parameters to the geophysical properties, this study treats the intermediate geophysical properties as latent (unobservable) variables. To perform inversion in such a latent variable model, the intractable likelihood function of the (hydro)geological parameters given the geophysical data needs to be estimated. This can be achieved by approximation with a Gaussian probability density function based on local linearization of the geophysical forward operator, thereby, accounting for the noise in the petrophysical relationship by a corresponding addition to the data covariance matrix. The new approximate method is compared against the general correlated pseudo-marginal method, which estimates the likelihood by Monte Carlo averaging over samples of the latent variable. First, the performances of the two methods are tested on a synthetic test example, in which a multivariate Gaussian porosity field is inferred using crosshole ground-penetrating radar first-arrival travel times. For this example with rather small petrophysical uncertainty, the two methods provide near-identical estimates, while an inversion that ignores petrophysical uncertainty leads to biased estimates. The results of a sensitivity analysis are then used to suggest that the linearized Gaussian approach, while attractive due to its relative computational speed, suffers from a decreasing accuracy with increasing scatter in the petrophysical relationship. The computationally more expensive correlated pseudo-marginal method performs very well even for settings with high petrophysical uncertainty.

A new approach to model geomaterials with heterogeneous properties in thermo-hydro-mechanical coupled problems

The main objective of this article is to present a new approach to model coupled thermo-hydro-mechanical problems considering geomaterials with heterogeneous properties. This approach has been implemented in the software CODE_BRIGHT and it provides the possibility of considering geomaterials with a spatially correlated heterogeneous field of porosity, following a normal distribution. This spatial correlation can be isotropic or anisotropic. An important feature of this approach is that material properties such as intrinsic permeability, thermal conductivity, diffusivity, retention curve, elastic modulus or cohesion are defined as a function of porosity and, thus, they become heterogeneous with spatial correlation and, eventually, anisotropic. A validation exercise and other basic numerical examples have been carried out to illustrate the possibilities of the proposed approach. The results, which have been compared with a homogeneous case, show that considering heterogeneous fields can be relevant in different modelling problems, especially coupled thermo-hydro-mechanical problems.

An adjoint‐based solver with adaptive mesh refinement for efficient design of coupled thermal‐fluid systems

AbstractA multi‐objective continuous adjoint strategy based on the superposition of boundary functions for topology optimization of problems where the heat transfer must be enhanced and the dissipated mechanical power controlled at the same time, has been here implemented in a finite volume (FV), incompressible, steady flow solver supporting a dynamic adaptive mesh refinement (AMR) strategy. The solver models the transition from fluid to solid by a porosity field, that appears in the form of penalization in the momentum equation; the material distribution is optimized by the method of moving asymptotes (MMA). AMR is based on a hierarchical nonconforming h‐refinement strategy and is applied together with a flux correction to enforce conservation across topology changes. It is shown that a proper choice of the refinement criterium favors a mesh‐independent solution. Finally, a Pareto front built from the components of the objective function is used to find the best combination of the weights in the optimization cycle. Numerical experiments on two‐ and three‐dimensional test cases, including the aero‐thermal optimization of a simplified layout of a cooling system, have been used to validate the implemented methodology.

Analytic circuit model for thermal drying behavior of electronic inks

Understanding the sintering process of conductive inks is a fundamental step in the development of sensors. The intrinsic properties (such as thermal conductivity, resistivity, thermal coefficient, among others) of the printed devices do not correspond to those of the bulk materials. In the field of biosensors porosity plays a predominant role, since it defines the difference between the geometric area of the working electrode and its electrochemical surface area. The analysis reported so far in the literature on the sintering of inks are based on their DC characterization. In this work, the shape and distribution of the nanoparticles that make up the silver ink have been studied employing a transmission electron microscopy. Images of the printed traces have been obtained through a scanning electron microscope at different sintering times, allowing to observe how the material decreases its porosity over time. These structural changes were supported through electrical measurements of the change in the trace impedance as a function of drying time. The resistivity and thermal coefficient of the printed tracks were analyzed and compared with the values of bulk silver. Finally, this work proposes an analytical circuit model of the drying behavior of the ink based on AC characterization at different frequencies. The characterization considers an initial time when the spheric nanoparticles are still surrounded by the capping agent until the conductive trace is obtained. This model can estimate the characteristics that the printed devices would have, whether they are used as biosensors (porous material) or as interconnections (compact material) in printed electronics.

The role of different sources of uncertainty on the stochastic quantification of subsurface discharges in heterogeneous aquifers

The quantification of subsurface discharges may significantly be affected by multiple sources of uncertainty, especially in highly heterogeneous aquifers. In this work, subsurface discharges within alluvial aquifers located in the province of Lecco (Lombardy, Northern Italy) are quantified relying upon a Monte Carlo (MC) framework under geological and conditioning data uncertainty sources. These are respectively employed through different geological conceptual models, whose facies spatial distributions are simulated through Sequential Indicator Simulations methodology (SISIM), and by conditioning stochastic facies simulations with different data constraints. Stochastic porosity fields are then achieved upon facies simulated fields by associating average porosity values to each local simulated facies. Then, stochastic subsurface discharge fluxes are quantified under steady hydraulic conditions upon saturated hydraulic fields, inferred from porosity ones. Then, the study of subsurface discharges is pursued via the Generalized Extreme Values (GEV) approach. This method enables both to interpret peak-over-thresholds (POTs) discharges spatial distribution, and to study the role geological and data uncertainty sources on their quantification, for example studying discharges conveyed towards a lake through the surrounding aquifer. These uncertainty-based analyses show that subsurface discharges clearly convey in correspondence of main alluvial aquifers, as well as the tendency to be collected along preferential pathways. Furthermore, the inspection of POTs empirical probability spatial distributions highlights that the geological uncertainty plays as a primary uncertainty source rather than the conditioning data. The latter plays a minor role, e.g. between different employed models, only in presence of a relevant geological difference between analyzed ones. These results are consistent with other kinds of information, such as hydrogeological sections and available piezometric maps, providing novel contributions to understanding how to quantify subsurface discharges in poorly monitored aquifers or how to treat various sources of uncertainty.

Spectral stochastic isogeometric analysis of bending and free vibration of porous functionally graded plates

In this paper, a systematic spectral stochastic isogeometric analysis (SSIGA) process is presented for the static bending and free vibration analyses of functionally graded (FG) plates with three-dimensional (3D) random porosity. The porosity is modeled as a Beta random field, represented compactly via the Karhunen-Loève expansion. A novel hierarchical locking-free quasi-3D shear deformation theory, called spectral displacement formulation (SDF), is proposed to approach exact 3D solutions and reflect more realistic effects of the random porosity field (RPF). Isogeometric analysis is utilized to meet the C1-continuity requirement of the SDF. The response surfaces of the porous FG plates are constructed non-invasively by the spectral collocation method. A new spectral stochastic post-processing process is developed to evaluate the probability characteristics of the responses and exclude the adverse convergence-in-probability property, which typically exists in sampling-based methods. Numerical examples illustrate the SSIGA process and demonstrate its effectiveness. The influences of the RPF parameters and the gradient index on the response statistics are investigated.

On the parameterisation of dual-continuum models for reactive transport modelling in fractured media

In sparsely fractured rocks, the rock matrix is an important geochemical buffer and provides significant retardation to contaminants advected through the flowing fractures. Accounting for geochemical reactions and mass exchange between these two regions is key to properly capture the overall buffering capacity and the related hydrogeochemical evolution of a fractured medium. Reactive transport modelling in these kinds of fractured media is routinely performed using Equivalent Continuous Porous Media (ECPM) models: i.e. continuum models based on permeability and porosity fields that somehow preserve the underlying fracture properties, which are in turn described by companion Discrete Fracture Network (DFN) models. However, the proper parameterisation of these models, in terms of mass exchange between fractures and the bordering matrix, is still a largely unresolved issue. Here, we leverage the Dual Continuum Disconnected Matrix Model (DCDMM) formulation included in the massively parallel code PFLOTRAN to propose a novel parameterisation approach that honours the local volumetric fracture density (P32). The proposed approach is first benchmarked against a semi-analytical solution with a problem that entails flow and transport along two different and consecutive fractures. Two demonstrative large-scale reactive transport problems are also presented and discussed: the first is related to the generation and migration of radiogenic helium and the second assesses the buffering capacity of a realistic fractured medium against the infiltration of acidic water. The latter simulation, which includes more than three hundred million transport degrees of freedom is one of the largest subsurface reactive transport models ever formulated and solved. This simulation was made possible by the highly efficient implementation of the DCDMM in PFLOTRAN, which makes the solution of the secondary continuum equations embarrassingly parallel.