Streamline-based simulation of nanoparticle transport in field-scale heterogeneous subsurface systems

Nanoparticle (NP) transport is increasingly relevant to subsurface engineering applications such as aquifer characterization, fracture electromagnetic imaging and environmental remediation. An efficient field-scale simulation framework is critical for predicting NP performance and designing subsurface applications. In this work, for the first time, a streamline-based model is presented to simulate NP transport in field-scale subsurface systems. It considers a series of behaviors exhibited by engineered nanoparticles (NPs), including time-triggered encapsulation, retention, formation damage effects and variable nanofluid viscosity. The key methods employed by the algorithm are streamline-based simulation (SLS) and an operator-splitting (OS) technique for modeling NP transport. SLS has proven to be efficient for solving transport in large and heterogeneous systems, where the pressure and velocity fields are firstly solved on underlying grids using finite-difference (FD) methods. After tracing streamlines, one-dimensional (1D) NP transport is solved independently along each streamline. The adoption of OS enhances flexibility for the entire solution procedure by allowing different numerical schemes to solve different governing equations efficiently and accurately. For the NP transport model, an explicit FD scheme is used to solve the advection term, an implicit FD scheme is used for the diffusion term and an adaptive numerical integration is used to solve the retention terms. The model is implemented in an in-house streamline-based code, which is verified against analytical solutions, a commercial FD reservoir simulaPreprint submitted to Advances in Water Resources September 13, 2020 Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 13 September 2020 doi:10.20944/preprints202009.0295.v1 © 2020 by the author(s). Distributed under a Creative Commons CC BY license. tor (ECLIPSE) and an academic FD colloid transport code (MNMs). For a 1D homogeneous case, the effluent breakthrough curves (BTC) produced by the in-house simulator are in good agreement with the analytical solution and MNMs, respectively. For a two-dimensional (2D) heterogeneous case, the BTC and concentration pattern of the in-house simulator all match well with the solution produced by commercial simulator. Simulations on a synthetic three-dimensional (3D) nanocapsule application engineering design case, are performed to investigate the effect of fluid and NP properties on the displacement pattern of an existing subsurface fluid.


