Detailed numerical evaluation of diffusion convection equation in layered reservoirs during tracer injection

Characterization of heterogeneous reservoirs such as multilayered or fractured systems is an important issue in different disciplines such as hydrology, petroleum and geothermal systems. One of the popular methods that can be used for this purpose is tracer tests. Better understanding of the mechanisms of mass transfer (convection–diffusion process) is essential for having a proper test interpretation. In this study, the solutions of different scenarios of tracer flow in a pair of high and low-permeable layered reservoirs including convection and diffusion mechanisms are discussed. Although analytical solutions generally provided exact solutions, they involve several assumptions and might be hard to use for complex problems. As a result, numerical methods are selected for the investigation of different scenarios and addressing cases that are beyond access of analytical methods. In this study, several scenarios of considering diffusion and convection in low and high permeable zones and effective parameters on tracer concentration are investigated. According to the results of this study, the higher the porosity ratio of low to high permeable layer, the more time is needed to get the final concentration value. Also, by increasing the value of the dispersivity coefficient, the time needed to increase the concentration decreases. In other words, the sharp increase in concentration for lower times is seen in higher dispersivity values. The concentration profile variation is affected by Peclet number. The difference among concentration profiles in different cases is considerable, especially in low Peclet numbers where the diffusion mechanism is dominant. This behavior is more common in low permeable mediums such as multilayered tight or shale reservoirs.


