# Visualization Method for Porous Groundwater Seepage Flow Field Based on Particle Flow: Case of Yancheng City in the East Coast of China

Authors: Liang He; Ling Chen; Suozhong Chen
Journal: Geofluids
 (2022)

---

## Abstract

Porous groundwater seepage flow field reflects the spatiotemporal dynamic characteristics of porous groundwater movement. The geometric point and vector line methods are the traditionally used visualization methods for groundwater seepage flow field; however, they exhibit difficulty in describing the characteristics and attributes of a seepage flow field. Therefore, a new particle flow method in the field of computer graphics is introduced in this study to visualize seepage flow field, including the storage conditions, transport mechanism, and seepage field properties of porous groundwater. This visualization method is discussed from four aspects: the generation of seed points, the life cycle of particles, the movement velocity of particles, and the formation of particle trajectory and smooth processing. The seepage field of confined aquifer III groundwater expressed via particle flow and the contour of critical groundwater level are superimposed and analyzed in Yancheng City in the east coast of China. Results show that particle flow can be effectively applied to the visualization of groundwater seepage flow field and provide spatial auxiliary decision-making support for the sustainable exploitation of groundwater and the formulation of geological environment protection measures.


---

## Body

## 1. Introduction

Groundwater seepage flow field is a mixed field composed of groundwater head scalar field and groundwater flow velocity field [1, 2]. In practical applications, groundwater seepage flow field is frequently obtained from the dynamic monitoring data of groundwater level or the numerical simulation of groundwater flow [3]. The precision of the numerical simulation of a groundwater flow field is extremely important to the research of hydrogeology as it directly affects the calculation precision of solute transport, multiphase flow, and the visualization of groundwater flow field [4–6]. In addition, the visualization of groundwater flow field is an important research content in the field of hydrogeology that shows the flow law of groundwater in front of people [7]. The visualization of groundwater flow field can be read by professional researchers, engineering practitioners, and even nonprofessionals and provide a scientific basis for related research fields and practical applications [8].Flow visualization can be divided into four categories in accordance with visual expression: geometric visualization [9, 10], volume visualization [11, 12], texture-based visualization [13], and feature-based visualization [14]. Geometric shape visualization is relatively simple, intuitive, and easy to understand; thus, it is the most widely used visualization method in engineering [15]. Geometric shape visualization can display and express scalar and vector field data. Scalar field data are divided into isoline and isosurface in accordance with different dimensions. Methods for visualizing the geometric shape of vector field data include the geometric point, time series point [16], vector line [17], dynamic graphical [18], vector surface, and vector volume methods [19].Geometric point and vector line methods are the most widely used methods for visualizing groundwater seepage flow field. The two methods have the simplest data mapping process with fast rendering speed, and they can better reveal the geometry of groundwater seepage flow field than other methods [20]. However, both methods are not ideal for describing the dynamics of groundwater seepage flow field; moreover, most of the given methods are static or false dynamic simulation under an instantaneous state [21]. In summary, describing the characteristics and dynamic properties of a flow field is difficult because of the superposition of successive pictures at adjacent moments [22]. As one of the effective methods for simulating irregular dynamic objects or scenes, the particle flow method is an important component of computer graphics. Particle systems are widely used in virtual reality, 3D simulation, game development, special film effects, visualization, and other fields [23, 24]. Particle systems are also used to simulate the dynamic effects of waves, flames, and smoke; the special effects of light, shadow, and lightning; the falling process of rain and snow; and the flowing process of liquid or gas [25–27]. Additionally, a particle system is an indispensable modeling tool for constructing 3D scenes [28]. The basic idea of this method is to gather as many simple particles as basic elements to form an irregular object and then a closed system [29]. This method fully reflects the dynamic, irregular, and random nature of fuzzy scenery. Any object can be regarded as composed of the simplest particles regardless of whether it is in the solid, liquid, or gaseous state [30]. In addition, particle flow is an indispensable modeling tool for constructing 3D scenes [31, 32].Moreover, flow field is difficult to visualize directly because of the complexity of groundwater storage conditions and the space–time dynamics of groundwater movement [33]. However, particle flow is a dynamic simulation of a closed system formed by a large number of particles with attributes. In addition, particle flow has a large number of properties, and researchers can modify and set these properties through effective human–computer interaction [34, 35]. Therefore, the dynamic visual simulation of groundwater seepage flow field based on the particle flow method is feasible and has certain theoretical importance and application value.


## 2. Materials and Methods

### 2.1. Study Area

The study was conducted in Yancheng City, which is located in the east coast of China and the central part of the Northern Jiangsu Plain. Its geographical coordinates are 32° 34′–34° 28′ N and 119° 27′–120° 54′ E, and its total land area is about 17,000 km2. Figures 1(a)–1(c) show the geographical location and hydrogeological profile of the study area. This area is densely covered by a water network that belongs to the coastal water network plain landform type [36]. The terrain in the area is flat, and it tilts slowly from southeast to northwest. At present, the drinking water source of residents in the rural areas of Yancheng City nearly depends entirely on groundwater. A certain scale of groundwater depression funnel has been formed in many key cities and towns, such as Yandu District and Dafeng City, because of the excessive exploitation of groundwater with the rapid development of the economy. The settlement of Dafeng City is the most serious, with an accumulative settlement above 600 mm, while that of Yancheng City exceeds 400 mm. Binhai County has a land subsidence rate of above 10 mm/a in the past 5 years, and it is facing serious ecological and environmental problems.Figure 1
Map of the study area and its geographical location. (a) Jiangsu Province, east coast of China; (b) Yancheng City; (c) spatial location of wells and hydrogeologic drill holes and (d) hydrogeological section.The major types of groundwater in the study area are porous phreatic and confined water, and the depth of the confined aquifer is 50–350 m. The lithology of the aquifer is mostly silt, fine sand, coarse sand, and gravel; water content is abundant and stable. Clay or silt clay layers between different aquifers form aquicludes [37]. Confined aquifers I, II, and III are distributed in this area (Figure 1(d)). Confined aquifer III is composed of lower Pleistocene medium-fine and medium-coarse sand with a thickness of 20.00–35.00 m, and its roof elevation ranges from −144.00 m to −294.60 m, indicating a good water yield. Confined aquifer III is the primary exploitation layer of groundwater. Therefore, confined aquifer III was selected as the research object in the current study.


### 2.2. Sample Collection

The study area has 133 groundwater dynamic monitoring wells, which are distributed in Yancheng County, Funing County, Dongtai City, and other areas. Among these wells, 55 are confined aquifer III, and the frequency of groundwater level dynamic monitoring is two times a month. The current study selected 10 years (2005–2014) of groundwater level monitoring data. The monitoring data of four time points, namely, early January, early April, early July, and early October, were chosen for analysis in view of the small change range of groundwater level within a short period and the close monitoring data of adjacent time intervals. The groundwater level changes of eight representative monitoring wells in the past 10 years are shown in Figure2. Overall, the average level of confined groundwater in the study area exhibited a downward trend from 2005 to 2009, and the decline in groundwater level tended to be stable from 2009 to 2014.Figure 2
Main monitoring well groundwater level changes.


### 2.3. Particle Flow Field

#### 2.3.1. Property Analysis of Particle Flow Field

The initial properties of the particle flow of porous groundwater seepage flow field include the number of initial particles, initial position, particle size, initial velocity and direction of movement, initial color, initial transparency, initial shape, and life cycle. The most important step in generating a flow field is to place tracking particles. The initial position and number of particles directly affect the final visualization quality. Improper placement or low number of particles can lead to failure in tracking important groundwater flow features. Meanwhile, too many particles can also block the line of sight and affect the depiction of major features. Therefore, appropriate initial position and initial number of particles are the keys to flow field visualization. Many researchers have proposed optimal solutions, such as particle placement guidance strategy, spatial average distribution, preferential streamline generation, double streamline, and entropy guidance, to solve the problem of particle placement in a flow field. In the current study, the initial position of particles in groundwater seepage flow field was at the grid node, and the spatial average distribution method was used to optimize the placement of streamline particles.Particle flow was constructed in the form of “point particles,” which can be approximately regarded as pixels with a certain size. Texture mapping technology was used to map the “point particles” to form visual fluid particles. The velocity of particles in porous groundwater seepage flow field has two attributes: size and direction. The contour of groundwater level was obtained by interpolating existing data of groundwater monitoring wells. The velocity magnitude and direction of a flow field at the grid nodes and the initial velocity of fluid particles can be obtained in accordance with Darcy’s law [38]. The life cycle of porous groundwater seepage flow field particles is reflected in the time scale of fluid particles appearing in a given display area, which is typically represented by number of frames. Life cycle determines the forward distance of fluid particles when finding the sample velocity.


#### 2.3.2. Calculation of Particle Flow Trajectory