Introduction
The transport and retention of nanoparticles (NPs) in porous media is relevant to a variety of applications. In groundwater, NPs can appear after release into the environment and subsequently cause facilitated transport of contaminants or may be environmental hazards themselves ( [1]). For filtration, separations, or microfluidic operations with NPs, the behavior of the NPs inside these porous domains dictates the effectiveness of the processes. In oil and gas and environmental cleanup operations, NPs are being developed and tested for use in site characterization, enhanced oil recovery (EOR), and conformance/mobility control ( [2] and [3]). For these applications, a wide range of functionalities are being designed, including singlecomponent particles that enhance contrast during subsurface imaging ( [4] and [5]), nano-capsules that release a payload at prescribed times or locations in a reservoir ( [6]), and smart NPs with chemical switches or memory to record reservoir conditions encountered between their injection and production/detection ( [7]). Additionally, in situations where larger particles can be used (e.g., in a propped fracture), micro-fabricated devices are being designed to collect and report data ( [8]).
A challenge for modeling NPs in any of these applications is that their flow and retention mechanisms make them behave differently than fluid or solutes that traditional models are designed for. If they are much smaller than pore sizes (a common situation), individual particles will not be susceptible to mechanical straining. However, they are highly sensitive to electrostatic interactions with surfaces, fluid interfaces, and other particles. These sur-face interactions, in turn, can result in reversible or irreversible capture ( [9]). Additionally, density differences can make NPs respond to gravity or inertial flows in ways that solutes do not ( [10]). Finally, engineered behaviors such as payload release or triggered responses add a more complex reaction component to the modeling, compared to traditional advective-dispersive-reactive transport ( [11,12]  Our understanding of basic NP behavior in porous media is founded on a long history of experiments and theory on colloids and interfaces, particularly examples such as ( [13]; [14]; [15]; [16]; [17] and [18]). Colloid Filtration Theory (CFT) provides a good point of reference for modeling the capture mechanisms ( Fig. 1) of small particles in porous media. Traditional CFT is derived from a single spherical collector (representing a particle in a granular porous medium) combined with modified Happel flow around a sphere ( [19]). CFT uses this flow field to predict the expected trajectory of finite-sized particles, which combined with surface forces leads to probabilistic values for capture of the particles. When mechanistic particle dynamics are upscaled and tied to the porosity of a material being modeled, the results can be used to estimate a continuum-scale capture efficiency representing the ability of the porous media to retain colloidal particles ( [20]). [21] gives a comprehensive review on transport and retention mechanisms and CFT models of NPs in porous media.
Based on mechanistic studies, various continuum models have been pro-posed to fit laboratory results by considering different transport mechanisms. Most continuum-based models were calibrated and validated against onedimensional column experiments by interpreting and fitting measured breakthrough curves (BTCs) from effluent NP concentrations. In classic CFT, a continuum model that only considers irreversible adsorption is used for colloids transport in porous media ( [13]). Because desorption of NPs is usually observed in experiments, continuum models that incorporate reversible adsorption are commonly used. For example, [22] shows that the experimental BTCs of engineered NPs (silver and carbon nanotubes) matches well with a continuum model that accounts for adsorption and desorption. [14] proposed a continuum model considering adsorption, desorption and site blocking, with adsorption limited by a maximum attachment capacity, which successfully fit fullerene NP breakthrough curves in a sand column. Recent experimental data indicate that NP capture in heterogeneous porous media is controlled by multiple mechanisms. For instance, [16] proposed a twosite continuum model considering both reversible and irreversible retention mechanisms. This model was successfully validated with five different sets of experimental data for silica NPs. Results also show that capture is greatly affected by fluid salinity. [18] developed a single-site reversible retention CFT model that accounts for the effect of different ionic strengths, so that the proposed model is able to reproduce the experimental observations as well as breakthrough curves under different values of fluid salinity. Continuum models with more complex mechanisms are summarized in the review paper by [9]. While various mechanistic and continuum models have been successful in reproducing and understanding NP transport behavior from 1D laboratory columns, studies of NP transport in realistic large-scale systems are scarce and need exists for new 3D modeling tools to take more of the NP transport mechanisms into account ( [9]). In addition, modeling NP transport is challenging from a computational standpoint in highly heterogeneous environments (both structural and chemical heterogeneity) such as those found in the subsurface, and with the complex reaction behavior that may be exhibited by smart NPs. Fast NP transport solutions are crucial for industrial application of engineered NPs with these complex situations.
In general terms, streamline-based simulation (SLS) is known to be an efficient method with low numerical dispersion. An important feature of modern SLS is that complex 3D flow and transport problems can be simplified to a series of 1D solutions along each streamline. This feature makes SLS particularly effective when dealing with large, highly heterogeneous systems particularly for advection dominant flow. Since the 1990s, SLS has been extended to simulate non-isothermal, multiphase, multi-compositional flow and mass transport in a wide range of reservoir scenarios ( [23]; [24]; [25]; [26]; [27] and [28]). It has been applied to many field-scale petroleum engineering applications including waterflooding, polymer flooding, carbon dioxide storage and the history matching of field production data ( [25]). Of particular note here, [24] successfully applied SLS to model reactive solute transport in porous media, where the solution of the SLS matches well with a particle tracking code and a commercial simulator. Results show that SLS has low numerical dispersion and high computational efficiency for advectiondominated transport.
In this paper, we present a new SLS model to simulate multi-component NP transport in heterogeneous subsurface systems. The model includes essential physics that affect nanoparticle transport and it incorporates a reaction matrix to accommodate the much broader range of behaviors one can expect from smart engineered nanoparticles compared to traditional reactive systems. An operator-splitting approach is used to develop a flexible and accurate numerical scheme for the NP transport along streamlines, including a new framework that implements a stable explicit FD scheme, a classic implicit FD scheme and an adaptive Runge-Kutta scheme to solve for advection, diffusion and nonlinear reactions, respectively. Additionally, the time-of-flight (TOF) coordinate transformation in SLS, is used to model time-triggered NP behaviors such as encapsulated release.
The model is verified with a series of analytical solutions under the limiting conditions for which these solutions are relevant. It is also compared for general accuracy and efficiency against a commercial FD reservoir simulator (ECLIPSE) and compared for NP transport against a specialized onedimensional finite difference code (MNMs, [29]). Finally, it is demonstrated on the more complex scenarios for which it was designed by simulating a series of hypothetical situations where the functionality of the reaction matrix is used in conjunction with a 3D heterogeneous reservoir model.
The remainder of the paper is organized as follows. In Section 2, the assumptions and governing equations of flow and transport for NPs are introduced, together with the streamline-based transport equations. In Section 3, the algorithm and numerical implementation are presented. In Section 4, several verification examples are presented and compared to analytical solutions and other simulators. In Section 5, a proposed field test for NP-based conformance control is simulated and analyzed for a reasonably large reservoir model.

Mathematical Formulation
In this section, the governing equations for NP transport and reaction are introduced. Transport mechanisms include advection and dispersion. Behaviors generalized into the reaction term include retention, time-triggered encapsulation, permeability loss, and concentration-dependent viscosity. Several assumptions are made in the current development, but could be relaxed if necessary: 1) Single-phase incompressible fluid flow, solved in a regular Cartesian grid system; 2) NP transport is dominated by advection and heterogeneity, meaning molecular diffusion, transverse dispersion and gravity are neglected. We ignore molecular diffusion becuase we expect to be modeling high-Pecletnumber conditions, and ignore gravity both for numerical simplicity and because current simulations do not exhibit significant density differences. Dispersion is more complicated; we have chosen to neglect transverse dispersion partly for numerical efficiency and partly because for mechanical dispersion, the longitudinal component is approximately an order of magnitude larger (than transverse) for the materials we are considering. Well-established operator splitting methods that could be used to relax these assumptions include ( [25], [30]).