Porosity ratio of low to high permeable layers
Essentially reservoir heterogeneity is a variation in reservoir properties as a function of space and this challengeable factor exists in a lot of underground reservoirs.Therefore, identifying reservoir properties and also the mechanisms of fluid flow in heterogeneous systems is a vital goal of reservoir engineers for better characterization and reservoir management 1 .Any prediction of reservoir performance in these types of reservoirs based on limited information may not be true in some cases and should only serve as a starting point for analyzing the mediums.Multilayered and fractured reservoirs are two common types of heterogeneous systems.In these reservoirs, the fluid flow will be done in two or more different environments.Fractured reservoirs are composed of two different medium, including a matrix system with high porosity and low permeability, as well as fracture network with low porosity and high permeability in which reservoir fluid may be flow between these two mediums.On the other hand multi-layer reservoirs can also include different environments with porosity and permeability distribution, where their fluid flow calculation may be somewhat similar to fractured reservoirs.Therefore, a better understanding of the mechanism of mass transfer (convection-diffusion process) in these reservoirs and also the identification of affected parameters to flow is an important topic in fluids flow phenomena.Besides the different methods of reservoir characterization, the tracer flow tests are the methods for the better investigation of the mass transfer process in heterogeneous reservoirs.Different investigations have been done on tracer flow in heterogeneous reservoirs such as naturally fractured and multilayer reservoirs [2][3][4][5][6][7] .
To better evaluation of multilayered reservoirs as a type of heterogeneous system, Brigham and Smith proposed a tracer flow model in a water flooding process.This model was able to describe permeability heterogeneity in layered reservoirs 8 .Yuen et al. presented a computer algorithm to determine the degree of heterogeneity in these reservoirs using Brigham equations 9 .Abbaszadeh-Dehghani and Brigham (1984) reported analytical expressions to analyze tracer flow in layered reservoirs using a nonlinear optimization technique 10 .Sato and Abbaszadeh used tracer pulse test response as a function of layer number to evaluate the maximum produced concentration from homogeneous and multilayered reservoirs 11 .Sometimes later Samaniego et al. presented new models to evaluate tracer flow behavior in heterogeneous reservoirs using dimensionless parameters 12 .Shook used a flow-storage capacity diagram to get an idea of flow geometry, swept volume, and saturation distributions in the layered reservoir 13 .
Suarsana et al. used the comparison of tracer test results and analysis of connectivity injector and producer to investigate pilot water flooding performance in a multilayered reservoir 14 .Jason et al. show that changes in tracer concentration gradients can be good indicators of changes in porosity (or water saturation) between layers in multilayered reservoirs 15 .Shen et al. presented a new formulation to evaluate multilayered reservoir behavior regarding limited crossflow 16 .Bahamon and Mora described a method to evaluate water channeling and assess sweep efficiency improvements in a multilayered high water cut field using an inter-well tracer test to better perform the water flooding process 17 .Davidescu et al. used a tracer test in the polymer and water injection process to evaluate the sweep efficiency of a multilayered reservoir 18 .
Some investigations also worked on tracer flow in the matrix-fracture mediums in fractured reservoirs as other types of heterogeneous systems.In 1986 a mathematical model was developed to estimate fracture width using tracer test results.The results showed good agreement in comparison with laboratory data 19 .After that Raven et al. investigated the tracer flow in a fractured geothermal reservoir by introducing the Peclet dimensionless number to evaluate the time of tracer breakthrough in matrix fracture media 20 .Ramirez et al. investigated the response of the tracer test in matrix fracture medium and tried to estimate fracture opening and also matrix diffusion coefficient in different geometries 7,21 .Samaniego and Rodriguez used the tracer data in an observation well to estimate some reservoir parameters in a fractured reservoir using a trial method 22 .Qasem et al. showed that parameters such as porosity and permeability and their relationship significantly affect tracer results in fracture-matrix systems 23 .After that Samaniego et al. in 2005 presented a model to investigate the tracer flow in the homogenous and heterogeneous reservoir.They reported that some reservoir parameters such as block size and also dispersion coefficient can be estimated using tracer analysis 12 .Kocabas and Maier used analytical and numerical models to investigate the tracer flow in a fractured reservoir 24 .Haddad et al. used analytical methods to estimate the dispersion coefficient in different geometries of matrix-fracture media by solving the convection-diffusion equation 25 .Aymen et al. used the tracer response in analytical and numerical models to estimate water cut in fractured reservoirs 26 .Abbasi et al. used tracer data and discussed the analytical solution of the mass transfer equation to evaluate the shape factor in a fractured reservoir 27 .Jing et al. estimated the future performance of an ultralow-permeability fractured reservoir using tracer test and production data 28 .Kumar et al. reported that the combination of tracer data, completion data and also stimulation data in a machine-learning framework could show better dynamic behavior in a fracture network system 29 .Also a brief history of activities performed in the analysis of tracers is presented in Appendix C 11,12,16,18,24, .
In this study, first, an introduction to the general convection-diffusion equation is presented.Then, tracer flow in a two-layered system as a heterogeneous reservoir is evaluated in three scenarios.In the first scenario (Base Case), convection in higher permeable and diffusion from higher to lower permeable layers are dominant mechanisms.In the second case, the diffusion toward the lower permeable layer (L2) is also included in the evaluation.Finally, the convection mechanism toward L2 is added to the calculations.In this work, three different cases of mass transfer in tracer injection with different boundary conditions are investigated.Because only the analytical solution of the simplest case is available in the literature, the numerical methods are considered to solve mass transfer equations in these cases.Accordingly in the next section the results of analytical and numerical methods are compared for base simple case to check the validity of numerical methods.After that, the sensitivity analysis of affecting parameters in the convection-diffusion equation is done.Finally, a comparison among different cases of tracer flow in these layers is investigated.

