Computational Transport Phenomena Lab
For more information visit: https://ctpl.kaust.edu.sa/Pages/Home.aspx
Recent Submissions

Base level changes based on Basin Filling Modelling: a case study from the Paleocene Lishui Sag, East China Sea Basin(Petroleum Science, Springer Science and Business Media LLC, 20200717) [Article]Estimation of base level changes in geological records is an important topic for petroleum geologists. Taking the Paleocene Upper Lingfeng Member of Lishui Sag as an example, this paper conducted a base level reconstruction based on Basin Filling Modelling (BFM). The reconstruction was processed on the ground of a previously interpreted seismic stratigraphic framework with several assumptions and simplification. The BFM is implemented with a nonlinear diffusion equation solver written in R coding that excels in shallow marine stratigraphic simulation. The modeled results fit the original stratigraphy very well. The BFM is a powerful tool for reconstructing the base level, and is an effective way to check the reasonableness of previous interpretations. Although simulation solutions may not be unique, the BFM still provides us a chance to gain some insights into the mechanism and dynamic details of the stratigraphy of sedimentary basins.

Fully discrete energy stable scheme for a phasefield moving contact line model with variable densities and viscosities(Applied Mathematical Modelling, Elsevier BV, 20200304) [Article]In this study, we propose a fully discrete energy stable scheme for the phasefield moving contact line model with variable densities and viscosities. The mathematical model comprises a Cahn–Hilliard equation, Navier–Stokes equation, and the generalized Navier boundary condition for the moving contact line. A scalar auxiliary variable is employed to transform the governing system into an equivalent form, thereby allowing the double well potential to be treated semiexplicitly. A stabilization term is added to balance the explicit nonlinear term originating from the surface energy at the fluid–solid interface. A pressure stabilization method is used to decouple the velocity and pressure computations. Some subtle implicit–explicit treatments are employed to deal with convention and stress terms. We establish a rigorous proof of the energy stability for the proposed timemarching scheme. A finite difference method based on staggered grids is then used to spatially discretize the constructed timemarching scheme. We also prove that the fully discrete scheme satisfies the discrete energy dissipation law. Our numerical results demonstrate the accuracy and energy stability of the proposed scheme. Using our numerical scheme, we analyze the contact line dynamics based on a shear flowdriven droplet sliding case. Threedimensional droplet spreading is also investigated based on a chemically patterned surface. Our numerical simulation accurately predicts the expected energy evolution and it successfully reproduces the expected phenomena where an oil droplet contracts inward on a hydrophobic zone and then spreads outward rapidly on a hydrophilic zone.

Sorption and Diffusion of Methane and Carbon Dioxide in Amorphous Poly(alkyl acrylates): A Molecular Simulation Study(The Journal of Physical Chemistry B, American Chemical Society (ACS), 20200129) [Article]Molecular simulations were carried out to understand the structural features and the sorption and diffusion behavior of methane and carbon dioxide in amorphous poly(alkyl acrylates) in the temperature range of 300600 K. The hybrid Monte Carlo/molecular dynamics approach was employed to address the effects of polymer swelling and framework flexibility on the gas sorption. Simulations show that the glasstransition temperature decreases with the sidechain length of poly(alkyl acrylate), consistent with experiments. This is due to the fact that the shielding of the polar ester groups increases with the sidechain length. The simulated sorption isotherms for methane and carbon dioxide were in agreement with the experimental data. The polymer swelling becomes more pronounced, especially in the case of sorption of carbon dioxide. A significant swelling occurs, possibly because of the greater interaction between carbon dioxide and the polar ester groups in the polymers. The uptake of methane and carbon dioxide by poly(alkyl acrylates) generally increases with the sidechain length. Our simulations confirm the experimental findings that the diffusion coefficients of methane and carbon dioxide in poly(alkyl acrylates) increase with the sidechain length. Interestingly, the activation energies of gas diffusion decrease with the sidechain length. The diffusion coefficients of the penetrants have an exponential relationship with the void fraction, which is in agreement with the free volume theory.

