Process-based upscaling of reactive flow in geological formations

Abstract Recently, there is an increased interest in reactive flow in porous media, in groundwater, agricultural and fuel recovery applications. Reactive flow modeling involves vastly different reaction rates, i.e., differing by many orders of magnitude. Solving the ensuing model equations can be computationally intensive. Categorizing reactions according to their speeds makes it possible to greatly simplify the relevant model equations. Indeed some reactions proceed so slow that they can be disregarded. Other reactions occur so fast that they are well described by thermodynamic equilibrium in the time and spatial region of interest. At intermediate rates kinetics needs to be taken into account. In this paper, we categorize selected reactions as slow, fast or intermediate. We model 2D radially symmetric reactive flow with a reaction-convection-diffusion equation. We show that we can subdivide the PeDaII phasespace in three regions. Region I (slow reaction); reaction can be ignored, region II (intermediate reaction); initially kinetics need to be taken into account, region III (fast reaction); all reaction takes places in a very narrow region around the injection point. We investigate these aspects for a few specific examples. We compute the location in phase space of a few selected minerals depending on salinity and temperature. We note that the conditions, e.g., salinity and temperature may be essential for assigning the reaction to the correct region in phase space. The methodology described can be applied to any mineral precipitation/decomposition problem and consequently greatly simplifies reactive flow modeling in porous media.


Introduction
The efficiency of improved or enhanced oil recovery (IOR/EOR) processes is often influenced by the composition of the flowing aqueous phase and is thus affected by the mass exchange between fluid and solid phases [1] . For example, the rheology of watersoluble polymers or the magnitude of the interfacial-tension reduction by surfactants strongly depends on the ionic strength and hardness (concentration of divalent cations) of the aqueous phase and strongly influence the displacement efficiency [2] . In recent years, the additional oil extracted by tuning the composition of the injected water has led to more detailed investigation of the nature of the interactions between the rock and the fluids residing in the pore [3] .
Additionally, dissolution and precipitation of minerals have a considerable effect on the permeability of the reservoir [4] . The permeability impairment can lead to high circulation costs and loss of injectivity [5] . Permeability impairment occurs both due to par-ticle release and subsequent retention [5] and precipitation of minerals [4] .
When a fluid is injected into an oil reservoir, the chemical species in the solution will interact with the species in the reservoir fluids contained in the pores and/or the substances on the rock grains. As these species are subject to transport in the reservoir, in the numerical modelling of the process it is essential to know the final composition of the solution at the end of each timestep [6][7][8] . This is determined by how far the reactions are completed within each grid cell during each timestep [7,8] . For fast reactions, the front develops over a short distance, while for slow reactions fronts develop over a longer distance from the injection point.
Subsurface formations are composed of several minerals, which adds another complexity in modelling of geochemical reactions [9,10] . The time required to reach equilibrium can vary from seconds to (a few thousand) years for different minerals [11,12] . Therefore, in the numerical simulations, during a single time step, the concentration of some species could be over-or underestimated. There are two important methods in geo-chemical modelling, viz., the Local Equilibrium Approach (LEA) and the Kinetic Approach (KA) [13] . The LEA, which is the most common ap-proach in reactive transport modelling, assumes that the reactions take place instantaneously and neglects the duration of the reaction. This approach is preferred because of its simplicity and the fact that it is computationally less expensive. However, the local equilibrium approach may over-or underestimate the concentration of the species in the solution, e.g., due to heterogeneity [14][15][16] .
Indeed, the reactions in nature do not always reach equilibrium within the time and space of interest. The Kinetic Approach describes the time-dependency of the reactions. Some reactions complete faster than others and therefore the assumption of local equilibrium does not hold in some practical cases. For example, in IOR/EOR applications, the injected fluid (which might be in equilibrium itself) is usually out of equilibrium with the rock (this is referred to as partial equilibrium) [1] . The reservoir rocks consist of several components whose reactions with the injected fluid have different rates. Therefore, the kinetics of the reactions in this case determines the composition of the fluid (and consequently the efficiency of the injected chemical). This poses challenges when the results of the lab-scale experiments are translated to large-scale field applications, because different time and length scales are involved. For this purpose, when the experiments are modelled, special attention should be given to correctly calculate the equilibrium time ( T eq ) and the equilibrium distance ( R eq ) of the relevant minerals. These quantities are used as guidelines for selection of the temporal and spatial resolution of the numerical simulations.
It is difficult (if not impossible) to develop a general upscaling technique to account for the geochemical reactions in different IOR/EOR processes. Apart from the geological uncertainties and the intrinsic heterogeneity of the formations, the information on the mineralogical composition of the rocks is scarce. Even if the mineralogy of the formation is partially known, the surface area available for the reactions is hard to measure. Notably, the bulk mineralogical composition determined from X-Ray Diffraction (XRD) data can significantly differ from the surface composition obtained by X-Ray Photoelectron Spectroscopic (XPS) [17] .
Therefore, it is suggested to follow a simpler approach, which is process specific. In such an approach, only relevant minerals affecting the efficiency of the process are listed. The main minerals considered are Na (K, Ca) -silicate (carbonate) minerals (see [9][10][11][12] for weathering rates). The particular interest is then to evaluate the time and distance that is required for the injected fluid to reach the equilibrium state. The equilibrium time ( T eq ) and the equilibrium distance ( R eq ) are indications of the temporal and spatial resolutions of the simulations.
It is therefore our objective to include mineral reactions in the modeling of reactive flow in porous media. This problem can be studied on two characteristic length scales R char : the well radius and the reservoir size. This choice defines also a characteristic time T char , see Section 2.3 . This gives the opportunity to study the equilibrium length ( R eq ) and equilibrium time ( T eq ) of these reactions with respect to the characteristic length and time. We distinguish three possibilities (I) R eq R char , the whole reservoir will be at the injection concentration (II) R eq R char , the reservoir will be at the initial concentration and (III) R eq ~R char , the concentration is varying throughout the whole reservoir. Note that in cases (I) and (II) the reservoir can be considered at equilibrium and kinetics do not have to be taken into account.
The paper is organized as follows. Section 2 gives the physical model and model equations. In Section 2.6 the model equations are reduced to a second order ODE with boundary conditions. Section 3 shows and compares analytical results for different values of the physical parameters. Section 4 computes R eq , and we determine whether we are in case (I, II, III) for a given set of minerals. We end the paper with some conclusions in Section 5 .