Governing Equations for Flow and Transport
The mass balance equation for single-phase incompressible flow in porous media is: where u [L/T] is the Darcy velocity and q [1/T] is a volumetric source or sink per unit volume. For 3D flow (not considering gravitational forces), the differential form of Darcy's law is: Introducing Darcy's law into the mass balance gives the governing equation for fluid flow: where k [1/L 2 ] is the permeability tensor, µ [M/L/T] is the fluid viscosity (dependent on NP behavior or concentration in some situations), p [M/L/T 2 ] is the pressure, and [dimensionless] is a mobility correction factor. This mobility correction term is dependent on the concentration of deposited or attached particles in the rock matrix and is used to account for structure-based permeability changes attributed to the presences of NPs (e.g., mechanical plugging by the particles themselves, stimulation by acid released from NPs, etc.). The advection-dispersion-reaction equation (ADRE) is used to describe the transport of NPs in porous media ( [16] and [17]). The ADRE for the i-th NP type can be expressed as follows: where c i [M/L 3 ] is concentration in the fluid, s i [M/L 3 ] is the concentration of the adsorbed NPs on the rock surface, D [L 2 /T] is the hydrodynamic dispersion tensor, v [L/T] is the interstitial velocity (which can be related to Darcy velocity as v = u/φ, and R is the reaction term used to represent all manner of changes to the NPs or caused by the NPs. We refer to these changes collectively as reaction behaviors, and examples include surface interactions that transform NPs from free to adsorbed or vice versa (R ret ), burst of a capsule that eliminates NPs of one type but releases other molecular or particle species (R cap ), and conformational changes in response to time or environment that transforms them to a different species in the model (R r ). Retention is often the most important reaction behavior ([9]).

Hydrodynamic dispersion
Based on the assumptions listed above, only longitudinal dispersion is considered and can be directly calculated along streamlines. The simplified dispersion tensor is calculated in terms of longitudinal dispersivity α L [L] as follows ( [31]): Considering a Cartesian coordinate system, |v| [L/T] is the magnitude of the interstitial velocity, and v x , v y , and v z [L/T] are the respective components of the velocity vector.

NPs retention
One of the most important reaction behaviors in many scenarios is adsorption/desorption. Based on [16]'s (2016) research, a two-site retention model captures many experimental observations, including both irreversible and reversible attachment. The retention reaction term Rret(c,s) can be expressed by a general two-site retention model as in Eq. 6 ( [16]):

Other NP reaction behaviors
When envisioning types of smart nanoparticles that might be developed in the future, the possibilities for complex reaction behaviors are almost endless. The model is designed to accommodate new functionalities, with operator splitting (described below) ensuring that new reaction terms can be added in a modular way. However, to keep a reasonable focus in this paper, we describe one specific example: a nanocapsule containing cross-linker as its payload, set to burst with a time-based trigger. This reaction mechanism requires monitoring burst time for one or more batches of the capsules that would be injected, converting from capsule concentration to payload-species concentration at the time the payload is released, and for the payload species to react with other components in solution (e.g., a cross-linker reacting with polymer that was injected separately).
Nanocapsules ( [11,12]) are a promising new technology for controlled delivery in the subsurface. A species of interest (e.g., acid, crosslinker, surfactant, etc.) is encapsulated into a nanocapsule, which then dissolves/bursts at a desired time, upon encountering an external trigger, or from an internal trigger. The idea is to release its payload at an optimal location in the reservoir, to protect the completion or near-wellbore region from the species in the nanocapsule, and/or to minimize loss of the species before it reaches the intended target. For the sake of brevity, we focus on modeling a time-triggered release from a nanocapsule.
This particular reaction behavior assumes that an encapsulated NP type i will release an encapsulated payload chemical j when its time since injection (or time since mixing) reaches the designed burst time. Assuming the payload (c j ) release follows a simple first-order reaction of c i -¿ c j and s i -¿ c j , the reaction model is expressed as follows: For encapsulated NPs attached on porous media (s i ), similar reaction models can be written as: where w i,j [M/M] is the releasing mass fraction defined as the mass fraction between each encapsulated payload j and the total mass of a nanocapsule i ( [11]), k re [1/T] is the empirical release-rate constant, and t i,burst is the designed burst time of a nanocapsule. If the capsule is designed to release based on a different trigger than burst time (e.g., a local concentration or temperature), the reaction term would be modified accordingly and tied to the relevant variable in the model. After internal species are released from the nanocapsule, they continue to follow their own relevant reaction models.
A complicating factor is that nanocapsules of the same type, but having different clock times toward their triggered release (e.g., because they were created or injected at different times on the surface), can end up mixed together in the subsurface. Modeling the release behavior for this situation is nontrivial, and is addressed using a Lagrangian splitting approach discussed below.

Viscosity
Nanofluids are colloidal suspensions of NPs in a base fluid. The viscosity behavior of these nanofluids is often experimentally determined as a function of volumetric concentration ( [32]). Various empirical models for different types of NPs have been proposed and two classes of empirical models can be summarized as follows ( [33]): where, a 0 , a 1 , a 2 , a 3 , a 4 and a 5 [dimensionless] are empirical coefficients, µ brine is the viscosity of the base fluid (brine), and . For a multicomponent nanofluid with nNP nanparticles, the volumetric concentration can be calculated as follows ( [34]): where c i and ρ i is the NP concentration and density of i-th component in a nanofluid mixture.

Permeability Alteration
NP injection has the potential to alter permeability, for instance damage from mechanical straining of the NPs or stimulation from injection of acidencapsulated nanocapsules. In the model, effects of permeability change are kept separate from viscosity in Eq. 3 (i.e., instead of lumping them as a single mobility term), which allows the NPs to cause permanent permeability change at a specific location, independently of viscosity changes that can continue to move with the flow.
Experimental studies on permeability loss due to straining have been performed, including [35] and [36]. Permeability reduction is modeled using a linear flow resistance factor expressed as follows ( [35]): where β adsorption is the formation damage coefficient, dimensionless. A situation that causes permeability change based on a different mechanism (e.g. acid dissolution of the porous medium) would need a different constitutive equation for k 0 /k, likely fitted empirically or tied to a simple permeabilityporosity relationship if an appropriate version is available.

Streamline Formulation
Streamline simulation includes the important feature that a 3D transport problem (Eq. 4) can be reformulated to a series of 1D sub-problems along each streamline in terms of the time-of-flight (TOF) ( [25]). Following the definitions and notations from [25], the TOF (τ ) is the time it takes for a neutral particle for travel to a specific distance L along a streamline. Mathematically, TOF can be expressed as: Based on Eq. 12, the following transform operator can be defined as Following the derivations for multiphase flow described in [23], the 3D transport equation of Eq. 4 can be transformed into multiple 1D equations using the above transform operator along streamlines and can be rewritten as: Detailed derivations are presented in Appendix A. The first, second and third terms on the right-hand-side of Eq. 14 describe the advection, longitudinal dispersion, and generalized NP reaction terms, respectively.

Numerical Implementation
Streamline-based simulation has become increasingly popular due to its advantage in computational efficiency for large, geologically complex and heterogeneous systems. Operator Splitting (OS) is a key aspect of streamlinebased simulation. As shown in Fig. 2, OS is used to decouple the flow (Eq. 3) and transport (Eq. 14) equations into an underlying 3D spatial grid and a 1D streamline "time-of-flight" grid ( [25]). In this scheme, the flow solution on the underling grid allows one to have different timesteps along a streamline. Additionally, the off-streamline physics vs streamline-oriented physics can be solved independently on the 3D grids and streamline grids, respectively. OS is also used to handle reaction term, which we view as essential in order to handle the nearly infinite type of reaction behaviors that smart nanoparticles may exhibit in the future. Instead of solving ADRE in a coupled manner, as shown in Fig. 3, the solution of advection (Eq. 16), dispersion (Eq. 18) and various reaction terms for NP transport (such as, retention in Eq. 6 and encapsulation in Eq. 7) can be solved independently with different numerical methods.
The use of operator splitting in the NP streamline-based simulations is summarized by the following steps: 1) Solve the pressure equation and compute face velocities on the underlying 3D grid. This step is repeated periodically depending on changes to fluid and rock properties, denoted by global pressure update timesteps ∆T .
2) Trace streamlines and construct 1D TOF grids for each streamline.
3) Solve the advection, dispersion, and various reactions independently along each streamline using local streamline timesteps ∆t. 4) Map transport solution variables (c i and s i ) from 1D TOF grids to the underlying 3D grids. 5) Update the concentration-dependent fluid and porous-media properties, such as concentration-dependent nanofluid viscosity (Eq. 9) and permeability loss (Eq. 11), in the underlying 3D grids If necessary, repeat Steps 1 & 2 to re-compute velocity field and obtain new streamlines. 6) Map transport solution variables (c i and s i ) back from the underlying 3D grids to 1D TOF grids. 7) Repeat Step 3 -6 until the end of simulation time Operator splitting allows use of a relatively large pressure-update timestep (global timesteps ∆T ) for the pressure solution in the 3D underlying grids, with a smaller transport timestep (local streamline timesteps ∆t) for the transport solution in 1D streamline grids. The flowchart for the numerical implementation of the streamline-based simulation is shown in Fig. 2 and the computational procedures for each step are described in detail in the following sub-sections. As mentioned above, the off-streamline physics, such as molecular diffusion, transverse diffusion and gravity, can be included by solving the diffusivity equation for NP concentration in the underlying 3D grids after Step 5 [30]. Step 2) Streamline Step 1) Underlying Grid Step 3-4) Streamline<-->Underlying Grid Step 2) Streamline Step 1) Underlying Grid Step 3-4) Streamline<-->Underlying Grid