Groundwater flow velocity (seepage velocity) refers to a type of simulated velocity that occurs when water flows through a section with skeletons and pores [38]. In 1856, French hydraulic scientist Henry Darcy found that porous groundwater flows from high to low groundwater level and concluded that
(1)q=KJ=−K∇h,where q is the seepage velocity and K (unit: meter/day) is the conductivity coefficient of the aquifer, which is a measure of an aquifer’s ability to transmit groundwater. The K conductivity coefficient is the core characteristic data of a groundwater environment. J is the hydraulic gradient that refers to the average head loss per unit length of a flow path. ∇h is the gradient of the function h; the direction points to the direction of water head reduction along the normal line of the water head isoline. In three dimensions, the seepage velocity at any point (x, y, and z) in the seepage field can be written as
(2)q=ui⟶+vj⟶+wk⟶.The component of a hydraulic gradientJx, Jy, and Jz in directions x, y, and z can be expressed as
(3)J=−∇h=−dhds,Jx=dhdx,Jy=dhdy,Jz=dhdz.A connection exists between the gradient at any point on the surface and the isoline that passes through that point [39]. The direction of the gradient of the functionz=fx,yatPx,yis the same as the normal direction of the isolinefx,y=Cpassing throughPx,yat this point following the law from the low value to the high value of the isoline. In addition, the module of the gradient is the derivative of the function in the normal direction. The relationship between the gradient at a point on the surface and the isoline through this point is illustrated in Figure 3.Figure 3
The relationship between the gradient at a point on the surface and the isoline through this point.The surface isoline gradient value was obtained, and then, the hydraulic gradient value was determined in accordance with the hydraulic gradient and the relationship between gradients. The method is described as follows.The variation characteristics of the equivalent surface (curved surface) of groundwater level in 3D space are projected onto a plane in 2D space. The (x, y) coordinate of a 2D plane is regarded as the independent variable, and the groundwater level H is considered the dependent variable [21]. The function H=fx,y reflects the regional distribution of groundwater level, and the statistical velocity must involve at least three groundwater level data. Thus, the local trend surface method was adopted to fit groundwater level data in the current study to improve the fitting degree of the trend surface and avoid the situation in which the number of the present groundwater level points near the calculated point is on the low side. In this work, a regular grid that can cover all the study areas was constructed and then combined with the groundwater level of monitoring wells. Kriging interpolation was adopted to calculate groundwater level at the grid intersection point in the study area.The groundwater level of grid points and original monitoring wells obtained via Kriging interpolation statistics was used to establish the trend surface of groundwater level dynamic change [40]. If the point is Px,y, then the radius R is defined and a circle with Px,y as the center and R as the radius is drawn. All the sampling points covered in the circle were used to establish groundwater level trend surface. Considering the global trend and fitting effect of Px,y, the trend surface model with the time of two was selected for the statistics of groundwater level surface. At this moment, the gradient vector of the point on the trend surface is
(4)gradfx,y=∂H∂x,∂H∂y=−a1+2a3x+a4y,a2+a4x+2a5y.Groundwater flows from the higher point to the lower point of mechanical energy; thus, the direction of hydraulic slope and gradient vector is opposite.(5)J=−a1−2a3x−a4y,−a2−a4x−2a5y.Then, the velocity of the pointPx,y is calculated in accordance with Equation (1).The advancing front technique (AFT) algorithm was used on the basis of the boundary data of the study area [41]. The AFT algorithm must be adjacent to each other to generate a regular grid for expressing the distribution of the water head through the monitoring well node objects with attached groundwater level within the scope of the study area. The groundwater level data provided by the dynamic monitoring wells for groundwater level were inserted into the calculation grid as constraint points.A groundwater dynamic monitoring well provides the actual monitoring groundwater level data in discrete space, and groundwater level is continuously distributed in space. Therefore, the spatial interpolation method is required to obtain the groundwater level of each calculation grid node at the initial time to form the initial groundwater level field. Then, this value is distributed to the discrete grid nodes to complete the assignment of the initial groundwater level [42]. The piezometric head value hi of each grid node position (xi and yi ) was extracted, and the digital elevation model of the water head was constructed by xi, yi, and hi to present the spatial distribution characteristics of aquifer groundwater level. We can obtain the movement position of the fluid proton after unit time step in accordance with the calculated velocity and direction of the fluid proton. We can also obtain the change rules for the movement position of groundwater fluid and then the dynamic simulation of porous groundwater seepage flow field by repeating the calculation process.


#### 2.3.3. Particle Flow Field Visualization

The visualization simulation of porous groundwater seepage flow field based on the particle flow method must consider the real-time effect and authenticity of flow field simulation [43]. The real-time effect of the simulation will become better, but its authenticity will become worse with a reduction in the overall number of particles. By contrast, the real-time effect of the simulation will become worse, but its authenticity will become better with an increase in the overall number of particles. Therefore, the overall number of particles in particle flow should be determined in accordance with the actual simulation situation. In addition, the choice of particle model will also influence the real-time effect and authenticity of the visual expression. The “point particle” graph element structure was selected to simulate groundwater fluid particles, while the graph element structure was used as the drawing unit of particle flow. The texture image of a “water droplet” was mapped onto the point particle rendering unit of particle flow by using the visualization toolkit (VTK) technology to increase the realistic effect of the simulation. This section focuses on the following issues: the generation of seed points, the life cycle of particles, and the formation and smoothing of particle trajectories.(1) Generation of Seed Points. The two methods for generating and distributing the seed points of particle flow are the average new particle number and object area methods [44]. The groundwater level of grid nodes can be obtained via interpolation through the groundwater level of monitoring wells. The generation distribution principle of the object area measurement method was adopted, and the seed points of particle flow were generated at the grid nodes to obtain the velocity of the flow field at the initial time. In this method, the specified average value and maximum variation range are the average value and variance of the generated particles per unit area of the screen. Thus, the number of generated particles depends on the size of the screen occupied by an object. This method can produce different degrees of details and better simulate groundwater seepage flow field than other methods.(2) Life Cycle of Particles. In accordance with the property analysis of groundwater flow particles, the process of particles from production to extinction exhibits a life cycle. The life cycle of a particle determines how long a particle exists in the simulation. In general, the number of frames is used to describe the life cycle of a particle. Therefore, the life cycle of a particle in the i frame can be expressed by the following equation:
(6)Finallifetimet=initiallifetime−t×attenuationperframe,where the initiallifetime is the initial life cycle of a particle, t is the time since the particle was generated, and attenuationperframe is the life decay rate of particles per frame, indicating the number of lives lost by particles per frame. The advantage of using number of frames as the life cycle unit of particles is that it is simple and easy to implement. The disadvantage is that the complexity in generating graphics will directly affect the life cycle and result in the uncertainty of time units. Corresponding new particles are generated at the location of particle extinction to maintain the balance of the total number of particles in the flow field simulation process.(3) Formation and Smoothing of Particle Trajectories. On the basis of the calculation method for particle flow trajectory introduced in Section 2.3.2, the position and velocity of particles at each time frame can be calculated in real time in accordance with existing data conditions. Thus, the movement trajectory of fluid particles in their life cycle was calculated in accordance with the hydrodynamic conditions, and the flow simulation of porous groundwater particle flow field was realized [45]. The real effect of particle flow trajectory is enhanced to the reality of a particle flow field. In the current study, the trajectory from the head position to the tail position was set, and the width of the trajectory was gradually reduced to form the “trailing” effect of particle flow. The life cycle of particles will end when the velocity of fluid particles reaches zero or the number of times of movement approaches five. The transparency of a particle is reduced when it is about to die to enhance the gradual effect of dying particles and reduce the sudden feeling of particle death.


#### 2.3.4. Groundwater Level Control Management

