Influence of Grain Size Transition on Flow and Solute Transport through 3D Layered Porous Media

Soils and other geologic porous media often have contrasting grain size layers associated with a grain size transition zone between layers. However, this transition zone is generally simpli ﬁ ed to a plane of zero thickness for modeling assumption, and its in ﬂ uence has always been neglected in previous studies. In this study, an approach combining a deposition process and a random packing process was developed to generate 3D porous media without and with consideration of the transition zone. The direct numerical models for solving the ﬂ ow and concentration ﬁ elds were implemented to investigate the in ﬂ uence of the grain size transition on ﬂ ow and solute transport. Our results showed that although the transition zone occupied 13.6% of the entire layered porous medium, it had little in ﬂ uence on the distribution of ﬂ ow velocity at the scale of the entire layered porous medium. However, the transition zone had a signi ﬁ cant in ﬂ uence on the local ﬂ ow ﬁ eld, which was associated with the increased spatial variability of velocity and the varied distribution of ﬂ ow velocity. This varied local ﬂ ow ﬁ eld could increase the solute residence time and delay the breakthrough time for solute transport. Although using both the advection-dispersion equation (ADE) and the mobile and immobile (MIM) models to ﬁ t the breakthrough curves (BTCs) for solute transport through layered porous media resulted in trivial errors, the ADE model failed to capture the in ﬂ uence induced by the local ﬂ ow ﬁ eld, especially the in ﬂ uence of the transition zone. In contrast, the MIM model was shown to be able to capture the in ﬂ uence of the transition on solute transport. It was found that the mass transfer rate α , a parameter of the MIM model, was signi ﬁ cantly improved by the presence of the transition zone, while it decreased as the transition zone fraction increased. Our study emphasized that the transition zone can vary the local ﬂ ow ﬁ eld at the pore scale, while it has little in ﬂ uence on the hydraulic properties (e.g., hydraulic conductivity) of the macroscale ﬂ ow ﬁ eld. However, the local ﬂ ow ﬁ eld varied by the transition zone has a signi ﬁ cant in ﬂ uence on solute transport.


