Analytical evaluation of steady-state solute distribution in throughdiffusion and membrane behavior test under non-perfectly flushing boundary conditions

Guannian Chen ,Yuchao Li ,Kristin M.Sample-Lor ,Shan Tong,,*

a School of Civil &Environmental Engineering and Geography Science,Ningbo University,Ningbo,315211,China

b MOE Key Laboratory of Soft Soils and Geoenvironmental Engineering,College of Civil Engineering and Architecture,Zhejiang University,Hangzhou,310058,China

c Ningbo ZCONE High-tech Company Limited,Ningbo,315000,China

d Department of Civil and Environmental Engineering,Villanova University,Villanova,PA,USA

Keywords: Diffusion testing Membrane behavior Coupled transport Clay barrier Transport modeling

ABSTRACT The through-diffusion and membrane behavior testing procedure using a closed-system apparatus has been widely used for concurrent measurement of diffusion and membrane efficiency coefficients of lowpermeability clay-based barrier materials.However,the common assumption of perfectly flushing conditions at the specimen boundaries could induce errors in analyses of the diffusion coefficients and membrane efficiencies.In this study,an innovative pseudo three-dimensional (3D) analytical method was proposed to evaluate solute distribution along the boundary surfaces of the soil-porous disks system,considering the non-perfectly flushing conditions.The results were consistent with numerical models under two scenarios considering different inflow/outflow positions.The proposed model has been demonstrated to be an accurate and reliable method to estimate solute distributions along the boundaries.The calculated membrane efficiency coefficient and diffusion coefficient based on the proposed analytical method are more accurate,resulting in up to 50% less relative error than the traditional approach that adopts the arithmetic mean value of the influent and effluent concentrations.The retardation factor of the clay specimen also can be calculated with a revised cumulative mass approach.Finally,the simulated transient solute transport matched with experimental data from a multi-stage through-diffusion and membrane behavior test,validating the accuracy of the proposed method.

1.Introduction

For low-permeability clay barrier materials that are typically used in containment systems,e.g.bentonite with hydraulic conductivity less than 5×10-11m/s,solute diffusion and the existence of membrane behaviors in the clay have an important impact on the long-term barrier performance (e.g.Rowe,2012;Malusis et al.,2020).The through-diffusion and membrane behavior test,among other testing methods,is widely accepted as one of the most accurate and reliable approaches that allow measurement of the effective diffusion coefficient(D*)and retardation factor(Rd)of the solute,as well as membrane efficiency (ω) of the specimen (e.g.Fritz,1986;Shackelford et al.,2016;Tang et al.,2014,2015;Musso et al.,2017;Fritz et al.,2020).

The through-diffusion and membrane behavior test usually is conducted in a closed-system apparatus,in which a cylindrical soil specimen is enclosed in a rigid-wall cell between two pieces of porous plastic disks that connect to two reservoirs containing different solute concentrations(i.e.with perfect flushing condition)to establish a chemical gradient across the soil (e.g.Malusis et al.,2001;Bohnhoff and Shackelford,2015;Meier and Shackelford,2017;Sample-Lord and Shackelford,2017,2018;Dominijanni et al.,2018;Fu et al.,2021;Tong and Sample-Lord,2022).During the test,the testing solutions (e.g.a chemical solution with two different concentrations) flush across the porous plastic disks at a constant flow rate that is precisely controlled by a closed flowpump system.Using this setup,the testing solutions only flush along the boundary with advective flow through the specimen prohibited (Malusis et al.,2001).The through-diffusion and membrane behavior test typically is performed until both steady-state conditions for diffusion coefficients and membrane efficiency are achieved.The cumulative flux of the solute that diffuses across the specimen from the high-concentration boundary to the lowconcentration boundary is monitored over time and used to calculate theD* (e.g.Shackelford,1991) and ω (e.g.Malusis et al.,2001).

In a through-diffusion and membrane behavior test,the top and bottom boundaries of the specimen are circulated with chemical solutions with different concentrations (Δc) that can generate a constant chemical gradient across the specimen.During the test,solute diffuses from the top boundary(i.e.the higher-concentration side) to the bottom boundary (the lower-concentration side),resulting in a decreased concentration along the top boundary(the top influent concentration,cot,reduces to the top effluent concentration,ct,eff),and an increased concentration along the bottom boundary (the bottom influent concentration,cob,increases to the bottom effluent concentration,cb,eff).Under a certain concentration gradient,the diffusion and adsorption properties of the solute in the soil specimen can be determined using the steady-state(timelag) method,whereby the cumulative diffusive mass flux of the solute(Q′t,mg/m2)from the bottom outflow is calculated based on measured solute concentration in the effluents and plotted as a function of cumulative elapsed time (t′) (Muurinen,1990;Shackelford,1991).D* andRdcan then be calculated by using a best-fit linear regression to the steady-state portion of the-t′curve(Shackelford and Lee,2003):