A Novel Energy Factorization Approach for the DiffuseInterface Model with PengRobinson Equation of State(SIAM Journal on Scientific Computing, Society for Industrial & Applied Mathematics (SIAM), 20200107) [Article]The PengRobinson equation of state (PREoS) has become one of the most extensively applied equations of state in chemical engineering and the petroleum industry due to its excellent accuracy in predicting the thermodynamic properties of a wide variety of materials, especially hydrocarbons. Although great effort has been made to construct efficient numerical methods for the diffuse interface models with PREoS, there is still not a linear numerical scheme that can be proved to preserve the original energy dissipation law. In order to pursue such a numerical scheme, we propose a novel energy factorization (EF) approach, which first factorizes an energy function into a product of several factors and then treats the factors using their properties to obtain the semiimplicit linear schemes. We apply the EF approach to deal with the Helmholtz free energy density determined by PREoS, and then propose a linear semiimplicit numerical scheme that inherits the original energy dissipation law. Moreover, the proposed scheme is proved to satisfy the maximum principle in both time semidiscrete form and cellcentered finite difference fully discrete form under certain conditions. Numerical results are presented to demonstrate the stability and efficiency of the proposed scheme.

Special Issue: Advanced numerical modeling and algorithms for multiphase flow and transport(Journal of Computational Physics, Elsevier BV, 20200106) [Article]

Bulk and interfacial properties of decane in the presence of carbon dioxide, methane, and their mixture.(Scientific reports, Springer Science and Business Media LLC, 20191224) [Article]Molecular dynamics simulations were performed to study the bulk and interfacial properties of methane + ndecane, carbon dioxide + ndecane, and methane + carbon dioxide + ndecane systems under geological conditions. In addition, theoretical calculations using the predictive PengRobinson equation of state and density gradient theory are carried out to compare with the simulation data. A key finding is the preferential dissolution in the decanerich phase and adsorption at the interface for carbon dioxide from the methane/carbon dioxide mixture. In general, both the gas solubility and the swelling factor increase with increasing pressure and decreasing temperature. Interestingly, the methane solubility and the swelling of the methane + ndecane system are not strongly influenced by temperature. Our results also show that the presence of methane increases the interfacial tension (IFT) of the carbon dioxide + ndecane system. Typically, the IFT of the studied systems decreases with increasing pressure and temperature. The relatively higher surface excess of the carbon dioxide + ndecane system results in a steeper decrease in its IFT as a function of pressure. Such systematic investigations may help to understand the behavior of the carbon dioxideoil system in the presence of impurities such as methane for the design and operation of carbon capture and storage and enhanced oil recovery processes.

A locally and globally phasewise mass conservative numerical algorithm for the twophase immiscible flow problems in porous media(Computers and Geotechnics, Elsevier BV, 20191214) [Article]In this work, we introduce a novel numerical method to solve the problem of twophase immiscible flow in porous media that is conservative to both phases. In the widely used implicit pressure, explicit saturation (IMPES) scheme, the conservation of mass of both the two phases are summed to form an equation involving the total Darcy's velocity. In the discretization of such an equation it becomes difficult to enforce the conservation of mass of each phase. To guarantee the conservation of mass of both phases locally and hence globally, we introduce a scheme in which the time discretization of the mass conservation equations is considered separately. Cellcentered finite difference (CCFD) methods are adopted for spatial discretization, where the variables of fluid properties (i.e. relative permeability and mobility) are upwinded separately according to the velocity of each phase and not according to the total velocity. Furthermore, this new scheme updates all phase velocities and uses them to update the corresponding phase saturation. In addition, a twoscale of timesplitting methods are adopted for pressure equation and saturation equations to improve the computational efficiency. For the sake of simplicity, we show a number of examples of twophase system in twodimensional geometry solved using the new scheme. It is shown that the new scheme is more embracing the physics and it can be more accurate than traditional IMPES scheme, particularly for the cases in which the phase velocities are in opposite direction, and conventional IMPES schemes fail.

A phasefield moving contact line model with soluble surfactants(Journal of Computational Physics, Elsevier BV, 20191204) [Article]A phasefield moving contact line model is presented for a twophase system with soluble surfactants. With the introduction of some scalar auxiliary variables, the original free energy functional is transformed into an equivalent form, and then a new governing system is obtained. The resulting model consists of two CahnHilliardtype equations and incompressible NavierStokes equation with variable densities, together with the generalized Navier boundary condition for the moving contact line. We prove that the proposed model satisfies the total energy dissipation with time. To numerically solve such a complex system, we develop a nonlinearly coupled scheme with unconditional energy stability. A splitting method based on pressure stabilization is used to solve the NavierStokes equation. Some subtle implicitexplicit treatments are adopted to discretize convection and stress terms. A stabilization term is artificially added to balance the explicit nonlinear term associated with the surface energy at the fluidsolid interface. We rigorously prove that the proposed scheme can preserve the discrete energy dissipation. An efficient finite difference method on staggered grids is used for the spatial discretization. Numerical results in both two and three dimensions demonstrate the accuracy and energy stability of the proposed scheme. Using our model and numerical scheme, we investigate the wetting behavior of droplets on a solid wall. Numerical results indicate that surfactants can affect the wetting properties of droplet by altering the value of contact angles.

