A Semi-Analytical Model of Contaminant Transport in Barrier Systems with Arbitrary Numbers of Layers

: In regions with sanitary landﬁlls, unsuitable liner designs can result in signiﬁcant soil and groundwater contamination, leading to substantial environmental remediation costs. Addressing this challenge, we propose a semi-analytical model for solute transport that uses the advection–dispersion– reaction equation in a multi-layered liner system. A distinctive feature of our model is its ability to account for inﬁltration velocity, arbitrary numbers of layers, thin layers such as geomembranes, and mass ﬂow. We validated our model against existing published models and applied it to a case study of a real sanitary landﬁll in the capital of Brazil. Through parametric analyses, we simulated contaminant transport across various layers, including the geomembrane (GM), geosynthetic clay liner (GCL), soil liner (SL), and compacted clay liner (CCL). The analyses showed the importance of choosing the most appropriate construction system based on the location and availability of materials. Considering toluene contamination, a GM molecular diffusion coefﬁcient ( D GM ) greater than 10 − 13 m 2 s − 1 exhibited similar efﬁciency when compared with CCL (60 cm thick). In addition, the results showed that the liner system may have the same efﬁciency in changing SL (60 cm thick) for a GCL (1 cm thick).


Introduction
Escalating loads of municipal solid waste (MSW) end up in open dumps and landfills, producing continuous flows of landfill leachate.The risk of incorporating highly toxic landfill leachate into the environment is important to be evaluated and measured in order to facilitate decision making for landfill leachate management and treatment [1].In most countries, sanitary landfill is still the primary approach to the disposal of MSW [2], which can generate pollution in soil, air, surface water, and groundwater because of the presence of leachate [3].The contamination of areas around solid urban waste dumps is a global challenge for the maintenance of environmental quality in large urban centers in developing countries [4].
Derived from rainfall, surface drainage, and the decomposition of MSW, landfill leachate contains high concentrations of dissolved organic matter (DOM), inorganic salts, heavy metals, and xenobiotic organic compounds [5].If the landfill is not well lined with impermeable material or the lining fails due to geological conditions, landfill leachate will enter the groundwater as a contaminant front, which is one of the most significant environmental problems of landfills [6].
One solution typically adopted to minimize and prevent groundwater contamination involves the construction of a liner composed of several materials [7,8].The bottom cover layers (liners) in landfills are designed to make the leachate front in the groundwater compatible with potability standards.In all cases, the materials used must have certain technical characteristics, such as durability, adequate resistance, and low hydraulic conductivity.
Hydraulic conductivity (k) describes the capacity of a porous medium to transmit a given liquid and is a function of both the medium and the liquid [9].
One way to construct these liners involves low-permeability soil layers overlaid with a geomembrane (GM).The geomembrane in a single-or composite-liner barrier system inhibits the diffusive and advective migration of contaminants (in gas or liquid state) through the base liners and cover layer and to the surrounding environment [10][11][12].A geomembrane is a polymeric sheet material that is impervious to liquid as long as it maintains its integrity [13].
Liners formed from other materials include soil liners (SLs), compacted clay liners (CCLs), and geosynthetic clay liners (GCLs).The choice of the most appropriate material depends on the characteristics of the waste, the operation of the landfill, and the cost of materials available in the region.A GCL is used in a landfill because of its relatively low cost compared with other materials, simple operational process, and low permeability [14,15].CCLs have advantages such as good damping capacity and cheapness, as well as disadvantages such as high swelling and shrinkage [16].According to the United States Environmental Protection Agency (US EPA), the minimum thickness of CCLs should be 60 cm and the hydraulic conductivity of compacted liners should not be more than 1 × 10 −7 cm/s [13].
A GCL is a relatively thin layer of processed clay (typically bentonite) either bonded to a geomembrane or fixed between two sheets of geotextile.A geotextile is a woven or nonwoven sheet material less impervious to liquid than a geomembrane, but more resistant to penetration damage [17].GCLs are 5 to 10 mm thick and their hydraulic conductivity typically ranges from 5 × 10 −12 to 5 × 10 −11 m/s [18], and the reason for such low hydraulic conductivity is the bentonite [19].
The liner materials should have an optimal cost/benefit ratio for each project.For this reason, detailed studies should be carried out for the retention of contaminants and thus for the prevention of contamination of the environment.Currently, analytical and numerical methods are essential for predicting the transport of leachate and contaminants in porous media.Analytical solutions, the basic tools for the initial design of liners, often involving simple assumptions [20,21], are preferred because they are continuous in space and time [22].
In clay barriers, leachate contamination is mainly transported via a molecular diffusion mechanism rather than any advection mechanism because of the low hydraulic conductivity of the clays and the low leachate head (less than 30 cm) in sanitary landfills [23].However, the waste piles of landfills may generate a very high leachate head [24].Several solutions described in the literature disregard the advective term [25][26][27][28][29], which is necessary in cases where the hydraulic head is significant.In such cases, the advective term should be considered in the transport process of the contaminant.
In addition to the advective term, the first-order term (decay, biodegradation) is also essential because of its effect on the mobility of organic contaminants [30].Firstorder kinetics is the most-used model because of its simple mathematical assumption in transport models [31].Similarly, considerations involving the generation or consumption of substances (zero-order term) have also interested researchers.
Analytical solutions are extensively used to simulate the transport of contaminants in a porous medium [32,33].However, landfill cover layers commonly comprise multilayered systems.Analytical solutions regarding multi-layered systems are described in the literature, but they include many simplifications [34,35].For example, Carr [22] proposed a semi-analytical solution to solving the advection-dispersion-reaction equation in a multilayered system.Despite the solution's novelty, the author did not consider the presence of thin layers (e.g., geomembrane), the partition coefficient, or mass flow calculation.
The semi-analytical model in this research may serve as a fundamental design tool that helps evaluate the potential contamination of groundwater and environmental damage caused to leachate in sanitary landfills.The model described in this paper is not only semianalytical, but it also incorporates the complexity of multi-layered systems, recognizing the importance of thin layers as geomembranes, the influence of partition coefficients, and the criticality of mass flow calculations.As shown within the research, the model agrees well with other models in the literature, which includes those by Chen et al. [36] and Rowe and Booker [37].
In the simulations, a low-permeability material was overlaid with a geomembrane, because in the installation of the geomembrane, it is recommended to prepare the foundation area.The prepared foundation must be smooth and free of projections preventing damage to the liner; thus, even with the execution of the geomembrane, it is common to have a CCL or SL to prevent any damage.Moreover, the simulations considering the liner system with GM and molecular diffusion coefficient (D GM ) greater than 10 −13 m 2 s −1 exhibited similar efficiency when compared with CCL (60 cm thick) for toluene contamination.Meanwhile, D GM with values less than 10 −14 m 2 s −1 proved to be an efficient liner system containing the toluene front.Hence, obtaining D GM with adequate precision is necessary to achieve the non-pollution of soil and groundwater.Considering a real case scenario in a sanitary landfill in Brazil, the simulations for the current liner system of the structure showed that benzene takes 50 years to reach the bottom of the liner.However, the life span of the landfill is just 13 years, showing the efficiency of the liner system and the applicability of the model.
Given the evolving complexities and challenges in environmental geotechnics, particularly concerning sanitary landfills, the need for more refined and versatile modeling techniques is urgent.Thus, the proposed model provides a more integrated understanding and prediction of any contaminant (dissolved organic matter, inorganic salts, heavy metals, and xenobiotic organic compounds) behavior as an essential tool for contemporary geotechnical engineers.This original perspective aims to elevate the standards of landfill liner design, making it more adaptable and effective in safeguarding our environment.