Groundwater level and its change trend are important criteria for measuring whether the development of groundwater resources is reasonable [46]. The controlling critical level of groundwater is the critical groundwater level that is determined to prevent various problems or disasters. In accordance with hydrogeological conditions, ecological functions, and geological environment functions caused by the change in groundwater level in different hydrogeological zones [47], the controlling critical level of groundwater is a set of static critical thresholds. The controlling critical level of groundwater can be divided into upper and lower critical groundwater levels. The upper and lower limits of critical groundwater level refer to the maximum and minimum groundwater levels, respectively, that do not cause various problems or disasters. The upper critical groundwater level is caused by excessive groundwater supply or insufficient exploitation. It is the critical groundwater level for preventing soil secondary salinization, road frost boiling, and the antifloating and waterproof design of underground structures, such as subways. The lower critical groundwater level is caused by insufficient groundwater recharge or overexploitation. It is the critical groundwater level used to prevent land subsidence, ground fissures, natural vegetation, and seawater intrusion.The visualization of groundwater seepage flow field was conducted on the basis of particle flow and the contour of groundwater level. On the one hand, the isosurface of groundwater level can reflect its spatial distribution. On the other hand, particle flow direction can indicate the movement trend of groundwater and the recharge relationship among groundwater aquifer group [48]. The specific application methods for groundwater particle flow visualization in the management of the critical groundwater level are as follows.
(1)
Visualization analysis method for groundwater’s controlling critical level (upper and lower limits of critical groundwater level) [49]: the critical groundwater level of an aquifer is consistent within the same critical groundwater level management range. Therefore, the critical water table surface can also be regarded as a constant groundwater level surface. The visual analysis of the study range was performed using the method of generating groundwater level isosurface by considering the geometric characteristics of the controlling critical groundwater level. In addition, the range of groundwater level, which is higher than the lower critical threshold and lower than the upper critical threshold, was rendered in orange by using the spatial data symbolization technology of a geographic information system (GIS) [50]. The orange area indicates the restricted exploitation site in this region. The range of groundwater level that is lower than the lower critical threshold value was rendered in red, and the red area represents the prohibited mining site. The remaining regions are reasonable mining areas. On the basis of this method, controlled critical groundwater level analysis can be performed for groundwater level in the study area(2)
Method for quantifying and visualizing the analytical results of critical groundwater level: in accordance with the analysis results of the controlling critical groundwater level in the study area, we can determine which specific area of the groundwater level is lower than the controlling critical groundwater level and the geographical location of these areas. The spatial topological analysis technique of GIS is used to analyze the topological relationship of the contour of critical groundwater level. The isosurface of groundwater level was generated from the groundwater level isoline and filled with color sequence. The contour of groundwater level has two forms: the closed and unclosed states. Thus, achieving optimization is necessary before the topological groundwater level isoline is establishedUnclosed contour: in accordance with the GIS topological data structure, a surface is composed of unclosed line segments because the existence of a boundary will make part of the isoline unable to close. The unclosed isopleth lines should be connected to a closed line by using the boundary line of the study area to obtain the closed isopleth line.Closed contour: the isosurface is constructed through the topological relationship between lines after closing the unclosed isolines. Then, inclusion analysis is performed on the constructed equivalent surface. Surfaces A and B are overlapping, and surface C is deleted if it exhibits an inclusion relationship with surfaces A and B. Surface C is not deleted if it has no containment relationship with the other surfaces. The presence and absence of containment relationships are depicted in Figures 4 and 5.Figure 4
Containment relationship.Figure 5
No containment relationship.By analogy, all the isosurfaces in the study area are traversed. Finally, the size of the processed face was calculated, and the area that was not “0” was rendered with color. The color sequence from blue to red indicates that the groundwater level is getting deeper.


## 2.1. Study Area

The study was conducted in Yancheng City, which is located in the east coast of China and the central part of the Northern Jiangsu Plain. Its geographical coordinates are 32° 34′–34° 28′ N and 119° 27′–120° 54′ E, and its total land area is about 17,000 km2. Figures 1(a)–1(c) show the geographical location and hydrogeological profile of the study area. This area is densely covered by a water network that belongs to the coastal water network plain landform type [36]. The terrain in the area is flat, and it tilts slowly from southeast to northwest. At present, the drinking water source of residents in the rural areas of Yancheng City nearly depends entirely on groundwater. A certain scale of groundwater depression funnel has been formed in many key cities and towns, such as Yandu District and Dafeng City, because of the excessive exploitation of groundwater with the rapid development of the economy. The settlement of Dafeng City is the most serious, with an accumulative settlement above 600 mm, while that of Yancheng City exceeds 400 mm. Binhai County has a land subsidence rate of above 10 mm/a in the past 5 years, and it is facing serious ecological and environmental problems.Figure 1
Map of the study area and its geographical location. (a) Jiangsu Province, east coast of China; (b) Yancheng City; (c) spatial location of wells and hydrogeologic drill holes and (d) hydrogeological section.The major types of groundwater in the study area are porous phreatic and confined water, and the depth of the confined aquifer is 50–350 m. The lithology of the aquifer is mostly silt, fine sand, coarse sand, and gravel; water content is abundant and stable. Clay or silt clay layers between different aquifers form aquicludes [37]. Confined aquifers I, II, and III are distributed in this area (Figure 1(d)). Confined aquifer III is composed of lower Pleistocene medium-fine and medium-coarse sand with a thickness of 20.00–35.00 m, and its roof elevation ranges from −144.00 m to −294.60 m, indicating a good water yield. Confined aquifer III is the primary exploitation layer of groundwater. Therefore, confined aquifer III was selected as the research object in the current study.


## 2.2. Sample Collection

The study area has 133 groundwater dynamic monitoring wells, which are distributed in Yancheng County, Funing County, Dongtai City, and other areas. Among these wells, 55 are confined aquifer III, and the frequency of groundwater level dynamic monitoring is two times a month. The current study selected 10 years (2005–2014) of groundwater level monitoring data. The monitoring data of four time points, namely, early January, early April, early July, and early October, were chosen for analysis in view of the small change range of groundwater level within a short period and the close monitoring data of adjacent time intervals. The groundwater level changes of eight representative monitoring wells in the past 10 years are shown in Figure2. Overall, the average level of confined groundwater in the study area exhibited a downward trend from 2005 to 2009, and the decline in groundwater level tended to be stable from 2009 to 2014.Figure 2
Main monitoring well groundwater level changes.


## 2.3. Particle Flow Field

### 2.3.1. Property Analysis of Particle Flow Field

The initial properties of the particle flow of porous groundwater seepage flow field include the number of initial particles, initial position, particle size, initial velocity and direction of movement, initial color, initial transparency, initial shape, and life cycle. The most important step in generating a flow field is to place tracking particles. The initial position and number of particles directly affect the final visualization quality. Improper placement or low number of particles can lead to failure in tracking important groundwater flow features. Meanwhile, too many particles can also block the line of sight and affect the depiction of major features. Therefore, appropriate initial position and initial number of particles are the keys to flow field visualization. Many researchers have proposed optimal solutions, such as particle placement guidance strategy, spatial average distribution, preferential streamline generation, double streamline, and entropy guidance, to solve the problem of particle placement in a flow field. In the current study, the initial position of particles in groundwater seepage flow field was at the grid node, and the spatial average distribution method was used to optimize the placement of streamline particles.Particle flow was constructed in the form of “point particles,” which can be approximately regarded as pixels with a certain size. Texture mapping technology was used to map the “point particles” to form visual fluid particles. The velocity of particles in porous groundwater seepage flow field has two attributes: size and direction. The contour of groundwater level was obtained by interpolating existing data of groundwater monitoring wells. The velocity magnitude and direction of a flow field at the grid nodes and the initial velocity of fluid particles can be obtained in accordance with Darcy’s law [38]. The life cycle of porous groundwater seepage flow field particles is reflected in the time scale of fluid particles appearing in a given display area, which is typically represented by number of frames. Life cycle determines the forward distance of fluid particles when finding the sample velocity.


### 2.3.2. Calculation of Particle Flow Trajectory