Generalized multiscale approximation of mixed finite elements with velocity elimination for subsurface flow(Journal of Computational Physics, Elsevier BV, 20191125) [Article]A frame work of the mixed generalized multiscale finite element method (GMsFEM) for solving Darcy's law in heterogeneous media is studied in this paper. Our approach approximates pressure in multiscale function space that is between finegrid space and coarsegrid space and solves velocity directly in the finegrid space. To construct multiscale basis functions for each coarsegrid element, three types of snapshot space are raised. The first one is taken as the finegrid space for pressure and the other two cases need to solve a local problem on each coarsegrid element. We describe a spectral decomposition in the snapshot space motivated by the analysis to further reduce the dimension of the space that is used to approximate the pressure. Since the velocity is directly solved in the finegrid space, in the linear system for the mixed finite elements, the velocity matrix can be approximated by a diagonal matrix without losing any accuracy. Thus it can be inverted easily. This reduces computational cost greatly and makes our scheme simple and easy for application. Comparing to our previous work of mixed generalized multiscale finite element method [E. T. Chung, Y. Efendiev, and C. S. Lee. Mixed generalized multiscale finite element methods and applications. Multiscale Modeling and Simulation, 13(1):338366, 2015.], both the pressure and velocity space in this approach are bigger. As a consequence, this method enjoys better accuracy. While the computational cost does not increase because of the good property of velocity matrix. Moreover, the proposed method preserves the local mass conservation property that is important for subsurface problems. Numerical examples are presented to illustrate the good properties of the proposed approach. If offline spaces are appropriately selected, one can achieve good accuracy with only a few basis functions per coarse element according to the numerical results.

Recent Progress on Phase Equilibrium Calculation in Subsurface Reservoirs Using Diffuse Interface Models(Springer International Publishing, 20191117) [Book Chapter]Compositional multiphase flow in subsurface porous media is becoming increasingly attractive due to issues related with enhanced oil recovery, greenhouse effect and global warming, and the urgent need for development in unconventional oil/gas reservoirs. One key effort prior to construct the mathematical model governing the compositional multiphase flow is to determine the phase compositions of the fluid mixture, and then calculate other related physical properties. In this paper, recent progress on phase equilibrium calculations in subsurface reservoirs have been reviewed and concluded with authors’ own analysis. Phase equilibrium calculation is the main approach to perform such calculation, which could be conducted using two different types of flash calculation algorithms: The NPT flash and NVT flash. NPT flash calculations are proposed early, well developed within the last few decades and now become the most commonly used method. However, it fails to remain the physical meanings in the solution as a cubic equation, derived from equation of state, is often needed to solve. Alternatively, NVT flash can handle the phase equilibrium calculations as well, without the pressure known a priori. Recently, Diffuse Interface Models, which were proved to keep a high consistency with thermodynamic laws, have been introduced in the phase calculation, incorporating the realistic equation of state (EOS), e.g. PengRobinson EOS. In NVT flash, Helmholtz free energy is minimized instead of Gibbs free energy used in NPT flash, and this energy density is treated with convexconcave splitting technique. A semiimplicit numerical scheme is designed to process the dynamic model, which ensures the thermodynamic stability and then preserve the fast convergence property. A positive definite coefficient matrix is designed to meet the Onsager Reciprocal Principle so as to keep the entropy increasing property in the presence of capillary pressure, which is required by the thermodynamic laws. The robustness of the proposed algorithm is verified via two numerical examples, one of which has up to seven components. In the complex fluid mixture, special phenomena could be capture from the global minimum of TPD functions as well as the phase envelope resulted from the phase equilibrium calculations. It can be found that the boundary between the singlephase and vapor–liquid phase regions will move in the presence of capillary pressure, and then the area of each region will change accordingly. Some remarks have been concluded at the end, as well as suggestions on potential topics for future studies.