Mathematical Modeling
This study employed the advection-dispersion-reaction equation tailored for heterogeneous soil profiles commonly found in landfill settings.The model explicitly incorporates various barrier systems, including GMs, GCLs, CCLs, and SLs; it reflects the typical configurations observed in modern sanitary landfills.
The assumptions of the model proposed are (i) initial concentrations within layers may vary in each material; (ii) the top surface concentration is constant or an arbitrary function of time; (iii) interface concentrations and mass fluxes are continuous between layers; and (iv) the bottom boundary condition implies that the concentration derivative is zero (Neumann boundary condition).
The parameters considered in the model are as follows: the hydrodynamic dispersion coefficient D i of each layer from i = 1, 2, . .., m; the molecular diffusion of geomembrane (D GM ); the seepage velocity v; the first-order decay constant λ; the zero-order constant γ; and the partition coefficients of the geomembrane S L,GM and S 0,GM (Figure 1).The governing equations and boundary conditions used in this paper are as follows The transport of contaminants in layer can be expressed using the advection-dispersion reaction equation [22]: The governing equations and boundary conditions used in this paper are as follows.The transport of contaminants in layer can be expressed using the advection-dispersionreaction equation [22]: where t is time (T); z is the direction of the transport of the contaminant (L); D h,i is the hydrodynamic dispersion coefficient of layer i (L 2 T −1 ); v i is the seepage velocity (LT −1 ); λ i is the first-order decay constant (T −1 ); γ i is the zero-order production/consumption constant (M L −3 T −1 ); and R i is the retardation factor for linear sorption determined using Equation (2): where ρ di is the dry density of layer I (ML −3 ); K di is the sorption distribution coefficient of layer i (L 3 M −1 ); and θ i is the volumetric water content of layer i (L 3 L −3 ).In geotechnical engineering, absorption and adsorption (referred to as sorption) can affect the concentration and mobility of solutes [38].
The contaminant biodegradation may be expressed by the first-order decay constant [7].The decay constant λ i is described as a function of the half-life (t 1/2,i ) of the contaminant [36]: where t 1/2 is the half-life of a contaminant of layer i (T), defined as the time required for the quantity of the contaminant to decay to half of its initial value.
Seepage velocity in a multi-layer system considering the thickness of the contaminant above the m layers is represented by the following Equation (3) [23]: where h w is the contaminant head (L); l i is the thickness of each layer in the system (L); and k e is the effective hydraulic conductivity of the system (LT −1 ), which is equivalent to the harmonic mean of the hydraulic conductivities of the layers [18]: where k i is the saturated hydraulic conductivity of layer i (LT −1 ).
The hydrodynamic dispersion coefficient can be expressed as the sum of the effective coefficient of molecular diffusion D m,I (L 2 T −1 ) and the dispersion coefficient D D,I (L 2 T −1 ): The solution to the multi-layer contaminant transport equation depends on the initial and boundary conditions.The initial condition can be expressed as follows [18,23]: where c i is the initial concentration of the contaminant in layer i (M L −3 ).
The top boundary condition of the imposed concentration can be approximated by a first-order rate equation [3,36,37]: where c 0 is the peak concentration of the contaminant (M L −3 ) and λ is the decay constant (T −1 ) described in Equation ( 3).When the concentration of the contaminant in the overlying medium is not available, a constant concentration c 0 is often applied.Thus, Equation ( 8) becomes [18,23,24] c u (0, t) = c 0 (9) In the transport of contaminants in multi-layers (Figure 1), we need to account for the continuity between layers.Hence, the contaminant concentration and dispersive flux are considered continuous at the interfaces between adjacent layers [22]: When considering the seepage velocity, we obtain [23] If the first layer is composed of the GM (Figure 1a), then the general transport equation, Equation (1), becomes [18,23,24] where D GM is the molecular diffusion coefficient of GM (L 2 T −1 ).At the interface between the GM and contaminant, the inlet boundary condition is [23] c GM (0, t) = c 0 (0, t)S 0,GM (14) where S 0,GM is the partition coefficient that indicates the ratio of the concentration of the contaminant to the concentration in GM.
The interface boundary condition between the GM and underlying SL is defined by the following equation [21]: where S L,GM is the partition coefficient indicating the ratio of the concentration of the GM contaminant to that in the soil's pore water at the interface of the two layers.We now assume that the interface boundary conditions can be described as [23] The top boundary condition at the inlet (z = 0) becomes Equation ( 17) into Equations ( 15) and ( 16), we obtain The GM is not at the top of the system but between two layers (Figure 1b), then Equations ( 17), (19), and ( 20) are still valid.
The bottom boundary condition at the last layer is of the following type [18,23,24]: Based on the initial and boundary conditions and the continuity of the interfaces, the solution of the equation for a multi-layer system can be proposed.In Appendix A, the solution is presented.
Based on the advection and dispersion mass balance, the mass flux J of an arbitrary layer i is [37,39] as follows:

Model Validation
The model was validated using the simulation made by Chen et al. [36].The authors analyzed data from existing numerical solutions [21,37] to compare with their analytical solutions.The parameters of the CCL were taken from Foose et al. [21], and the parameters of the SL were obtained from Rowe and Booker [37].All parameters are listed in Table 1 and are related to the toluene contaminant.The contaminant concentration profiles of the two-layered soil were compared to the semi-analytical solution (Figure 2) and the solution presented by Chen et al. [36].The boundary conditions included a zero initial concentration (c i = 0), constant top-surface concentration (c u [0,t] = 100 µg L −1 ), and a Neumann bottom boundary (Equation ( 21)).The times (10 and 20 years) are simulation periods often used to design liner systems.The triangle in the lower left corner represents the contaminant head (e.g., leachate).
As shown in Figure 2, the semi-analytical solution used in this research agrees well with other models in the literature, which includes those by Chen et al. [36] and Rowe and Booker [37] (the latter of which was validated by the authors of the former).Benzene and acetone were selected to represent hydrophobic and hydrophilic organic compounds, respectively, as these contaminants are usually found in leachate compounds [36].Differently from Chen et al. [36], the model in the present study also accounts for the presence of a GM. Figure 3 shows the validation with a GM (1.5 mm) and a CCL (0.6 m).Table 2 presents the simulation parameters along with the values for the high-density polyethylene (HDPE) geomembrane [36].The results confirm the applicability of the model.
The contaminant concentration profiles of the two-layered soil were compared to the semi-analytical solution (Figure 2) and the solution presented by Chen et al. [36].The boundary conditions included a zero initial concentration (ci = 0), constant top-surface concentration (cu[0,t] = 100 µg L −1 ), and a Neumann bottom boundary (Equation ( 21)).The times (10 and 20 years) are simulation periods often used to design liner systems.The triangle in the lower left corner represents the contaminant head (e.g., leachate).As shown in Figure 2, the semi-analytical solution used in this research agrees well with other models in the literature, which includes those by Chen et al. [36] and Rowe and Booker [37] (the latter of which was validated by the authors of the former).Benzene and acetone were selected to represent hydrophobic and hydrophilic organic compounds, respectively, as these contaminants are usually found in leachate compounds [36].Differently from Chen et al. [36], the model in the present study also accounts for the presence of a GM. Figure 3 shows the validation with a GM (1.5 mm) and a CCL (0.6 m).Table 2 presents the simulation parameters along with the values for the high-density polyethylene (HDPE) geomembrane [36].The results confirm the applicability of the model.As shown in Figure 2, the semi-analytical solution used in this research agrees well with other models in the literature, which includes those by Chen et al. [36] and Rowe and Booker [37] (the latter of which was validated by the authors of the former).Benzene and acetone were selected to represent hydrophobic and hydrophilic organic compounds, respectively, as these contaminants are usually found in leachate compounds [36].Differently from Chen et al. [36], the model in the present study also accounts for the presence of a GM. Figure 3 shows the validation with a GM (1.5 mm) and a CCL (0.6 m).Table 2 presents the simulation parameters along with the values for the high-density polyethylene (HDPE) geomembrane [36].The results confirm the applicability of the model.Table 2. Parameters of materials for the simulation [36].