Pressure solution
As shown in Fig. 2, the first step in streamline-based simulation is to calculate the pressure field and face velocities for each block on the underlying grid. The steady-state single-phase flow equation (Eq. 3) is solved based on the standard finite difference method with a 7-point stencil ( [37]). A no-flow boundary condition is prescribed on the reservoir boundaries of the entire model. The boundary conditions for wells are either constrained by constant flow rate or constant bottom-hole pressure (BHP). The classic Peaceman well model ( [38]) is employed in the modeling. Details of the implementation of the pressure solution can be found in [37]. Once the pressures in each gridblock are obtained, the Darcy velocityu on each of the six faces of a gridblock can be calculated using Eq. 2 between two neighboring grids. In addition, the permeability and viscosity field are frequently updated (e.g. every ∆T ), in effect linearizing the nonlinear pressure equation, for effects related to formation damage or viscosity change from the NP behavior or other factors such as fluid mixing. The frequency of global pressure-updates ∆T is influenced by non-linearity of simulation model. In this paper, the ∆T determined by a simple trial-and-error method, where several simulation runs are performed with decreasing ∆T until the solution variables converge. This process is demonstrated in the numerical example Case 2. More details of advanced time-step control scheme can be found in [30]'s work.

Streamline tracing
Knowing the interstitial velocity field on the underlying Cartesian grid, [39]'s semi-analytical method is used to trace a streamline grid-by-grid from injector to producer. The velocity components in the x, y and z-directions are assumed to vary linearly in each gridblock. Therefore, the exit coordinates and TOF can be analytically calculated within a particular gridblock based on given inlet coordinates.
Launching streamlines from gridblocks containing wells cannot guarantee every underlying grid cell will contain as least one streamline, especially for large heterogeneous models. The gridblocks containing no streamlines are marked as missed gridblocks, for which fluid properties must still be defined ( [23]). For example, mapping NP concentrations from streamlines back to underlying grids requires that each grid cell has at least one streamline passing through it. After tracing streamlines from injector(s) to producer(s), new streamlines can be launched from the center of each missed gridblock and traced backward to an injector and forward to a producer. This process is repeated until there are no missed grids ( [23]). One other note is that Pollock's algorithm cannot be directly applied to a grid cell containing an injector or a producer due to the existence of a sink or source ( [39]). The fill-grid method (Eq. 15) can be adopted to estimate the average travel time for grid cells containing wells ( [40]).