Layer Charge Effects on Adsorption and Diffusion of Water and Ions in Interlayers and on External Surfaces of Montmorillonite(ACS Earth and Space Chemistry, American Chemical Society (ACS), 20191028) [Article]Molecular simulations using classical force fields were performed to study the adsorption and diffusion properties of water and ions in the highcharge (Arizonatype) montmorillonite clays with varying relative humidity (RH) at 298.15 K. The simulation results of basal distances derived from swelling free energy curves and of water uptake are in good agreement with experiments. Overall, the simulated selfdiffusion coefficients of the interlayer species are in reasonable agreement with experiments and lower than those estimated for the external surfaces. Influence of the magnitude of the layer charge was studied by comparing these simulation results with those obtained for the lowcharge (Wyomingtype) montmorillonite clays. Most importantly, these comparisons confirm the experimental finding that the highcharge clay generally shifts swelling transitions toward lower RH values. Therefore, the adsorption and dynamics of water and ions are significantly different in the low A nd highcharge clays near the transition RH values. We find that the amount of water in an interlayer hydration state is mostly independent of the layer charge, probably due to steric reasons. In contrast, the externally adsorbed water content increases with increasing layer charge. Furthermore, the mobility of water and ions is generally lower in the highcharge montmorillonite than in the corresponding lowcharge system. However, mobility of cations in the mesopores of highcharge montmorillonite is typically higher than that of the corresponding lowcharge system which might be attributed to the presence of the tetrahedral substitutions in the latter case.

Low salinity effect on the recovery of oil trapped by nanopores: A molecular dynamics study(Fuel, Elsevier BV, 20191025) [Article]Low salinity waterflooding (LSW) is an effective method for enhancing the oil recovery from many reservoirs, and its success has been traced to a host of low salinity effects. In this work, we perform molecular dynamics simulations to study the feasibility of recovering oil trapped by nanopores by lowering the reservoir salinity. The oil is initially trapped by a slit nanopore, with a portion of the oil protruding from the pore entrance. After the reservoir salinity is lowered, the thin brine films that separate the oil and pore walls become thicker to drive some of the trapped oil out of the pore. We quantify the free energy profile of this process and clarify the underlying molecular mechanisms. Interestingly, the brine film growth is dominated by the water transport from the brine reservoir into the pore rather than by the depletion of ions from the brine film. These results provide molecular evidence that low salinity brines benefit the recovery of the oil trapped by nanopores. They highlight that when ion depletion from thin brine films is suppressed, the osmosis of water can play a fundamental role in the expansion of the brine films; thus, the enhanced oil recovery. The slow osmosis of water through thin brine films and thus the slow displacement of oil from the pore may help explain the anomalously slow oil recovery reported in micromodeling experiments of LSW.

A full multigrid multilevel Monte Carlo method for the single phase subsurface flow with random coefficients(arXiv, 20191010) [Preprint]The subsurface flow is usually subject to uncertain porous media structures. In most cases, however, we only have partial knowledge about the porous media properties. A common approach is to model the uncertain parameters of porous media as random fields, then the statistical moments (e.g. expectation) of the Quantity of Interest(QoI) can be evaluated by the Monte Carlo method. In this study, we develop a full multigridmultilevel Monte Carlo (FMGMLMC) method to speed up the evaluation of random parameters effects on singlephase porous flows. In general, MLMC method applies a series of discretization with increasing resolution and computes the QoI on each of them. The effective variance reduction is the success of the method. We exploit the similar hierarchies of MLMC and multigrid methods and obtain the solution on coarse mesh $Q^c_l$ as a byproduct of the full multigrid solution on fine mesh $Q^f_l$ on each level $l$. In the cases considered in this work, the computational saving due to the coarse mesh samples saving is $20\%$ asymptotically. Besides, a comparison of Monte Carlo and QuasiMonte Carlo (QMC) methods reveals a smaller estimator variance and a faster convergence rate of the latter approach in this study.