Introduction
Characterization of solute transport through layered porous media plays an essential role in many environmental, industrial, and engineering processes, including exploitation of deep geothermal energy, groundwater contamination remediation [1], geothermal energy recovery [2], enhanced oil recovery [3], and geological repository [4]. Typical examples of layered porous media systems are aquifer-aquitard systems [5][6][7][8] and fracture-rock matrix systems [9][10][11][12]. Layers exhibit different pore space structures and hydraulic properties due to grain size contrast. In the zone near the layer interface, the grain size changes from one side to the other [13,14]. The transition of grain size between layers is due to the complex geological process [15][16][17][18]. The deposited grains in different geological periods may mix, resulting in a grain size transition zone with a small thickness. However, the grain size transition zone in layered porous media is generally neglected and roughly assumed to be a "zero thickness" plane due to its small thickness. This rough assumption may be reasonable for estimating the hydraulic properties of the layered porous media at the field or laboratory scale [19]. However, it is inappropriate for characterizing the solute transport process which is largely controlled by flow behavior at the pore scale.
Pore space structure defines how pore space is distributed over an area. The heterogeneity of geological formations is ubiquitous and is controlled by the geometric properties of the pore space structure. The layered porous medium is the simplest type of heterogeneity. The heterogeneity causes multiscale velocity field distributions [20], which induces anomalous diffusion transport in the context of transport in many systems [21,22]. For example, a typical representative of this phenomenon is the preferential flow path formed by a series of continuous pores as it provides a high-velocity path for fluid flow and solute transport [23]. Moreover, the low-velocity domain controls the contaminant accumulation and release [24,25], and the mass transfer between a stagnant domain and flow domain could lead to long-tailed breakthrough curves (BTCs). In the zone near the interface between layers, especially in layered porous media where the grain size varies greatly in different layers, a complex flow field and anomalous mass transfer have been observed in several studies [26][27][28]. Therefore, any slight variation in the pore space structure of the transition zone could significantly affect the solute transport process. However, little attention has been paid to the transition zone and its effects on solute transport.
The relationship between the pore space structure, fluid flow, and solute transport process in heterogeneous porous media is essential [29][30][31][32]. The mathematical models developed based on the macroscopic scales do not consider the effects of pore space structure. Accurately measuring or describing the small-scale pore space structure (e.g., connectivity, pore size, and pore shape) of natural formations remains a challenge. With the development of computer performance, direct simulations, such as standard Computational Fluid Dynamics (CFD) [33] and Lattice Boltzmann Method (LBM) [34,35], have been proposed. These methods are commonly used to investigate the flow and transport processes in porous media reconstructed from actual samples [36,37]. Direct numerical simulation is able to provide information about the internal state of the porous medium, such as the fluid flow and the basic solute transport processes occurring at the pore scale [38]. This information is beneficial to analyze the solute transport process and relate it to the pore space structure [39]. Willingham et al. [40] investigated the effects of grain size, grain orientation, and intraparticle porosity upon the extent of a fast bimolecular reaction using LBM. They found that consideration of grain size alone was not sufficient to explain mixing at the pore scale. Voller [41] illustrated the relationship between the nature of heterogeneity in a given system and the emergence of anomalous transport through direct numerical simulations. Swanson et al. [42] estimated transport model parameters via calibration of a transport model to column-scale solute tracer tests and compared them to results from geophysical or material characterization methods that measure the physical properties of porous media.
The characteristics of non-Fickian transport behavior (e.g., early breakthrough and long tail) in porous media could be quantified or accurately estimated by fitting observed BTCs with various models [43]. Traditional Fickian advection-dispersion equations (ADEs) do not account for heterogeneity in pore structure and thus cannot capture anomalous transport characteristics caused by the influence of the transition zone. To capture the non-Fickian transport behavior, many types of time nonlocal transport models have been proposed. These models include the Continuous-Time Random Walk (CTRW) framework [44], the Tempered Time Fractional Advection-Dispersion Equation (tt-fADE) model [45], and the nonequilibrium mobile and immobile (MIM) model [46]. The nonequilibrium mobile and immobile (MIM) model has been commonly used to describe solute transport through layered porous media [47][48][49]. Recently, some studies [12,50] have attempted to link the MIM model parameters and pore space structure to explain non-Fickian transport behavior. de Vries et al. [51] investigated the effects of a wide range of key parameters (e.g., porosity and permeability) obtained by varying the pore throat sizes within the porous medium on the main parameters of the MIM model. Their results showed that the micropore structure of the porous media can affect the mass transfer coefficients, but the mechanism behind these effects needs further investigation.
The aim of this study was to investigate the difference in flow and solute transport through layered porous media under the condition with and without considering the transition zone. First, an approach combining a deposition process and a random packing process was developed to generate three-dimensional (3D) layered porous media without and with considering the transition zone at the pore scale. The geometric properties of the pore space in developed layered porous media were estimated to support the analysis. Two parameters (i.e., pore throat and pore size) relating to pore space structure were proposed and calculated. The flow and solute transport processes were simulated directly by coupling the Navier-Stokes (N-S) equation and the advection-diffusion equation. The flow fields in both the transition zone and the entire layered porous media were analyzed. The spreading of the solute plume through the layered porous media was measured using the spatial moment method. The BTCs were fitted with the MIM model and the ADE model to estimate different parameters reflecting the characteristics of the solute transport process.