Parametric Evaluation
Parametric studies were conducted to validate the versatility of the model under varying conditions, encompassing combinations of compacted clay liners (CCLs), SLs, geomembranes (GMs), and GCLs.Table 3 lists the material parameters for toluene, as documented by Xie et al. [40].
Table 4 shows scenarios altering the diffusion coefficient and GM thickness in the presence of a CCL (parameters in Table 3).In the simulation, the values of D GM were 10 −11 m 2 s −1 for a 0.75 mm HDPE GM [41], 10 −13 m 2 s −1 for a 1.5 mm HDPE GM [40], and 10 −14 m 2 s −1 for a 2.5 mm HDPE GM [41].The values of S 0,GM and S L,GM were set to 100.In addition, the dispersion was ignored, as in the studies from where the parameters were obtained.Figure 4 shows values for the relative concentration (Figure 4a) and the flow at the base (Figure 4b) in the four scenarios.The first three scenarios (SI, SII, and SIII) involved constant heads of toluene (0.3 m) and CCL (0.6 m) but used different diffusion coefficients for GM (D GM ).For a comparative analysis, scenario SIV, having only a CCL (0.6 m), was considered.The results show that D GM values greater than 10 −13 m 2 s −1 (SI and SII) exhibited similar efficacies when compared with the presence of CCL alone (SIV).However, D GM = 10 −14 m 2 s −1 (SIII) proved more efficient than SIV.These analyses show the importance of choosing the most appropriate construction system based on the location and availability of materials.Furthermore, they show the need to obtain D GM with adequate precision to achieve the non-pollution of groundwater.Figure 4b shows the flow at the base of the system.
Sustainability 2023, 15, x FOR PEER REVIEW 9 of 19 achieve the non-pollution of groundwater.Figure 4b shows the flow at the base of the system.
Figure 5 shows a two-layer system composed of a GM and CCL, simulating scenario SIII in Table 4 but with varying toluene head values (hw = 0.3, 1.5, and 10 m).The initial concentration was considered null, and the top boundary condition equaled 100 mg L -1 .Figure 5a,b illustrate how the toluene breakthrough curve and the flow were affected, respectively, depending on the head of the contaminant (Equation ( 4)).Similarly, Figure 6 shows results for only a 0.6 m thick CCL.A comparison of Figures 5 and 6 explains that the GM causes the greater mitigation of contaminant transport.However, depending on the maximum value allowed for the contaminant (e.g., groundwater, soil) and the landfill's lifespan, sometimes, only CCL may be used.Figure 5 shows a two-layer system composed of a GM and CCL, simulating scenario SIII in Table 4 but with varying toluene head values (h w = 0.3, 1.5, and 10 m).The initial concentration was considered null, and the top boundary condition equaled 100 mg L −1 .Figure 5a,b illustrate how the toluene breakthrough curve and the flow were affected, respectively, depending on the head of the contaminant (Equation ( 4)).Similarly, Figure 6 shows results for only a 0.6 m thick CCL.A comparison of Figures 5 and 6 explains that the GM causes the greater mitigation of contaminant transport.However, depending on the maximum value allowed for the contaminant (e.g., groundwater, soil) and the landfill's lifespan, sometimes, only CCL may be used.Figure 7 shows the results for a three-layer system using the parameters for CCL and SL presented in Table 3.The GM parameters matched those in scenario SII in Table 4 (DGM = 10 −13 m 2 s −1 and a thickness of 1.5 mm).In addition, the initial concentration was null, and the upper boundary condition was considered constant at 100 mg L −1 with a toluene head of 0.3 m. Figure 7a illustrates a type of material distribution (CCL, GM, and SL) often used in landfills because of the intercalation of the GM with different soils.Figure 7b shows the intercalation of GM with GCL and a layer of CCL at the base.In both analyses of Figure 7, the contaminant first-order term (decay/biodegradation) was considered only in soil layers (CCL and SL) using Equation (3).Thus, the higher the contaminant half-life t1/2, the greater the flow at the system's base.In addition, the system formed by CCL, GM, and SL provides a lower flow of toluene, mainly because of the thickness and hydraulic conductivity values of the total system (Figure 7a).However, even with the different thicknesses of GCL (0.01 m) and SL (0.6 m), the change in flow is not significant.Figure 7 shows the results for a three-layer system using the parameters for CCL and SL presented in Table 3.The GM parameters matched those in scenario SII in Table 4 (D GM = 10 −13 m 2 s −1 and a thickness of 1.5 mm).In addition, the initial concentration was null, and the upper boundary condition was considered constant at 100 mg L −1 with a toluene head of 0.3 m. Figure 7a illustrates a type of material distribution (CCL, GM, and SL) often used in landfills because of the intercalation of the GM with different soils.Figure 7b shows the intercalation of GM with GCL and a layer of CCL at the base.In both analyses of Figure 7, the contaminant first-order term (decay/biodegradation) was considered only in soil layers (CCL and SL) using Equation (3).Thus, the higher the contaminant half-life t 1/2 , the greater the flow at the system's base.In addition, the system formed by CCL, GM, and SL provides a lower flow of toluene, mainly because of the thickness and hydraulic conductivity values of the total system (Figure 7a).However, even with the different thicknesses of GCL (0.01 m) and SL (0.6 m), the change in flow is not significant.
To analyze the difference in GM parameters and their effect on the transport of toluene, relative concentrations were compared, as shown in Figure 8, for a three-layer system (CCL, GM, and SL) with a toluene head of 0.3 m.The parameters used for CCL and SL in Figure 8a,b are found in Table 3.The GM parameters from scenarios SII and SIII (Table 4) are used in Figure 8a,b, respectively.The imposed concentration was 100 mg L −1 .Thus, we see that a factor of 10 in D GM causes a significant delay in the contamination front.
intercalation of GM with GCL and a layer of CCL at the base.In both analyses of Figure 7, the contaminant first-order term (decay/biodegradation) was considered only in soil layers (CCL and SL) using Equation (3).Thus, the higher the contaminant half-life t1/2, the greater the flow at the system's base.In addition, the system formed by CCL, GM, and SL provides a lower flow of toluene, mainly because of the thickness and hydraulic conductivity values of the total system (Figure 7a).However, even with the different thicknesses of GCL (0.01 m) and SL (0.6 m), the change in flow is not significant.To analyze the difference in GM parameters and their effect on the transport of toluene, relative concentrations were compared, as shown in Figure 8, for a three-layer system (CCL, GM, and SL) with a toluene head of 0.3 m.The parameters used for CCL and SL in Figure 8a,b are found in Table 3.The GM parameters from scenarios SII and SIII (Table 4) are used in Figure 8a,b, respectively.The imposed concentration was 100 mg L −1 .Thus, we see that a factor of 10 in DGM causes a significant delay in the contamination front.Figure 9 shows a situation in which the upper boundary condition is equal to c0(0, t) = 100 exp (−λt) in a two-layer system (GM + SL) with a different half-life.The first-order term in this analysis was only considered for the boundary condition and not for the SL (Equation ( 8)).The toluene head was 1 m, the SL parameters are from Table 3, and the GM parameters were from scenario SII (Table 4).For smaller values of t1/2, the decrease was more significant in the toluene breakthrough curve.Therefore, analyses considering this boundary condition must be well evaluated for safety reasons.If t1/2 tends to infinity (solid line in Figure 9), the situation becomes similar to the analyses in Figure 5. Figure 9 shows a situation in which the upper boundary condition is equal to c 0 (0, t) = 100 exp (−λt) in a two-layer system (GM + SL) with a different half-life.The first-order term in this analysis was only considered for the boundary condition and not for the SL (Equation ( 8)).The toluene head was 1 m, the SL parameters are from Table 3, and the GM parameters were from scenario SII (Table 4).For smaller values of t 1/2 , the decrease was more significant in the toluene breakthrough curve.Therefore, analyses considering this boundary condition must be well evaluated for safety reasons.If t 1/2 tends to infinity (solid line in Figure 9), the situation becomes similar to the analyses in Figure 5.
The parameters in Table 5 were used to simulate a four-layer system considering the generation of a contaminant (zero-order term).Figure 10 shows the results of an analysis incorporating a protective layer (PL), drainage layer (DL), geomembrane (GM), and a compacted clay liner (CCL).A contaminant was present at the top of the PL (c 0 = 100 µL −1 ), and a zero initial condition was used for all layers.The seepage velocity and the dispersive terms were ignored in all layers.However, a contaminant (γ DL > 0) was generated at the DL, indicating problems with this layer.As shown in Figure 10, the flux changed by a factor of ten, comparing the non-generation (γ DL = 0) with the highest generation considered (γ = 10 −2 ).This result shows the value of knowing the presence of a contaminant that can alter the flux in the system significantly.Figure 9 shows a situation in which the upper boundary condition is equal to c0(0, t) = 100 exp (−λt) in a two-layer system (GM + SL) with a different half-life.The first-order term in this analysis was only considered for the boundary condition and not for the SL (Equation ( 8)).The toluene head was 1 m, the SL parameters are from Table 3, and the GM parameters were from scenario SII (Table 4).For smaller values of t1/2, the decrease was more significant in the toluene breakthrough curve.Therefore, analyses considering this boundary condition must be well evaluated for safety reasons.If t1/2 tends to infinity (solid line in Figure 9), the situation becomes similar to the analyses in Figure 5.The parameters in Table 5 were used to simulate a four-layer system considering the generation of a contaminant (zero-order term).Figure 10 shows the results of an analysis incorporating a protective layer (PL), drainage layer (DL), geomembrane (GM), and a compacted clay liner (CCL).A contaminant was present at the top of the PL (c0 = 100 µL −1 ), and a zero initial condition was used for all layers.The seepage velocity and the dispersive terms were ignored in all layers.However, a contaminant (γDL > 0) was generated at the DL, indicating problems with this layer.As shown in Figure 10, the flux changed by a factor of ten, comparing the non-generation (γDL = 0) with the highest generation considered (γ = 10 −2 ).This result shows the value of knowing the presence of a contaminant that can alter the flux in the system significantly.