Thermodynamically consistent modelling of twophase flows with moving contact line and soluble surfactants(Journal of Fluid Mechanics, Cambridge University Press (CUP), 20190927) [Article]Droplet dynamics on a solid substrate is significantly influenced by surfactants. It remains a challenging task to model and simulate the moving contact line dynamics with soluble surfactants. In this work, we present a derivation of the phasefield moving contact line model with soluble surfactants through the first law of thermodynamics, associated thermodynamic relations and the Onsager variational principle. The derived thermodynamically consistent model consists of two CahnHilliard type of equations governing the evolution of interface and surfactant concentration, the incompressible NavierStokes equations and the generalized Navier boundary condition for the moving contact line. With chemical potentials derived from the free energy functional, we analytically obtain certain equilibrium properties of surfactant adsorption, including equilibrium profiles for phasefield variables, the Langmuir isotherm and the equilibrium equation of state. A classical droplet spread case is used to numerically validate the moving contact line model and equilibrium properties of surfactant adsorption. The influence of surfactants on the contact line dynamics observed in our simulations is consistent with the results obtained using sharp interface models. Using the proposed model, we investigate the droplet dynamics with soluble surfactants on a chemically patterned surface. It is observed that droplets will form three typical flow states as a result of different surfactant bulk concentrations and defect strengths, specifically the coalescence mode, the noncoalescence mode and the detachment mode. In addition, a phase diagram for the three flow states is presented. Finally, we study the unbalanced Young stress acting on triplephase contact points. The unbalanced Young stress could be a driving or resistance force, which is determined by the critical defect strength.

Darcyscale phase equilibrium modeling with gravity and capillarity(Journal of Computational Physics, Elsevier BV, 20190905) [Article]The modeling of multiphase fluid mixture and its flow in porous media is of great interest in the field of reservoir simulation. In this paper, we formulate a novel energybased framework to model multicomponent twophase fluid systems at equilibrium. PengRobinson equation of state (EOS) is used to model the bulk properties of each phase, though our framework works well also with other equations of state. Our model reduces to the conventional compositional grading if restricted to one spatial vertical dimension together with the assumption of monodisperse poresize distribution (all pores being one size). However, our model can be combined with a general distribution of pore size, which can generate interesting behaviors of capillarity in porous media. In particular, the model can be used to predict the capillary pressure of twophase fluid as a function of saturation, with a given poresize distribution. This model is the quantitative study of the first time in the literature for the capillarity of a twophase fluid with partial miscibility. We proposed an unconditionalstable energydecay numerical algorithm based on convexconcave splitting, which has been demonstrated to be both robust and efficient using numerical examples. To verify our model, we simulate the compositional grading of a binary fluid mixture consisting of carbon dioxide and normal decane. To demonstrate powerful features of our model, we provide an interesting example of fluid mixture in a porous medium with wide pore size distribution, where the competition of capillarity and gravity is observed. This work represents the first effort in the literature that rigorously incorporates capillarity and gravity effects into EOSbased phase equilibrium modeling.

Heat and mass transfer in a viscous nanofluid containing a gyrotactic microorganism over a stretching cylinder(Symmetry, MDPI AG, 20190905) [Article]This work consists of a theoretical boundary layer analysis of heat and mass transport in a viscous fluidembracing gyrotactic microorganism over a cylinder. The flow governing equations are modeled through boundary layer approximations. The governing nonlinear partial differential equations are lessened to a set of nonlinear ordinary differential equations using similitude transformation. The boundary layer equations are elucidated numerically, applying the spectral relaxation method with the aid of the computational software MATLAB. The impact of several pertinent parameters on flow convective characteristic phenomena are explored through the use of graphs and tables and are discussed with indepth physical descriptions. In addition, the friction factor, the rate of heat transfer, rate of mass transfer, and the density number of the motile microorganism are also presented with respect to the above controlled parameters. It is noticed that for the increasing values of the magnetic parameter with reductions and enhancements, the density of the motile microorganism is a declining function of, and the concentration field enhances with the strengthening of, whereas it reduces with the rise of. Furthermore, the streamline patters are emphasized for the impact of controlled flow variables. Current outcomes are compared with the available results from previous cases and are observed to be in agreement.