whereQis the diffusive mass flux of the solute as measured in the bottom reservoir;Δcis the concentration difference between the average concentrations of the two porous plastic disks;neis the effective porosity of the soil specimen;2HandAare the thickness and cross-sectional area of the soil specimen,respectively;andLis the inflow-to-outflow distance.

Membrane behavior of the bentonite is quantified by the membrane efficiency,ω (0 ＜ω ＜1),with unity representing an ideal semi-permeable membrane with selective restriction of a certain type of solute (e.g.anion exclusion in bentonite) and zero representing no selective restriction (Mitchell,1993).If the specimen exhibits membrane behavior,a differential pressure(ΔP)will develop due to the applied concentration gradient (Malusis et al.,2001).During the test,water pressures at the top and bottom boundaries are monitored to determine ΔP(=Pt-Pb),and the steady-state values of boundary water pressures are used to calculate ω (Malusis et al.,2001):

where Δπ is the theoretical maximum chemico-osmotic pressure difference across the specimen in accordance with the Van’t Hoff equation (Barbour and Fredlund,1989);ΔPeis the measured differential boundary water pressure across the specimen;ν is the number of ions per molecule of the solute (e.g.2 for potassium chloride(KCl));is the universal gas constant(8.314 J/mol K);andTis the absolute temperature (assumed to be 293 K for room temperature in this study).

In both diffusion and membrane behavior evaluations,an accurate measurement of Δcis essential.In previous studies,both midpoint and average solute concentrations are generally calculated by taking the arithmetic mean of that of the influent and effluent solutions at each boundary (ct_ave=(c0+ct,eff)/2;cb_ave=(cDIW+cb,eff)/2)(e.g.Shackelford,2013a,b;Malusis and Daniyarov,2016;Shackelford et al.,2016).However,errors may arise due to the actual solute concentration profiles within the porous plastic disks being unknown,resulting in inaccurate interpretations of ω andD*.

In this paper,a pseudo three-dimensional(3D)analytical model of steady-state solute flux is proposed in consideration of nonperfectly flushing boundary conditions by simulating cylindricalshaped apparatus.The relative errors for the proposed methods are investigated numerically in comparison with a previously established method.A transient simulation is also conducted to measure breakthrough curves of an example using the determined parameter values from the proposed method.

2.Theory

2.1.Geometry of 3D models

A 3D model is established to illustrate the solute distribution in a soil-porous disks system under steady-state conditions by a typical through-diffusion and membrane behavior testing apparatus.The geometry of the model comprises a cylindrical soil specimen sandwiched between two identical cylindrical porous disks with an identical diameter of 2R.As shown in Fig.1,the thickness of the soil specimen and the thin porous disks are 2Handh,respectively.

Fig.1.Planarized diagram of the model of a soil specimen and porous disks in a typical through-diffusion and membrane behavior test.

To validate the analytical model and evaluate the impact of the inflow and outflow locations on the solute distribution in the system,a 3D numerical model is established with the COMSOL Multiphysics®software.The origin of the Cartesian coordinate system is located at the centroid of the cylindrical soil specimen.In addition,two identical cylinders are added to the external surface of each porous disk that represents the inflow and outflow ports connecting the porous disk to the influent and effluent reservoirs.The diameter(rs)of each connecting port is set to be 1/100 to that of the soil specimen,to minimize the size impact (that is,3.5×10-4m),and a depth of 3.2 mm.The axes of all cylinders are parallel to thez-axis,and the ports that connect the porous disk to the influent and effluent reservoirs are set on thex-zsurface.