Generation of Layered Porous Media.
In order to investigate the influence of grain size transition zone on fluid flow and solute transport through layered porous media, the 3D layered porous media were generated. The generation approach of layered porous media combines a deposition process [52,53] under the force of gravity and a random packing process [54]. A sedimentation algorithm was adopted to simulate the deposition of solid grains. This algorithm contains three main steps. First, a spherical grain represents the coarse solid grain that is generated and randomly placed on the top of the computational domain. It should be noted that the computational domain for grain deposition is larger than the computational domain for fluid flow. Second, 2 Lithosphere under the force of gravity, the coarse solid grain continuously sinks until it reaches the bottom or contacts the other three settled solid grains. When the deposition process ends, the sinking solid grains reach a local minimum of their potential energy. The above steps were performed circularly until the coarse grains exceed the middle plane of the computational domain. Third, the sinking coarse grain is switched to a fine grain and continually sinks until the computational domain is filled up. Finally, the original layered porous media were generated. Figure 1 shows the original layered porous media, which is the layered porous media without considering the transition zone and hereafter referred to as PM0. In PM0, the coarse grain layer and the fine grain layer were separated by a "zero thickness" middle plane directly. The probability distributions of the grain radius are shown in Figure 1(b). The radii of the coarse or fine grains follow the truncated lognormal distribution, which agree with the natural grains [55]. The layered porous medium was heterogeneous and had tortuous flow paths. In PM0, there are 1873 spherical grains in a two-size range constituting the coarse grain layer and fine grain layer. The average radius r (±δ, standard deviation) of coarse grains and fine grains was set to 1:0 ± 0:20 mm and 0:60 ± 0:12 mm, respectively. The coefficient of variation (COV) (COV = δ/ r) was 0.2, and the porosity for each layer was 0.44 [56,57].
Based on PM0, three layered porous media with different fractions of the grain size transition zone were generated and referred to as PM1-3. The generation process is shown in Figure 2 and described as follows: the coarse grains and fine grains in the middle zone of PM0 were selected. The middle zone corresponds to a range of 1 mm on both sides of the "zero thickness" middle plane. The selected grains were mixed and relocated through a random packing procedure. This procedure, reflecting the mixed sedimentation process of two-size range solid grains, produced the transition zone. As a result, the layered porous medium PM1 with the same domain size as PM0 was obtained. It should be noted that except for the transition zone, there was no change in the two layers to eliminate the uncertainty caused by the different pore structures. The presence of the mixed zone of the coarse and fine grains between the layers is the only difference between PM0 and PM1.
To obtain a layered porous medium with a larger fraction of the transition zone, the transition zone of PM1 was duplicated once and reinserted into PM1 to produce PM2. The transition zone of PM1 was duplicated twice and reinserted into PM1 to create PM3. The larger fraction of the transition zone means that the mixed sedimentation process of two-size range solid grains had happened over a longer period. For PM1-3, the transition zone fraction varies from 5.0%, 9.1%, and 13.6%, respectively, across the layered porous media. It should be mentioned that the transition zone can be considered a layer rather than a "nonzero thickness" transition zone between the layers if the transition zone fraction is large enough. Several parameters were calculated to evaluate the variation of the pore space structure on the transition zone (see Table 1).
As shown in Table 1, the overall pore space structure of PM0-3 showed a slight change. The r of grains and the porosity for PM0-3 were maintained at about 0.72 mm and 0.44, respectively. In Table 1, the pore throat τ represents the small gap between every two grains and the pore size d is defined as the equivalent diameter of the pore space surrounded by nearby grains and their pore throats [56]. τ and d were introduced as two metrics to access the pore space structure, since these two parameters were found to be related to the velocity probability density functions (PDFs) at low velocities and induce anomalous solute transport [58]. We first calculated τ based on the size and center of the grains in the sections. Then, the small gaps between each grain were plotted and translated into a binary image to determine the pore size using a watershed algorithm. As shown in Figure 3, the dominant colors on the two grain size layers are significantly different and represent two types of pore sizes. Near the transition zone, the pore size suddenly contracts while the grain size changes. For the entire layered porous media, the τ and the d did not show much variation. Since the transition zone fraction in the porous media was relatively small, the pore space structure parameters calculated based on the entire layered porous media could not effectively reflect this local structure change.