Application to Real Sanitary Landfill
The model developed was applied in a real sanitary landfill to verify its performance.The case study was the Brasília Sanitary Landfill (BSL), located in the capital of Brazil.The BSL was opened on 17 January 2017.As of June 2023, the BSL has already received 4,212,142 tons of waste [42].The BSL's liner system comprises a triple layer of 0.5 m topsoil CCL (CCL1) above an HDPE geomembrane that overlays a 1.5 m thick CCL (CCL2) layer.The current total clay thickness is 2 m and the life span of the BSL is 13 years.
All parameters of the landfill materials for the benzene contaminant are shown in Table 6, as mentioned in Baran [43].To assess the performance of the geomembrane with an extra layer of clay, the parameters shown for CCL3 were defined.

Application to Real Sanitary Landfill
The model developed was applied in a real sanitary landfill to verify its performance.The case study was the Brasília Sanitary Landfill (BSL), located in the capital of Brazil.The BSL was opened on 17 January 2017.As of June 2023, the BSL has already received 4,212,142 tons of waste [42].The BSL's liner system comprises a triple layer of 0.5 m topsoil CCL (CCL1) above an HDPE geomembrane that overlays a 1.5 m thick CCL (CCL2) layer.The current total clay thickness is 2 m and the life span of the BSL is 13 years.
All parameters of the landfill materials for the benzene contaminant are shown in Table 6, as mentioned in Baran [43].To assess the performance of the geomembrane with an extra layer of clay, the parameters shown for CCL3 were defined.Breakthrough and flow curves for different benzene heads (0.5, 1, and 2 m) were evaluated in a three-layer system, as shown in Figure 11.The initial concentration of the contaminant was considered null.The upper boundary condition was constant at 5 µg L −1 , the maximum allowed value (MAV) for water potability according to Brazilian standards [44].The benzene concentration reached the MAV after 100 years of simulation for all benzene heads.The landfill was designed for a lifespan of approximately 13 years.Therefore, the designed liner is effective against benzene transport, exhibiting good performance.The model can be applied for the most critical contaminants in the leachate of the site to assess safety against groundwater contamination.
Using the values in Table 6, we compared the current construction system of BSL (CCL1, GM, and CCL2) with a system including only one CCL layer of 2.5 m (CCL3) and a benzene head of 2 m.The comparative layer did not contain GM. Figure 12 shows the breakthrough curve of the current landfill system (dashed line) and the modified system (continuous line).The result shows that adding 0.5 m of soil thickness results in a more effective system-without GM-for a period longer than 50 years, which is significantly greater than the expected lifespan of the landfill.However, a cost-benefit analysis must be carried out before the optimal construction system can be defined.The benzene concentration reached the MAV after 100 years of simulation for all benzene heads.The landfill was designed for a lifespan of approximately 13 years.Therefore, the designed liner is effective against benzene transport, exhibiting good performance.The model can be applied for the most critical contaminants in the leachate of the site to assess safety against groundwater contamination.
Using the values in Table 6, we compared the current construction system of BSL (CCL1, GM, and CCL2) with a system including only one CCL layer of 2.5 m (CCL3) and a benzene head of 2 m.The comparative layer did not contain GM. Figure 12 shows the breakthrough curve of the current landfill system (dashed line) and the modified system (continuous line).The result shows that adding 0.5 m of soil thickness results in a more effective system-without GM-for a period longer than 50 years, which is significantly greater than the expected lifespan of the landfill.However, a cost-benefit analysis must be carried out before the optimal construction system can be defined.
(CCL1, GM, and CCL2) with a system including only one CCL layer of 2.5 m (CCL3) and a benzene head of 2 m.The comparative layer did not contain GM. Figure 12 shows the breakthrough curve of the current landfill system (dashed line) and the modified system (continuous line).The result shows that adding 0.5 m of soil thickness results in a more effective system-without GM-for a period longer than 50 years, which is significantly greater than the expected lifespan of the landfill.However, a cost-benefit analysis must be carried out before the optimal construction system can be defined.