NP transport solution
From the constructed 1D streamlines, the NP transport equations can be solved along each streamline by discretizing a streamline into 1D grids in terms of TOF. In order to simulate the potentially complex behaviors exhibited by nanoparticles (e.g., including nanosensors and nanocapsules), a flexible and robust numerical scheme is needed. One example discussed below is core-shell encapsulation conformance control technology (see Section 5), which involves four species and five reaction behaviors. Operator splitting helps make this and similar situations tractable, as shown in other applications with complex reactive transport behaviors with an arbitrary number of species ( [41]). For the current model, a numerical scheme based on OS is developed to solve NP transport along each streamline as shown in Fig. 3: In our approach, explicit finite-difference, implicit finite-difference and adaptive Runge-Kutta (RK) methods are implemented sequentially to solve the advection, dispersion and reaction terms, respectively. This scheme provides a stable solution for both advection and dispersion by using the TOF grids along each streamline, which results in an unconditionally stable NP transport solution by setting the uniform TOF grid size ∆τ equal to the local streamline timestep ∆t ( [42]). The detailed implementation and derivation are described below.

Advection
For each NP species i, the advection term in Eq. 14 can be solved using explicit finite difference method as follows: The local streamline time-steps ∆t should be selected based on Courant-Fredrich-Lewy (CFL) conditions to guarantee numerical stability, ∆t/∆τ ≤ 1. If a maximum CFL of 1 is adopted, the spacing of 1D regular streamline grids, ∆τ is equal to the local streamline solution time step, ∆t. Substituting CFL condition of ∆t/∆τ = 1 into Eq. 16 gives: Eq. 17 gives the exact and unconditionally stable solution ( [42]) for advection, which implies that advection of the concentration from the previous node will move forward to the next node in one timestep. This solution removes the timestep restriction for stability in the conventional treatment for advection in streamline simulation by equating spatial and temporal step size ( [24] and [43]).

Longitudinal dispersion
For each NP species i, the longitudinal dispersion term in Eq. 14 can be solved using an implicit finite difference method as follows: The discretization scheme for advection is also adopted for dispersion. Substituting ∆τ = ∆t into Eq. 18 gives: whereD L,k+1/2 andD L,k−1/2 are the dispersion coefficients at each interface of two gridblocks, and which can be directly estimated by interpolating the local velocity at the interfaces: At any position on the 1D regular streamline grids (Fig. 4), the local interstitial velocity |v| can be represented using linear interpolation based on Pollock's algorithm as follows: where A x , A y and A z [1/T] are the components of the linear velocity slope within an underlying grid. For example, the velocity gradient in the x direc-

NP retention
NP retention is solved on each grid at each local timestep. For NP species i, the retention part of the reaction term in Eq. 14 is a system of nonlinear ordinary differential equations (such as, Eq. 6), which can be calculated using a numerical integration method. In this paper, the adaptive Cash-Karp RK algorithm ( [44]) is implemented, which will adaptively subdivide the local timestep ∆t into several sub-solution steps to control the relative error within a given tolerance (which in this work is set to 1x10 −6 ).

Time-triggered encapsulation behavior
If nanocapsules are injected in a way other than a single pulse, then we may expect to encounter gridblocks with a mixture of time-triggered nanocapsules having different burst times. This problem is not restricted to nanocapsules, but is relevant for modeling any time-triggered behavior in which the start time for the trigger is not uniform. Reasons might include particles with a chemical timer that are created in sequential batches, or a technology where injection is continuous, but each particle's internal clock begins when it enters the wellbore. (Particles with a range of different burst times is a related issue, which is handled by treating these as separate NP species in the model.) Numerically solving for the burst-and-release reaction behavior for such mixtures is nontrivial. We propose a Lagrangian splitting approach to solve the reaction behavior described in Eqs. 7-8, where injected nanocapsules are split into different batches based on their injection time (Fig. .19). For each batch, the burst time is tracked, and the reactive transport equation is solved independently. The brute-force implementation of the Lagrangian splitting approach is to split the nanocapsules into several n batches where n=injection time/timestep. Obviously, the computational cost of the bruteforce implementation will grow significantly for a simulation with a long period of injection. For example, a plan to inject nanocapsules for 100 days with a timestep of 1 day will require solution of 100 ADR equations at each simulation timestep. Instead, we developed an optimized implementation where the ADR equations are grouped into several grouped batches (n batch ). The size of the batches (n batch ) is determined from the ratio between the burst time of the nanocapsule and simulation timestep. The details of the optimized implementation for Lagrangian splitting approach can be found in the Appendix. B.