Methodology General convection-dispersion equation
According to physical phenomena, a solute (such as tracer) molecular diffusion occurs from higher to lower concentration.Diffusion in porous media can be explained using Fick's first law.The mass of diffused fluid is related to the concentration gradient as below: in which F, D d and ( dC dx ) are mass flux, diffusion coefficient and concentration gradient, respectively.On the other hand, if concentration changes with time, Fick's second law is used to show the diffusion variation as a function of space and time.This equation in one dimension may be expressed as: Besides the diffusion phenomena, convection is another important mechanism for solute transport in porous media.The one-dimensional advection transport equation is in the following form where v x is the average linear velocity.
According to studies done by Ogata and Bear, the convection-dispersion equations in representative elemental volume (REV) of a porous medium are in the following forms 58,59 : where n e , dA, and i are effective porosity, the cross-sectional area of the element, and the normal direction to cross-section, respectively.The total mass of solute per unit cross-sectional area transported in the i direction per unit time, Fi, is the sum of the advective and the dispersive transports are given by: The amount of solute entering and leaving the representative elemental volume is: The mass of solute variation in representative elemental volume also can be expressed as: Using the mass conservation law, the rate of mass change in the representative elemental volume must be equal to the difference in the mass of the solute entering and leaving: Substituting Eq. (6) into Eq.(10) and canceling ne from both sides yields the general 3D convection-diffusion equation in Cartesian coordinate: The above relation is a 3D mass transfer equation in porous media without any chemical reaction.In a homogeneous medium, D x , D y , and D z do not vary in space.However, because the coefficient of the hydrodynamic dispersion is a function of the flow direction, even in an isotropic, homogeneous medium, then D x #D y #D z should be considered in calculations.
(1) www.nature.com/scientificreports/If dilution of solute at the advancing edge of flow occurs, the mixing will be mechanical dispersion.In this condition, the mixing which occurs along the advancing edge of flow is longitudinal dispersion and the mixing in the direction normal to the flow path is called transverse dispersion.
If the mechanical dispersion is described by Fick's law for diffusion and it is assumed as a function of the average linear velocity, then the coefficient of mechanical dispersion will be calculated by introducing the dispersivity.Dispersivity denoted by α is a property of porous media.Dispersivity times the average linear velocity i.e. α.v will be called the coefficients of mechanical dispersion The combination of this parameter and molecular diffusion (D*) is called hydrodynamic dispersion coefficient which is represented by the following formulas; where D L and D T are hydrodynamic dispersion coefficient parallel (longitudinal) and perpendicular (transverse) to the principal direction of flow.
For those domains where the average linear velocity v i is uniform in space, Eq. ( 11) for one-dimensional flow in a homogeneous isotropic media will be presented in the following form where only longitudinal dispersion coefficient (D L ) will be used: Also for two-dimensional flow with the direction of flow parallel to the x-axis, Eq. ( 11) will be presented in the following form using longitudinal and transverse dispersion coefficients: On the other hand, the radial flow from a well in polar coordinate can be written as bellow 59 : where r and V are the radial distance of the well and average pore velocity respectively.In this study, the above equations will be discretized and solved numerically to investigate the convection-diffusion mechanisms in heterogeneous reservoirs.These types of reservoirs may be presented by high and lowpermeable systems such as Matrix-Fracture or multilayered media.

Numerical modeling of convection-diffusion flow in layered reservoirs
Numerical modeling is a powerful tool to evaluate physical phenomena such as convection diffusion in heterogeneous reservoirs.Although the analytic solution of an equation or expression is the exact solution and guarantee that the resulting quantitative predictions are accurate, due to some complicated boundary conditions and limitations of analytical methods, numerical calculations are more useful.In the following sections of this study, numerical methods are used to evaluate physical phenomena such as convection diffusion in heterogeneous reservoirs.This study investigates the convection-diffusion mechanism of two-layered reservoirs with different degrees of contrast between layer's permeability.Different cases of flow in this reservoir are introduced to solve by numerical modeling.
First of all, it is assumed that tracer flow is occurred in two-layered systems from the wellbore toward the higher permeable layer (L1) by convection mechanism.Also, diffusion in vertical direction is the dominant mechanism for mass transfer from the higher permeable layer to the lower permeable layer (L2).This case may be similar to a dual porosity system in fractured reservoirs.In this case (base case), diffusion exist in L1 but it is not the dominant mechanism.
In the second case, the diffusion toward the lower permeable layer (L2) is also included in the evaluation.Finally, the convection mechanism toward L2 will be added to the calculations.
To solve the convection-diffusion equation, first of all, the dominant mechanism in L2 will be assumed to be the diffusion from L1 in the vertical direction using the following formula where subscript numbers 1 and 2 represent the higher and permeable layer properties respectively: where D2 is the diffusion coefficient in L2.Also, the general convection-diffusion equation in L1 will be presented in the following form and is similar for all cases.The equation is in the one-dimensional system in radial coordinates.The terms on the right-hand side of the equation are velocity-dependent.( 12) www.nature.com/scientificreports/ In the above equation the terms φ 2 φ 1 , h b/2 , α, q o , h and D 1 represent the porosity ratio of L 2 to L 1 , thickness of layer L2, dispersivity, production rate, thickness of L1 and diffusion coefficient in L1 respectively.The derivation of this equation is presented in Appendix A.
If the diffusion coefficient term is assumed to equal zero in high permeability medium which is not unusual, then the equations will be presented in the following short form: To better solve the tracer flow in porous media according to convection-diffusion equations, the following dimensionless parameters were introduced: The terms, C 1D , C D2 , and t D represent dimensionless concentrations of L1 and L2 layers and dimensionless time respectively.Also, h R is the dimensionless ratio of L2 layer thickness to wellbore radius and Z D represents the dimensionless ratio of the vertical coordinate of the L2 layer to wellbore radius.The terms α D and D 2_1D show dimensionless dispersivity and dimensionless diffusion coefficients of L1 to L2 respectively.On the other hand, to define the relation between convection and diffusion mechanisms, dimensionless Peclet number is introduced.Peclet number (P e ) is the ratio of advective velocity to the molecular diffusion coefficient.Fundamentally the value of this number is lower in regions of lower permeable than higher permeable layers.
According to the literature, five dispersion flow regimes can be described in different mediums.These regimes may be divided according to the value of the Peclet number and its dependency on the experimentally estimated dispersion ratio.For example for sand-pack medium, Fig. 1 shows the five dispersion regimes as follows: 1. Pe < 0.3 presents the diffusion regime in which convection is not the case and diffusion is dominant.2. The regime 0.3 < Pe < 5 is the transition zone in which the effect of convection is important but the diffusion effect is strong.3.In the range 5 < Pe < 300 the convection dominates dispersion, but the effect of molecular diffusion cannot be neglected.4. When 300 < Pe < 10 5 the flow regime is the purely convective regime.5.For Pe > 10 5 dispersion is in the turbulent regime in which the Peclet number is no longer the only correlating parameter, as the Reynolds number should also be considered.In this condition or flow through porous media, this regime is not of interest.
In this study, the relation between the velocity of the high permeable layer (L1) and diffusion between high and low permeable layers (L2) will be presented with Pe 2_1 as defined in Eq. (28).On the other hand, the relation between convection and diffusion in the low permeable layer will be introduced with Pe 2 in Eq. (46).According