Direct Numerical Model for Fluid Flow and Solute
Transport. In this study, the velocity field u was obtained by solving the continuity and Navier-Stokes equations. The fluid in the flow field is assumed to be incompressible and not affected by gravity.
where ρ is the density of the fluid (ML -3 ), u = ½u, v, w is the velocity vector (LT -1 ), p is the fluid pressure (ML -1 T -2 ), and μ is the fluid dynamic viscosity (ML -1 T -1 ). The standard water properties at 293.15 K were used, for example, ρ = 998:2 kg/m 3 and μ = 1:002 × 10 −3 Pa · s. The boundary conditions were set as follows: the left and right boundaries, which represented inlet and outlet boundaries, respectively, were set as the first-type boundary. The steady flow was driven from left to right by the given pressure gradient dp. The no-slip boundary condition was applied at the solid-liquid interfaces. There was no flux at the surface of the cylinder. The transient conservative solute transport in porous media was simulated by solving the advection-diffusion equation: where c is the solute concentration (ML -3 ), t is the time (T), and D m is the molecular diffusion coefficient (L 2 T -1 ). The velocity vector u is from the flow field based on the solution of equations (1) and (2). The same molecular diffusion coefficient was adopted for all cases and set to 1:0 × 10 −9 m 2 /s. This D m value relates to the dissolved organic contaminant (e.g., dichloroethene, toluene, MTBE ethylbenzene, and 3 Lithosphere xylene) relevant to groundwater quality problems [59]. A dimensionless parameter (i.e., pore volume, Pv) related to the time is introduced as follows: where Q is the flow rate (L 3 T -1 ) per length and V is the total volume of the pore space domain (L 3 ). In this study, solute was introduced by a pulse injection method, with a flux weight distribution along the inlet boundary. The whole space in the porous medium was initially saturated with pure water. A total of 0.05 Pv of a solute, with an initial concentration c 0 of 1 mol/m 3 , was injected. Flow and solute transport simulations were performed in pore-scale regions by imposing pressure gradients resulting in certain dimensionless Peclet number (Pe) defined by which is the ratio of the characteristic diffusion time ðτ D = ð2 rÞ 2 /D m Þ to the characteristic advection time (τ a = 2 * r/u av ),    The classic ADE, which is based on Fick's law, was used to fit the BTCs, and the analytical solution for the pulse injection condition is given by where the D f ,cde is the fitted longitudinal dispersion coeffi-  5 Lithosphere dispersivity (L), u cde is the average velocity (LT -1 ) in the longitudinal direction, A is the total area of the injection cross section that is perpendicular to the flow direction, m 0 is the total amount of solute contained in the pulse, and n is the effective porosity (−). Unless the flow velocity is extremely low (or the groundwater is nearly stagnant), molecular diffusion is usually neglected when studying solute transport in porous media. In this study, the Peclet numbers were set to 10 and 100, indicating that advection was dominant. The velocity was large enough to neglect the effects of molecular diffusion.
where C m and C im represent the solute concentration in the mobile and immobile domains, respectively; θ m and θ im are the water contents in the mobile and immobile domains (or porosities), respectively, and the sum of them is equal to the porosity of the porous medium; D f ,mim is the fitted dispersion coefficient; and α is the rate coefficient associated with mass transfer between the mobile and immobile domains. A partitioning coefficient β = θ m /ðθ m + θ im Þ is introduced as the mobile water fraction. In this study, to investigate the influence of the transition zone and its fraction on solute transport, the average velocity u m in flow direction and dispersion coefficient D f ,mim were obtained by fitting the BTCs from the direct simulations to the ADE model. For the MIM inverse model, three best-fit parameters D f ,mim , α, and β were calculated. And the solute dispersivity α L ðα L = D f ,mim /u m Þ was estimated for all cases. The initial estimate of parameters has a great influence on the accuracy of the inverse parameters. An accurate initial estimate can improve the convergence to a5 minimum error. The initial value of u was set to be the value calculated by the direct simulation. The STANMOD V2.08 based on the CXTFIT code was used for the inversion of both models. The R 2 value and the globe errors E i , for assessing the fitting goodness, were calculated for both models. The E i is given by where c fit,i is the fitted result and c s,i is the BTCs from numerical simulations and N is the total amount of the data in BTCs.

