Importance of Inﬁltration Rates for Fate and Transport of Benzene in High-Tiered Risk-Based Assessment Considering Korean Site-Speciﬁc Factors at Contaminated Sites

: Widely used conservative approaches for risk-based assessments of the subsurface transport processes have been calculated using simple analytical equations or general default values. Higher-tier risk assessment of contaminated sites requires the numerical models or additional site-speciﬁc values for input parameters. Previous studies have focused on the development of sophisticated models ﬁt into risk-based frameworks. Our study mainly aims to explore the applicability of site-speciﬁc parameters and to modify the risk-based fate and transport model according to the types of the site-speciﬁc parameters. To apply the modiﬁed fate and transport equation and the site-speciﬁc default inﬁltration range, this study assessed the source depletion, leachate concentrations, and exposure concentration of benzene, which is a representative organic contaminant. The numerical models consist of two continuous processes, the fate and transport of contaminants from (1) the soil to the groundwater table in the vadose zone and subsequently (2) from the groundwater table to exposure wells in the saturated zone. Spatially varied Korean domestic recharge data were successfully incorporated into site-speciﬁc inﬁltration parameters in the models. The numerical simulation results were expressed as transient time series of concentrations over time. Results presented the narrow range of predicted concentrations at the groundwater table when site-speciﬁc inﬁltration was applied, and the dilution–attenuation factors for the unsaturated zone (DAF unsat ) were derived based on the prediction. When a contaminant travels to the longest path length of 10 m with a source depth of 1 m in the vadoze zone, the simulated DAF unsat ranged from 3 to 4. The highest DAF unsat simulation results are close to 1 when contaminants travel to a source depth of 5 m and the shortest path length of 1 m. In the saturated aquifer below the contaminated sites, the variation in exposure concentration with time at monitoring wells is detected differently depending on the depth of the saturated zone. investigated a test bed of Cheongju area, where the groundwater cultivation area was precisely calculated. Figure 7a shows the monthly temporal variation of precipitation and inﬁltration assessed by Chang and Chung [43] from January 2011 to January 2014 using a hydrological analysis tool, the SWAT model [42] developed by the USDA Agricultural Research Service (ARS), to simulate the groundwater stasis state. As input data for the SWAT model, meteorological data such as precipitation, temperature, wind speed, insola-tion, and relative humidity were used at the Cheongju Meteorological Observatory located in the Musimcheon Basin. The average amount of groundwater recharge in this area was estimated to be ≈ 13.9% of precipitation.