Introduction to different cases of flow in numerical modeling
In this section, the convection-diffusion equation for two layered defined systems will be presented in three following cases named I, II, and III.In all cases, it is assumed that the convection from L1 to L2 is ignored.Also, the diffusion in L1 wouldn't include in calculations due to its low effect on mass transfer.It is considered that Eq. ( 31) will be solved numerically as the mass transfer equation of L1 for all of the following cases.Table 1 describes introduced mass transfer mechanisms in defined cases. ( Dependence of the dispersion coefficient ratio on Peclet number Pe and the various dispersion regimes 60,61 . Table 1.Description of mass transfer mechanisms in defined cases.The schematic of diffusion-convection flow in pours media, in this case, is shown in Fig. 2A where h high_k and h low_k are the half length of L1 and L2 layers respectively.In this figure, the convection exists in L1 and the diffusion is considered to exist form L1 to L2.As is presented in Fig. 2B after discretization and introduction of the coefficient and the known matrixes, the unknown matrix will be solved using numerical methods.The figure shows a sample of grids discretization in radial (r) and vertical (z) directions respectively.The number of total grids will be calculated using the following formula:

I
where N low_k and N high_k are the number of grids in vertical and radial directions for lower and the high permeable layers respectively.The total number of grids in high permeable layer equals N high_k and the total number of grids in the low permeable layer equals NN-N high_k .It is assumed that L1 and L2 are divided by 10 grids for low permeable layer and 20 grids for high permeable layer in Fig. 2B and a total of 220 grids according to above formula.In this condition the location of each grid will be described by (r,z).It is assumed that the numerical calculations will be done for a pair of higher-lower permeability layers.In this case, Eq. ( 30) will be used to solve diffusion from L1 to L2 in which the initial and boundary conditions are in the following forms: www.nature.com/scientificreports/Also to solve the mass transfer equation in L1 according to Eq. ( 28), the initial and boundary conditions will be presented in the following form: The diffusion term may be included in the mass transfer of the higher permeability layer (L1), but it is not so important because the convection is the dominant mechanism in this layer.