Groundwater flow velocity (seepage velocity) refers to a type of simulated velocity that occurs when water flows through a section with skeletons and pores [38]. In 1856, French hydraulic scientist Henry Darcy found that porous groundwater flows from high to low groundwater level and concluded that
(1)q=KJ=−K∇h,where q is the seepage velocity and K (unit: meter/day) is the conductivity coefficient of the aquifer, which is a measure of an aquifer’s ability to transmit groundwater. The K conductivity coefficient is the core characteristic data of a groundwater environment. J is the hydraulic gradient that refers to the average head loss per unit length of a flow path. ∇h is the gradient of the function h; the direction points to the direction of water head reduction along the normal line of the water head isoline. In three dimensions, the seepage velocity at any point (x, y, and z) in the seepage field can be written as
(2)q=ui⟶+vj⟶+wk⟶.The component of a hydraulic gradientJx, Jy, and Jz in directions x, y, and z can be expressed as
(3)J=−∇h=−dhds,Jx=dhdx,Jy=dhdy,Jz=dhdz.A connection exists between the gradient at any point on the surface and the isoline that passes through that point [39]. The direction of the gradient of the functionz=fx,yatPx,yis the same as the normal direction of the isolinefx,y=Cpassing throughPx,yat this point following the law from the low value to the high value of the isoline. In addition, the module of the gradient is the derivative of the function in the normal direction. The relationship between the gradient at a point on the surface and the isoline through this point is illustrated in Figure 3.Figure 3
The relationship between the gradient at a point on the surface and the isoline through this point.The surface isoline gradient value was obtained, and then, the hydraulic gradient value was determined in accordance with the hydraulic gradient and the relationship between gradients. The method is described as follows.The variation characteristics of the equivalent surface (curved surface) of groundwater level in 3D space are projected onto a plane in 2D space. The (x, y) coordinate of a 2D plane is regarded as the independent variable, and the groundwater level H is considered the dependent variable [21]. The function H=fx,y reflects the regional distribution of groundwater level, and the statistical velocity must involve at least three groundwater level data. Thus, the local trend surface method was adopted to fit groundwater level data in the current study to improve the fitting degree of the trend surface and avoid the situation in which the number of the present groundwater level points near the calculated point is on the low side. In this work, a regular grid that can cover all the study areas was constructed and then combined with the groundwater level of monitoring wells. Kriging interpolation was adopted to calculate groundwater level at the grid intersection point in the study area.The groundwater level of grid points and original monitoring wells obtained via Kriging interpolation statistics was used to establish the trend surface of groundwater level dynamic change [40]. If the point is Px,y, then the radius R is defined and a circle with Px,y as the center and R as the radius is drawn. All the sampling points covered in the circle were used to establish groundwater level trend surface. Considering the global trend and fitting effect of Px,y, the trend surface model with the time of two was selected for the statistics of groundwater level surface. At this moment, the gradient vector of the point on the trend surface is
(4)gradfx,y=∂H∂x,∂H∂y=−a1+2a3x+a4y,a2+a4x+2a5y.Groundwater flows from the higher point to the lower point of mechanical energy; thus, the direction of hydraulic slope and gradient vector is opposite.(5)J=−a1−2a3x−a4y,−a2−a4x−2a5y.Then, the velocity of the pointPx,y is calculated in accordance with Equation (1).The advancing front technique (AFT) algorithm was used on the basis of the boundary data of the study area [41]. The AFT algorithm must be adjacent to each other to generate a regular grid for expressing the distribution of the water head through the monitoring well node objects with attached groundwater level within the scope of the study area. The groundwater level data provided by the dynamic monitoring wells for groundwater level were inserted into the calculation grid as constraint points.A groundwater dynamic monitoring well provides the actual monitoring groundwater level data in discrete space, and groundwater level is continuously distributed in space. Therefore, the spatial interpolation method is required to obtain the groundwater level of each calculation grid node at the initial time to form the initial groundwater level field. Then, this value is distributed to the discrete grid nodes to complete the assignment of the initial groundwater level [42]. The piezometric head value hi of each grid node position (xi and yi ) was extracted, and the digital elevation model of the water head was constructed by xi, yi, and hi to present the spatial distribution characteristics of aquifer groundwater level. We can obtain the movement position of the fluid proton after unit time step in accordance with the calculated velocity and direction of the fluid proton. We can also obtain the change rules for the movement position of groundwater fluid and then the dynamic simulation of porous groundwater seepage flow field by repeating the calculation process.


### 2.3.3. Particle Flow Field Visualization

The visualization simulation of porous groundwater seepage flow field based on the particle flow method must consider the real-time effect and authenticity of flow field simulation [43]. The real-time effect of the simulation will become better, but its authenticity will become worse with a reduction in the overall number of particles. By contrast, the real-time effect of the simulation will become worse, but its authenticity will become better with an increase in the overall number of particles. Therefore, the overall number of particles in particle flow should be determined in accordance with the actual simulation situation. In addition, the choice of particle model will also influence the real-time effect and authenticity of the visual expression. The “point particle” graph element structure was selected to simulate groundwater fluid particles, while the graph element structure was used as the drawing unit of particle flow. The texture image of a “water droplet” was mapped onto the point particle rendering unit of particle flow by using the visualization toolkit (VTK) technology to increase the realistic effect of the simulation. This section focuses on the following issues: the generation of seed points, the life cycle of particles, and the formation and smoothing of particle trajectories.(1) Generation of Seed Points. The two methods for generating and distributing the seed points of particle flow are the average new particle number and object area methods [44]. The groundwater level of grid nodes can be obtained via interpolation through the groundwater level of monitoring wells. The generation distribution principle of the object area measurement method was adopted, and the seed points of particle flow were generated at the grid nodes to obtain the velocity of the flow field at the initial time. In this method, the specified average value and maximum variation range are the average value and variance of the generated particles per unit area of the screen. Thus, the number of generated particles depends on the size of the screen occupied by an object. This method can produce different degrees of details and better simulate groundwater seepage flow field than other methods.(2) Life Cycle of Particles. In accordance with the property analysis of groundwater flow particles, the process of particles from production to extinction exhibits a life cycle. The life cycle of a particle determines how long a particle exists in the simulation. In general, the number of frames is used to describe the life cycle of a particle. Therefore, the life cycle of a particle in the i frame can be expressed by the following equation:
(6)Finallifetimet=initiallifetime−t×attenuationperframe,where the initiallifetime is the initial life cycle of a particle, t is the time since the particle was generated, and attenuationperframe is the life decay rate of particles per frame, indicating the number of lives lost by particles per frame. The advantage of using number of frames as the life cycle unit of particles is that it is simple and easy to implement. The disadvantage is that the complexity in generating graphics will directly affect the life cycle and result in the uncertainty of time units. Corresponding new particles are generated at the location of particle extinction to maintain the balance of the total number of particles in the flow field simulation process.(3) Formation and Smoothing of Particle Trajectories. On the basis of the calculation method for particle flow trajectory introduced in Section 2.3.2, the position and velocity of particles at each time frame can be calculated in real time in accordance with existing data conditions. Thus, the movement trajectory of fluid particles in their life cycle was calculated in accordance with the hydrodynamic conditions, and the flow simulation of porous groundwater particle flow field was realized [45]. The real effect of particle flow trajectory is enhanced to the reality of a particle flow field. In the current study, the trajectory from the head position to the tail position was set, and the width of the trajectory was gradually reduced to form the “trailing” effect of particle flow. The life cycle of particles will end when the velocity of fluid particles reaches zero or the number of times of movement approaches five. The transparency of a particle is reduced when it is about to die to enhance the gradual effect of dying particles and reduce the sudden feeling of particle death.


### 2.3.4. Groundwater Level Control Management

Groundwater level and its change trend are important criteria for measuring whether the development of groundwater resources is reasonable [46]. The controlling critical level of groundwater is the critical groundwater level that is determined to prevent various problems or disasters. In accordance with hydrogeological conditions, ecological functions, and geological environment functions caused by the change in groundwater level in different hydrogeological zones [47], the controlling critical level of groundwater is a set of static critical thresholds. The controlling critical level of groundwater can be divided into upper and lower critical groundwater levels. The upper and lower limits of critical groundwater level refer to the maximum and minimum groundwater levels, respectively, that do not cause various problems or disasters. The upper critical groundwater level is caused by excessive groundwater supply or insufficient exploitation. It is the critical groundwater level for preventing soil secondary salinization, road frost boiling, and the antifloating and waterproof design of underground structures, such as subways. The lower critical groundwater level is caused by insufficient groundwater recharge or overexploitation. It is the critical groundwater level used to prevent land subsidence, ground fissures, natural vegetation, and seawater intrusion.The visualization of groundwater seepage flow field was conducted on the basis of particle flow and the contour of groundwater level. On the one hand, the isosurface of groundwater level can reflect its spatial distribution. On the other hand, particle flow direction can indicate the movement trend of groundwater and the recharge relationship among groundwater aquifer group [48]. The specific application methods for groundwater particle flow visualization in the management of the critical groundwater level are as follows.
(1)
Visualization analysis method for groundwater’s controlling critical level (upper and lower limits of critical groundwater level) [49]: the critical groundwater level of an aquifer is consistent within the same critical groundwater level management range. Therefore, the critical water table surface can also be regarded as a constant groundwater level surface. The visual analysis of the study range was performed using the method of generating groundwater level isosurface by considering the geometric characteristics of the controlling critical groundwater level. In addition, the range of groundwater level, which is higher than the lower critical threshold and lower than the upper critical threshold, was rendered in orange by using the spatial data symbolization technology of a geographic information system (GIS) [50]. The orange area indicates the restricted exploitation site in this region. The range of groundwater level that is lower than the lower critical threshold value was rendered in red, and the red area represents the prohibited mining site. The remaining regions are reasonable mining areas. On the basis of this method, controlled critical groundwater level analysis can be performed for groundwater level in the study area(2)
Method for quantifying and visualizing the analytical results of critical groundwater level: in accordance with the analysis results of the controlling critical groundwater level in the study area, we can determine which specific area of the groundwater level is lower than the controlling critical groundwater level and the geographical location of these areas. The spatial topological analysis technique of GIS is used to analyze the topological relationship of the contour of critical groundwater level. The isosurface of groundwater level was generated from the groundwater level isoline and filled with color sequence. The contour of groundwater level has two forms: the closed and unclosed states. Thus, achieving optimization is necessary before the topological groundwater level isoline is establishedUnclosed contour: in accordance with the GIS topological data structure, a surface is composed of unclosed line segments because the existence of a boundary will make part of the isoline unable to close. The unclosed isopleth lines should be connected to a closed line by using the boundary line of the study area to obtain the closed isopleth line.Closed contour: the isosurface is constructed through the topological relationship between lines after closing the unclosed isolines. Then, inclusion analysis is performed on the constructed equivalent surface. Surfaces A and B are overlapping, and surface C is deleted if it exhibits an inclusion relationship with surfaces A and B. Surface C is not deleted if it has no containment relationship with the other surfaces. The presence and absence of containment relationships are depicted in Figures 4 and 5.Figure 4
Containment relationship.Figure 5
No containment relationship.By analogy, all the isosurfaces in the study area are traversed. Finally, the size of the processed face was calculated, and the area that was not “0” was rendered with color. The color sequence from blue to red indicates that the groundwater level is getting deeper.


