Multiscale modelling of dual-porosity porous media; a computational pore-scale study for flow and solute transport

Many environmental and agricultural applications involve the transport of water and dissolved constituents through aggregated soil profiles, or porous media that are structured, fractured or macroporous in other ways. During the past several decades, various process-based macroscopic models have been used to simulate contaminant transport in such media. Many of these models consider advectivedispersive transport through relatively large inter-aggregate pore domains, while exchange with the smaller intra-aggregate pores is assumed to be controlled by diffusion. Exchange of solute between the two domains is often represented using a first-order mass transfer coefficient, which is commonly obtained by fitting to observed data. This study aims to understand and quantify the solute exchange term by applying a dual-porosity pore-scale network model to relatively large domains, and analysing the porescale results in terms of the classical dual-porosity (mobile-immobile) transport formulation. We examined the effects of key parameters (notably aggregate porosity and aggregate permeability) on the main dual-porosity model parameters, i.e., the mobile water fraction ( φm ) and the mass transfer coefficient ( α). Results were obtained for a wide range of aggregate porosities (between 0.082 and 0.700). The effect of aggregate permeability was explored by varying pore throat sizes within the aggregates. Solute breakthrough curves (BTCs) obtained with the pore-scale network model at several locations along the domain were analysed using analytical solutions of the dual-porosity model to obtain estimates of φm and α. An increase in aggregate porosity was found to decrease φm and increase α, leading to considerable tailing in the BTCs. Changes in the aggregate pore throat size affected the relative flow velocity between the intraand inter-aggregate domains. Higher flow velocities within the aggregates caused a change in the transport regime from diffusion dominated to more advection dominated. This change increased the exchange rate of solutes between the mobile and immobile domains, with a related increase in the value of the mass transfer coefficient and less tailing in the BTCs. © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Soil and groundwater pollution by a broad range of industrial and agricultural contaminants is an ever-increasing problem worldwide. One issue exacerbating effective management of the subsurface is the problem of preferential flow of surface applied chemicals such as fertilizers, pesticides, trace elements and pathogenic microorganisms. Much evidence exists that preferential flow through especially the vadose zone is contributing to surface and subsurface pollution problems (e.g., Flury et al., 1994;Abbaspour et al., 2001;Hendrickx and Flury, 2001;Allaire et al., 2009, Vogel et al., 2010Zhang et al., 2013; term is used, to lead to mobile-immobile (MIM) type dual-porosity models of the form ( Coats and Smith, 1964;van Genuchten and Wierenga, 1976 ): where the subscripts m and im refer to the mobile and immobile region, c is the average concentration, θ is the volumetric water content, φ m is the mobile water fraction defined as θ m / θ , D m is the dispersion coefficient, q is the volumetric fluid flux, and α is a first-order mass transfer coefficient.
The dual-porosity model given by Eqs. (1) and (2) involves several assumptions, the most important being that advective transport in the smaller intra-aggregate pores can be neglected. This implies that the overall pore water velocity distribution within the porous medium is approximated by the step function, with one part of the medium having an average pore water velocity equal to v m = q/ θ m , while water in the other part is completely stagnant.
This assumption is generally not met since the intra-aggregate (soil matrix) region often has some non-zero permeability, even if small compared to the inter-aggregate (or fracture) region.
Another simplifying assumption of Eqs. (1) and (2) is that solute exchange between the mobile and immobile regions can be described using a quasi-empirical first-order rate term proportional to the difference between the average concentrations of the mobile and immobile regions. Various attempts have been made to obtain a more physical basis of the mass transfer coefficient, α, in terms of such parameters as the diffusion coefficient and the shape and size of the aggregates or soil matrix. This has led to a number of analytical models that explicitly considered diffusion from the inter-aggregate region into immobile intra-aggregate regions of various shapes (e.g., Rasmuson and Neretnieks, 1980;Tang et al., 1981;Sudicky and Frind, 1982;van Genuchten et al., 1984;van Genuchten, 1985 ). These and related studies also allowed derivation of approximate relationship for the mass transfer coefficient, including through the use of Laplace transforms or moment analyses, as exemplified in studies by van Genuchten and Dalton (1986), Bolt (1979), Barker (1985), van Genuchten and Dalton (1985), Parker and Valocchi (1986), Goltz and Roberts (1987) , and Hantush and Marino (1988) .
Another approach for estimating the parameters φ m and α in Eqs. (1) and (2) is by direct measurement (e.g., Clothier et al., 1992;Jaynes and Shao, 1999 ), or by analyzing a large number of previously published data such as shown by Maraqa (2001) . A latter study revealed an approximately linear relationship between the mass transfer coefficient and the residence time of the solute in the transport domain ( Pontedeiro et al., 2010 ).
Both of the above assumptions (i.e., negligible advection within the aggregates, and the use of a first-order exchange term) require more research. This includes how best to account for aggregate shape and size, which are known to vary widely and may involve various spherical, blocky, columnar, and prismatic geometries, or mixtures thereof (e.g., Tisdall and Oades, 1982;Hillel, 2003 ). Imaging techniques can be of great value for determining the size and shapes of aggregates, and hence can reduce the dependency on using idealized shapes ( Bultreys et al., 2015( Bultreys et al., , 2016. The mass transfer coefficient is, in actuality, a more complex integrated parameter whose value depends on many porous media characteristics including pore and aggregate geometry, the solute diffusion coefficient, the intra-aggregate permeability, the relative magnitude of the mobile and immobile region, as well as the dynamics of the overall system such as the concentration gradients and the applied flow velocity. At relatively low pore water velocities, the time scale of solute diffusion into aggregates may be comparable with the transport time scale within the macropore domain, in which case a larger fraction of the solutes may diffuse into or out of the aggregates, thus limiting any tailing in observed BTCs. At the other extreme, at relatively high pore water velocities, macroand micro-pores may become essentially disconnected, leading to negligible solute exchange and possibility of dual peaks in the BTCs ( Zhou et al., 2014 ).
The value of the mass transfer coefficient, α, is commonly obtained by fitting macroscopic models to observed solute BTCs ( Toride et al., 1995 ). Since BTCs often show tailing, this approach may be very time-consuming in terms of getting appropriate resolution in the data. Moreover, the BTC data then provide information only of the macroscopic concentrations at selected observation points within the medium, or from column outflow experiments, thus providing little insight into the internal concentration distributions and interactions between the mobile and immobile zones. This makes it difficult to estimate contribution of different transport processes into the BTCs and extrapolate such BTC data to other transport regimes. Methods are hence needed to provide information in a systematic manner about the internal state of aggregated media, including quantification of the basic transport processes operating at the microscopic level. An alternative would be to obtain the BTCs using pore scale modelling . For transport in unsaturated non-aggregated media,  used pore network modelling to obtain several BTCs in this manner, which would be difficult and time consuming to obtain experimentally.
Advanced 3D X-ray microtomography and related imaging techniques are now being increasingly used to obtain non-destructive visualizations of pore structures ( Allaire et al., 2009;Dal Ferro et al., 2013;Mangalassery et al., 2013;Zhou et al., 2013;Martínez et al., 2015) . This includes studies of dynamic processes such as fluid flow and structural dynamic processes ( Cnudde and Boone, 2013;Wildenschild and Sheppard, 2013 ). Soil aggregates for such studies could be imaged and analysed individually to collect data on pore morphology and connectivity ( Zhou et al., 2013;Dal Ferro et al., 2013 ), or direct fluid flow experiments could be performed on a dual-porosity medium. Resulting information can then be utilized to construct pore network structures needed for pore-scale flow and/or transport models.
Various pore-scale modelling approaches are now being pursued, with differences pertaining to the specific mathematical formulation such as the use of direct numerical solutions ( Bijeljic et al., 2013;Fathi et al., 2017aFathi et al., , 2017, Lattice Boltzmann methods ( Jafari et al., 2011 ) or pore network modelling , including differences in the invoked modelling resolution and required computational time ( Sahimi, 2011;Bultreys et al., 2016 ). Multi-scale pore network modelling have been applied for simulating two-phase flow ( Bultreys et al., 2015;Jiang et al., 2013;Prodanovi ć et al., 2015 ) as well as solute transport ( Bijeljic et al., 2013 ). and also using grain filling and pore filling method ( Mehmani and Prodanovi ć, 2014 ).
The aim of our study was to use a pore-scale network model to simulate flow and transport in dual-porosity domains containing a large number of aggregates. The pore network modelling approach assumes that the porous medium continuum can be divided into pore elements made up of pore bodies representing the larger voids in the medium, and pore throats representing narrow openings connecting the pore bodies Bultreys et al., 2016 ). Applying mass balance equations, flow and transport is simulated within each individual pore. Averaging over a large number of pores will then allow estimation of the macroscopic transport properties for porous media containing unimodal pore sizes , or for multi-scale media with multimodal pore size distributions ( Bultreys et al., 2015;Mehmani and Prodanovi ć, 2014;Mehmani et al., 2015 ). The obtained information at the pore-scale makes it possible to relate macroscopic transport properties to the underlying physical pore-scale processes and pore size distributions, and their connectivities . Using such multi-scale porous media requires description of the connectivity between the macro-and micro-porosity domains, which computationally can be particularly challenging for heterogeneous multimodal pore scale domains.
Our objective is to develop a dual-porosity pore structure to represent an aggregated medium, and to use this structure to simulate flow and transport within the composite medium in order to obtain pore-scale distributions of the solute concentration. The resulting solute breakthrough curves are then compared with solutions of the macroscopic equations given by Eqs (1) and (2) to obtain estimates of the fraction of mobile water ( φ m ), and the mass transfer coefficient ( α). A large number of simulations will be carried out to obtain dependencies of φ m and α on such soil aggregate properties as porosity and internal permeability. Resulting insight may be useful for field scale models to improve predictions needed at the larger scales for evaluating alternative management or contaminant remediation strategies, as well as industrial porous media applications.

Theoretical development
In this section we describe the pore-scale network model that was used to simulate flow and transport in the aggregated dualporosity pore structure. After a brief review of the pore-scale network modelling approach in general, details are provided on how the pore-scale model was used to generate pore structures and conduct simulations for an aggregated medium.

Pore network model
Pore network modelling considers a porous medium as a system of pore elements composed of pore bodies and pore throats. While pore bodies, due to their larger size, mainly control the porosity of the medium, the smaller pore throats determine the hydraulic conductivity. The complete pore structure of an aggregated soil requires the use of two types of pore spaces with different sizes. The aggregated medium in our study was obtained by first generating a macropore domain. Since pore connectivity is an important topological property of the medium, we used a pore network model with random pore connectivity ( Raoof & Hassanizadeh, 2012 ) to represent the randomness of real porous media. Once the macropore domain was generated, several aggregates were superimposed within the overall domain. Each aggregate of certain size was assumed to contain a large number of micropores. Algorithms were developed to connect the collection of pores located at the outer surface of the aggregates to their adjacent macropores. We generated for this purpose several domains with different aggregate fractions, with the aggregates themselves randomly placed within the medium. An example is shown in Fig. 1 . Once the pore structures were constructed, flow and transport processes were simulated within the pore network domain using the PoreFlow package developed by  , which is capable of simulating saturated and variably-saturated flow, as well as multi-component reactive transport, within arbitrary pore structures.
Since the adopted pore-scale model has been explained in detail elsewhere Hassanizadeh, 2010, 2012;, only a brief overview is provided here of the main equations governing flow and transport within the porous medium. Assuming laminar flow, the volumetric flow rate within each pore was obtained using the Hagen-Poiseulle equation: where q ij is the volumetric volume rate through the pore throat between two adjacent connected pore bodies i and j, p i and p j are pressures of the two adjacent pore bodies, and g ij is the conductance of the pore throat which, for a pore throat with a cylindrical cross section, is given by: where R ij is the pore throat radius, l ij is the throat length, and μ is the fluid dynamic viscosity. Considering incompressible flow, the continuity equation for flow is applied at pore junctions: where q ij is the volumetric flow rate within a pore throat from pore body i to pore body j , and z i is the pore coordination number of pore i . Combining Eqs.
(3) -(5) for all pores results in a linear system of equations with a sparse, symmetric and positive-definite coefficient matrix, which is to be solved for the pore body pressures. Considering the sample as a representative elementary volume (REV), the average pore water velocity, v , can be calculated as :   Khan et al. (2012) 1-5 1-85 0.40 Synchrotron-based nano- CT Peth et al. (2008) 3-5 5-50 0.13 X-ray microtomography Rabbi et al. (2015) 2-4 5.2-60 0.048 X-ray microtomography Zhou et al. (2013) 1.5-2 10-500 0.14 X-ray microtomography where Q tot is the total discharge rate through the pore network, which can be determined at the inlet or outlet of the pore network as the sum of all fluxes, L is the length of the pore network, and V f is the total volume of the fluid phase within the pore network.
The intrinsic permeability, κ, of the sample can be determined using Darcy's law: where μ is the dynamic viscosity, P is the pressure difference between the inlet and outlet pores, and A is the cross-sectional area of the pore network. Solute transport through the pore network was described for the general case involving both advective transport and diffusion ( Vasilyev et al., 2012 ). Calculations were done by considering each pore element (i.e., a pore body or a pore throat) as the control volume. We used a backward Euler method for the temporal discretization and first-order upwind and central schemes for spatial discretization of the advection and diffusion terms, respectively. For a given pore body i , one can write the following mass balance equation: where V i is the volume of the pore body, Q i is the total volumetric rate going out of the pore body, A ij is the cross-sectional area of the pore throat, D 0 is the ionic or molecular diffusion coefficient, c i is the concentration in the pore body, c ij is the concentration in the pore throat, l ij is again the length of the throat, while N throat in represents the number of pore throats flowing into pore body i . Similarly, the mass balance equation for solutes in a pore throat ij may be written as: which assumes that pore body j is the upstream node. Eqs. (8) and (10) were solved using a fully implicit scheme available in the PoreFlow software package . Flux-averaged breakthrough curves at selected points were obtained by averaging concentrations over the network cross-section at successive times as follows: where c ( x,t ) is the normalized average concentration at location x and time t , N x t refers to the total number of pores at location x , and c 0 is the input concentration of the solute.

Application to aggregated media
Imaging techniques such as X-ray microtomography are being increasingly used to visualize the inner structures of aggregates and to provide information on pore size distribution, aggregate porosity, and aggregate size and shape ( Mangalassery et al., 2013;Martínez et al., 2015;Zhou et al., 2013;Dal Ferro et al., 2013 ). We note that studies using only X-ray microtomography often show lower aggregate porosities as compared to studies using multiple measuring techniques ( Czachor et al., 2015;Dal Ferro et al., 2012. The resolution of X-ray tomography has been shown to be a limiting factor since it tends to neglect the presence of belowresolution pores. Additional techniques such as mercury intrusion hence are often used to obtain a more realistic value for the porosity of the aggregates ( Dal Ferro et al., 2013 ). Focused ion beamscanning electron microscopy (FIB-SEM) can provide much higher resolution images, as compared to X-ray tomography, for analysing microporosity ( Hemes et al., 2015 ). Using a combination of multiple techniques applicable to different scales will be most valuable for multi-scale pore topology applications such as in this study. Table 1 shows aggregate sizes, pore sizes and aggregate porosities as obtained from the literature. We used these data to establish as realistic ranges of properties as possible for the aggregates used in our calculations.
Transport processes within the aggregated media we considered were studied by first creating a reference model. The reference model was then modified to consider i) domains with different aggregate porosities by changing pore body sizes, and ii) domains with different aggregate permeabilities by changing pore throat sizes to explore the effects of flow velocity differences between the macro-and aggregate domains. The reference model had a physical size of 100 × 20 × 20 mm 3 and contained a total of 10 0 0 aggregates. The mean distance between pore bodies in the aggregates was taken to be 0.1 of the separation distance between the macropores, while the mean coordination number was taken to be 5. The minimum and maximum aggregate diameters were set to 1.0 and 3.0 mm, respectively. Pore body sizes and aggregate throat sizes were assigned from a truncated log-normal distribution given by ( Raoof & Hassanizadeh, 2012 ).
where R min , R max and R m are the minimum, maximum and mean of the distribution, and σ 2 is the variance of the distribution. The macropore and micropore bodies each had their own size distribution. The macropore throat sizes were set as the smallest of the adjacent pore body sizes. The throat lengths were calculated by determining the length between the center points of two adjacent pore bodies of the throat, while subtracting the radii of the two pore bodies. Table 2 provides the pore size statistics used for the calculations, with Cases I-3 and II-2 serving as the reference models. Case I calculations were meant to show the effects of aggregate porosity, and Case II calculations the effects of aggregate permeability. Based on the data in Table 1 , the single aggregate porosity (Case I in Table 2 ) was varied between 0.082 and 0.700 in 6 steps by changing the aggregate pore body size distributions. The Table 2 Minimum ( R min ), maximum ( R max ) and mean ( R m ) pore body radii (Case I) and pore throat radii (Case II), as well as standard deviations ( σ ), of the pore-size distributions used in the calculations.

Porosity effects
Permeability effects aggregate throat radii as well as the macropore sizes were kept constant in these simulations to examine the effects of varying aggregate porosities on flow and transport. For the Case II simulations ( Table 2 ), the aggregate permeability was varied by changing the aggregate throat sizes in 4 steps, while the aggregate body and macropore body/throat sizes were kept constant. Simulation II-1 had smaller throats, and simulations II-3 and II-4 larger throats, compared to the reference model (Case II-2). The presence of macropores and aggregates with distinctly different micropore sizes created various pore systems. To explore the contribution of the different pore systems, porosities were calculated separately for the macropores and aggregates. Pore bodies were presumed to be spherically shaped, and pore throats to be cylindrical capillaries. Using the pore volumes of each domain, the domain porosities could be obtained, with the total porosity of the dual-porosity medium being simply the sum of the macropore and micropore domain porosities. For all simulations we kept the mean pore water velocity constant at a value of 1.50 × 10 −5 m/s. The molecular diffusion coefficient, D 0 , was set at 1.6 × 10 −9 m 2 /s, the viscosity at 0.001 Pa s, and the fluid density at 1 g/cm 3 .
For all simulations we injected into the pore network a pulse of one pore volume (equivalent to 66 min) of a tracer having a relative concentration of 1.0. BTCs of the average concentrations as obtained with the pore network model were analysed using the CXTFIT program ( Toride et al., 1995 ) within the STANMOD software of Šimunek et al. (20 0 0) , leading to estimates of the solute dispersivity, λ = D m/ v , the fraction of mobile water, φ m , and the mass transfer coefficient, α, in Eqs. (1) and (2) . This in turn also produced values of the mobile water content, θ m , and the immobile water content, θ im . The intrinsic permeability was determined for the domain as a whole, as well as for the macropore domain separately (i.e., excluding the presence of the aggregates). In addition, simulations were performed to obtain the intrinsic permeability of the single aggregates. For the Case II simulations, based on the simulated pore velocities, the characteristic time scales of diffusion and advection were calculated to explore the effects of the porewater velocity ratio between the macropore and aggregate domains on the calculations.

Results and discussion
Before presenting results of the various pore network calculations, we first summarize the pore network that was used for the reference simulations.

Reference model
The structured macropore network consisted of 12,929 pore throats connected to each other at 6850 junctions serving as the pore bodies, while the mean pore coordination number was 3.7. A total of 10 0 0 aggregates were randomly placed within the macropore network. Statistics of the generated pore network are given in Table 3 . The total number of pores within the created computational domain was 3,357,025. The macropore domain porosity was 0.198, and the mean porosity of the single aggregates 0.204. Fig. 2 shows the generated pore body and pore throat size distributions of the reference model, while Fig. 3 shows the generated dual-porosity pore network.

Case I: porosity effects
Based on the information in Table 1 , six different pore body size distributions were generated as displayed in Fig. 4 . Table 4 summarizes the porosities of the micropore and macropore domains for the Case I simulations. The porosity of the single aggregates ranged from very small (0.082) to very large (0.700). The aggregate domain porosity ranged from 0.013 to 0.114, and the total domain  Table 4 . porosity (micropore and macropore domains combined) from 0.214 to 0.315 (the porosity of the macropore domain was kept constant at 0.198). Fig. 5 shows the simulated solute concentration distributions for Cases I-1, I-3, and I-6 after injecting 1.0 pore volume of solute into the medium. The distributions are given at T = 1.1 pore volume, with pore volume defined as T = q t /( φ t L ), where L is the length of the medium (100 mm). Concentration distributions at other times are provided in the Supplementary Material, as well a video clip showing the effect of aggregate porosity. A pulse injection was used to show both solute diffusion into the aggregate domain and subsequent back diffusion from the aggregate domain into the macropore domain. The plots in Fig. 5 show increased solute spreading as the porosity of single aggregates increased. Case I-6 with its very high single aggregate porosity ( φ sa = 0.700) was included as a mere limiting case. The other extreme would be if no aggregate porosity were to be present ( φ sa = 0), in which case no intra-aggregate diffusion would take place and the standard advection-dispersion equation presumably would apply. Fig. 6 shows calculated BTCs (indicated by symbols) obtained with the pore network model at 30, 60 and 90 mm distances from the inlet for all Case I scenarios. Also included are the fitted curves obtained with the macroscopic mobile-immobile (MIM) dual-porosity model given by Eqs. (1) and (2) . We optimized for this purpose the parameters ( λ m , φ m and α) simultaneously to the three curves (i.e., at x = 30, 60 and 90 mm) obtained for each domain with a particular porosity of single aggregates. R 2 values for all Case I calculations, as well for the various Case II examples to be discussed later, were always higher than 0.9995, thus reflecting excellent descriptions of the pore-scale modelling results using the macroscopic formulation. Fig. 7 shows the fitted BTCs at 90 mm for all six scenarios. Similar BTCs at 30 and 60 mm are provided in the Supplementary Materials.
The results in Fig. 6 indicate increased spreading of the BTCs obtained at larger distances. Also, while the BTCs of Case I-1 were nearly symmetrical, Case I-6 with its highest porosity of single aggregates produced far more asymmetry and tailing in the distributions. Higher porosities of single aggregates (and hence more relatively immobile water) clearly enhanced the extent of preferential flow as exhibited by early breakthrough and late-time tailing of the curves in Fig. 7 . This agrees with previous study done on highly heterogenous rocks by Mehmani et al. (2015) who showed that an increasing fraction of microporosity causes more tailing of the observed solute plume. Our findings are also consistent with several other previous studies (e.g., van Genuchten and Wierenga, 1976;Brusseau, 1993;Shukla, 2013 ). Table 5 provides the transport parameters obtained for the six Case I calculations. Results indicate that a higher aggregate porosity caused a decrease in the mobile water fraction, φ m , and an Table 4 Micro-and macro-porosities for the Case I calculations with variable aggregate porosities. increase in the mass transfer coefficient, α. The mobile water content, θ m , was found to be approximately 0.21 for all Case I simulations, which is slightly higher than the macroporosity of 0.198. For all of the simulations, the mobile water content was slightly larger than the macro domain porosity. This is because the calculated fraction of mobile water ( φ m ) at the macroscopic level is determined by the flow lines affected by the presence of both macropores and aggregates. Fig. 8 shows that a clear relationship exists between the mobile water content, φ m , and the fraction of macropores. As compared with the porosity of macropores, the fraction of mobile water content using the fitting resulted in 5-10% larger values. Fig. 9 shows the relationship we found between the mass transfer coefficient, α, and the porosity of single aggregates, φ sa . The mass transfer coefficient increased with aggregate porosity as more solute mass could be transported into aggregates having higher porosity values and higher permeabilities as discussed next.
To further analyse and interpret the pore scale modelling results, we calculated the permeability ( K ) of the macro domain and the aggregate domain, as well as the permeability ratio between the aggregates and the macropore region ( K agg /K macro ). The perme-ability of the macropore domain was found to be 2.24 × 10 −4 mm 2 , with the results for K agg and K ratio given in Table 6 . When the porosity of the single aggregates increased, pores come to be closer to each other, which resulted in higher aggregate permeability values and larger permeability ratios. With a larger permeability, more fluid may flow into and through the aggregates, leading to higher α and lower φ m values ( Fig. 10 ).

Case II: aggregate permeability effects
Utilizing the input parameters from Table 1 , pore body size distributions were generated for the various examples. We used domains with four different pore-throat size distributions for the aggregates ( Fig. 11 ). Flow and solute transport simulations were carried out for each pore size distribution. Note that the reference model is Case II-2, which corresponds to Case I-3 of the first set of simulations. Table 7 provides the porosity values for the Case II simulations. The data show that an increase in the pore throat radii increased the single aggregate porosity as well as the total porosity.
The progress of one pore volume of solute moving through and dispersing into the domains of Cases II-1, II-2, and II-4 at 1.1 pore  Table 5 Optimized MIM transport parameters for the Case I simulations.
Porosity effects domain shows that solutes moved much faster into and through the macropore domain for Case II-4, while the residence time of solutes in the aggregates was higher for Case II-1. These results indicate that the exchange rate of the solute tracer between the aggregate and macropore domains for Case II-1 was lower than for Cases II-2 and II-4. A higher solute exchange rate was observed especially for Case II-4 due to the fact that solutes in the aggregates were almost instantly leached back into the macropore domain (also see Supplementary Material). = 5.52 10 -5 sa + 9.59 10 -6 R 2 = 0.996 Fig. 9. Plot of the mass transfer coefficient ( α) against the porosity of single aggregates ( φsa ) for Case I.

Table 6
Calculated permeability values for the aggregates, and the permeability ratio ( K agg /K macro ) for Case I.
show that reduced pore throat sizes lead to earlier breakthrough of solute, increased attenuation of the peak concentration, and more positive skewness of the BTCs in the form of tailing (e.g., Case II-1). Increasing pore throat sizes, on the other hand, resulted in delayed breakthrough curves with less tailing and more symmetrical BTCs. This is consistent with the visualizations in Fig. 12 and the Supplementary Materials, which show a more uniform solute front for Case II-4 as compared to Case II-1. The latter case exhibits more preferential flow and more pronounced tailing.
Utilizing the BTCs, transport parameters were again obtained for each simulation using the CXTFIT program. The results are given in Table 8 . The obtained mobile fraction, φ m , was found to be higher than the volume fraction of macropores for Cases II-1,  Table 7 . II-2, and II-3, while φ m was lower than the volume fraction of macropores for Cases II-4. These results are due to the hydrodynamic effects caused by the presence of aggregates in the domain. Smaller pore throats create relatively impermeable aggregates, which causes the overall fluid flow path through the macropores to become more tortious, leading to higher immobile fractions at the macroscopic scale. On the other hand, the presence of aggregates with larger pore throats causes a less tortuous flow path as aggregates provide a noticeable contribution to the total flow, leading to lower immobile water fractions at the macroscopic scale. The mass transfer rate, α, was lowest for Case II-1 and highest for Case II-4. This is in agreement with our earlier observations about the visualizations. Table 9 shows calculated pore water velocities of the macropore and micropore (aggregate) domains. The data indicate that larger pore throat sizes lead to higher pore water velocities. Fig. 15 shows plots φ m and α against the velocity ratio, v agg/ v ma . The mass transfer rate, α, decreased with a lower velocity ratio. In contrast with Case I, an increase in α did not lead to longer tailing but produced more symmetric BTCs. Increasing the aggregate pore water velocities caused the dominant transport process in the micropore domains to change from diffusion dominated to advection dominated. Interactions between the aggregate and macropore domains consequently increased, leading to higher values of the mass transfer coefficient, α. We note here that the two-region mobile-immobile model (MIM) considers the immobile domain to be completely Table 7 Micro-and macro-porosities for the Case II calculations with variable aggregate permeability. When solute-free water is injected and the transport regime in the aggregate domain is diffusion dominated, then the characteristic time for solutes to diffuse back into the macropore domain is given by t dif = l 2 / D m , where l is the effective length of the aggregates and D m the molecular diffusion coefficient. For an advection dominated transport regime, on the other hand, the characteristic transport time is described by t adv = l / v , where l is now the domain length. These time scales are characteristic for the extent of BTC tailing when the amount of mass stored in the aggregates is noticeable, and when the velocity ratio is low. The advection and diffusion time scales for the aggregate and macropore domains are listed in Table 10 . For these calculations we used the mean fluid velocity in the pore bodies and throats, while for the characteristic length of the aggregates we used the mean aggregate radius.
We used the transport time scales to identify the dominant transport mechanism in the aggregate domain. When t adv,agg < t dif,agg , transport in the aggregate domain is primarily advection dominated, but when t adv,agg > t dif,agg the aggregated domain becomes more dominated by diffusion. The data in Table 10 indicate that Case II-1 was very much diffusion dominated, while the experiments towards Case II-4 became increasingly more advection dominated. Tailing will be observed when the characteristic time  for diffusive transport in the aggregates is similar or becomes less in comparison to the advection time scale in the macropore domain. For Case II-1, when the transport domain is leached, the amount of solute tracer remains relatively high in the micropore domain since the time of solute tracer to be transported back into the macropore domain is much lower than the front of the solute-free water passing through the domain. This effect becomes less visible for Case II-2, and is almost absent for Case II-4. Extensive tailing is observed when the transport regime in the aggregates is diffusion dominated and the transport regime in the macro domain is advection dominated. For flow conditions when the velocity ratio between the two domains is close to 1, fluid passes through the aggregates as well as through the macropore domain, leading to relatively symmetrical BTCs. Fluid flow at lower velocity ratios is relatively much slower in the aggregates, leading to long tails in the BTC. These findings agree with other studies, such as those by Brusseau (1993) and Bijeljic et al. (2013) , who found that the degree of non-equilibrium and dominance of the velocity variations between the inter-aggregate and intra-aggregate pores influence the severity of tailing.
Values of the calculated permeabilities are listed in Table 11 . With increasing pore throat sizes, the permeability increased, with the permeability ratio increasing from a value of 7.24 × 10 −5 for Case II-1 to 0.56 for Case II-4. This means that for Case II-4 the fluid penetrated into aggregates almost as easy as it passed through the macropores. For Case II-1 the ratio is very low, meaning that fluid flow now takes place much less through the aggregate pores as compared to the macropores. This is consistent with the fact that considerable tailing was observed in the BTCs of Case II-1, while Case II-4 produced more symmetrical BTCs.

Table 8
Optimized MIM transport parameters for the Case II simulations.  Fig. 15. Values of the mobile water fraction ( φm ) and the mass transfer coefficient ( α) plotted against the velocity ratio ( v agg/ v ma ) for Case II.

Table 10
Characteristic advection and diffusion times (in minutes) for the Case II simulations.

Table 11
Calculated permeability values for the aggregates, and the permeability ratio ( K agg /K macro ) for Case II.

Conclusions
We investigated the effects of multi-scale pore sizes in a dualporosity porous medium on fluid flow and solute transport processes. Several dual-porosity pore networks were constructed by placing a large number of aggregates into the macropore domain. After construction of the aggregate domains, fluid flow and solute transport were simulated using the PoreFlow  pore network model to obtain three-dimensional flow fields and solute distributions, as well as the BTCs at several locations along the flow path. The resulting BTCs were fitted analysed in terms of the conventional macroscopic dual-porosity (MIM) transport model to estimate several key macroscopic transport parameters. Use of pore scale modelling provided information to link the transport properties to the underlying pore scale processes, which would be very difficult and time-consuming to do using only experimental methods.
A pore network model consisting of 10 0 0 aggregates, each having an average of 966 intra-aggregate pores, was generated to produce a domain with total of 3,357,025 pores. By varying the pore sizes, different aggregates with porosity values between 0.082 and 0.700 were obtained. We showed that increasing aggregate porosities did lead to higher values of the mass transfer coefficient and more tailing in the BTCs. Higher fluid velocities within the aggregate domain also increased the mass transfer coefficient, but produced less tailing in the calculated BTCs. On the other hand, lower pore water velocities in the micropore domain produced lower mass transfer coefficients but more extensive tailing in the BTCs. The velocity ratio between the aggregated and macro domain could explain the magnitude of tailing observed in the BTCs.
We showed that dual-porosity pore network models are attractive tools for investigating the basic processes governing solute interactions between the macropore domain (containing interaggregate pores) and the micropore or aggregate domain (containing intra-aggregate pores). The developed pore network model provides a flexible means of analysing the effects of pore scale transport properties such as aggregate porosity and permeability on macroscopic flow and transport parameters of dual-porosity porous media.
While in this study we focused on the effects of aggregate porosity and permeability on macroscopic transport, further studies are needed to better understand all of the pore scale properties of aggregates. In our work we varied the aggregate porosity by varying the size of the pore bodies. The aggregate porosity can vary also as a result of changes in the number of pore bodies, the pore coordination number, and the size of the throats. The influence of the range, variance, skewness, and curtosis of the pore size distribution on the mass transfer coefficients must also be investigated further. We note also that aggregate domain porosity can be influenced by the size and number of aggregates, which have been kept constant in this study. We furthermore did not consider any direct connections between the aggregates. These connections will have an influence on the residence time of the solute in the aggregated domain. They very likely will influence the magnitude of transferred mass, and how solutes are exchanged between the micropore and macropore domains. The use of suitable imaging techniques remains indispensable in all of this to most accurately obtain the exact micropore structures of the aggregates, and the porous medium as a whole.