Physical and mathematical model
In this section we will present single phase transport and reaction (dissolution and precipitation of the omni-present minerals), using a reaction convection diffusion model. In the model equations we will use capital letters for dimensional variables, variables with subscript c for characteristic scales and small letters for dimensionless variables.

Reaction rates
The methodology presented in this work can be applied to general dissolution/precipitation reactions. As illustration we consider the following explicit reactions: K-feldspar: For the reaction rates S 1 we apply the often adopted simplifying assumption (see, however, [18] ) of using the law of mass action and obtain (see [10] ) where IAP is the ionic activity product, K is the equilibrium constant, is the reactive surface area per cubic meter water; κ [ mol/ (m 2 s ) ] is a temperature dependent rate constant. In practice the value of n is chosen to mimic experimentally determined reaction rates, and here chosen to be equal to one. The rate constant can be considered in its Arrhenius form where k 10 is the frequency factor, E a is the activation energy, R is the gas constant and T is temperature in K .
As concentration of interest, C , we use the dissolved ions, (e.g., Al 3+ ) or dissolved molecules (e.g. H 4 SiO 4 ). Moreover we assume that the reactions occur independent of each other and that we can single out one specific ion or molecule. This approach cannot describe the decomposition of the mineral, but does allow to classify reactions into slow, intermediate and fast reactions.
We assume furthermore isothermal flow at constant pH , which means that we assume that the ion activity product is related to the dissolved mineral concentration C , "independent" of other concentrations, which are present in excess so that they can be assumed to be constant. Various temperatures and pH values can be considered. Given its approximate nature it is not warranted to assume non-linear dependences, and we assume that IAP is proportional to C , i.e., We note that similar simplifying assumptions are made in the literature (see [7] ). We then lump the proportionality constant and K together into the constant C 0 and obtain

Transport equations
We assume incompressible flow, which means that we have the following equation for the interstitial velocity V Transport is described with a reaction convection diffusion equation where the reaction rate S is given in Section 2.1 . The concentration initially present in the reservoir is denoted by C 0 [mol / (kgwater)]; we assume that the reservoir initially is in equilibrium, meaning that we disregard slow naturally occurring variations that take place after times longer than our timeframe of interest. This means that, from a practical point of view, the reservoir is initially at the equilibrium concentration (at the initially present environmental conditions). This means that reaction will take place if we perturb the initial concentration C 0 to C . We are looking for radially symmetric solutions V = V (r, t ) ˆ r , C = C(r, t ) , which means that equation (6) reduces to and equations (7) -(5) yield