Two scenarios considering different inflow-to-outflow distances(L),denoted as M1 and M2,are shown in Fig.2.In the M1 scenario,the connecting ports are located at the two opposite ends of the surface of each porous disk (i.e.L=2R),representing an ideal scenario in which the testing solution flushes evenly through the entire surface area that results in a relatively lower concentration gradient across the porous disk.In the M2 scenario,the connecting ports are located at 1/3 and 2/3 of the diameter on the surface of each porous disk (i.e.L=2R/3),representing the actual testing apparatus used in most of the experimental studies (e.g.Malusis et al.,2001;Shackelford et al.,2016;Tong and Sample-Lord,2022).The geometrical parameters of the two scenarios are generally consistent with the specimen dimensions reported by Kang and Shackelford(2009,2010),as summarized in Table 1.

Table 1Parameters of the 3D numerical model with two different scenarios.

Fig.2.Geometries of 3D numerical model with (a) M1,and (b) M2 scenarios(dimension in m).

2.2.Flow in porous disks

To simulate the testing solution flushing across porous disks,flow circulation is established at the top and bottom surfaces of each porous disk.In this study,the influent at the top porous disk is a chemical solution with a source concentration ofc0,whereas that at the bottom porous disk is deionized water(DIW).The advection within the soil-porous disks system is expressed in the form of Darcy’s law (Javandel et al.,1984;Bear and Verruijt,1987):

whereuis the pore-water pressure;k is the permeability assumed to be 1 ×10-1m/s and 1 ×10-12m/s for the porous disk (e.g.a porous disk commonly used in the experimental studies;TO-6,GenPore,PA,USA) and soil specimen (e.g.a high-quality sodiumbentonite) adopted in the models,respectively;and vais the apparent velocity (i.e.Darcy velocity) of the flow through the porous disks.In the numerical model,the advective flow is driven by the pressure head of the influent,u0.The boundary conditions are set as follows:

where n is the normal vector to the corresponding surface and should be distinguished from porosity.Definitions of the boundaries can be found in Fig.1.

2.2.1.Flow distribution

In the analytical model,the advection is allowed to occur in the porous disks,whereas that in the soil specimen is inactivated due to the low permeability as well as representing the closed-system setup in the through-diffusion and membrane behavior tests.To simplify the simulation,the flow area within the porous disks is set to be infinite to avoid Bessel equations in the analytical solutions,assuming negligible impact from the geometrical boundary of the porous disk to the water pressure distribution.Thus,the distribution of water pressure can be calculated by

For each half of the porous disk,the seepage lines appear to be circles that center on they-axis (i.e.the centerline between the influent and effluent ports),according to Eq.(7):

where υ is they-axis coordinate of the center of each circular seepage line.Since any given point on thisx-yplane,except for the influent and effluent source,only belongs to one seepage line,the υ value is unique and can be determined by Eq.(8).The intersection of each seepage line and they-axis isy0.Note that they×y0will always be positive since the seepage lines will never cross thexaxis.The flow nets of the proposed analytical model and numerical models are shown in Fig.3.

Fig.3.Flow nets of the(a)Analytical(infinite)model,and numerical models of the(b)M1,and (c) M2 scenarios (dimension in m).

In the numerical model,the flow nets of the two scenarios are significantly different,as shown in Fig.3b and c.For example,the flow net of the M1 scenario is relatively uniform,similar to the flow nets in the central area of the M2 scenario and the analytical model(with infinite flow area).The outer flow nets of the analytical model and the M2 model are much wider and sparsely distributed,corresponding to the longer flow path and larger cross-sectional width,as shown in Fig.3a and c.However,the outer seepage lines of these two models are also significantly different,since the circular flow lines in the M2 model are compressed by the boundary of the porous disk,resulting in a higher flow density in the modeled area.

2.2.2.Effectiveflow rate

The centerline of the porous disk (i.e.they-axis) is regarded as the critical axis for the flow rate determination.Since no geometrical boundary is defined for the analytical model,the area within the range ofy0=[-R,R]is regarded as the effective flow area.Based on this definition,the effective flow area of the M2 model is much smaller than the actual surface area of the porous disk,as shown in Fig.3a.Then,the effective flow rate of the analytical model,q,which is the total flow rate passingy0=[-R,R]is equivalent to the effective influent and effluent fluxes due to mass balance:

The values ofqin the M1 and M2 scenarios based on the analytical model are 9.5×10-5u0m3/s and 1.91×10-4u0m3/s withk=0.1 m/s andh=0.0032 m,respectively,whereas that based on numerical models are 9.49 ×10-5u0m3/s and 2.27 ×10-4u0m3/s,respectively.For the M1 scenario,the simulated flow rate matches perfectly with that calculated in the analytical model,since the boundaries of the former are identical to the seepage lines of the latter.In this case,the assumption of the infinite flow area for the analytical model has no impact on the results.However,the calculated flow rates based on the analytical models are significantly different from that based on the numerical M2 model.For the M2 scenario,some of the outer seepage lines are compressed into the surface area of the porous disk,resulting in a relatively higher total flow rate.In this case,a correcting factor of 0.84 is needed to match the actual flow rate to the proposed analytical model in M2.