Convection in L1-diffusion from L1 to L2-diffusion in L2
In this case, the diffusion toward the L2 is introduced as another existing mechanism of mass transfer in addition to diffusion from L1 to L2.The schematic of diffusion-convection flow, in this case, is shown in Fig. 3.
In this condition the mass transfer equation in L2 can be shown in the following forms: Also, the initial and boundary conditions (BC) in this case are presented below: In addition to boundary conditions in the previous case, two new B.C's (Eqs.44 and 45) will be added to numerically solve the mass transfer equation.According to new B.C., the tracer concentration near the wellbore (injection point, r D = 1) is equal to one and its value in far distances from the injection point will be zero.

Convection in L1-convection in L2-diffusion from L1 to L2-diffusion in L2
In this case, it is considered that the convection also exists in L2 as a tracer flow mechanism in addition to before mentioned ones.In the discretization process, the boundary condition of tracer flow from the wellbore to its adjacent grid is different from other grids.In this case, the rates of tracer flow from the wellbore toward the L1 and L2 will be different.Due to this difference, a new dimensionless Peclet number (P em ) will be introduced to describe the convection flow of the tracer from the wellbore toward the lower permeability layer (L2).This case is similar to the dual permeability mechanism in fractured reservoirs.The initial and boundary conditions are similar to the previous case.The schematic of flow in case III is illustrated in Fig. 4.
The definition of the new Peclet number and also general convection-diffusion equation in L2 will be presented in the following forms: To solve the mass transfer equation numerically in the three mentioned cases, some assumptions were considered in the modeling.The time of simulation (dimensionless time) is assumed to be 10 4 .The value of the porosity ratio (the porosity of L2 to L1 ratio) is considered to equal five.The value of time difference (dt) is assumed to www.nature.com/scientificreports/be equal to 0.5, the value of horizontal changes (dz) equals to 0.5 and the value of the radial change (dr) equals to 5 in numerical modeling.Therock and fluid properties and also the porosity ratio of the low permeable to the high permeable layer is assumed to be constant; therefore, the difference in porosity ratio in different grids is ignored.On the other hand the properties of the grids are assumed to be constant with time, except for the changes in concentration.The dimensions of the grids are assumed unchanged in numerical calculations.The possible existence of diffusion mechanism from the low permeable layer to the high permeable layer has been neglected.The other assumed parameters are presented in Table 2 which they may be changed for sensitivity analysis in the following section.
On the other hand, other assumptions will be introduced individually in each case.The total number of discretization grids from the injection well to a production well is assumed to be 220 grids for L1 and L2 in which N low_k and N high_k are equal to 10 and 20 respectively according to Eq. ( 32).In our study, according to (r, z + 1) nomination in the schematic of Fig. 2B, the first grid number of the higher permeable layer will be described by grid (1,11) and the final grid of permeable layer will be shown by (20,11).In another word, the value of r changes from 1 to N high_k, and the value of z + 1 will change from 1 to N low_k + 1 in grid numbering.
The detailed procedure to solve the convection-diffusion problem is presented in Appendix B for a twelve (12) grid system in which (N low_k + 1 × N high_k) equals 12 grids according to Fig. 2.

Results and discussion
In this section, the result of convection-diffusion equations in the above-mentioned cases will be discussed.First of all, a comparison between analytical and numerical calculations will be evaluated to check the accuracy of numerical modeling.Then sensitivity analysis on affected parameters in the concentration profile and also the comparison between cases will be investigated in more details.

Comparison between analytical and numerical modeling of tracer flow in a multilayered reservoir
Although analytical solutions are generally considered to be more rigorous than numerical methods due to providing exact solutions, sometimes these methods might not be able to handle complex problems with specific assumptions.As it mentioned in previous section, three different cases with different boundary condition will be discussed in this study in which only the first case is solved using analytical methods in literature 25 .
It is obvious that in order to compare the numerical and analytical solutions, the values of φ 2 φ 1 , hr, αD, rD, and also simulation time are assumed to be equal in both solution methods.
The analytical solution of a mass transfer equation is solved in base case (Eqs.29 and 30) which is simple case, using the following formula in Laplace (S) domain 25 : In which the c D1 and c D2 represent the concentration values of higher and lower permeability layers in Laplace domain respectively.The above equations can be solved by introducing the Ai as Airy function.www.nature.com/scientificreports/Finally the above equations will be inversed in time domain using the Stehfest algorithm.Due to complexity of assumptions in other cases, in this study, the numerical methods are selected to use for the investigation and evaluation of mass transfer equations in the three mentioned cases, but due to check the accuracy, the first case (Convection in L1-diffusion from L1 to L2) is solved by both the analytical and numerical methods.Figure 5 shows the comparison of the result where the average difference between the analytical and numerical values of calculated concentration is about 0.13%.Also, the maximum difference between the two methods is about 9.4% which shows good agreement and a reasonable range of difference.