∂C ∂T
which reduces to
Equations (8) and (10) are made dimensionless using the characteristic length R c , characteristic time T c (discussed below) and the initial equilibrium concentration C 0 as follows: i.e., c = 0 corresponds to the equilibrium concentration.
and equation (10) reads where we used the Peclet number which denotes the ratio between convection and diffusion. We also used the first Damkohler number which denotes the ratio between reaction and convection. Furthermore we will also use the Damkohler number II which denotes the ratio between reaction and diffusion.

Discussion of the scales
Two scales have to be chosen: the characteristic length R c and a characteristic interstitial velocity V c (or a characteristic time T c ).We are mainly interested in clogging around the inlet, which means that we will set R c = R w (well radius) in our computations for the field scale. For reasons of flexibility we will maintain a general R c , to allow other choices, like R c = R res (reservoir radius).
Solving equations (14) yields this means that V varies throughout the reservoir and multiple reasonable choices for V c can be made. However, choosing the character velocity as the velocity at R w , i.e., V c = V (R w ) implies that the constant K 1 = 1 which simplifies our expressions; in this work we will make this choice.
Using the characteristic velocity V c and the characteristic length R c we have the characteristic time Note that the characteristic time depends on whether we choose the reservoir radius or the well radius as our characteristic length. In Section 4 we will discuss the behaviour of four minerals (anhydrite, calcite, kaolinite, K-feldspar) under two conditions (low salinity and alkaline flood) both at low injection velocity/low Pe (named "field") and higher injection velocity/higher Pe (named "experiment"). The field and experimental parameters are given in Tables 2 and 3 ; the parameters concerning the four minerals are given in Table 4 .