2.3.Steady-state solute distribution in the soil specimen and porous disks

The steady-state solute distribution in the soil specimen-porous disks system is established based on the steady-state advectiondispersion equations(Javandel et al.,1984;Bear and Verruijt,1987)and the continuity equation of the solute flux at the interfaces between the soil specimen and porous disks:

wherecis the solute concentration in pore water;andnpand Dpare the porosity and diffusion coefficient of the porous disk,respectively.

Eqs.(10) and (11) represent the solute distribution in each porous disk and the soil specimen,respectively.Eq.(12)represents the flux continuity between porous disks and the soil specimen.Dirichlet boundaries are adopted for the inflow boundaries of the porous disks where the testing solutions(i.e.salt solutions or DIW)are injected at constant concentrations(Sin,tfor the top porous disk andSin,bfor the bottom porous disk).Neumann boundaries are adopted for the outflow boundaries (Sout) with an open boundary with zero diffusive flux.The rest of the boundaries(Sim,pfor porous disks andSim,sfor the soil specimen) are set as no-flux (impermeable) and Neumann boundaries due to zero-advection.The boundary conditions are expressed as follows:

Since the solute distribution in the soil-porous disks system is antisymmetric,the sum of the solute concentration of any two mirrored points with respect to the central cross-sectional plane of the specimen can be considered as the source concentrationc0.Thus,only the upper half of the model (i.e.the upper porous disk and the upper half of the soil specimen) is analyzed herein.The central cross-sectional plane of the soil specimen is considered as the new boundary(i.e.z=0),on which the solute concentration is constantly equivalent to the average of the top and bottom influent concentrations:

To simplify the derivation of analytical solutions of the solute concentration profile in the soil specimen-porous disks system,the horizontal diffusion of the solute in the soil specimen is also ignored.A detailed discussion of the negligible impact of the horizontal diffusion is provided in the Supplement.The governing equation of solute distribution in the soil specimen is simplified as

The solute distribution in the porous disk is critical in the entire transport process through the soil-porous disk system:

wherect(x,y) is the function of the solute distribution in the top porous disk that needs to be solved.

Horizontal diffusion within porous disks is neglected due to the relatively high flow rate that is commonly used in the throughdiffusion and membrane behavior tests (i.e.21.6 mL/d).Considering the relatively small thinness of porous disks used in typical through-diffusion and membrane behavior tests,(i.e.3.2 mm,Kang and Shackelford 2009,2010),the solute transport in porous disks can be further reduced to be a two-dimensional problem(Zou et al.,2016).The testing solution is assumed to be sufficiently mixed along the vertical direction (z-axis) in porous disks so that no concentration gradient exists alongz-direction:

Using boundary conditions shown in Eqs.(13)-(15),Eq.(19)can be written as follows:

whereIXis the part of seepage lines that connect to the influent port,as shown in Fig.3a.The integral in Eq.(20)represents the time required for the advective flow to reach an objective point from the influent port.Combining water pressure distribution (Eq.(7)) and Darcy’s law (Eq.(3)),Eq.(20) can be estimated as

where τ is the dimensionless arriving time of the solute at an objective point,which can be written as

Note that φ in Eq.(22)is the angle from the positivey-axis to the characteristic radius of an objective point on the seepage lines,as illustrated in Fig.3a.By this definition,φ is negative when the characteristic radius rotates clockwise to the positivey-axis.φinis the φ value of the inflow point with the same υ to the objective point to maintain the same seepage line.Thus,the expressions of φ are proposed as

The value of φ ranges from-π to π,and may exceed the typical range of-π/2 to π/2 as the arctan function returns when 2Ris larger thanL(i.e.M2).Finally,the hypothetical source radiusrscan be eliminated with the substitution of Eq.(9) into Eq.(21):

The 3D concentration distribution in the soil-porous disk system can be obtained by combing Eq.(24) with Eq.(18).

3.Analysis and discussion