Results and Discussion
3.1. Flow Field and Concentration Field. The flow field and its heterogeneity highly control the solute transport process. Figure 4 shows the flow fields in PM0-3. There was little difference between the flow fields for the coarse grain layer and the fine grain layer. This is because the pore space structure of the coarse grain layer and the fine grain layer in PM0-3 was kept constant. As can be seen in the zoomed snapshots of the flow field, the color contrast in the transition zone of PM1-3 was smaller than that in the corresponding middle zone of PM0. This indicates that the transition zone changed the flow field and caused a relatively uniform velocity distribution.
In Figure 4(a), the high local flow velocity in PM0 was captured at narrow pore throats due to the preferential flow path parallel to the flow direction. This preferential flow path was formed by the strong contraction of the pore space between the layers (Figure 3(b)). However, the presence of the transition zone reduced this strong contraction of the pore space and increased the pore throat and pore size in PM1-3. The peak velocity in the transition zone in PM1 was 4:87 × 10 −5 m/s and was reduced by 21.6% compared to that in the middle zone of PM0.
The probability density functions (PDFs) of the flow velocity in the entire layered porous media and the transition zone are obtained as shown in Figure 5. From Figure 5(a), it can be seen that the PDFs in PM0-3 were highly consistent on the scale of the entire layered porous media. This indicates that the transition zone, even though it occupies 13.6% of the entire layered porous media, has little influence on the distribution of flow velocity at the scale of the entire layered porous media. However, for the transition zone, the PDFs in PM1-3 are significantly different from those in PM0 (see Figure 5(b)). The PDFs of high flow velocity in the transition zone of PM1-3 showed a decrease compared to those in the middle zone of PM0, indicating that the distribution range of flow velocity in PM1-3 was reduced by the presence of the transition zone. This is due to the fact that as the transition zone fraction increased, the strong contraction of the pore space between the layers was reduced. The preferential flow path with high velocity decreased.
The coefficient of flow velocity variation, a global measure of the velocity spatial variability [30,54], was calculated in PM0-3 to quantitatively evaluate the velocity heterogeneity in the flow fields: where the σ u is the standard variance of the flow velocity, u av is the average velocity magnitude, and V is the volume of the pore space occupied by fluid.
In Table 2, it can be seen that the CV u for the entire layered porous media in PM0-3 was almost the same, indicating that the transition zone has little influence on the velocity spatial variability at the scale of the entire layered porous media. This was consistent with the result of the   PDFs of velocity (see Figure 5(a)) for the entire layered porous media, where the distribution range of flow velocity was almost the same in PM0-3. However, the CV u for the middle zone in PM0 was different from that for the transition zone in PM1-3. The difference between CV u in the transition zones of PM0 and PM1 was 16.64%. Compared with PM1, the CV u in PM2 and PM3 was increased by 5.24% and 11.12%, respectively. As the transition zone fraction increased, the velocity spatial variability in the transition zone increased.
In addition, hydraulic conductivity was slightly increased in PM1-3 compared to PM0 due to the presence of the transition zone. This variation in hydraulic conductivity was trivial, since the difference in hydraulic conductivity in different layers can naturally be several orders of magnitude [60]. This also shows that the influence of the transition zone on the flow field cannot be reflected by the hydraulic conductivity. However, in this study, the properties of the pore space structure between PM0-3, such as porosity, pore size, and throat, were almost the same at the scale of the entire porous media domain (see Table 1). There was an actual change in the pore structure of the transition zone. This change has significant influence on the local flow fields. For solute transport, the changed local flow field could significantly affect the concentration distribution, which in turn affects the evolution of the solute transport process. As shown in Figures 6(a)-6(d), the front of the solute plume passing through the transition zone in PM1 was different from the process in PM0 because the transition zone changed the pore throat and pore size. In PM0, the fronts of solute plume spread mainly along one preferential flow path.
However, in PM1-3, the fronts of solute plume spread along multiple paths in the transition zone and had a broader size in the transverse direction. Due to the preferential flow path, the concentration gradient was more pronounced in PM0 than in PM1-3.