Conclusions
Analytical models of contaminant transport are robust tools for parametric and numerical validation.Notably, both analytical and semi-analytical approaches maintain continuity across spatial and temporal domains.This research extends a semi-analytical model of the advection-dispersion-reaction equation to encapsulate these benefits.

Conclusions
Analytical models of contaminant transport are robust tools for parametric and numerical validation.Notably, both analytical and semi-analytical approaches maintain continuity across spatial and temporal domains.This research extends a semi-analytical model of the advection-dispersion-reaction equation to encapsulate these benefits.
Model validation, employing materials such as compacted clay liners (CCLs), SLs, GCLs, and geomembranes (GMs), effectively demonstrated the model's capability to simulate multi-layer systems.The inherent heterogeneity of these materials indicates their critical importance when designing geostructures, especially landfill liners.The accurate comprehension of the transport of contaminants (like leachate) is pivotal for refining these geostructural designs.
The use of analytical and semi-analytical models eases parametric analysis greatly.Our research indicates the value of determining the geomembrane's diffusion coefficient through lab analyses, facilitating an accurate prediction of contaminant transport and highlighting the trade-offs between opting for a geomembrane (GM) or augmenting the thickness of traditional soil liners (e.g., CCLs and SLs).Regardless, the geomembrane's essential role as an advective flow barrier remains undebatable, necessitating rigorous checks to ensure its integrity.Moreover, incorporating first-order (decay/biodegradation) and zero-order (species source) terms is critical for a holistic prediction of contaminant transport dynamics.
The open-source model introduced here may be regarded as a versatile tool for simulating multi-layer systems to aid in optimizing landfill liner designs and thus contributing directly to minimizing soil contamination and mitigating groundwater pollution.Aiming for more comprehensive application, future research should investigate the applicability of the solution comparing with laboratory results analyzing multilayers of materials.Furthermore, for a holistic consideration of the contaminant transport phenomenon, the implementation of the material saturation (CCL and SL) and temperature in the model is recommended.