Implementation
The model has been implemented in an in-house three-dimensional streamline code, which referred to internally and in this paper as NPSL (NanoParticle StreamLine). The current user interface is either through text input files or a Microsoft Excel front end that allows more intuitive control but then launches the same code. The reaction matrix is an essential part of the functionality. It allows a user to define complex NP scenarios by assigning relevant reaction behaviors to each species in the simulation design by activating them on a grid. For instance, if species A is a nanocapsule, the user might check "time-triggered burst" and "reversible adsorption" as two of the behaviors. If species B is acid contained inside a nanocapsule, the user might check "comes from a capsule" and "permeability change" as two relevant behaviors. Each reaction behavior activated by the user requires associated parameters (e.g., burst time, adsorption/desorption rate parameters, reaction stoichiometry, and dissolution kinetics for the behaviors listed above).

Model Verification
In this section we present case studies on 1D and 2D scenarios that have been performed to verify specific aspects of the model against exact solutions, the academic colloid transport code MNMs ( [29]), and the commercial reservoir simulator ECLIPSE ( [45]).
The NP transport solution is first verified against an analytic solution and MNMs simulation for a 1D homogeneous model. In the second case, the functionalities for solute transport and viscosity change are compared to a commercial simulator on a 2D heterogeneous model. Finally, an error analysis is performed to quantify the numerical errors of the presented solution algorithm. The base model and simulation parameters for all cases are shown in Table 1

Case 1: Verification for advection, dispersion and retention
To evaluate our implementation of the streamline model for solving the advection-dispersion equation, the NP transport solution is verified against an analytical solution and an established numerical solution (MNMs). As shown in Fig. 5, a 1D reservoir with a length of 298.70 m has been discretized into 50 gridblocks in the x direction. An injection well is located at the center of the first grid block (1,1,1) and given a constant injection rate of 79.49 m3/day. A production well is located at the center of the last grid block (50,1,1) and given a constant bottom hole pressure of 27.58 MPa. A NP slug is injected at a concentration of 1 kg/m3 for a period of 10 days.
The detailed model and simulation parameters are given in Table 2. Other parameters are the same as shown in Table 1 above.

Injector
Producer Grid: 50x1x1   For the analytic comparison, we set the irreversible adsorption coefficient to zero and make the surface adsorption capacity infinite, which allows direct comparison with the model described by [46]: The analytical solution of the above problem with pulse injection can be derived based on Laplace transformation and superposition as follows ( [47]):  Fig. 6 shows the effluent breakthrough curve (BTC) with different Peclet numbers and Damköhler numbers (Fig. 7), comparing the NPSL and analytical solutions of Eq. 23. The numerical solutions from NPSL match well with analytical solutions for a wide range of Da and Pe. In particular, we note that the sharp front is accurately captured without any oscillations, that high-Da behavior is indicated by the lower NP concentrations in the production well, and that low-Pe behavior is indicated by the dispersive spread of the effluent concentration and an earlier breakthrough time.
Solution of the nonlinear NP transport equations (Eq. 6) is also validated with an academic colloid transport code, MNMs, which has successfully reproduced nanoparticle column experiments ( [29,48]). This second comparison allows modeling of an adsorptive capacity of the porous medium quantified by s 2,max (see Eq. 6). In this case, parameters remain the same as before except the retention capacity s 2,max is set equal to 0.1 mg/g and reversible attachment rate coefficient are varied as k ra = 0, 0.078, 0.39 day 1. As shown in Fig. 7, NPSL shows excellent agreement with MNMs for various Da and Pe numbers. Note that to obtain these matching results, the MNMs system required 500 grid cells in the x-direction and smaller timesteps (∆t = 0.01 days) compared to the streamline approach.

Case 2: Variable fluid properties
The NPSL simulator was also compared against the commercial reservoir simulator ECLIPSE in a heterogeneous 2D model, with a focus on simulating variable fluid properties. The scenario used in this comparison is a simple polymer injection scenario, without formation damage, retention, or reaction. By injecting polymer, the fluid viscosity will be linearly increased with increasing of polymer concentration. The implemented fluid viscosity model can be expressed as follows ( [45]) where µ 0 is the initial polymer viscosity at the injection concentration c 0 . Eq. 25 describes a simple linear relationship between concentration and viscosity. The minimum and maximum values are the brine viscosity and initial polymer viscosity, respectively. Fig. 8 shows the permeability field from a portion of the 8th layer of the SPE 10 Comparative Solution Project -dataset 2 ( [49]). An injector is located at the center of the model and four producers are located at the four corners. A polymer slug is injected at a concentration of 1 kg/m 3 for a period of 400 days. Additional model parameters are shown in Table 3. All other parameters are the same as in Table 1. The concentration map and BTC for NPSL and ECLIPSE are shown in Figs. 9-10, respectively. Qualitatively, the NPSL solution matches well with ECLIPSE, but the streamline-based method shows less numerical dispersion  and presents a much sharper concentration front than the conventional fully implicit FD simulator. As expected, the BTC for ECLIPSE shows earlier breakthrough time than NPSL because of the higher numerical dispersion. These differences between finite-difference and streamline solutions are well documented ( [25]), and are pointed out here simply to explain the differences between the NPSL solution and a well-accepted numerical simulation tool. The nonlinearity introduced by variable fluid properties must be addressed by periodically re-solving the pressure field and regenerating the streamlines (at each global pressure update timestep ∆T ). These procedures require the polymer concentration to be mapped from the streamlines to the underlying grid, which allows the in-situ fluid viscosity to be updated by Eq. 25. As mentioned earlier, the update frequency (∆T ) is determined by some form of a trial-and-error method. As shown in Fig. 11, the bottom-holepressure (BHP) of the injector converges when ∆T is reduced from 200 days to 40 days. Further decreasing ∆T to 20 days introduces more numerical mapping errors ( [26]) and the BHP shows abnormal fluctuations. Thus, ∆T is set at 40 days with 50 pressure updates for all SPE10 simulations in this paper. (The mapping error issue of SLS is discussed by [26] and can be addressed by advanced timestep control technologies if needed ( [30]).)