## 2.3.1. Property Analysis of Particle Flow Field

The initial properties of the particle flow of porous groundwater seepage flow field include the number of initial particles, initial position, particle size, initial velocity and direction of movement, initial color, initial transparency, initial shape, and life cycle. The most important step in generating a flow field is to place tracking particles. The initial position and number of particles directly affect the final visualization quality. Improper placement or low number of particles can lead to failure in tracking important groundwater flow features. Meanwhile, too many particles can also block the line of sight and affect the depiction of major features. Therefore, appropriate initial position and initial number of particles are the keys to flow field visualization. Many researchers have proposed optimal solutions, such as particle placement guidance strategy, spatial average distribution, preferential streamline generation, double streamline, and entropy guidance, to solve the problem of particle placement in a flow field. In the current study, the initial position of particles in groundwater seepage flow field was at the grid node, and the spatial average distribution method was used to optimize the placement of streamline particles.Particle flow was constructed in the form of “point particles,” which can be approximately regarded as pixels with a certain size. Texture mapping technology was used to map the “point particles” to form visual fluid particles. The velocity of particles in porous groundwater seepage flow field has two attributes: size and direction. The contour of groundwater level was obtained by interpolating existing data of groundwater monitoring wells. The velocity magnitude and direction of a flow field at the grid nodes and the initial velocity of fluid particles can be obtained in accordance with Darcy’s law [38]. The life cycle of porous groundwater seepage flow field particles is reflected in the time scale of fluid particles appearing in a given display area, which is typically represented by number of frames. Life cycle determines the forward distance of fluid particles when finding the sample velocity.


## 2.3.2. Calculation of Particle Flow Trajectory

Groundwater flow velocity (seepage velocity) refers to a type of simulated velocity that occurs when water flows through a section with skeletons and pores [38]. In 1856, French hydraulic scientist Henry Darcy found that porous groundwater flows from high to low groundwater level and concluded that
(1)q=KJ=−K∇h,where q is the seepage velocity and K (unit: meter/day) is the conductivity coefficient of the aquifer, which is a measure of an aquifer’s ability to transmit groundwater. The K conductivity coefficient is the core characteristic data of a groundwater environment. J is the hydraulic gradient that refers to the average head loss per unit length of a flow path. ∇h is the gradient of the function h; the direction points to the direction of water head reduction along the normal line of the water head isoline. In three dimensions, the seepage velocity at any point (x, y, and z) in the seepage field can be written as
(2)q=ui⟶+vj⟶+wk⟶.The component of a hydraulic gradientJx, Jy, and Jz in directions x, y, and z can be expressed as
(3)J=−∇h=−dhds,Jx=dhdx,Jy=dhdy,Jz=dhdz.A connection exists between the gradient at any point on the surface and the isoline that passes through that point [39]. The direction of the gradient of the functionz=fx,yatPx,yis the same as the normal direction of the isolinefx,y=Cpassing throughPx,yat this point following the law from the low value to the high value of the isoline. In addition, the module of the gradient is the derivative of the function in the normal direction. The relationship between the gradient at a point on the surface and the isoline through this point is illustrated in Figure 3.Figure 3
The relationship between the gradient at a point on the surface and the isoline through this point.The surface isoline gradient value was obtained, and then, the hydraulic gradient value was determined in accordance with the hydraulic gradient and the relationship between gradients. The method is described as follows.The variation characteristics of the equivalent surface (curved surface) of groundwater level in 3D space are projected onto a plane in 2D space. The (x, y) coordinate of a 2D plane is regarded as the independent variable, and the groundwater level H is considered the dependent variable [21]. The function H=fx,y reflects the regional distribution of groundwater level, and the statistical velocity must involve at least three groundwater level data. Thus, the local trend surface method was adopted to fit groundwater level data in the current study to improve the fitting degree of the trend surface and avoid the situation in which the number of the present groundwater level points near the calculated point is on the low side. In this work, a regular grid that can cover all the study areas was constructed and then combined with the groundwater level of monitoring wells. Kriging interpolation was adopted to calculate groundwater level at the grid intersection point in the study area.The groundwater level of grid points and original monitoring wells obtained via Kriging interpolation statistics was used to establish the trend surface of groundwater level dynamic change [40]. If the point is Px,y, then the radius R is defined and a circle with Px,y as the center and R as the radius is drawn. All the sampling points covered in the circle were used to establish groundwater level trend surface. Considering the global trend and fitting effect of Px,y, the trend surface model with the time of two was selected for the statistics of groundwater level surface. At this moment, the gradient vector of the point on the trend surface is
(4)gradfx,y=∂H∂x,∂H∂y=−a1+2a3x+a4y,a2+a4x+2a5y.Groundwater flows from the higher point to the lower point of mechanical energy; thus, the direction of hydraulic slope and gradient vector is opposite.(5)J=−a1−2a3x−a4y,−a2−a4x−2a5y.Then, the velocity of the pointPx,y is calculated in accordance with Equation (1).The advancing front technique (AFT) algorithm was used on the basis of the boundary data of the study area [41]. The AFT algorithm must be adjacent to each other to generate a regular grid for expressing the distribution of the water head through the monitoring well node objects with attached groundwater level within the scope of the study area. The groundwater level data provided by the dynamic monitoring wells for groundwater level were inserted into the calculation grid as constraint points.A groundwater dynamic monitoring well provides the actual monitoring groundwater level data in discrete space, and groundwater level is continuously distributed in space. Therefore, the spatial interpolation method is required to obtain the groundwater level of each calculation grid node at the initial time to form the initial groundwater level field. Then, this value is distributed to the discrete grid nodes to complete the assignment of the initial groundwater level [42]. The piezometric head value hi of each grid node position (xi and yi ) was extracted, and the digital elevation model of the water head was constructed by xi, yi, and hi to present the spatial distribution characteristics of aquifer groundwater level. We can obtain the movement position of the fluid proton after unit time step in accordance with the calculated velocity and direction of the fluid proton. We can also obtain the change rules for the movement position of groundwater fluid and then the dynamic simulation of porous groundwater seepage flow field by repeating the calculation process.


## 2.3.3. Particle Flow Field Visualization

The visualization simulation of porous groundwater seepage flow field based on the particle flow method must consider the real-time effect and authenticity of flow field simulation [43]. The real-time effect of the simulation will become better, but its authenticity will become worse with a reduction in the overall number of particles. By contrast, the real-time effect of the simulation will become worse, but its authenticity will become better with an increase in the overall number of particles. Therefore, the overall number of particles in particle flow should be determined in accordance with the actual simulation situation. In addition, the choice of particle model will also influence the real-time effect and authenticity of the visual expression. The “point particle” graph element structure was selected to simulate groundwater fluid particles, while the graph element structure was used as the drawing unit of particle flow. The texture image of a “water droplet” was mapped onto the point particle rendering unit of particle flow by using the visualization toolkit (VTK) technology to increase the realistic effect of the simulation. This section focuses on the following issues: the generation of seed points, the life cycle of particles, and the formation and smoothing of particle trajectories.(1) Generation of Seed Points. The two methods for generating and distributing the seed points of particle flow are the average new particle number and object area methods [44]. The groundwater level of grid nodes can be obtained via interpolation through the groundwater level of monitoring wells. The generation distribution principle of the object area measurement method was adopted, and the seed points of particle flow were generated at the grid nodes to obtain the velocity of the flow field at the initial time. In this method, the specified average value and maximum variation range are the average value and variance of the generated particles per unit area of the screen. Thus, the number of generated particles depends on the size of the screen occupied by an object. This method can produce different degrees of details and better simulate groundwater seepage flow field than other methods.(2) Life Cycle of Particles. In accordance with the property analysis of groundwater flow particles, the process of particles from production to extinction exhibits a life cycle. The life cycle of a particle determines how long a particle exists in the simulation. In general, the number of frames is used to describe the life cycle of a particle. Therefore, the life cycle of a particle in the i frame can be expressed by the following equation:
(6)Finallifetimet=initiallifetime−t×attenuationperframe,where the initiallifetime is the initial life cycle of a particle, t is the time since the particle was generated, and attenuationperframe is the life decay rate of particles per frame, indicating the number of lives lost by particles per frame. The advantage of using number of frames as the life cycle unit of particles is that it is simple and easy to implement. The disadvantage is that the complexity in generating graphics will directly affect the life cycle and result in the uncertainty of time units. Corresponding new particles are generated at the location of particle extinction to maintain the balance of the total number of particles in the flow field simulation process.(3) Formation and Smoothing of Particle Trajectories. On the basis of the calculation method for particle flow trajectory introduced in Section 2.3.2, the position and velocity of particles at each time frame can be calculated in real time in accordance with existing data conditions. Thus, the movement trajectory of fluid particles in their life cycle was calculated in accordance with the hydrodynamic conditions, and the flow simulation of porous groundwater particle flow field was realized [45]. The real effect of particle flow trajectory is enhanced to the reality of a particle flow field. In the current study, the trajectory from the head position to the tail position was set, and the width of the trajectory was gradually reduced to form the “trailing” effect of particle flow. The life cycle of particles will end when the velocity of fluid particles reaches zero or the number of times of movement approaches five. The transparency of a particle is reduced when it is about to die to enhance the gradual effect of dying particles and reduce the sudden feeling of particle death.