Figure 1 .
Figure 1.Schematic diagram of the transport of contaminant through multiple layers: (a) GM as first layer and (b) GM as intermediate layer.

Figure 1 .
Figure 1.Schematic diagram of the transport of contaminant through multiple layers: (a) GM as first layer and (b) GM as intermediate layer.

Figure 2 .
Figure 2. Comparison of the semi-analytical model and the numerical model of Chen et al. [36].

Figure 2 .
Figure 2. Comparison of the semi-analytical model and the numerical model of Chen et al. [36].

Figure 2 .
Figure 2. Comparison of the semi-analytical model and the numerical model of Chen et al. [36].

Figure 3 .
Figure 3. Validation of the proposed model with the contaminants (a) acetone and (b) benzene [36].

Figure 4 .
Figure 4. Parametric analysis of a two-layer system: (a) relative concentration and (b) flux at the base.Figure 4. Parametric analysis of a two-layer system: (a) relative concentration and (b) flux at the base.

Figure 4 .
Figure 4. Parametric analysis of a two-layer system: (a) relative concentration and (b) flux at the base.Figure 4. Parametric analysis of a two-layer system: (a) relative concentration and (b) flux at the base.

Figure 4 .Figure 5 .
Figure 4. Parametric analysis of a two-layer system: (a) relative concentration and (b) flux at the base.