Introduction
The Soil Screening Guidance for risk-based soil screening levels (SSLs) of the United States Environmental Protection Agency (USEPA) provides a framework for screening contaminated soils [1]. For instance, the EPA selected a default dilution and attenuation factor (DAF), which is defined as the ratio of the source zone concentration to the receptor point concentration, of 20 for U.S. alluvial aquifers based on the EPACMTP model [2]. The five organic contaminants in the unsaturated zone in various soil conditions using several analytical solutions.
One of the general assumptions in the lower-tiered approach is infinite sources without losses over a period. This assumption may be oversimplified or too conservative for easily volatile and degradable chemicals, such as benzene and naphthalene [25,29]. During field studies, it is also important to characterize the source zone [20,30]. Feenstra et al. [31] described a method for calculating the pore water concentration of soil using residual NAPL. Garg et al. [32] showed that the reduction of contaminants via biodegradation differs depending on the composition of hydrocarbons. Mackay et al. [33] showed that the depletion of sources based on decomposition in the unsaturated zone is relatively higher than that decomposition caused by transport from the groundwater at oil leak sites.
Default values are used for the chemical and physical properties of contaminants and some properties of hydrogeological factors [2]. Site-specific input variables consider hydraulic conductivity values reflecting soil characteristics, the soil pH, various biological mechanisms, and seasonal changes from the unsaturated zone to the groundwater. In Korea, the groundwater flow model in the unsaturated zone is actively used to understand groundwater level fluctuations and predict the amount of groundwater recharge [34].
Infiltration rate is a parameter that affects the fate and transport of contaminants in the unsaturated zone [24,25,35]. Verginelli and Baciocchi [24] considered infiltration factors when simulating source depletion, and the leaching concentration eventually approaches zero. Ryu [25] analytically conducted a sensitivity test by varying the infiltration rate from 10 −6 to 1 m. As a result, the concentration of benzene at the water table of the unsaturated zone was diluted and attenuated by 10 4 to 10. Chang et al. [35] considered the range of infiltration from 0.05 to 1 m/year.
The determination process of the groundwater exposure concentration and the following input default values illustrated in the Korean soil risk assessment guidelines [36] is mostly based on the Soil Screening Guidance for risk-based SSLs of the USEPA [1] based on RBCA's Tier 1 approaches [3]. The specific objectives of this study are to apply Korean site-specific parameters to a risk-based framework. Based on previous studies that have attempted to quantify the exposure concentration of the subsurface pathway in the risk assessment framework [24,25,27,37], we explore high-tiered approaches of determination of the benzene pathway in the subsurface at the contaminated sites by applying site-specific parameters, especially infiltration data. The method considers (1) sophisticated temporal evaluation in the numerical model and (2) site-specific input parameters using the Korean domestic hydrological database. In the numerical simulations, two site-specific factors in the governing equation were considered: (1) a leaching process considering source depletion in the vadose zone at the contaminated sites and (2) the temporally varying infiltration rate.

