## 4 Results and Discussion  
### Hydrogeological Model  
The top soil layer of the treatment area was excavated on site to form a 2.3 m deep pit, and the excavated portion is shown in the left solid red line area in Figure 1. Based on the geologic profile information obtained from the on-site exploratory wells, layering can be performed using COMSOL's parametric surface function. Based on the engineering geological units of the foundation bedrock, the study area, from a depth of 2.3 m to 30 m below the pit, is divided into five layers. The stratigraphy is summarized in Figure 2. From top to bottom, these layers are 1--silty clay, 2--rounded gravel, 3--clay, 4--rounded gravel, and 5--aquifier with relative confining properties. Layers 1 to 4 constitute the aquifer, while the remaining layers below are conceptualized as aquifers with relative confining properties. The dashed circles in Figure 2 represent well locations, with Well 7 and Well 8 serving as injection wells, while the others are extraction wells. The wells have a diameter of 0.2 m and a depth of 14 m.  
In this study, we employed a method involving homogeneous extraction wells and injection wells for the remediation of various categories of pollutants present in the soil and soil pores. Pollution elimination was achieved by injecting water into pressure-inlet wells and pumping water from eight surrounding wells. We utilized Darcy's law interface, with the central two wells designated as pressure inlets, maintaining a water column pressure of 5 m at the wellheads. We established an analytical function using COMSOL to account for boundary pressures in wells of varying heights. Non-flow boundary conditions were set to simulate impermeable membranes and low-permeability rocks to prevent flow to the surroundings. The top interface represented the boundary between soil and water, and as there was no groundwater seeping to the surface, it was also set as a non-flow boundary. Initial conditions included the initialization of pore pressures, where locations with zero depth were considered to be in contact with the atmosphere, thus having pore pressure equal to atmospheric pressure. At different depths, pore pressure was approximately equal to static water pressure, and we used analytical functions to set these initial pressure values. Since the testing period occurred during winter with lower precipitation and a lower groundwater Table, we set the soils in contact with the atmosphere in the excavation pit to have 95% variable saturation flow. Thus, the top two layers are described using [PERSON]'s equation, while the deeper ones are below the water level and are described using [PERSON]'s law. Based on on-site testing data, we adjusted the porosity and permeability of different soil layers to account for seasonal variations and soil properties. The porosity and hydraulic conductivity of each layer of the site were provided by the Zhejiang Provincial Exploration Institute through exploratory wells, taking cores, and indoor experiments. For more detailed information, please refer to Table 1.  
For the setting of initial values, an analysis was conducted summarizing the types and concentrations of pollutants measured in the on-site soil. It was found that the majority of pollutants consisted of volatile organic compound s(VOC) and semi-volatile organic compounds (SVOC), as well as extractable petroleum hydrocarbons (EPH). By organizing this information, the distribution of pollutants with different diffusion properties in the well field coordinate system was obtained. Interpolation functions for pollutant initial values were defined using COMSOL. The primary substance exchange in the process of pollutant removal occurs between groundwater and subsurface soil, as well as the diffusion of pollutants between ground