Figure 5 .Figure 6 .
Figure 5. Toluene transport for different heads: (a) relative concentration and (b) flux at the base.Figure 5. Toluene transport for different heads: (a) relative concentration and (b) flux at the base.Sustainability 2023, 15, x FOR PEER REVIEW 10 of 19

Figure 6 .
Figure 6.Toluene transport for different heads: (a) relative concentration and (b) flux at the base.

Figure 7 .
Figure 7. Toluene flow at the base of three-layer systems: (a) CCL, GM, and SL system; (b) GM, GCL, and CCL system.

Figure 7 .
Figure 7. Toluene flow at the base of three-layer systems: (a) CCL, GM, and SL system; (b) GM, GCL, and CCL system.

Figure 9 .
Figure 9. Imposed concentration with first-order rate equation for different half-lives: (a) relative concentration; (b) flux at the base of the liner system.

Figure 9 .
Figure 9. Imposed concentration with first-order rate equation for different half-lives: (a) relative concentration; (b) flux at the base of the liner system.

Figure 11 .
Figure 11.Contaminant transport in possible multi-liner system for the BSL: (a) concentration profile; (b) flux in a three-layer system.

Figure 11 .
Figure 11.Contaminant transport in possible multi-liner system for the BSL: (a) concentration profile; (b) flux in a three-layer system.

Figure 12 .
Figure 12.Comparison of constructive models of landfill liners.

Figure 12 .
Figure 12.Comparison of constructive models of landfill liners.