Numerical investigation of convection-diffusion mechanisms
Case I: (Convection in L1-diffusion from L1 to L2) Figure 6 shows the dimensionless concentration-time profile after numerical solving the mass transfer equation in the higher permeable layer (L1) regarding the convection mechanism.According to this figure, the value of concentration changes from its initial value (i.e.zero) to the final value (i.e.one) in different discretization grids.The figure shows that in the first grid (near the wellbore) the concentration rapidly yields the final value while its value in other grids gradually increases.It can be said that in a special time, the more the distance from the wellbore, the lower the concentration value of discredited grids will yield.The concentration value in all grids (51)    becomes unity sometimes later which these times are sensitive to affected parameters in the convection-diffusion equation.
As it is shown in Fig. 7, the higher the porosity ratio, the more time is needed to get the final concentration value.It can be said that in cases of higher porosity ratio, due to higher pore volume in L2, it will need more time to yield the final value of concentration in L1 and L2.
Also, this figure shows that the time to reach the final value of concentration in L1 will be increased in higher numbers of the discretized grid for all values of porosity ratio.As an example, it is shown that the concentration profile will yield faster to the final value in the numerical grid (5, 11) in comparison to the numerical grid (20,  11) because of the lower distance to the source of tracer injection.
On the other hand, Fig. 8 shows that by increasing the value of the dispersivity coefficient, the time needed to gradually increase the concentration will decrease i.e. the higher the dispersivity coefficient, the sharp and fast increase in concentration occurs in lower times.
The value of the Peclet number is another important affected parameter in concentration-time profile evaluation.By increasing the value of the Peclet number, the concentration will be faster received to its final value.Also due to less effect of the convection mechanism in lower values of Peclet number, the time to reach the final concentration value is increased in low Peclet numbers.This behavior is depicted in Fig. 9.
The variation of concentration versus location at different times for the first layer also is shown in Fig. 10.According to the curves illustrated in this figure, for a special location in the radial direction of tracer flow (discretized grids), the value of concentration will change at different dimensionless times.It can be said that the more time consume, the higher the value of concentration in a special location will occur to meet the final concentration value.
After discussing the concentration variation in L1, the profile of concentration time in L2 will be evaluated in the following figures.In this condition, only diffusion from L1 to L2 is included in calculations for this layer.
As is shown, the shape of the concentration profile in grids near layer L1 (Fig. 11a) is different in comparison to the grids far from this layer (Fig. 11b).It is also obvious that the trend of increasing concentration is faster in grids near L1 than in the grids for this layer.
Although the trend of concentration value approaches to unity in L2 like L1, because of the diffusion mechanism and also lower permeability values, the final concentration will be yielded at later times.

Case-II: (convection in L1-diffusion from L1 to L2-diffusion in L2)
In this case, the diffusion mechanism in L2 will be included in addition to diffusion from L1 to L2. Therefore two new boundary conditions will be used in numerical calculations of mass transfer for this layer.The concentration profile for the higher permeability layer (L1) is similar to the previous case.
In this case, adding a new diffusion mechanism may not have so significant effect on mass transfer calculation.It is because, in higher values of the Peclet number, the convection mechanism is dominant.In this condition, the concentration-time profile wouldn't be so different in comparison to case 1 as shown in Fig. 12.To better understand the importance of the diffusion mechanism, especially in L2, a sensitivity analysis is done on Peclet number values.To compare the results, two grid numbers are selected and the concentration-time profile will be evaluated for different Peclet numbers in cases I and II.To do this, the grid number (10, 1) will be selected.Because the effect of diffusion is so important in low Peclet number values, the values of 10, 100, and 1000 for this dimensionless number were considered for comparison.
A comparison of the results of the concentration-time profile in Case I and Case II is shown in Fig. 13.According to mentioned figure, the difference between concentration profiles in the two cases is considerable in a low Peclet number of 10 where the diffusion mechanism is dominant.This behavior is more common in low permeable mediums such as single/multilayered tight or shale reservoirs.
The above figure also shows that the effect of tracer diffusion in L2 is not considered when values of the Peclet number increase in the range of convection dominance.