Breakthrough Curves and Spatial Moment.
In order to investigate the influence of the transition zone on the BTCs, the flux-averaged BTCs were calculated in PM0-3 at three different cross sections. Three cross sections were located at the coarse grain layer, the fine grain layer, and the outlet boundary, respectively. The results are summarized in Figure 7. The flux-averaged concentration was normalized as In Figure 7, the BTCs of PM1-3 showed differences from the BTCs of PM0. All measured BTCs exhibited the long-tail behavior and the deviation from the Gaussian profile. As the solute plume was diluted during the transport, the peak concentration of BTCs decreased at three distinct cross sections (see Figure 7(a)). From Figure 7(b), it can be seen that although the average velocities of the entire layered porous media were the same, the breakthrough time for the BTCs of the coarse grain layer in PM0-3 was slightly different. Since the pore structure of the coarse grain layer in PM0-3 was the same (see Table 1), this difference in the breakthrough time for the BTCs of the coarse grain layer in PM0-3 is caused by increasing the average velocity of the coarse grain layer. This result is in accordance with previous 8 Lithosphere results [61]. However, when the solute plume was transported across the transition zone, the breakthrough time of the fine grain layer in PM0-3 was inversely different from that of the coarse grain layer in PM0-3 (see Figures 7(b)  and 7(d)). This indicates that the transition zone can increase the solute residence time and delay the breakthrough time. The presence of the transition zone reduced the geometric mutation leading to a widening of the pore throat and an increase in pore size. As shown in Figure 4, the advection in the preferential flow path of PM0 was stronger than that in PM1-3, which is related to the fact that the solute plume spread rapidly in the preferential flow path. Compared with PM1-3, the front of the solute plume could travel a longer distance with the same Pv after passing through the middle zone in PM0. Consequently, an earlier first breakthrough time was observed in PM0. Another interesting finding was the influence of the transition zone on the peak concentration of BTCs. The peak concentration of BTCs at the outlet boundary was higher in PM1-3 than in PM0. For Pe = 10, the peak concentration in PM1-3 increased by 1.99%, 4.27%, and 3.89%, respectively. For Pe = 100, the peak concentration in PM1-3 increased by 1.94%, 3.79%, and 3.99%, respectively. These differences indicated that the maximum concentration of the solute center was greater in PM1-3 than in PM0 when it left the outlet boundary. As the transition zone fraction increased, the solute center yielded higher peak concentration. This means that the dilution process of the solute center was weakened due to the presence of the transition zone.
Spatial moments are commonly used to characterize the distribution of solute concentration. The zeroth moment represents the total mass in the volume: The first spatial moment M i,l ðtÞ represents the location of the plume centroid in the main flow direction, i.e., the longitudinal direction: The second central spatial moment M 2,l,l ðtÞ represents the mean square displacement from the plume centroid in the longitudinal direction: In fact, M 2,l,l ðtÞ is a measure of spreading and dispersion of the plume. The first and second spatial moments of the plume distributions for PM0-3 were determined and are shown in Figures 8 and 9, respectively.
In order to analyze the influence of the transition zone on the distribution of solute concentration, both the first and second spatial moments were normalized by the 10 Lithosphere length of layered porous media L and the zeroth moment, respectively. From Figure 8(a), it can be seen that the overall trend of the normalized first spatial moment in PM0-3 was almost the same for a given Pe. However, it was found that the normalized first spatial moment at Pv = 0:5 was larger in PM0 than in PM1-3 (see Figures 8(b) and 8(c)). This means that the solute plume is transported faster in PM0 than in PM1-3 due to the presence of the transition zone. In Figure 9, the normalized M 2,l,l showed the mass spread around its plume center in the main flow direction. The results showed that the influence of the transition zone on solute transport was dependent on the transition zone fraction. There was less effective longitudinal spread in PM1 than in PM0 for Pe = 10 and Pe = 100. As the transition zone fraction increased, M 2,l,l in PM3 kept rising and exceeded the value in PM0. The coefficient of determination R 2 and the global error E i were determined as shown in Table 3. From Table 3, it can be seen that all R 2 were more than 0.99 and all E i were trivial. These results demonstrated that the robustness of both MIM and ADE models performed well for fitting the solute transport across layered porous media. However, these results reflect that the errors from both MIM and ADE models are not sensitive to the transition zone. In other words, the errors from both MIM and ADE models are insufficient to capture the influence of the transition zone 11 Lithosphere on the solute transport, even though the transition zone has a significant influence on the characteristics of the BTCs (see Figure 7). Table 4 shows the estimated parameters from both ADE and MIM models. The ratio of u cde to u was around 1 for all cases, indicating that the u cde estimated from the ADE  For the MIM model, all β were around 0.5 because the pore structure of the entire layered porous media was the same. However, the α in PM0-3 was significantly different, indicating that the transition zone has a significant influence on the mass transfer between the mobile and immobile domains. As shown in Table 4, the α in PM1 increased by 163.7% for Pe = 10 compared to that in PM0. As the transition zone fraction increased, the α in PM2-3 decreased. Similarly, for Pe = 100, the α showed a significant increase in PM1 while the α decreased in PM2-3. On the one hand, the pore space structure in the transition zone caused the variation of α. When the pore throat and pore size widened and increased, respectively, the flow field became more homogeneous. There was a larger contact area between the mobile and immobile regions in the widened pore space, which led to the increase in α.
For both ADE and MIM models, there was a smaller α L in PM1-3 than in PM0. For the ADE model, the α L was mainly due to the heterogeneity of the flow field. The transition zone caused a smaller CV u of the flow field in PM1-3 than in PM0. However, for the MIM model, the α L not only was controlled by the mass transfer between the mobile and immobile domains but also is a monotonically increasing function of the residence times in immobile zones [62]. For the same Pe, the α in PM0 was small, indicating the weak mass transfer between the mobile and immobile domains.