Fig.4 shows the distribution of steady-state solute concentration(ct)for both the M1 and M2 scenarios,with the upper half and lower half in each plot representing the results from the analytical and numerical models.Simplifications of the diffusion process were kept in the numerical model used in this section,and further validations using original 3D numerical models are provided in the Supplement (Figs.S1 and S2).In the simulations,the advection process is driven by a constant water head,and the values of the pressure head of the influent (u0) are set to be 2 × 10-6m and 2 × 10-7m to provide flow rates similar to that in the actual experimental tests;the value ofneD* is 1 × 10-10m2/s in the simulation.Thectvalue gradually decreases from the inflow port to the effluent port in all cases,while the decreasing rate is relatively faster in the lower flow rate scenarios (i.e.Fig.4b and d).

Fig.4.Distribution of solute concentration(ct)in the porous disks from the analytical and numerical models for(a)M1 with u0=2×10-6 m,(b)M1 with u0=2×10-7 m,(c)M2 with u0=2 ×10-6 m,and (d) M2 with u0=2 ×10-7 m.

As illustrated in Fig.4a and b,the results from the analytical and numerical models for the M1 scenario are identical,indicating that Eqs.(21)-(24) are adequately reliable for describing the solute concentration distribution on the boundary surface of the specimen in a through-diffusion and membrane behavior test when the flow net is accurate.Thus,no further modification is needed for the analytical model in the determination of effluent and mid-point concentrations if the testing setup is identical to the M1 scenario.

However,the results from the analytical and numerical models for the M2 scenario are significantly different,as shown in Fig.4c and d.The solute concentration distributions calculated based on the two models are similar in the inner circular area (with a diameter connecting inflow and outflow ports),indicating that the advection term calculated based on Eq.(7) is reliable in this area.However,in the outer area,taking the central line (i.e.y-axis) as a reference,the solute concentrations calculated by the analytical model are always higher,mainly due to the further underestimated length of seepage lines when the flow rate has already been underestimated.For the M2 scenario,a factor of 0.84 is needed to modify the calculation ofct,midfor the actual flow rate to recover the central flow rate when using Eq.(24):

The effluent concentration can further be estimated by averaging the total effluent flux.The intersection of the seepage lines and they-axis (y0) is used for the characterization of flux and terminal effluent concentration of the seepage lines.Since the coordinate of the effluent port is the same for each seepage line and cannot be used to distinguish the seepage lines,the dimensionless arriving time τ (Eq.(22)) of solute at the effluent port can be considered as twice that at they-axis,according to the symmetric characteristic of the seepage lines.A modification is also needed for the M2 scenario as the solute concentration in the outer area in the analytical model could have been overestimated.Combining Eq.(7)and Eq.(24),the effluent concentration can be expressed as follows:

Modification factors in Eqs.(26) and (28) need to be calibrated separately for test apparatus withL/Rvalues other than that of M1 or M2.Thect,midandct,effcalculated by the analytical model (Eqs.(25)-(28)) compared to that from the numerical model are shown in Fig.5.For the M1 scenario,the values ofct,midandct,efffrom the analytical model and the numerical model are identical,as shown in Fig.5a.However,ct,midcalculated based on the traditional approach (i.e.using the arithmetic mean of the influent and effluent concentrations at each boundary to represent thect,mid)may result in a relative error of a few percentages even when the flow rate is extremely high(e.g.1×10-9m3/s withh=0.0032 m).For the M2 scenario,ct,midis more insensitive to the decrease ofct,effdue to the relatively higher velocity ratio through the mid-point compared to the total flow rate.Thus,using the arithmetic mean of the influent and effluent concentration to calculate the diffusion coefficient may not be adequately accurate.The analytical model,on the contrary,fits the numerical results well for bothct,mid,andct,effwith proper flow rate modifications.

Fig.5.Comparison of ct,eff and ct,mid values of analytical and numerical models for (a)the M1,and (b) the M2 scenarios.

A dimensionless factor,qH/(neD*L2),representing the relative degree of the advection across the porous disk over the vertical diffusion through the soil specimen,is shown as the topx-axis in Fig.5.For a given scenario,qH/(neD*L2) is the merely parameter that influences the steady-state concentration,according to Eq.(24)and the nondimensionalized governing equations.Thus,Fig.5 can be used as design charts for determining effective diffusion coefficientD*based on measured flow rateqand effluent concentration at the top boundaryct.Fig.5 and the proposed analytical solutions are applicable for any closed-system diffusion and membrane behavior test that is applied to low-permeability,high radius-toheight ratio soil specimens,and thin porous disks.Note that the flow rate (q) values in Fig.5 are the measured values and need no modification when Fig.5a and b are used as design charts.