Controlling Factors of Degassing in Crosslinked Polyethylene Insulated Cables(Polymers, MDPI AG, 20190902) [Article]Here, we analyze the degassing process of a byproduct (methane) formed during the peroxideinduced crosslinking of polyethylene. A diffusion model based on Fick’s law is used to obtain the controlling factors of degassing in a crosslinked polyethylene (XLPE) insulated power cable (132 kV with 18 mm of insulation). We quantitatively analyze different scenarios of the diffusion of methane through the XLPE insulation and two semiconductor layers under various in situ degassing conditions. The analyzed degassing conditions include heat transfer and its effect on the diffusion properties, the different transport and boundary conditions due to the free spaces within the cable conductor, and the nonuniform distribution of methane concentrations within the insulation layers. Our simulation results clearly demonstrate that the free spaces between the copper strands in the cable conductor significantly affect the degassing efficiency. However, the temperaturediffusion coupling has a relatively minor effect on the overall degassing efficiency due to the rapid temperature increase of the polymer layers during the initial stages of degassing. Moreover, we find that the nonuniform distribution of methane in the initial stages also plays an important role in degassing in the cable, but this effect varies significantly during the degassing process.

Advances in Gaussian random field generation: a review(Computational Geosciences, Springer Science and Business Media LLC, 20190805) [Article]Gaussian (normal) distribution is a basic continuous probability distribution in statistics, it plays a substantial role in scientific and engineering problems that related to stochastic phenomena. This paper aims to review stateoftheart of Gaussian random field generation methods, their applications in scientific and engineering issues of interest, and opensource software/packages for Gaussian random field generation. To this end, first, we briefly introduce basic mathematical concepts and theories in the Gaussian random field, then seven commonly used Gaussian random field generation methods are systematically presented. The basic idea, mathematical framework of each generation method are introduced in detail and comparisons of these methods are summarized. Then, representative applications of the Gaussian random field in various areas, especially of engineering interest in recent two decades, are reviewed. For readers’ convenience, four representative example codes are provided, and several relevant uptodate opensource software and packages that freely available from the Internet are introduced.

Structure, Thermodynamics, and Dynamics of Thin Brine Films in OilBrineRock Systems.(Langmuir : the ACS journal of surfaces and colloids, American Chemical Society (ACS), 20190802) [Article]Thin brine films are ubiquitous in oilbrinerock systems such as oil reservoirs and play a crucial role in applications such as enhanced oil recovery. We report the results of molecular simulations of brine films that are confined between model oil (ndecane) and rock (neutral or negatively charged quartz slabs), with a focus on their structure, electrical double layers (EDLs), disjoining pressure, and dynamics. As brine films are squeezed to ∼0.7 nm (∼3 water molecule layers), the structures of the waterrock and wateroil interfaces change only marginally, except that the oil surface above the brine film becomes less diffuse. As the film is thinned from ∼1.0 to ∼0.7 nm, ions are enriched (depleted) near the rock (oil) surface, especially at a bath ion concentration of 0.1 M. These changes are caused primarily by the reduced dielectric screening of water and the weakened ion hydration near wateroil interfaces and, to a smaller extent, by the increased confinement. When the brine film is ∼1.0 nm thick, hydration and EDL forces contribute to the disjoining pressure between the charged rock and the oil. The EDL forces are reduced substantially as the ion concentration increases from 0.1 to 1.0 M, and the magnitude of the reduction is close to that predicted by the PoissonBoltzmann equation. When the brine film is thinned from ∼1.0 to ∼0.7 nm, the disjoining pressure increases by ∼10 MPa, which is mostly due to an increase in the hydration forces. The first layer of water on the rock surface is nearly stagnant, even in 0.74 nmthick brine films, whereas the viscosity of water beyond the first layer is bulklike, and the slip coefficient of oilwater interfaces is close to that under unconfined conditions. The insights that are obtained here help lay a foundation for the rational application of technologies such as lowsalinity waterflooding.

Parallel reservoir simulators for fully implicit complementarity formulation of multicomponent compressible flows(Computer Physics Communications, Elsevier B.V., 20190730) [Article]The numerical simulation of multicomponent compressible flow in porous media is an important research topic in reservoir modeling. Traditional semiimplicit methods for such problems are usually conditionally stable, suffer from large splitting errors, and may accompany with violations of the boundedness requirement of the numerical solution. In this study we reformulate the original multicomponent equations into a nonlinear complementarity problem and discretize it using a fully implicit finite element method. We solve the resultant nonsmooth nonlinear system of equations arising at each time step by a parallel, scalable, and nonlinearly preconditioned semismooth Newton algorithm, which is able to preserve the boundedness of the solution and meanwhile treats the possibly imbalanced nonlinearity of the system. Some numerical results are presented to demonstrate the robustness and efficiency of the proposed algorithm on the Tianhe2 supercomputer for both standard benchmarks as well as realistic problems in highly heterogeneous media.