To provide theoretical support for the protection of dispersed drinking water sources of groundwater, we need to accurately evaluate the time and scope of groundwater pollution hazards to human health. This helps the decision-making process for remediation of polluted soil and groundwater in service stations. In this study, we conducted such an evaluation by coupling numerical modeling with a health risk assessment. During the research, soil and groundwater samples were collected and analyzed for 20 pollutants. Fifty-six percent of the heavy contaminants and 100% of the organic contaminants exhibited maximum values at the location of the oil depot. Gray correlation analysis showed that the correlation between background samples and soil underlying the depot was 0.375–0.567 (barely significant to insignificant). The correlation between the reference sequence of other points was 0.950–0.990 (excellent correlation). The correlation of environmental impact after oil depot leakage followed the order: organic pollutants > heavy metals > inorganic pollutants. The groundwater simulation status and predictions indicated that non-carcinogenic health risks covered 25,462 m2 at the time of investigation, and were predicted to extend to 29,593 m2 after five years and to 39,873 m2 after 10 years. Carcinogenic health risks covered 21,390 m2 at the time of investigation, and were predicted to extend to 40,093 m2 after five years and to 53,488 m2 after 10 years. This study provides theoretical support for the protection of a dispersed drinking water source such as groundwater, and also helps the decision-making process for groundwater and soil environment improvement.

Local head loss caused by fracture intersection is often ignored because there has not been a simple method to calculate it until now. Relevant research shows that neglecting the local flow resistance leads to inaccurate results, especially when the velocity and cross angle are large. Therefore, it is necessary to find a portable method for calculation. Physical experiments of single fracture with different apertures (e = 0.77, 1.18, 1.97, 2.73 mm) were set up first to study the flow characteristics, showing obvious non-Darcian flow, which can be depicted by the Forchheimer equation when the flow velocity is sufficiently large. The computational fluid dynamics (CFD) software ANSYS FLUENT was used to build numeric simulation models. A good correlation between CFD simulation results and physical experiment results was found (Pearson's correlation coefficient > 0.99). Then, the CFD models of flexural crack with different angles from 30 • to 150 • were established to compute the pressure drop of flexural crack at different velocity. It was found that the local head loss of the flexural crack varied with the bending angle, and its coefficient was expressed by the deformation of the logistic equation. By using this model, as well as a frictional head loss equation fitted by Forchheimer equation, the head loss of crossed fissures with fixed fracture aperture could be easily calculated. Water 2020, 12, 232 2 of 15 applications, a quadratic equation (J = Av + Bv 2 ) and an exponential equation (J = Cv −m ) are used to under limited conditions [8-11]. The Forchheimer equation can express the characteristics of non-Darcian flow in fracture [12-15]: − ∇P = aQ + bQ 2 (1) where Q is volumetric flow velocity, a and b are model coefficients, and −∇P is the pressure gradient. The factors affecting the flow of fissure fluid are complex and include the shape of the fissure, degree of contact, roughness, fluid viscosity, and external pressure [16-19]. Liu et al. [20] studied the law of fluid movement in a fracture network and found that the coefficient (a, b) of the Forchheimer equation decreases with an increase in the fracture aperture, and the error in model fitting can be reduced by expressing the coefficients of the Forchheimer equation as a power equation of fracture aperture [14]. Olson et al. [21] studied the effect of the ratio of the fracture opening to its height on the flow pattern, and Shu et al. [22] studied the variation in head loss with the width of the fissure and the velocity of flow when water flows through L-shaped fissures. Li et al. [23] studied the correlation between the head loss of a fluid in cross fissures and the width as well as roughness of the fissures.This study proposes a portable calculation model where the aperture and the shape of flexural crack are considered to predict head loss directly. It is difficult to consider all factors when studying the characteristics of flow in fractures. In order to investigate the influence of each one, the synergistic effects of each condition need to be reduced, ...