Model Applications
In this section, three examples of injection scenarios of smart NPs are presented to demonstrate the versatility and applicability of the presented streamline NP transport model.

Case 3: 1D laboratory column experiment
In the first example application, NPSL is used to perform inverse modeling of a silica NP column experiment ( [50]), in which a 3 PV slug of 5 wt% silica NPs was injected into a one-foot column of 90 wt% 250-297 µm ,QMHFWRU%+3SVL 3UHVXSGDWHV 3UHVXSGDWHV 3UHVXSGDWHV 3UHVXSGDWHV Figure 11: Convergence of BHP for NPSL with different pressure-update frequencies equal to 10, 25, 50 and 100 pressure updates, corresponding to ∆T =200, 80, 40 and 20 days, respectively.
Boise sandstone mixed with 10 wt% kaolinite. Using the NPSL code, the dispersivity of the column sample was fit against the measured BTC of a tracer. The retention parameters are then fitted against the BTC of the NP injection ( Table 4). As shown in Fig. 12, the BTC of the NPs has a delayed breakthrough curve and a tail, compared to the tracer. The NPSL solver matches well with the NP experiment (R 2 ¿ 0.97).

Properties
Units Value  1D

Case 4: Encapsulated NP tracer injection into a 2D reservoir
Recently, some advancing technology that utilize encapsulated tracers for subsurface aquifer or reservoir characterization ( [51]). This example illustrates how the model could be used to help optimize treatment design. It simulates NP encapsulation technology that can deliver a payload or tracer deeper into a reservoir by protecting it from the reservoir environment (which may cause loss from adsorption or premature reaction). A quarter five-spot problem is used with a 2D heterogeneous reservoir from the 85th layer of the SPE10B model ( [49]). Fig. 13 shows the permeability field and flow streamlines. The retention behavior parameters of the NP are taken from Table 4. All other simulation parameters are the same as Case 2. As shown in Table 5, three simulation scenarios are tested. In case 4a, 1 PV of tracer is injected in the reservoir as the reference solution; In case 4b, 1 PV of NP tracer is injected but is susceptible to retention behavior. In case 4c, the same amount of NPs are injected but they are nanocapsules with tracer inside. The mass fraction and burst time of the nanocapsules are 100% and 30 days, respectively.

Physics
Equation 4a 4b 4c Advection, Dispersion (4) X X X NP retention (6) -X X Encapsulation (7-8) --X Fig. 14 shows the concentration map of the NP tracer and the encapsulated NP tracer. By comparing the results between Case 4b and Case 4c, the encapsulated NPs (Case 4c) are released far from the injection well, creating less absorptive loss and the earlier breakthrough time. Fig. 15 is a plot of the produced NP concentration, which shows that retention behavior in Case 4b and Case 4c leads to delayed breakthrough and lower peak concentration (due to NPs being captured via adsorbing to the formation rock). However, a larger amount of the NPs are produced in Case 4c, where they are protected by encapsulation during the first 30 days.
Generally, time-trigger encapsulation technology can be designed to protect the payload flowing through the reservoir from being affected by the reservoir environment in the vicinity of injectors. The encapsulation technology may release the payload in a variety of ways (temperature, encountering a chemical, PH value, etc.).