## 2.3.4. Groundwater Level Control Management

Groundwater level and its change trend are important criteria for measuring whether the development of groundwater resources is reasonable [46]. The controlling critical level of groundwater is the critical groundwater level that is determined to prevent various problems or disasters. In accordance with hydrogeological conditions, ecological functions, and geological environment functions caused by the change in groundwater level in different hydrogeological zones [47], the controlling critical level of groundwater is a set of static critical thresholds. The controlling critical level of groundwater can be divided into upper and lower critical groundwater levels. The upper and lower limits of critical groundwater level refer to the maximum and minimum groundwater levels, respectively, that do not cause various problems or disasters. The upper critical groundwater level is caused by excessive groundwater supply or insufficient exploitation. It is the critical groundwater level for preventing soil secondary salinization, road frost boiling, and the antifloating and waterproof design of underground structures, such as subways. The lower critical groundwater level is caused by insufficient groundwater recharge or overexploitation. It is the critical groundwater level used to prevent land subsidence, ground fissures, natural vegetation, and seawater intrusion.The visualization of groundwater seepage flow field was conducted on the basis of particle flow and the contour of groundwater level. On the one hand, the isosurface of groundwater level can reflect its spatial distribution. On the other hand, particle flow direction can indicate the movement trend of groundwater and the recharge relationship among groundwater aquifer group [48]. The specific application methods for groundwater particle flow visualization in the management of the critical groundwater level are as follows.
(1)
Visualization analysis method for groundwater’s controlling critical level (upper and lower limits of critical groundwater level) [49]: the critical groundwater level of an aquifer is consistent within the same critical groundwater level management range. Therefore, the critical water table surface can also be regarded as a constant groundwater level surface. The visual analysis of the study range was performed using the method of generating groundwater level isosurface by considering the geometric characteristics of the controlling critical groundwater level. In addition, the range of groundwater level, which is higher than the lower critical threshold and lower than the upper critical threshold, was rendered in orange by using the spatial data symbolization technology of a geographic information system (GIS) [50]. The orange area indicates the restricted exploitation site in this region. The range of groundwater level that is lower than the lower critical threshold value was rendered in red, and the red area represents the prohibited mining site. The remaining regions are reasonable mining areas. On the basis of this method, controlled critical groundwater level analysis can be performed for groundwater level in the study area(2)
Method for quantifying and visualizing the analytical results of critical groundwater level: in accordance with the analysis results of the controlling critical groundwater level in the study area, we can determine which specific area of the groundwater level is lower than the controlling critical groundwater level and the geographical location of these areas. The spatial topological analysis technique of GIS is used to analyze the topological relationship of the contour of critical groundwater level. The isosurface of groundwater level was generated from the groundwater level isoline and filled with color sequence. The contour of groundwater level has two forms: the closed and unclosed states. Thus, achieving optimization is necessary before the topological groundwater level isoline is establishedUnclosed contour: in accordance with the GIS topological data structure, a surface is composed of unclosed line segments because the existence of a boundary will make part of the isoline unable to close. The unclosed isopleth lines should be connected to a closed line by using the boundary line of the study area to obtain the closed isopleth line.Closed contour: the isosurface is constructed through the topological relationship between lines after closing the unclosed isolines. Then, inclusion analysis is performed on the constructed equivalent surface. Surfaces A and B are overlapping, and surface C is deleted if it exhibits an inclusion relationship with surfaces A and B. Surface C is not deleted if it has no containment relationship with the other surfaces. The presence and absence of containment relationships are depicted in Figures 4 and 5.Figure 4
Containment relationship.Figure 5
No containment relationship.By analogy, all the isosurfaces in the study area are traversed. Finally, the size of the processed face was calculated, and the area that was not “0” was rendered with color. The color sequence from blue to red indicates that the groundwater level is getting deeper.


## 3. Results and Conclusion

### 3.1. Groundwater Level Distribution

The groundwater level of each node in the regional grid was obtained via interpolation calculation based on the known groundwater level data of monitoring wells in the study area [38]. On this basis, we can obtain groundwater level distribution and variation within the scope of the study area and draw the contour groundwater level map and its regional scope. Regional groundwater level distribution is a necessary information for groundwater flow mapping. The data required for flow field drawing can be calculated and processed on the basis of the results of groundwater level isoline and isosurface. Finally, mapping the flow field data and rendering and displaying porous groundwater seepage flow field were realized on the basis of different visualization methods. Figure 6 shows the groundwater level isosurface map in the study area. The colors in the figure represent the depth of the groundwater level, with blue indicating the shallowest groundwater level and red indicating the deepest groundwater level.Figure 6
Map of the groundwater level isosurface.


### 3.2. Comparison of Flow Field Visualization

The dynamic adaptive grid in the study area was generated via the finite element method and the AFT algorithm on the basis of existing boundary data of the study area [40]. The monitoring wells that are evenly distributed in the study area and have a large average variation of groundwater level were selected and inserted into the generated grid. Kriging interpolation simulation was performed to obtain groundwater level at each grid node [41]. The velocity of the flow field at the grid nodes was calculated in accordance with Darcy’s law. The contour of groundwater level was constructed, and the hydraulic gradient value was calculated on the basis of the obtained regional groundwater level distribution. Finally, the velocity and direction of groundwater movement at each grid node in the study area were obtained in accordance with Darcy’s law.


#### 3.2.1. Geometric Point Method

The geometric point method, which visualizes flow field in the form of arrows, can directly and vividly display groundwater movement trend and recharge relationship in a region. Its data mapping process is simple, easy, and fast; thus, this method is widely used. However, the point graph method is a static simulation in the instantaneous state. This method can be combined and played frame by frame through the simulation results under multiple step sizes, and thus, it presents an animation effect. However, the geometric point method is still a false dynamic simulation technique that cannot efficiently realize the dynamic visualization of groundwater seepage flow field. The length of the arrows is set to be the same, and the direction of the arrows is used to indicate the direction of the movement speed of groundwater seepage flow field at the current position. The color of the arrow is used to indicate the magnitude of groundwater seepage flow field velocity at the current location. The visualization results of porous groundwater seepage flow field based on the point plot method are presented in Figure7.Figure 7
The visualization results of porous groundwater seepage flow field based on point plot method.


#### 3.2.2. Vector Line Method

One of the advantages of the vector line method over the geometric point method is that it can achieve better continuity in flow field visualization simulation results. Therefore, the characteristics and details of a flow field can be displayed and revealed well. However, the visual simulation of the streamline form should be judged artificially in accordance with the groundwater level isoline or isosurface. Moreover, this method is also a static flow field simulation in instantaneous state. Streamline distribution can be calculated on the basis of the regional groundwater level isoline in accordance with the principle that the streamline is the normal direction of the groundwater level isoline. Then, the vector-based visual simulation results of porous groundwater seepage flow field can be drawn. The number of regional streamlines was set to 90 in this section considering the visual effect of flow field simulation. The visualization results of porous groundwater seepage flow field based on the vector line method are provided in Figure8.Figure 8
The visualization results of porous groundwater seepage flow field based on vector line method.


#### 3.2.3. Particle Flow Method