Case-3: Convection in L1-convection in L2-diffusion from L1 to L2-diffusion in L2
In this case, the convection mechanism in L2 is also included in addition to other defined mechanisms in previous cases.Although the concentration-time profile for L1 is not the case of change, due to including the convection mechanism in calculations of L2, the concentration profile will be affected by a considerable change in comparison to the other cases.The results of this case can simulate the convection-diffusion flow mechanisms in the real porous medium.
According to Eq. ( 44), a new Peclet number will introduce and used to include the convection in L2.In this condition, the value of the Peclet number for L2 will be selected as equal to 100. Figure 14 shows the comparison of the concentration-time profile for cases I, II, and III.The difference in profile near the injection well (tracer source) is significant in case 3 compared to the two before-mentioned cases (Fig. 14A,B), but the difference www.nature.com/scientificreports/becomes less in the more distant locations (Fig. 14C,D).It can be said that the effect of including convection in the low permeable layer is considerable at initial locations, but the effect of diffusion from higher to lower permeable layer will be more important and affect the results in farther locations.

Conclusions
In the current study, a detailed numerical evaluation of the convection-diffusion equation in layered reservoirs is discussed using tracer response.The dispersivity coefficient, the effect of some related parameters to Peclet number, as well as changes in the concentration profile of the tracer, especially at low Peclet numbers, have also been investigated.On the other hand, the time of importance and dominance of Convection and Diffusion mechanisms in the low permeable layer is also one of the main subjects that was evaluated in this study.The main results of this paper are summarized as follows: • The diffusion mechanism is dominant and has a more visible effect in low permeable mediums such as single/ multilayered tight or shale reservoirs.Also, diffusion mechanism is more considerable in lower values of Peclet number.• As the porosity ratio of low-to-high permeability layer increases, more time is needed to reach the final con- centration value.Also by increasing the dispersivity coefficient, the sharp and fast increase in concentration occurs in smaller times.• Although the trend of concentration in the low and high permeable layers are similar to each other, because of diffusion mechanism and also permeability values, the maximum concentration in low permeable layer will yield at later times.• The effect of including convection in the low permeable layer is considerable at initial locations near the tracer injection point, but the effect of diffusion from the higher to lower permeable layer will be more important and affect the results in farther locations.

( 32 )Figure 2 .
Figure 2. Schematic of case I (A)/Layer discretization for a sample reservoir N L1 = 20 and N L2 = 10.Total number of grids in L1 and L2 are 20 and 200 grids respectively (B).

Figure 5 .
Figure 5.Comparison between concentration profile in analytical and numerical methods in semilog scale (Case1: convection in L1-diffusion from L1 to L2).

Figure 6 .
Figure 6.The concentration-time profile for L1.Grid numbering according to (r,z) gridding.Total number of grid in this condition is calculated according to Eq. (32).

Figure 7 .Figure 8 .
Figure 7. Sensitivity analysis of porosity ratio in L1 for numerical grids 5 and 20.The higher the porosity ratio, the more time is needed to get final concentration value.

Figure 9 .
Figure 9. Sensitivity analysis of Peclet number values in L1 for numerical grids 10 and 20.The convection mechanism is dominant in higher Peclet values.

Figure 10 .
Figure 10.Concentration-distance profile for L1 in different times.The more the time consuming, the higher the value of concentration in special location will be occurs to meet final concentration value.

Figure 13 .
Figure 13.The comparison between profiles of concentration in cases I and II.The difference between concentration profiles in two cases is considerable in low Peclet numbers.

Figure 14 .
Figure 14.The comparison of concentration-time profile for cases I, II and III.The difference in profile near the injection well (A, B) is higher than this profile in locations farther from injection point (C, D).