Case 5: Encapsulation conformance control in a 3D reservoir
The purpose of this case is to demonstrate the applicability of NPSL in 3D field-scale application model with complex reaction behavoir. A pilot simulation of encapsulation conformance control ( [52] and [11,12]) is simulated on a quarter five-spot problem based on a portion of the SPE10B model (Fig. 16). An injector is located at one corner of a rectangular model and a producer is located at the opposite corner. The wells are completed in all layers of the model. The detailed model parameters are shown in Table 6.
As shown in Table 7 three simulation scenarios are shown: waterflooding and polymer flooding are simulated as reference solutions in Case 5a and Case 5b, respectively. In Case 5b, 0.67 PV of 10 cp polymer is injected in the reservoir. In Case 5c, 0.67 PV of time-triggered nanocapsules containing a cross-linker and a special chemical base fluid is injected into the reservoir. The burst time of the nanocapsules was set as 10 days. After the nanocapsules reach their burst time, the cross-linker is released and can react with   Figure 16: Permeability map (left) and streamline with TOF map (right) of case 5 the base fluid to produce a 10 cp polymer. In order to compare the results with those in Case 5b, a cross-linking reaction (cross-linker + base fluid = polymer) is used to produce a polymer with the same properties used in Case 5b. This system also assumes the reaction reaches to steady state within 1 simulation timestep (5 days). Similar to retention behavior, the cross-linking reaction can be solved using numerical integration where the reaction model of can be expressed as:  (26) --X Note that shear-thinning effects and adsorption are neglected for all polymers simulated. To evaluate the effect of different schemes on the conformance control, a pre-existing tracer is uniformly distributed in the domain. This allows the different displacement schemes (injected water, injected polymer, or in-situ-created polymer) to be compared through the production response of the tracer that was originally present. (The intent is to mimic the sweep of oil by injected fluid, even though the model is currently being run for single-phase conditions.) Retention and formation damage parameters vary with the different rock and NP encapsulation types and are usually determined by column experiments. In this case, experimental NP retention parameters are adopted from a large-scale lab experiment ( [10]). Formation damage is also considered by assuming the maximum permeability reduction is 50%, which means the formation damage coefficient, β adsorption , can be calculated based on Eq. 11. The viscosity of the nanofluids can be determined based on the concentration and density of the NPs. Because silica, iron and iron oxide NPs are often used in subsurface research ( [16]), an average nanoparticle density of 3,000 kg/m 3 is used to estimate the volume concentration and the viscosity of the nanofluid. The exponential function shown in Eq. 9 is used to compute fluid viscosity, with empirical coefficients of a 4 =12.5 and a 5 =6.356 for iron oxide nanoparticles ( [33]). The detailed NPs encapsulation and fluid properties are provided in Table 8.

Properties
Units   Fig. 17 shows the distribution of the unswept tracer and the recovery ratio for the three different scenarios, at the end of each simulation. A lower viscosity contrast (µ displacing /µ displaced = 1) between the displacing and the displaced fluid in Case5a will lead to a lower recovery ratio ( [37]). As shown in Fig. 17, in Case 5a, most of the displacing fluids are passing through the high permeability channels (especially in the upper layers of the model) leading to a large quantity of displaced tracer remaining in the lower permeability regions. The simulation cases with higher viscosity contrasts (Case 5b & 5c) show better sweep efficiency. The use of encapsulated crosslinker to create the polymer in situ (Case 5c) has a slightly increased the recovery ratio by sweeping the marginal regions towards the top left corner. In addition, as shown in Fig. 18, the injection pressure of Case 5c is significantly reduced compared to the conventional polymer treatment (Case 5b). While these results represent only one set of conditions, they illustrate how the NPSL could be used to optimize sweep efficiency and injection pressure by varying treatment design, including factors such as volume of the initial pad, total volume of capsules, burst time, concentrations of base polymer and crosslinker, or even sequencing the pad and crosslinker to enhance in-situ mixing.

Conclusions
The paper extends the efficient streamline-based method to simulate nanoparticle transport and complex reaction behavior in field-scale models. The model considers a number of important characteristics of NP transport, including the retention and fluid properties of nanofluids, encapsulation and formation damage. An improved numerical scheme for NP streamline-based simulation is presented and several conclusions can be drawn as follows.
1) The proposed method is a robust and accurate method for nanoparticle transport especially for large-scale problems. Based on operator splitting and streamline-based methods, the improved numerical scheme guarantees the transport solution is oscillation-free, bounded, and stable even when large timesteps and the conditions of high Peclet number (Pe =1000) are used. Taking a 1D case as an example, in order to match the analytical solutions, the conventional fully implicit finite different scheme requires 10 times finer grid and 50 times smaller timesteps than the streamline-based method.
2) A series of comprehensive verification cases were also conducted by comparing the proposed method to analytical solutions and solutions from conventional FD-based academic codes and commercial simulators. The solution for advection, diffusion, retention and variable fluid properties all match well with the reference solutions and experimental data.
3) A pilot engineering design of encapsulation conformance control is presented. Results indicate that the NPs encapsulation polymer flooding is a promising technology in improving the performance of conventional polymer flooding, not only for improved sweep but lower injection pressure.
4) The method's computational speed suggests it is a promising tool for optimizing NP treatment design for subsurface applications.
Applying the transform operation from Eq. 13 into the advection term of Eq. A-2, the advection in terms of TOF coordinates along a streamline can be expressed as follows: Similarly, applying Eq. 13 again into the dispersion term (only streamlineoriented longitudinal dispersion considered) of Eq. A-2 gives: Thus, the 1D NP transport equation can be expressed by substituting Eqs. A-3 and A-4 into Eq. A-2 as follows: Appendix B -Algorithm of Lagrangian splitting approach As shown in Fig. .19, assume a nanocapsule has a designed burst time of t burst , the injection is started at time t 0 , and continues for a period t inj . If the solution timestep is ∆t, the nanocapsule can be split into n batch = t burst /∆t batches, injected in ceil (t burst /∆t) rounds of batches. During each injection round m = 1, 2, . . . , ceil (t burst /∆t), the timestep nt inj,l to inject a batch l = 1, 2, . . . , n batch can be calculated by nt inj,l = l + n batch (m − 1) + t 0 /∆t − 1. The burst timestep nt burst,i for the batch can be determined as nt burst,i = nt inj,l + t burst /∆t . The NP transport equations for advection (Eq. 21), dispersion (Eq. 22), and reaction terms are solved for each batch individually. Only when the global timestep reaches the burst timestep for a batch will the payload releasing equations (Eqs. 9-10) be solved using an adaptive Cash-Karp RK algorithm.