The relative errors of the effective diffusion coefficient calculated based on the traditional method(Eq.(1),the arithmetic mean of influent and effluent concentrations),and the proposed analytical solution (Eq.(27)) compared to the values used in numerical analyses,are shown in Fig.6.The analytical solution functions perfectly with relative error ranges between ±0.01% for the M1 scenario and-20 to 10% for the M2 scenario.On the contrary,relative errors of the traditional method are generally several percentages off and increase rapidly with the decrease ofct,eff,as the M2 scenario is relatively more difficult to achieve perfect flushing due to the low concentration zone in the outer area.Overall,the performance of both the traditional method and the proposed analytical method are acceptable for the M1 scenario(with relative errors between-10% and 10% for most conditions) with the proposed solution (Eq.(27)) performing excellently.However,the relative error of the two methods for the M2 scenario is significantly larger.The relative error of the traditional method exceeds-10% whenct,eff/c0is less than 0.92 and increases to nearly-75%,indicating that the calculatedD* may only be 1/4 of the theoretical value.For the proposed analytical method,the relative error is controlled within-20%-10%,which is within an allowable relative error of ±10% whenct,eff/c0is higher than 0.7,demonstrating a significant improvement in the analysis.

Fig.6.Comparison of relative error in calculated effective diffusion coefficients based on the traditional method (Eq.(1)) and the analytical solution (Eq.(27)).

In a through-diffusion and membrane behavior test,the retardation factorRdof the soil specimen can be determined using a cumulative mass approach similar to that of the column tests(Shackelford,1995;Shackelford and Hong,2020),corresponding to the scenario when the average pore concentration of the entire testing specimen (i.e.the soil specimen and porous disks) being 0.5c0at steady-state.A detailed estimation example will be presented in the subsequent section.

4.An application example and model verification

To validate the analytical model,an application example is conducted by applying the model to an experimental data set from a multi-staged through-diffusion and membrane behavior test conducted by Tong and Sample-Lord (2022).The soil specimen used in the experiment was a GCL-grade polymer-enhanced bentonite (i.e.BPC-3.2).The testing solution was potassium chloride (KCl) with concentration ranges of 5-400 mmol/L.Chloride(Cl-)is selected as the indicator solute.The connecting ports in the testing setup are located at approximately 1/3 and 2/3 of the specimen diameter,similar to that in the M2 scenario in the proposed model.The porosity (ne) and the thickness (2H) of the BPC-3.2 specimen were 0.9 and 1.8 cm,respectively.The flow circulation rate was 21.6 mL/d for 5-100 mmol/L stages and then increased to 39.5 mL/d for the 200 mmol/L and 400 mmol/L stages to compensate for the long testing duration.Due to mass balance,the solute mass introduced to the system(c0)will be distributed in three ways:(i)diffuse through the soil specimen and appear in the bottom effluent solution (cb,eff),(ii) retained in the top boundary and appear in the top effluent solution (ct,eff),and (iii) retained in the soil specimen.Thus,the sum of the steady-state values ofcb,effandct,effis considered as the actualc0for each stage.

The effective diffusion coefficient (D*) is determined by using Eqs.(27)-(28) via MATLAB® programming.Total cumulative mass with respect to the flow rate and influent/effluent concentrations are plotted with the elapsed time.Linear regressions of the steadystate portion of theQL-tare conducted,in whichQLis defined to be the absolute value of the intercept with they-axis that represents the retained solute mass per cross-sectional area of the soil specimen during each stage:

Then,the retardation factor (Rd) of the solute in the soil specimen can be expressed as

whereciis the influent concentration at the top boundary from the previous testing stage.

Values of diffusion coefficient,membrane efficiency,and retardation factor that estimated based on the experimental data from Tong and Sample-Lord(2022)using two different boundary solute concentration methods are presented in Table 2.The values ofD*estimated based on the proposed method (Eq.(27)) increase with increasing pore-water concentration in the specimen,generally consistent but slightly higher by 25%-40% than that reported by Tong and Sample-Lord(2022).The values of ω calculated based on the proposed method are approximately 15%-35% lower than that based on the arithmetic mean method,and the difference between the two methods decreases with increasing pore-water concentration.The values ofRdgenerally range from 2 to 3.5 for all stages except for the 5 mmol/L stage,indicating a relatively small adsorption behavior of chloride onto the bentonite that is consistent with results reported by Shackelford and Redmond(1995)and Chen et al.(2020).For the 5 mM stage,theRdvalue is less than 1,which is possibly impacted by the remaining Cl-from the flushing period.