The particle flow method can accurately reflect and display the dynamic properties of groundwater seepage flow field compared with traditional flow visualization methods, such as the geometric point and vector line methods. Moreover, the particle flow method uses a large number of particles with relevant properties to replace fluid protons in the flow field and then gather to form a system as a whole to achieve flow field visualization [45].Many particles in the particle flow method are represented by the primitive structure of “point particles,” and texture mapping technology based on VTK enhances the authenticity of the simulation effect. In addition, the dynamics of groundwater seepage flow field simulation effect is enhanced to present a more intuitive image of the dynamic effect. The following visualization strategy, which is set in the life cycle of a single particle, is adopted when a particle is newly generated and its size is still small. The size of a particle increases gradually with time until its death at the end of its life cycle. Particles gradually become thicker in the process of motion by controlling the size of particles and presenting a sports effect of “trailing” form with a thicker head and a thinner tail. A particle that is about to die has reduced transparency to ease the sudden feeling when the particle disappears.The total number of particles in the area was set at 240 to consider the operation efficiency and visual effect of visual simulation. The specific settings of the initial and dynamic attributes of particle flow are provided in Table1. The number of particle position movements was artificially specified to be five times. The life cycle of a particle ends if the number of its position movement reaches five times when particle velocity is greater than 0 and the particle does not leave the study area.Table 1
Attribute assignment table for particle flow visualization.
Attribute typeAttribute assignmentInitial attributeNumber of initial particles240Seed point positionAt the location of the grid nodeInitial color and transparencyWhite, 1Particle size and shape1.2 cm, point particleInitial velocityCalculating the flow velocity at the node by Darcy’s lawLife cycleNumber ofmoves≤5Particlevelocity>0Particles are not out of rangeDynamic attributeGeneration of particlesCodetermination of initial properties of particle flowMotion of a particleAccording to the initial position, velocity, and Darcy’s lawExtinction of particlesSurvival conditions that do not meet the life cycleThe simulation results of porous groundwater particle flow field at three adjacent time nodes were intercepted to more intuitively reveal the dynamic visualization simulation results of porous groundwater seepage flow field based on the particle flow method. The results of the visualization of porous groundwater seepage flow field based on the particle flow method are presented in Figure9.Figure 9
The visualization results of porous groundwater seepage flow field based on particle flow method.
(a)(b)(c)The particle flow method can better express the dynamics and authenticity of groundwater seepage flow field in accordance with its unique system characteristics and simulation methods compared with the characteristics of instantaneous static or pseudodynamic simulations of other flow field visualization simulation methods. In contrast with the point chart method, the particle flow method experiences difficulty in describing the velocity of a flow field [24, 33]. Moreover, particle flow simulation has a relatively large calculation and thus has higher system performance requirements. The theory and method of particle flow were introduced into the visual simulation of porous groundwater seepage flow field to realize the dynamic visual expression of a flow field. Furthermore, particle flow was combined with groundwater numerical simulation. Particle attributes suitable for porous groundwater were designed in accordance with the dynamic calculation results of porous groundwater seepage flow field. These attributes can realize the detailed simulation and realistic effect of flow field visualization and can exhibit a variety of physical states, and thus, they exhibit certain practicability and relative advantages under specific conditions.The monitoring data of the groundwater level of the third confined aquifer in Yancheng City, Jiangsu Province, in July 2013 were selected, and the groundwater level isosurface and particle flow field in the study area are plotted in Figure8. Groundwater level in the eastern part of the city is evidently higher than that in the western part. In addition, groundwater level in the urban area of the city is generally low. The flow direction of groundwater and the supply relationship of groundwater in the study area can be more intuitively observed through the flow direction of particle flow by superposing the contour map of particle flow field and groundwater level. Simultaneously, we can clearly see which area is lower than the lower and upper critical threshold values by analyzing critical groundwater level. As shown in Figure 9, the orange range represents the restricted mining area, the red section represents the prohibited mining area, and the remaining regions are reasonable mining areas. We can count the areas where groundwater level is lower than the lower limit of groundwater level and in which areas should the exploitation of groundwater be controlled through the analysis of critical groundwater level in the study area.


## 3.1. Groundwater Level Distribution

The groundwater level of each node in the regional grid was obtained via interpolation calculation based on the known groundwater level data of monitoring wells in the study area [38]. On this basis, we can obtain groundwater level distribution and variation within the scope of the study area and draw the contour groundwater level map and its regional scope. Regional groundwater level distribution is a necessary information for groundwater flow mapping. The data required for flow field drawing can be calculated and processed on the basis of the results of groundwater level isoline and isosurface. Finally, mapping the flow field data and rendering and displaying porous groundwater seepage flow field were realized on the basis of different visualization methods. Figure 6 shows the groundwater level isosurface map in the study area. The colors in the figure represent the depth of the groundwater level, with blue indicating the shallowest groundwater level and red indicating the deepest groundwater level.Figure 6
Map of the groundwater level isosurface.


## 3.2. Comparison of Flow Field Visualization

The dynamic adaptive grid in the study area was generated via the finite element method and the AFT algorithm on the basis of existing boundary data of the study area [40]. The monitoring wells that are evenly distributed in the study area and have a large average variation of groundwater level were selected and inserted into the generated grid. Kriging interpolation simulation was performed to obtain groundwater level at each grid node [41]. The velocity of the flow field at the grid nodes was calculated in accordance with Darcy’s law. The contour of groundwater level was constructed, and the hydraulic gradient value was calculated on the basis of the obtained regional groundwater level distribution. Finally, the velocity and direction of groundwater movement at each grid node in the study area were obtained in accordance with Darcy’s law.


### 3.2.1. Geometric Point Method

The geometric point method, which visualizes flow field in the form of arrows, can directly and vividly display groundwater movement trend and recharge relationship in a region. Its data mapping process is simple, easy, and fast; thus, this method is widely used. However, the point graph method is a static simulation in the instantaneous state. This method can be combined and played frame by frame through the simulation results under multiple step sizes, and thus, it presents an animation effect. However, the geometric point method is still a false dynamic simulation technique that cannot efficiently realize the dynamic visualization of groundwater seepage flow field. The length of the arrows is set to be the same, and the direction of the arrows is used to indicate the direction of the movement speed of groundwater seepage flow field at the current position. The color of the arrow is used to indicate the magnitude of groundwater seepage flow field velocity at the current location. The visualization results of porous groundwater seepage flow field based on the point plot method are presented in Figure7.Figure 7
The visualization results of porous groundwater seepage flow field based on point plot method.


### 3.2.2. Vector Line Method

One of the advantages of the vector line method over the geometric point method is that it can achieve better continuity in flow field visualization simulation results. Therefore, the characteristics and details of a flow field can be displayed and revealed well. However, the visual simulation of the streamline form should be judged artificially in accordance with the groundwater level isoline or isosurface. Moreover, this method is also a static flow field simulation in instantaneous state. Streamline distribution can be calculated on the basis of the regional groundwater level isoline in accordance with the principle that the streamline is the normal direction of the groundwater level isoline. Then, the vector-based visual simulation results of porous groundwater seepage flow field can be drawn. The number of regional streamlines was set to 90 in this section considering the visual effect of flow field simulation. The visualization results of porous groundwater seepage flow field based on the vector line method are provided in Figure8.Figure 8
The visualization results of porous groundwater seepage flow field based on vector line method.


### 3.2.3. Particle Flow Method

The particle flow method can accurately reflect and display the dynamic properties of groundwater seepage flow field compared with traditional flow visualization methods, such as the geometric point and vector line methods. Moreover, the particle flow method uses a large number of particles with relevant properties to replace fluid protons in the flow field and then gather to form a system as a whole to achieve flow field visualization [45].Many particles in the particle flow method are represented by the primitive structure of “point particles,” and texture mapping technology based on VTK enhances the authenticity of the simulation effect. In addition, the dynamics of groundwater seepage flow field simulation effect is enhanced to present a more intuitive image of the dynamic effect. The following visualization strategy, which is set in the life cycle of a single particle, is adopted when a particle is newly generated and its size is still small. The size of a particle increases gradually with time until its death at the end of its life cycle. Particles gradually become thicker in the process of motion by controlling the size of particles and presenting a sports effect of “trailing” form with a thicker head and a thinner tail. A particle that is about to die has reduced transparency to ease the sudden feeling when the particle disappears.The total number of particles in the area was set at 240 to consider the operation efficiency and visual effect of visual simulation. The specific settings of the initial and dynamic attributes of particle flow are provided in Table1. The number of particle position movements was artificially specified to be five times. The life cycle of a particle ends if the number of its position movement reaches five times when particle velocity is greater than 0 and the particle does not leave the study area.Table 1
Attribute assignment table for particle flow visualization.
Attribute typeAttribute assignmentInitial attributeNumber of initial particles240Seed point positionAt the location of the grid nodeInitial color and transparencyWhite, 1Particle size and shape1.2 cm, point particleInitial velocityCalculating the flow velocity at the node by Darcy’s lawLife cycleNumber ofmoves≤5Particlevelocity>0Particles are not out of rangeDynamic attributeGeneration of particlesCodetermination of initial properties of particle flowMotion of a particleAccording to the initial position, velocity, and Darcy’s lawExtinction of particlesSurvival conditions that do not meet the life cycleThe simulation results of porous groundwater particle flow field at three adjacent time nodes were intercepted to more intuitively reveal the dynamic visualization simulation results of porous groundwater seepage flow field based on the particle flow method. The results of the visualization of porous groundwater seepage flow field based on the particle flow method are presented in Figure9.Figure 9
The visualization results of porous groundwater seepage flow field based on particle flow method.
(a)(b)(c)The particle flow method can better express the dynamics and authenticity of groundwater seepage flow field in accordance with its unique system characteristics and simulation methods compared with the characteristics of instantaneous static or pseudodynamic simulations of other flow field visualization simulation methods. In contrast with the point chart method, the particle flow method experiences difficulty in describing the velocity of a flow field [24, 33]. Moreover, particle flow simulation has a relatively large calculation and thus has higher system performance requirements. The theory and method of particle flow were introduced into the visual simulation of porous groundwater seepage flow field to realize the dynamic visual expression of a flow field. Furthermore, particle flow was combined with groundwater numerical simulation. Particle attributes suitable for porous groundwater were designed in accordance with the dynamic calculation results of porous groundwater seepage flow field. These attributes can realize the detailed simulation and realistic effect of flow field visualization and can exhibit a variety of physical states, and thus, they exhibit certain practicability and relative advantages under specific conditions.The monitoring data of the groundwater level of the third confined aquifer in Yancheng City, Jiangsu Province, in July 2013 were selected, and the groundwater level isosurface and particle flow field in the study area are plotted in Figure8. Groundwater level in the eastern part of the city is evidently higher than that in the western part. In addition, groundwater level in the urban area of the city is generally low. The flow direction of groundwater and the supply relationship of groundwater in the study area can be more intuitively observed through the flow direction of particle flow by superposing the contour map of particle flow field and groundwater level. Simultaneously, we can clearly see which area is lower than the lower and upper critical threshold values by analyzing critical groundwater level. As shown in Figure 9, the orange range represents the restricted mining area, the red section represents the prohibited mining area, and the remaining regions are reasonable mining areas. We can count the areas where groundwater level is lower than the lower limit of groundwater level and in which areas should the exploitation of groundwater be controlled through the analysis of critical groundwater level in the study area.