Conceptualized Numerical Model to Simulate Pathways in the Unsaturated and Saturated Zone
The key parameters for model input, simulation process, and expected output are summarized in Figure 1. The simulation was conducted using a simplified conceptual model with minimal input variables and various assumptions. Hydrogeological and chemical input variables are selected for the flow and transport equations. As presented in Figure 1, infiltration rate, source geometry (source depth and width), source mass and concentration, soil physical properties (porosity, volumetric water contents, volumetric air contents, fraction of organic carbon), and chemical properties (distribution coefficient, Henry's law constant) are selected as the key parameters in the unsaturated zone. The model is simulated using a one-dimensional finite-difference model, and results are expressed in the concentration at water level in Darcy flux. As a next step, modeling of the saturated aquifer is performed. Since the saturated model is a two-dimensional crosscut, the width and depth of the model are required. The temporal concentration profile obtained from the unsaturated zone is used as an input for the recharge in the saturated model. Then, assuming there is a monitoring well is installed, the concentration at the exposure point can be observed. Biodegradation was not considered in this analysis.
Water 2021, 13, x FOR PEER REVIEW 4 saturated aquifer is performed. Since the saturated model is a two-dimensional cross the width and depth of the model are required. The temporal concentration profile tained from the unsaturated zone is used as an input for the recharge in the satur model. Then, assuming there is a monitoring well is installed, the concentration at exposure point can be observed. Biodegradation was not considered in this analysis.  Figure 2 displays the conceptual model for the unsaturated zone and the satur aquifer. The conceptual model has been developed by previous studies [2,27]. We si lated the pathways of the dissolved concentration of benzene from the source to groundwater table and then the groundwater table to the exposure point. The satur aquifer below the contaminated site is often defined as a mixing zone because some servative approaches assume a mixed process in this area. We did not consider fur migration out of the mixing zone because geometry and complex boundary condition the flow system with additional site investigation should be considered. Mazzieri e [27] also evaluated the transport of contaminants from soil to the mixing zone to ana the dilution of the contaminants.    Figure 2 displays the conceptual model for the unsaturated zone and the saturated aquifer. The conceptual model has been developed by previous studies [2,27]. We simulated the pathways of the dissolved concentration of benzene from the source to the groundwater table and then the groundwater table to the exposure point. The saturated aquifer below the contaminated site is often defined as a mixing zone because some conservative approaches assume a mixed process in this area. We did not consider further migration out of the mixing zone because geometry and complex boundary conditions in the flow system with additional site investigation should be considered. Mazzieri et al. [27] also evaluated the transport of contaminants from soil to the mixing zone to analyze the dilution of the contaminants. saturated aquifer is performed. Since the saturated model is a two-dimensional crosscut, the width and depth of the model are required. The temporal concentration profile obtained from the unsaturated zone is used as an input for the recharge in the saturated model. Then, assuming there is a monitoring well is installed, the concentration at the exposure point can be observed. Biodegradation was not considered in this analysis.  [2,27]. We simulated the pathways of the dissolved concentration of benzene from the source to the groundwater table and then the groundwater table to the exposure point. The saturated aquifer below the contaminated site is often defined as a mixing zone because some conservative approaches assume a mixed process in this area. We did not consider further migration out of the mixing zone because geometry and complex boundary conditions in the flow system with additional site investigation should be considered. Mazzieri et al. [27] also evaluated the transport of contaminants from soil to the mixing zone to analyze the dilution of the contaminants.   Model 1 was applied for evaluating the unsaturated zone, and it assumes the occurrence of one-dimensional fate and transport processes. In the model, the z-direction was defined as the direction of gravity and the flow direction for a one-dimensional finitedifference model. The model assumed unit length for the x-direction and y-direction of the grids. The grid size of the z-direction was set to 0.1 m. The depth of the model was varied from 1 to 10 m for the estimation of the DAF of the unsaturated zone. When the case study was simulated, the length was then set to 1 m and 5 m, respectively.
In model 2, we simulated transport using a two-dimensional model for the saturated aquifer that does not consider the horizontal geometry of the contaminant source. The length of the contaminant source parallel to the groundwater was 50 m, and the depth of the aquifer was 10 m. Therefore, the model domain was set to 50 m horizontally and 10 m vertically. A grid of ∆x = 5 m, ∆y = 1 m, and ∆z = 1 was built in 11 rows: 1 column and 10 layers. Using the unit width of the column, the model can describe the crosscut of the contaminated site shown in two-dimensional system. We assumed that groundwater flowing in two directions flows through the model area and is discharged to the right boundary, and the change in concentrations over time was monitored at the right boundary. (1) Leaching from source Analytical Equations (1)-(5) were used to calculate the concentration at the water table. Verginelli and Baciocchi [24] established an integrated method to determine the entire pathway of the contaminants from leaching to entering the groundwater as follows: where C t is the soil exposure concentration (M/M), ρ b (M/L 3 ) is the volume density of the soil, H is the dimensionless Henry's constant, K d (L 3 /M) is the soil distribution coefficient, θ w is the soil moisture content, and θ a is the soil air content. If the calculated C L (0) value exceeds the saturation concentration, it is replaced with the saturation concentration C L (0) or the solubility value. The calculation for the soil distribution coefficient of organic contaminants is described below. The equation corresponding to the denominator of Equation (1) is called the soil-water partition coefficient, K SW . Assuming an environment without NAPL, the change in the concentration over time due to infiltration and biodegradation can be expressed by the following equation: where d s is the thickness of the contaminant source [L], R is the retardation coefficient [L/M], and λ source is the first-order biodegradation rate constant in the source [T −1 ]. The first term in Equation (3) reflects the source depletion due to leaching. The second term expresses the source depletion due to biodegradation, and µ is the source depletion coefficient expressed as the sum of the two values.
If C L (0) calculated by Equation (1) exceeds the saturation concentration and is replaced by the solubility value, Equation (2) can be rewritten as follows: where S is the solubility of the contaminant (M/L 3 ) and t* (T) is the point at which the soil discharge concentration decreases below the solubility due to the depletion of the source. As Equations (3) and (4) reflect the average infiltration rate, we modified the equation to incorporate the temporal change in infiltration: where µ i is the decomposition constant calculated on the precipitation in the i th term among stress periods of number n. In this study, the i th stress period refers to the i th year in the simulation time.
(2) Fate and transport in the unsaturated zone (Model 1)-Soil to groundwater pathway Models related to solute movement and behavior were used to predict advection, dispersion, and diffusion adsorption and are expressed as a one-dimensional partial differential equation as shown in the following equation. The model assumes that the vertical flow is dominant in the unsaturated subsurface. The equation is mainly based on the transport-diffusion equation with linear sorption. The contaminant pathways in the saturated zones are calculated based on the assumption of groundwater saturation flow in a homogeneous medium. The governing equation for flow in the saturated zone can be expressed as a partial differential equation as follows.
is the flow rate per unit volume flowing into or out of the groundwater system, [T −1 ] is the source/sink, and x, y, and z are Cartesian coordinates in which x and y are horizontal and z is vertical. D is the dispersion The left and boundary heads are fixed so that the head slope remains at 0.01. The concentration and flux discharged from model 1 were used as input to areal recharge.
The governing equation for transport in the saturated zone can be expressed as a partial differential equation as follows. The equation is based mainly on the transportdiffusion equation with linear sorption. Decay is considered only in a case study. The boundary condition and the initial condition of the saturated model is shown below: at t = 0, C = 0 at all nodes; at x = 0, C = 0; at x = L, ∂C ∂x = 0. The code used the hydrogeological default values presented in Table 1 as default input variable values, which are based on the Korean soil risk assessment guideline [34]. The distance between the bottom of the source and water table, L, varies from 1 to 10 m, and the thickness of the contaminant source, D, varies 1 to 5 m in this study. The first-order decay coefficient is set to zero. The guideline recommends a fraction of organic carbon content, f oc of 0.002, representing sandy soil types. The Korean Guideline by the Ministry of Environment (MOE) [36] follows the guideline by EPA [2] for a dimensionless Henry constant and distribution of coefficient values. The infiltration rate was applied by calculating the range of the domestic areal-recharge rate that was nationally investigated.

National Infiltration Data Acquisition for Model Input
To determine the range of effective infiltration rates that we can find in Korea, we used the recharge data from the Korean groundwater survey conducted by the Ministry of Land, Infrastructure, and Transport (MOLIT) and the MOE. This survey has been conducted for 100 areas since 1997 to investigate the aquifer characteristics and assess the water budget of groundwater resources by region. It also provides basic data on the groundwater flow direction, groundwater recharge rate, groundwater vulnerability map, and sustainable development [38]. In the processes related to the surveys, The Korea Water Resources Corporation, Korea Institute of Geoscience and Mineral Resources, Korea Institute of Construction Technology, Korea Groundwater, and the Geothermal Association calculated the recharge amount introduced into the groundwater to evaluate the potential groundwater development for each administrative district at the city and county levels. Multiple analysis methods, such as water budget analysis, Natural Resources Conservation Service (NRCS)-Curve Number (CN) [39], Water-Table Fluctuation (WTF) method [40], Hybrid WTF [34,41], and Soil and Water Assessment (SWAT) [42] were selected depending on the case. Figure 3 shows the spatial distribution of precipitation, recharge, and recharge rate surveyed since 1997 on the Korean peninsula. The graph shows data from basic surveys of 100 cities and counties conducted since 1997 that were collected by the Integrated Ground Water Information Service [38]. In the survey, the region's annual precipitation is described as the average value for the last 30 years from the investigation. The map shows that the spatial distribution of annual precipitation varies by region. The highest precipitation occurred on the southern coast, whereas the dry region had annual precipitation of 1100 mm or less and was formed around the middle-eastern area. The recharge values and recharge rate follow a pattern similar to the spatial distribution of precipitation by region. The white color represents the areas that have not been investigated yet. and recharge rate follow a pattern similar to the spatial distribution of precipitation by region. The white color represents the areas that have not been investigated yet.     and recharge rate follow a pattern similar to the spatial distribution of precipitation by region. The white color represents the areas that have not been investigated yet.   Chang et al. [35] conducted sensitivity tests by varying the hypothetical infiltration value from 0.1 to 1 m/year. This study used a range of infiltration that reflected actual soil conditions in Korea. We explored the distribution of the concentration at water level, C w , by varying the value within the range of 96-384 mm for the annual infiltration rate that was obtained from the statistical analysis shown in Figure 4, 1-10 m for the distance of the groundwater table from the source, and 1-5 m for the thickness of the contaminant source.

Results and Discussion
3.1. Estimation of the Probability Distribution of DAF unsat Figure 5 shows maximum concentrations and their arrival times at the water table along the unsaturated path when finite-difference one-dimensional transport was simulated. The spread of data points is illustrated when we conducted multiple calculations by varying the value within the range of 96-384 mm for the annual infiltration rate, 1-10 m for traveling distance, and 1-5 m for the source thickness. The rest of the input parameters used the default values shown in Table 1. This C W /C L at y-axis was also designated as the leachate attenuation factor by Verginelli and Baciocchi [24]. The results show that the maximum concentration is inversely proportional to the arrival time, given the combination of all input variables.
Water 2021, 13, x FOR PEER REVIEW 9 conditions in Korea. We explored the distribution of the concentration at water level by varying the value within the range of 96-384 mm for the annual infiltration rate was obtained from the statistical analysis shown in Figure 4, 1-10 m for the distance o groundwater table from the source, and 1-5 m for the thickness of the contaminant sou Figure 5 shows maximum concentrations and their arrival times at the water t along the unsaturated path when finite-difference one-dimensional transport was si lated. The spread of data points is illustrated when we conducted multiple calculat by varying the value within the range of 96-384 mm for the annual infiltration rate, 1 m for traveling distance, and 1-5 m for the source thickness. The rest of the input par eters used the default values shown in Table 1. This CW/CL at y-axis was also design as the leachate attenuation factor by Verginelli and Baciocchi [24]. The results show the maximum concentration is inversely proportional to the arrival time, given the c bination of all input variables. To express the reduction in contaminant concentration, DAFunsat has been adopte an indicator from the Ryu [25] defined as the ratio of the concentration of the source to concentration at the water table. Since we deal with time-varying concentrations in study, we redefined DAFunsat as the ratio of the initial source concentration to the hig concentration at the water level. Therefore, the y-axis of Figure 5 is equal to 1/DAFuns Figure 6 shows that the estimated range of concentrations became relatively nar and converged to a specific value, indicating that the prediction for the possible con tration at the water level may be indicated as a function of the thickness of the source the distance to the water table from the source. All DAFunsat values for all cases w greater than 1, indicating that the leaching concentration from the source was diluted attenuated. To express the reduction in contaminant concentration, DAF unsat has been adopted as an indicator from the Ryu [25] defined as the ratio of the concentration of the source to the concentration at the water table. Since we deal with time-varying concentrations in this study, we redefined DAF unsat as the ratio of the initial source concentration to the highest concentration at the water level. Therefore, the y-axis of Figure 5 is equal to 1/DAF unsat . Figure 6 shows that the estimated range of concentrations became relatively narrow and converged to a specific value, indicating that the prediction for the possible concentration at the water level may be indicated as a function of the thickness of the source and the distance to the water table from the source. All DAF unsat values for all cases were greater than 1, indicating that the leaching concentration from the source was diluted and attenuated. The distribution graph at the top left of Figure 6 is the simulated result for the case of the distance of the water table from the source of 1 m with a source thickness of 1 m. The DAF ranged between 1.6 and 2, meaning that the water level concentration is diluted and attenuated 1.6 to 2 times compared to the soil exposure concentration. Each plot pattern shows a gradual decrease in concentration as the contaminant travel distance increases to 2, 5, and 10 m. In addition, they show a gradual increase in concentration as the source depth increases to 2 and 5 m. When a contaminant travels to the longest path length of 10 m with a source depth of 1 m, the simulated DAFunsat ranged from 3 to 4. The highest DAFunsat simulation results are close to 1 when contaminants travel to a source depth of 5 m and the shortest path length of 1 m.

Case Study: Site Description and Simulation of the Benzene Pathway Using Hydrological Assessment
This study introduced the time-dependent varying infiltration rate to consider sitespecific infiltration. The time-series pattern of site-specific infiltration rate was determined using the groundwater content of the area, and the method for calculating source deple- The distribution graph at the top left of Figure 6 is the simulated result for the case of the distance of the water table from the source of 1 m with a source thickness of 1 m. The DAF ranged between 1.6 and 2, meaning that the water level concentration is diluted and attenuated 1.6 to 2 times compared to the soil exposure concentration. Each plot pattern shows a gradual decrease in concentration as the contaminant travel distance increases to 2, 5, and 10 m. In addition, they show a gradual increase in concentration as the source depth increases to 2 and 5 m. When a contaminant travels to the longest path length of 10 m with a source depth of 1 m, the simulated DAF unsat ranged from 3 to 4. The highest DAF unsat simulation results are close to 1 when contaminants travel to a source depth of 5 m and the shortest path length of 1 m.

Case Study: Site Description and Simulation of the Benzene Pathway Using Hydrological Assessment
This study introduced the time-dependent varying infiltration rate to consider sitespecific infiltration. The time-series pattern of site-specific infiltration rate was determined using the groundwater content of the area, and the method for calculating source depletion was explored. The study adopted the infiltration data from the previous study that investigated a test bed of Cheongju area, where the groundwater cultivation area was precisely calculated. Figure 7a shows the monthly temporal variation of precipitation and infiltration assessed by Chang and Chung [43] from January 2011 to January 2014 using a hydrological analysis tool, the SWAT model [42] developed by the USDA Agricultural Research Service (ARS), to simulate the groundwater stasis state. As input data for the SWAT model, meteorological data such as precipitation, temperature, wind speed, insolation, and relative humidity were used at the Cheongju Meteorological Observatory located in the Musimcheon Basin. The average amount of groundwater recharge in this area was estimated to be ≈13.9% of precipitation. lation, and relative humidity were used at the Cheongju Meteorological Obse cated in the Musimcheon Basin. The average amount of groundwater recharge i was estimated to be ≈13.9% of precipitation. Figure 7b shows that precipitation from the Representative Concentration 4.5 (RCP 4.5) scenario, which is based on the Intergovernmental Panel on Clima (IPCC)'s Fifth Assessment Report (AR5), was applied as an 80-years future pr scenario from 2021 to 2100 (http://www.climate.go.kr/home/CCS/contents_new eapoint_basic.php, accessed on 17 September 2021). Figure 7b also shows that effective infiltration was calculated by multiplying an 80-years RCP 4.5 precipi nario by the recharge rate, 13.9%. The impacts of climate change are not discuss are beyond the scope of this study. To examine the temporal change in concentrations at the exposure point th pathways, Equations (4) and (5) were applied, assuming that the initial leacha tration, (0), is below the saturation concentration or solubility of the con Based on Equation (4), we expected that the leachate concentration may expone crease over time. Therefore, the four periods were regarded as representative ment points of the decay curve considering the log scale. Figure 8 shows the model predicted results that simulated the plume of th inant in the saturation zone after reaching the water table through the unsatur The blue color indicates the concentration of contaminants as 0 and is shown in the relative concentration is 1. The graph shows the concentration distribution urated zone after 5, 15, and 30 years of simulation. We observed that the plum the right in the direction of the groundwater flow and simultaneously spreads The plume is in the form of a curved surface. It also shows a concentration le similar to the inflow concentration at the upper right boundary. In contrast, solute concentration is present at the left and bottom sections, where the gro   Figure 7b also shows that the future effective infiltration was calculated by multiplying an 80-years RCP 4.5 precipitation scenario by the recharge rate, 13.9%. The impacts of climate change are not discussed, as they are beyond the scope of this study.
To examine the temporal change in concentrations at the exposure point through the pathways, Equations (4) and (5) were applied, assuming that the initial leachate concentration, C L (0), is below the saturation concentration or solubility of the contaminant. Based on Equation (4), we expected that the leachate concentration may exponentially decrease over time. Therefore, the four periods were regarded as representative measurement points of the decay curve considering the log scale. Figure 8 shows the model predicted results that simulated the plume of the contaminant in the saturation zone after reaching the water table through the unsaturated zone. The blue color indicates the concentration of contaminants as 0 and is shown in red when the relative concentration is 1. The graph shows the concentration distribution in the saturated zone after 5, 15, and 30 years of simulation. We observed that the plume moves to the right in the direction of the groundwater flow and simultaneously spreads vertically. The plume is in the form of a curved surface. It also shows a concentration level almost similar to the inflow concentration at the upper right boundary. In contrast, almost no solute concentration is present at the left and bottom sections, where the groundwater flows through the hydraulic gradient. The color near the bottom was blue, indicating that the concentration was maintained at a remarkably low level. Figure 8b shows a relatively smaller concentration when a higher depth of source and a longer distance of unsaturated zone are applied.   Figure 9 shows the temporal profiles of source concentration, the concentration at the water table, and exposure point concentrations at different depths. The graph shows that the breakthrough curves are not smooth because of the seasonally varied infiltration that vertically flows toward the groundwater level. The top layer (layer 1) is most affected by precipitation fluctuation, and the curves get smoother approaching the bottom of the model. Source depletion decreases exponentially, and the concentration of the water table reaches its highest concentration about 10 years after leaching from the source. Contaminants introduced into saturated aquifer show a lower concentration while mixing with the existing freshwater, and the lower it goes, the lower the concentration. After 16 years, the concentration of layer 2 exceeds that of layer 1. After 20 years, layer 3 will exceed the concentrations of the upper layer. As observed in the plume movement of the model crosscut shown in Figure 8, the contaminant gradually moves down, and the recharge was introduced on the top layer at a low concentration after a certain time has elapsed. However, the lower layers after layer 6 show a gradually decreasing pattern as a whole because they are areas not affected by the contaminant's travel path. Figure 9b is the expected concentration when the mass of the initial source was assumed to be 1 mg/kg. At the end of the simulation, the bottom indicates a relatively high concentration because benzene concentration from recharge through the vadose zone decreases over time. Figure 9 shows the temporal profiles of source concentration, the concentration at the water table, and exposure point concentrations at different depths. The graph shows that the breakthrough curves are not smooth because of the seasonally varied infiltration that vertically flows toward the groundwater level. The top layer (layer 1) is most affected by precipitation fluctuation, and the curves get smoother approaching the bottom of the model. Source depletion decreases exponentially, and the concentration of the water table reaches its highest concentration about 10 years after leaching from the source. Contaminants introduced into saturated aquifer show a lower concentration while mixing with the existing freshwater, and the lower it goes, the lower the concentration. After 16 years, the concentration of layer 2 exceeds that of layer 1. After 20 years, layer 3 will exceed the concentrations of the upper layer. As observed in the plume movement of the model crosscut shown in Figure 8, the contaminant gradually moves down, and the recharge was introduced on the top layer at a low concentration after a certain time has elapsed. However, the lower layers after layer 6 show a gradually decreasing pattern as a whole because they are areas not affected by the contaminant's travel path. Figure 9b is the expected concentration when the mass of the initial source was assumed to be 1 mg/kg. Figure 10 showed that the delayed time and lower concentration after the contaminant traveled a longer distance from the top of the source to the water table in the unsaturated zone. cut shown in Figure 8, the contaminant gradually moves down, and the recharge was introduced on the top layer at a low concentration after a certain time has elapsed. However, the lower layers after layer 6 show a gradually decreasing pattern as a whole because they are areas not affected by the contaminant's travel path. Figure 9b is the expected concentration when the mass of the initial source was assumed to be 1 mg/kg.

Conclusions
Risk-based management and assessment of contaminated sites is a way to evaluate the potential hazards of contaminants based on considering linkages between pollution sources, pathways, and receptors (i.e., human body or ecosystem) though measuring the probability of adverse events [44]. In addition, it makes the decisions for the cleanup of contaminated soil and sediment.
Sophisticated computational methods and site-specific data are important for increasing the reliability of risk assessment results and reducing uncertainty. Site-specific databases can be used as a basis for new site development, enabling mid to long-term risk prediction using existing remediation strategies without large-scale investments, providing a quantitative evaluation tool for long-term monitoring management.
This study emphasizes the importance of on-site investigation and site-specific risk assessment of contaminated sites. Risk-based analysis of numerical simulations of contaminants through the soil to groundwater determined the DAFs based on the geometry of contaminated sites. Conservative models, proposed in the ASTM-RBCA framework, assume steady-state leaching in the unsaturated zone, a constant leaching rate without reduction of the source, and mixing based on a mass-balance approach underneath the contaminated sites. The method established in this study considered temporal changes for

Conclusions
Risk-based management and assessment of contaminated sites is a way to evaluate the potential hazards of contaminants based on considering linkages between pollution sources, pathways, and receptors (i.e., human body or ecosystem) though measuring the probability of adverse events [44]. In addition, it makes the decisions for the cleanup of contaminated soil and sediment.
Sophisticated computational methods and site-specific data are important for increasing the reliability of risk assessment results and reducing uncertainty. Site-specific databases can be used as a basis for new site development, enabling mid to long-term risk prediction using existing remediation strategies without large-scale investments, providing a quantitative evaluation tool for long-term monitoring management.
This study emphasizes the importance of on-site investigation and site-specific risk assessment of contaminated sites. Risk-based analysis of numerical simulations of contaminants through the soil to groundwater determined the DAFs based on the geometry of contaminated sites. Conservative models, proposed in the ASTM-RBCA framework, assume steady-state leaching in the unsaturated zone, a constant leaching rate without reduction of the source, and mixing based on a mass-balance approach underneath the contaminated sites. The method established in this study considered temporal changes for leaching in the unsaturated zone, source depletion, and numerical model simulation with sophisticated boundary conditions underneath the contaminated sites. We considered the data that estimated 20 years recharge by SWAT in a sub-basin in Korea as an infiltration input of the numerical model. In addition, the reported recharge rate data collected from the national survey for the entire Korean peninsula were selected for the risk-based estimation of groundwater flow and transport. This approach enabled the narrowing of the possible range of uncertainty for the predicted spatial and temporal contaminant concentration profiles at the groundwater level. This study identified the approximate distribution of the infiltration rate in Korea. The case study also confirmed that this computing method we applied could be suitable for regional-scale investigation.
The current guidelines for risk assessment are sometimes believed to overestimate the risk of contaminants due to the conservative and safe assumptions adopted. However, as observed in the simulation results, different from the conservative assumption to be entirely mixed and diluted in the aquifer, the expected concentration from the saturated models showed time-dependent varying concentration at each depth.
This study illustrated the application of numerical analysis considering high-tiered risk-based assessment. The proposed application is simple, and only the distribution of the Korean infiltration values was identified. The distribution characteristics for other indicators should be identified. In particular, the distribution characteristics of the source depth and the distance between the groundwater level and the water source need to be further investigated in future studies. Additionally, the approximate theoretical results derived from this study should be validated against extensive field data with other organic contaminants. These proposed additional studies, as USEPA used algorithms to implement Monte Carlo simulations, can determine DAFs for unsaturated regions in Korea based on the development of appropriate and reasonable statistical analysis methods.