Summary of the mathematical model
The governing equations are and ∂c ∂t + v ∂c ∂r We inject a solution at concentration C = C in j = C 0 at a constant injection velocity at r = r w , i.e., This means that we have a Danckwerts boundary condition for the concentration at the well r = r w where c is the dimensionless injection concentration For large values of Pe , the boundary condition (24) can be approximated by  Table 2 Parameters and scales (field).
well radius reservoir radius R res Pe 100 Table 3 Parameters and scales (experimental).
Pe 5000 which means that in this case the concentration is approximately the injection concentration at r = r w . At r = r res we have ∂c ∂r and at t = 0 we have the initial condition The dimensionless parameters Pe and Da I are given by Note that Pe is independent of the length scale chosen (because we , but that both Da I and Da II = Da I Pe depend on the length scale.

Final note on the length scale
Besides that the Da I depends on the length scale, the location of the boundary conditions depends on this choice as well. This changes the nature of the solutions of our problem.
If we use R c = R w as our characteristic length scale, our boundary conditions are at whereas if we take R c = R res as our characteristic length scale we have our boundary conditions at This means that because the functional form of the solution will change -the problems at these two length scales should be treated separately. In this work we will focus on the R c = R w choice.

Analytical steady-state solution and transients
In this section we will solve our problem (20) -(28) on the characteristic scale R c = R w .
The solution of the velocity profile is straightforward The solution of problem (20) -(28) is the superposition of a time independent steady state profile and a transient profile, i.e., which can be rewritten as The boundary conditions (24) Table 4 Mineral parameters low sal. Alaska at 60 C (field and experimental) and Alkali field, see [19] . The Damkohler numbers are calculated using equation (18) and Tables 2 and 3 .
and initial condition The solution of the steady state problem is already sufficient to answer our research question; computation of the transients is not necessary in this paper.

Analytical steady-state solutions of the model equations
Equations (35) -(36) determine c ss ( r ) for all values of Pe and Da II . Finding the complete explicit solution however is difficult due the large variation in the values of the Pe and Da II . We will subdivide the PeDa II − parameter space in three regimes and for each of the three regimes we will determine (an approximation to) c ss ( r ) using a different method.
Note that the assumption depends on the regime, i.e., we will have a different functional form of the approximation for each regime. The solution for Regime A is given in Section 3.1 , the solution for Regime B is given in Section 3.2 and the solution for Regime C is given in Section 3.3 . We then merge these results and discuss the full PeDa II -phasespace in Section 3.4 , using the appropriate approximations (A,B,C) for c ss in each point of the phasespace.

Regime a: Da II Pe 2
In Regime A the diffusion term d 2 c ss dr 2 can be neglected and equations (35) -(36) can now be simplified and solved (see Appendix A.1 for details) and the final result is as follows .
In Figs. 1 and 2 we plot c ss ( r ) for a fixed Pe = 10 5 and three different Da II numbers; Da II = 10 −4 , Da II = 10 2 and Da II = 10 8 . We observe three types of behaviour; for small Da II the whole reservoir is almost at the injection concentration, see the upper curve in Fig. 1 . for high Da II the whole reservoir remains at initial concentration, see Fig. 2 and note the different scales on the r -axis. For intermediate Da II we observe mixed behaviour, see the lower curve in Fig. 1 .
We will discuss this in more detail in Section 4 .

Regime b: Da II ~Pe 2
In regime B we use the following rescaling of the radial coordinate to solve equations (35) -(36) . We find the following equation for c ss (see Appendix A.2 for a detailed derivation) For small values of Da II we observe a non-horizontal profile (see Fig. 3 ); for larger values of Da II the entire reservoir remains at the initial concentration (see Fig. 4 ).

Regime c Da II Pe 2
In Regime C the convection term In this case we observe that the whole reservoir is at the initial concentration. Note that the width of the concentration layer ; this width may become in the order of the pore length, which means that application of our Darcy scale model is questionable in this part of regime C.

Phase space
In order to discuss the behaviour of c ss ( r ) we introduce the following two distances: r 0.99 and r 0.01 (see Fig. 6 ).
The distance r 0.99 is defined as follows: the reservoir is almost at the injection concentration (99% of the injection concentration or higher) until r = r 0 . 99 . This means that the whole reservoir is (eventually almost) at the injection concentration for large r = r 0 . 99 .
The distance r 0.01 is defined in a similar way: the reservoir is almost at the (zero) initial concentration (1% of the injection concentration or lower) after r = r 0 . 01 . Note that the parameter r 0.01 is also a measure of the width of the concentration profile, see Fig. 6 ). This means that the reservoir remains at the initial concentration if this width r 0 . 01 − 1 is very small, i.e., r 0.01 ≈ 1.
The distances r 0.99 and r 0.01 allow us to subdivide the full PeDa II -parameter space in four regions with distinct physical behaviour.
For some values of ( Pe, Da II ) the reaction is so slow that the whole reservoir is (eventually) at the injection concentration. We will call this subdomain of the ( Pe, Da II ) phase space "Region I"  ("no reaction"). Region I is characterized by high values of r 0.99 ; for our figures we will use r 0 . 99 = 100 as boundary for Region I.
High reaction rates lead to a situation where all reaction takes place at the inlet; almost no reaction takes place in the remainder of the reservoir, which remains at the initial concentration. We will call this subdomain of the ( Pe, Da II ) phase space "Region III" ("equilibrium"). Region III is characterized by values of r 0.01 close to one; for our figures we will use r 0 . 01 = 1 . 01 as boundary for Region III.
The range of values of Pe and Da II in between Regions I and III in the phasespace is called Region II ("kinetics"). In this part of the phasespace c ss ( r ) is neither at the initial nor at the injection concentration but varies between these values. In this Region II kinetics has to be taken into account.
Very high reaction rates may lead to a situation where the reaction takes places over distances over the order of a single pore R char r 0 . 01 ∼ 10 −6 m. In this case the Darcy-approximation may be questionable. We will call this subdomain of the ( Pe, Da II ) phase space Region IV ("Darcy-Brinkman").
In summary: Region I ("no reaction") The reaction rate is so low that the reaction can be neglected. The whole reservoir is eventually at the injection concentration.
Region II ("kinetics") Intermediate reaction rate; kinetics have to be taken into account in the whole reservoir. Region III ("equilibrium" High reaction rate; all reactions take place very close to the inlet, where the concentration varies from injection to initial concentration. Most of the reservoir remains at the initial concentration. Region IV ("Darcy-Brinkman") The width of the reaction zone becomes of the order of the poresize, which means that our Darcy-scale equations may not hold anymore.

Quantitative behaviour of minerals in three regions of the phase space
By way of example we choose three minerals that are representative for the three regions I, II, ad III. For all regions we introduce t eq , being the time needed to reach the steady state. For regions II and III, we also introduce the equilibrium distance r eq ≈ r 0.01 . This equilibrium distance r eq is the distance beyond which the reservoir is (almost) at the initial concentration, see Fig. 6 .
For an exact computation of the equilibrium time t eq the transient problem (37) -(39) has to be solved. We will not do so in this paper; instead we estimate t eq using the solution of the velocity profile (32) .
which means that we can associate a dimensionless equilibrium time t eq to an equilibrium distance as follows We will use equation (45) to estimate the equilibrium times in this section.
In regions I-III, we observe the following three types of behaviour Region I no reaction; the reservoir is (depending on r and t ) either at injection or at initial concentration. Region II r eq ≈ 1; kinetics have to be taken into account until t ≈ t eq . However, for t t eq the concentration profile may be approximated by the steady state profile c ss ( r ). Region III r eq 1; the whole reservoir remains at the initial concentration.

Discussion of K-Feldspar (regions i and II)
For K-Feldspar we have Da II = 1 . 0 · 10 −7 (low sal) or Da II = 2 . 0 · 10 −9 (alkaline pH > 7) (see Table 4 ). Using Fig. 7 (yellow dots)  Fig. 7. The full PeDa II -phase space is subdivided in four regions, viz., region I (reaction can be neglected), region II (kinetics has to be takes into account), region III (for the other two cases r 0.99 is even bigger). This corresponds to a dimensional distance R 0 . 99 = 450 m, which is of the order of the size of the reservoir: eventually the whole reservoir is at the injection concentration. We can approximate the full time dependent concentration profile (using equation (44) by This means that at, e.g., r = 10 (corresponding to 1 m) the concentration is equal to zero until t = 99 2 (corresponding to T = 99 2 T c ≈ 12 days) and equal to c afterwards.

Discussion of kaolinite (region II)
For Kaolinite we have Da II = 2 . 0 · 10 −3 (low sal) or Da II = 5 . 0 · 10 −6 (alkaline) (see Table 4 ). Using Fig. 7 (purple dots) we observe that Kaolinite is in Region II in the PeDa II -phasespace. This is confirmed by the computation of r 0.99 ; using equation (40) we find for, e.g., Pe =  is reached. The time required to reach the (spatially dependent) equilibrium solution can be estimated using equation (44) by which corresponds to T eq = 2 . 5 · 10 5 days. This means that kinetics has to be taken into account.
The time required to reach the (spatially dependent) equilibrium solution can be estimated using equation (44) which corresponds to T eq ≈ 0.05 days, i.e., roughly 1 h. This means that for 0 ≤ T ≤ T eq kinetics need to be taken into account. For T ≥ T eq the steady state saturation profile (42) may be used. For the alkaline case we find ( Pe = 10 2 and Da II = 2 . 0 · 10 8 ) using equation (43) r 0 . 01 = 1 . 0 0 03 , which means that Calcite remains at the initial concentration in this case.

Discussion of anhydrite (region III)
For Anhydrite we have Da II = 2 . 7 · 10 7 (low sal) or Da II = 1 . 0 · 10 7 (alkaline) (see Table 4 ). Using Fig. 7 (red dots) we observe that Anhydrite is in Region III in the PeDa II -phasespace. This is confirmed by the computation of r 0.01 ; using equation (43) we find for Pe = 10 2 , Da II = 1 . 0 · 10 7 that r 0 . 01 = 1 . 001 (and smaller for the other cases), which means that Anhydrite remains at the initial concentration in all cases.

Conclusion
Reactive flow modeling involves vastly different reaction rates, i.e., differing by many orders of magnitude. Consequently some reactions proceed so slow that in the time frame considered they can be disregarded. Other reactions occur so fast that they are well described by thermodynamic equilibrium in the time and spatial region of interest. Very fast and very slow reactions can be decoupled from the reactions that occur at intermediate rates; at intermediate rates kinetics needs to be taken into account. It is possible to categorize selected reactions as slow, fast or intermediate. We model 2D radially symmetric reactive flow with a reactionconvection-diffusion equation. Using the solution of the steady state profile we subdivide the PeDa II phasespace in three regions. Region I (slow reaction); reaction can be ignored and minerals are convected along with the flow. The minerals are either at the injection or initial concentration, depending on location and time. Region III (fast reaction); all reaction takes places in a very narrow region around the injection point. The reservoir remains at the initial concentration beyond this narrow region. Kinetics can be neglected in (almost) all of the computational domain. Finally, Region II (intermediate reaction); initially kinetics need to be taken into account. Eventually, beyond an equilibrium time a spatially dependent concentration profile is reached; after this equilibrium time, which can be much longer than the time of interest this concentration profile can be used; so in general kinetics has to be taken into account. We show in explicit examples, the location in phase space of a few selected minerals depending on salinity, pH and temperature. The computation of the Pe and Da II numbers for a (any) specific process makes it possible to determine in which of the Regimes I-III a process is, i.e., for which time window kinetics has to be taken into account. We note that the conditions, e.g., pH, salinity and temperature may be essential for assigning the reaction to the correct region in phase space. The methodology described can be applied to any mineral precipitation/decomposition problem and consequently greatly simplifies the model equations.

A1. Concentration profile in regime a
In Regime A the diffusion term 1 which means that this contribution can be neglected; solving for K yields K = c e α which means that we obtain .
We will now first show that the diffusion term can indeed be neglected. Using the solution (A.3) we find dc ss which is satisfied due to our restriction on Regime A.

A2. Concentration profile in regime b
In regime B we use the following rescaling of the radial coordinate Due to the second boundary condition (36) we have K 2 = 0 which means that we find Using the first boundary condition (36) we can determine K 1 and we find c ss = c

A3. Concentration profile in regime c
In Regime C the convection term , which is indeed satisfied due to our assumption on regime C.

Appendix B. Analysis of the PeDa II -phasespace
Using the notion of the distances r 0.99 and r 0.01 , we can subdivide the ( Pe, Da II )-phase space in four regions.
The first distance r = r 0 . 99 is associated with Region I and is defined as follows c ≥ 0 . 99 c for 1 ≤ r ≤ r 0 . 99 , (B.1) i.e., the reservoir is at injection concentration until r 0.99 . This means that if r 0.99 ≈ r res the entire reservoir is (eventually) almost at injection concentration. We define our Region I now as follows: where the dimensionless reservoir size r res = 10 3 with our choice of scales. The second distance r 0.01 is defined similarly i.e., the reservoir is at initial concentration beyond r 0.01 . This means that if r 0.01 ≈ 1 most of the reservoir is at the initial concentration. We define our Region III now as follows: The region in between Region I and Region III is Region II, i.e., Region II = { (P e, Da II ) | r 0 . 01 > 1 + 10 −2 , r 0 . 99 < r res } . Using the explicit analytical expressions for c ss derived in Section 3 we will subdivide the full PeDa II -phase space in each of the three regimes in (up to) four regions in the following subsections.

B1. Regime a
In regime A we have Da II Pe 2 ; for our numerical pictures we will use Da II ≤ 10 −2 Pe 2 . We will first compute the boundary of Region I. Using the explicit formula (A. 3) c ss = c e −α(r 2 −1) , α = Da II 2(P e − 1) and combining the definition of r 0.99 and the definition of the boundary of Region I as follows c ss (r = r 0 . 99 = 10 3 ) = 0 . 99 c , we can derive a relation between Da II and Pe on the boundary between regions I and II, i.e., Da II = −2 ln (0 . 99) This means that Region I is below the line Da II = 2 · 10 −8 Pe .
A similar computation yields the boundary between Regions II and III, where we have We finally compute the boundary between regions III and IV, Combining these three (blue) boundary lines we obtain the phase space in Fig. B.9 . The red line Da = 10 −2 Pe 2 denotes the boundary of Regime A; the equations in this subsection are only valid below the red line.

B2. Regime b
In regime B we have Da II ~Pe 2 ; for our numerical pictures we will use 10 −2 Pe 2 < Da II < 10 2 Pe 2 . We will first compute the boundary between Regions I and II. Using the explicit formula (A.12) :

B3. Phasespace in regime c
In regime C we have Da II Pe 2 ; for our numerical pictures we will use Da II > 10 2 Pe 2 . We will use the explicit formula (A. 17 to compute the boundaries between Regions II and III and between Regions III and IV; the boundary between Regions I and II is out of the range of our interest of ( Pe, Da II ) values. We observe that the boundaries are at a constant Da II number. For the boundary between Regions II and III we have ln (0 . 01) = − Da II · 10 −2 ⇒ Da II ≈ 2 · 10 5 and for the boundary between Regions III and IV we find similarly ln (0 . 01) = − Da II · 10 −5 ⇒ Da II ≈ 2 · 10 11 . These boundaries are denoted by the green lines in Fig. B.11 . Furthermore the red line Da = 10 2 Pe 2 denotes the boundary of Regime C. Note that the "Region I" in Regime B is out of the range of our interest of ( Pe, Da II ) values.

Supplementary material
Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ijheatmasstransfer.2020. 119969