## 3.2.1. Geometric Point Method

The geometric point method, which visualizes flow field in the form of arrows, can directly and vividly display groundwater movement trend and recharge relationship in a region. Its data mapping process is simple, easy, and fast; thus, this method is widely used. However, the point graph method is a static simulation in the instantaneous state. This method can be combined and played frame by frame through the simulation results under multiple step sizes, and thus, it presents an animation effect. However, the geometric point method is still a false dynamic simulation technique that cannot efficiently realize the dynamic visualization of groundwater seepage flow field. The length of the arrows is set to be the same, and the direction of the arrows is used to indicate the direction of the movement speed of groundwater seepage flow field at the current position. The color of the arrow is used to indicate the magnitude of groundwater seepage flow field velocity at the current location. The visualization results of porous groundwater seepage flow field based on the point plot method are presented in Figure7.Figure 7
The visualization results of porous groundwater seepage flow field based on point plot method.


## 3.2.2. Vector Line Method

One of the advantages of the vector line method over the geometric point method is that it can achieve better continuity in flow field visualization simulation results. Therefore, the characteristics and details of a flow field can be displayed and revealed well. However, the visual simulation of the streamline form should be judged artificially in accordance with the groundwater level isoline or isosurface. Moreover, this method is also a static flow field simulation in instantaneous state. Streamline distribution can be calculated on the basis of the regional groundwater level isoline in accordance with the principle that the streamline is the normal direction of the groundwater level isoline. Then, the vector-based visual simulation results of porous groundwater seepage flow field can be drawn. The number of regional streamlines was set to 90 in this section considering the visual effect of flow field simulation. The visualization results of porous groundwater seepage flow field based on the vector line method are provided in Figure8.Figure 8
The visualization results of porous groundwater seepage flow field based on vector line method.


## 3.2.3. Particle Flow Method

The particle flow method can accurately reflect and display the dynamic properties of groundwater seepage flow field compared with traditional flow visualization methods, such as the geometric point and vector line methods. Moreover, the particle flow method uses a large number of particles with relevant properties to replace fluid protons in the flow field and then gather to form a system as a whole to achieve flow field visualization [45].Many particles in the particle flow method are represented by the primitive structure of “point particles,” and texture mapping technology based on VTK enhances the authenticity of the simulation effect. In addition, the dynamics of groundwater seepage flow field simulation effect is enhanced to present a more intuitive image of the dynamic effect. The following visualization strategy, which is set in the life cycle of a single particle, is adopted when a particle is newly generated and its size is still small. The size of a particle increases gradually with time until its death at the end of its life cycle. Particles gradually become thicker in the process of motion by controlling the size of particles and presenting a sports effect of “trailing” form with a thicker head and a thinner tail. A particle that is about to die has reduced transparency to ease the sudden feeling when the particle disappears.The total number of particles in the area was set at 240 to consider the operation efficiency and visual effect of visual simulation. The specific settings of the initial and dynamic attributes of particle flow are provided in Table1. The number of particle position movements was artificially specified to be five times. The life cycle of a particle ends if the number of its position movement reaches five times when particle velocity is greater than 0 and the particle does not leave the study area.Table 1
Attribute assignment table for particle flow visualization.
Attribute typeAttribute assignmentInitial attributeNumber of initial particles240Seed point positionAt the location of the grid nodeInitial color and transparencyWhite, 1Particle size and shape1.2 cm, point particleInitial velocityCalculating the flow velocity at the node by Darcy’s lawLife cycleNumber ofmoves≤5Particlevelocity>0Particles are not out of rangeDynamic attributeGeneration of particlesCodetermination of initial properties of particle flowMotion of a particleAccording to the initial position, velocity, and Darcy’s lawExtinction of particlesSurvival conditions that do not meet the life cycleThe simulation results of porous groundwater particle flow field at three adjacent time nodes were intercepted to more intuitively reveal the dynamic visualization simulation results of porous groundwater seepage flow field based on the particle flow method. The results of the visualization of porous groundwater seepage flow field based on the particle flow method are presented in Figure9.Figure 9
The visualization results of porous groundwater seepage flow field based on particle flow method.
(a)(b)(c)The particle flow method can better express the dynamics and authenticity of groundwater seepage flow field in accordance with its unique system characteristics and simulation methods compared with the characteristics of instantaneous static or pseudodynamic simulations of other flow field visualization simulation methods. In contrast with the point chart method, the particle flow method experiences difficulty in describing the velocity of a flow field [24, 33]. Moreover, particle flow simulation has a relatively large calculation and thus has higher system performance requirements. The theory and method of particle flow were introduced into the visual simulation of porous groundwater seepage flow field to realize the dynamic visual expression of a flow field. Furthermore, particle flow was combined with groundwater numerical simulation. Particle attributes suitable for porous groundwater were designed in accordance with the dynamic calculation results of porous groundwater seepage flow field. These attributes can realize the detailed simulation and realistic effect of flow field visualization and can exhibit a variety of physical states, and thus, they exhibit certain practicability and relative advantages under specific conditions.The monitoring data of the groundwater level of the third confined aquifer in Yancheng City, Jiangsu Province, in July 2013 were selected, and the groundwater level isosurface and particle flow field in the study area are plotted in Figure8. Groundwater level in the eastern part of the city is evidently higher than that in the western part. In addition, groundwater level in the urban area of the city is generally low. The flow direction of groundwater and the supply relationship of groundwater in the study area can be more intuitively observed through the flow direction of particle flow by superposing the contour map of particle flow field and groundwater level. Simultaneously, we can clearly see which area is lower than the lower and upper critical threshold values by analyzing critical groundwater level. As shown in Figure 9, the orange range represents the restricted mining area, the red section represents the prohibited mining area, and the remaining regions are reasonable mining areas. We can count the areas where groundwater level is lower than the lower limit of groundwater level and in which areas should the exploitation of groundwater be controlled through the analysis of critical groundwater level in the study area.


## 4. Conclusions

This study explored the feasibility and applicability of the particle flow method in the visual simulation of porous groundwater flow dynamics and performed related experimental verification. The characteristics of existing conventional visualization methods were analyzed on the basis of the characteristics of the porous water movement mechanism and flow field. The particle flow method in computer graphics was introduced into the visual simulation of porous groundwater seepage flow field to fully represent the dynamic properties of a flow field. The simulation results of different visualization methods were compared and analyzed. The spatial differentiation of groundwater seepage field and the critical groundwater level of groundwater control were integrated and superimposed. The spatial distribution characteristics of groundwater exploitation intensity and critical groundwater level were revealed from the 2D perspective on the basis of the groundwater seepage field, and the regional scope of the critical groundwater level was automatically delineated.However, the spatiotemporal variation of porous groundwater seepage flow field is a complex and uncertain process influenced by many factors. Variations in groundwater seepage flow field and the complexity of visual simulation are caused by the difference in hydrogeological conditions. The focus of the current study is to introduce the idea of particle flow to visualize the expression of porous groundwater seepage flow field. Therefore, this study did not conduct more in-depth research and analysis of the complexity of porous groundwater seepage flow field, such as turbulence and eddy in groundwater seepage flow field. Hydrogeological conditions must be fully considered, and special conditions existing in a flow field must be properly dealt with to reflect the rules of porous groundwater seepage flow field movement more accurately.



---

*Source: 4850968-2022-11-15.xml*