Conclusion
In this study, the influence of the grain size transition zone on fluid flow and solute transport through layered porous media was investigated. An approach combining a deposition process and a random packing process was developed to generate 3D pore-scale layered porous media. The geometrical property parameters (e.g., pore throat and pore size) of the pore space structure in the generated 3D porescale layered porous media were evaluated. A direct numerical model for fluid flow and solute transport was solved to obtain the flow field and concentration field in 3D layered porous media. The coefficient of flow velocity variation and the PDFs of flow velocity were determined to characterize the velocity spatial variability and the velocity distribution, respectively. Both the ADE and MIM models were used inversely to fit the BTCs.
For the flow field, the CV u and the PDFs of the flow velocity at the entire layered porous media are insensitive to the presence of the transition zone. Our results showed that the transition zone, even though it occupies 13.6% of the entire layered porous media, has little influence on the distribution of flow velocity at the scale of the entire layered porous media. However, the presence of the transition zone reduces the strong contraction of the pore space between the coarse grain layer and the fine grain layer, resulting in an increase in the pore throat and pore size in the transition zone. The CV u and PDFs of flow velocity in the transition zone varied with increasing transition zone fraction. This variation of the flow field in the transition zone has little influence not only on the distribution of flow velocity and the velocity spatial variability but also on the hydraulic conductivity at the scale of the entire layered porous media. Moreover, it should be mentioned that this was the main reason why the transition zone was not treated as a new layer from the point of view of hydrogeology.
For the concentration field, the analysis of BTCs and spatial moments showed that the transition zone has a significant influence on the characteristics of solute transport. Although the transition zone did not affect the flow field at the entire layered porous media, the preferential flow path was changed at the local scale of the transition zone. This changed local flow field by the transition zone may increase the solute residence time and delay the breakthrough time. All measured BTCs exhibited the long-tail behavior and the deviation from the Gaussian profile. Both the ADE