To verify the accuracy of the proposed method,a transient model of solute transport in each stage is also simulated based on the M2 scenario,considering the full thickness of the soil specimen and all directions of solute diffusion,with all of the conditional parameters adopted from the steady-state portion of the testing data using the proposed method (Eqs.(27) and (30)),as shown in Table 2.In this case,transient solute transport is achieved by replacing the right-hand-side term in Eqs.(10)-(11),zero,byRd∂c/∂t.The seepage field in the simulation is assumed to be steadystate,and the boundary pressure is increased at the 230th day(i.e.within stage 250 mmol/L) to better simulate the realistic transport process.Values ofD* from Table 2 (i.e.Eq.(27)) are adopted for each stage.Rdis set to be 1 forc0=5 mmol/L and 2.5(equivalent to the average value of theRdfor all the other stages)for the rest of the stages.

The analytical results based on parameters calculated using the proposed method (Eqs.(27) and (30)) are compared with the experimental data(after Tong and Sample-Lord,2022)and plotted in Fig.7.The overall correlation indices(R2)between the simulated breakthrough curves and the measured data are 0.991 and 0.813(0.994 if neglecting the abnormal 4800 mg/Lcb,effdata at 214th day)for top and bottom effluent curves,respectively.The high correlation indices illustrate that the simulated breakthrough curves match well with the measured data,especially for the last two stages in which the flow rate is doubled.Effluent concentration from the bottom of the specimen slightly increased during the 228-230 day of the test,due to the increment of KCl concentration from 100 mM to 200 mM at the top of the specimen.Then rapidly decreased with the increased flow circulation rate comparing the previous stages,and increased again as the solute diffused through the soil specimen until reaching steady-state.The excellent simulation results captured for this process are shown in the enlarged plot in Fig.7.The aforementioned example demonstrated that the proposed method can be used to simulate the boundary concentrations and obtain more accurate transport properties out of through-diffusion and membrane behavior tests,and further provides more accurate predictions of barrier performance of lowpermeability containment materials.

5.Conclusions

In this study,an innovative pseudo-3D analytical method is proposed to evaluate the impact of non-perfectly flushing boundary conditions on the diffusive transport and membrane behavior of soil specimens in a through-diffusion and membrane behavior test.Two scenarios(M1 and M2 scenarios)with different positions of cylindrical inflow/outflow ports that connect the porous stones and reservoirs are considered in the analyses.The following conclusions can be drawn:

(1) The accuracy of the proposed advection model was demonstrated by comparing the analytical model to numerical models under two scenarios with two different influent-toeffluent distances.In the two scenarios,the flow distribution of the M1 scenario was consistent with that from numerical models,while that for the M2 scenario resulted in a relatively larger error in the outer area of the porous disk.

(2) The solute distribution of the M1 scenario based on the analytical model,including the midpoint concentration for the measurement of membrane efficiency (ω),fits well with that based on the numerical model,whereas that for the M2 scenario requires a correcting factor of 0.84 due to the more dispersed seepage lines.

(3) Diffusion coefficient (D*) values determined based on the proposed analytical method (with proper modification for M2)are highly accurate,reducing the relative error caused by the traditional method(i.e.using the arithmetic mean of the influent and effluent concentration to calculate concentration difference in Eq.(1)) by up to 50%.To validate the proposed model,the simulated boundary concentrations are compared with that measured from a multi-stage throughdiffusion and membrane behavior test and are used to calculateD*,Rd,and ω using parameters from the steadystate portion of each testing stage.

(4) A transient simulation using the established parameters fitted the measured breakthrough curves of the example excellently,further verifying the accuracy of the proposed analytical model.Thus,the proposed model can be used to obtain accurate transport properties from the throughdiffusion and membrane behavior tests.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment

The financial support received from the Ministry of Science and Technology of the People’s Republic of China (Grant No.2019YFC1806002) and National Natural Science Foundation of China(Grant Nos.42107174,42077241)is gratefully acknowledged.The authors would also like to acknowledge the financial support provided by the U.S.National Science Foundation (NSF) under Award 1812550 for the collection of the experimental data.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2023.04.010.