*Article* **Research on the Application of Typical Biological Chain for Algal Control in Lake Ecological Restoration—A Case Study of Lianshi Lake in Yongding River**

**Pengfei Zhang 1,2, Xiaoyu Cui 2, Huihuang Luo 2,\*, Wenqi Peng <sup>2</sup> and Yunxia Gao <sup>1</sup>**


**Abstract:** Maintaining the health of lake ecosystems is an urgent issue. However, eutrophication seriously affects lakes' ecological functions. Eutrophication is also the main target of lake ecological restoration. It is vital to carry out research on lake eutrophication control and energy flow evaluation in ecosystems scientifically. Based on in situ survey results for the aquatic life data for Lianshi Lake from 2018 to 2019, the Ecopath model was used to establish an evaluation index system for the typical biological chain to screen out the key species in the water ecosystem, and the fuzzy comprehensive evaluation (FCE) method was used to screen all the biological chains controlling algae. A combination of the FCE coupled with the Ecopath screening method for typical biological chains for algal control was applied to the Lianshi Lake area; the results show that the typical biological chain for algal control is phytoplankton (Phyt)–zooplankton (Zoop)–macrocrustaceans (Macc)–other piscivorous (OthP). Upon adjusting the biomass of Zoop and Macc in the typical biological chain for algal control to three times that of the current status, the ecological nutrition efficiency of Phyt was increased from 0.308 to 0.906. The material flow into the second trophic level from primary producers increased from 3043 to 8283 t/km2/year. The amount of detritus flowing into primary producers for sedimentation decreased from 7618 to 2378 t/km2/year. Finally, the total primary production/total respiratory volume (TPP/TR) decreased from 9.224 to 3.403, the Finn's cycle index (FCI) increased from 13.6% to 17.5%, and the Finn's average energy flow path length (FCL) increased from 2.854 to 3.410. The results suggest that the problem of eutrophication can be solved by introducing Zoop (an algal predator) and Macc to a large extent, resulting in improved ecosystem maturity. The research results can facilitate decision making for the restoration of urban lake water ecosystems.

**Keywords:** alga control; typical biological chain; Ecopath model; ecological restoration; Lianshi Lake

#### **1. Introduction**

The rapid expansion of the world's population has exacerbated the degradation of global lake ecosystems [1]. According to statistics, more than 60% of the lakes in the world are in different degrees of eutrophication [2]. The increasing eutrophication of lakes has become a global water environment problem [3] (e.g., for Lake Canyon in the United States [4], Lake Geneva in Switzerland [5], and Sugarloaf Lake in Australia [6]).

The interaction between lake water environments and water ecology is complex. It is important to ensure the integrity of the ecosystem while improving the quality of the water environment [7–9]. Therefore, how to deal with lake eutrophication and restore aquatic ecosystems has become an urgent problem to be solved in current limnology, environmental science, freshwater ecology and other disciplines, and has been closely followed by international scholars [10,11].

**Citation:** Zhang, P.; Cui, X.; Luo, H.; Peng, W.; Gao, Y. Research on the Application of Typical Biological Chain for Algal Control in Lake Ecological Restoration—A Case Study of Lianshi Lake in Yongding River. *Water* **2021**, *13*, 3079. https:// doi.org/10.3390/w13213079

Academic Editor: Antonio Zuorro

Received: 4 September 2021 Accepted: 24 October 2021 Published: 2 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

At present, the recognized theories of lake ecological restoration mainly include the nutrient salt concentration limit theory by Scheffer [12], multi-steady state theory by Lewontin [13] and biological manipulation theory by Hrbacek [14]. Among them, the biological manipulation theory has the most promising application prospects and has led to many successful cases of lake ecological restoration. For example, the biological manipulation process implemented by putting silver carp (Silc) and bighead carp (Bigc) in Donghu Lake in Wuhan is one of the main reasons for the disappearance of cyanobacteria blooms in the lake [15]. Olin et al. [16] reported that the removal of common carp (Comc) from 10 lakes in southern Finland effectively reduced the biomass of cyanobacteria and the degree of algal outbreaks. Shapiro et al. [17] adjusted the ratio of other piscivarious (Othp) and plankton-eating fish from 1:1.65 to 1:2.2 in the Round Lake, and the concentrations of total nitrogen (TN), total phosphorus (TP) and chlorophyl-a (Chl-*a*) in the lake were reduced to varying degrees. However, biological manipulation cannot enable all ecosystems to achieve the expected ecological functions and may lead to changes in the species diversity of the ecosystem, the decline of the average trophic level in the system, and the destruction of habitats [18,19]. For example, Razlutskij et al. [20], through 72 days of outdoor experiments, found that introducing *Carassius auratus* increased the biomass of Phyt. Recent ecological studies have revealed significant negative effects of crucian carp (Cruc) on the water quality and ecological states of shallow lakes, e.g., increasing nutrient levels, leading to reduced water clarity [21,22]. The current international guidance on the application of biological manipulation technology to the process of lake ecological restoration remains far from sufficient [23].

In light of this, this study intended to establish a screening method for typical biological chains of ecosystems and propose a biomass control strategy for the typical biological chain of alga control. The purpose of this study was to provide strong theoretical support for the management of urban lake eutrophication and water ecological restoration, and it consisted of the following: (1) Key species and target biological chains were screened, based on the Ecopath model; the key species of the ecosystem were analyzed, and the biological chain with phytoplankton (Phyt) as the primary producer and including the key species was selected. (2) Typical biological chains were screened, and an evaluation index that could reflect the efficiency and stability of the biological chain was constructed. The index weight was used to select typical biological chains based on the fuzzy comprehensive evaluation method. (3) The regulatory impact was analyzed, the biomass of key species in the typical biological chain was regulated, and then the impacts on the regulation target and the ecosystem were judged.

#### **2. Materials and Methods**

#### *2.1. Overview of the Study Area and Data Sources*

2.1.1. Overview of the Study Area

As a river-type lake in the urban section of the plain of the Yongding River, Lianshi Lake is located at the junction of Mentougou and Fengtai District in the southwest of Shijingshan District (116◦6 58~116◦9 34 E, 39◦56 3~39◦53 15 N), which belongs to the temperate continental monsoon climate, with high temperatures and rain in summer and cold and dry conditions in winter. The total length of Lianshi Lake is 5.8 km, the average width of the lake is 376 m, the widest point is about 500 m, the total water area is 106 hm2, and the average depth of the water body is 1.6 m. The phytoplankton in Lianshi Lake are dominated by *Pseudanabaena moniliformis* (Cyanophinyta) and *Limnothrix*. Lianshi Lake belongs to the northern freshwater city lake. The distribution of the study area and monitoring points are shown in Figure 1.

**Figure 1.** The location of Lianshi Lake and the distribution of monitoring points.

#### 2.1.2. Data Sources

The data sources included long-term field sampling survey data, formula calculations and reference comparisons. The main source of the biomass *Bi* of each function group was a field survey of the biomass at six points in Lianshi Lake from 2018 to 2019. The macroinvertebrates were collected using a Surber sampler with a mesh diameter of 40 mesh (500 μm) and a sampling area of 0.09 m2, and 70% alcohol was added to the test tube bottle to be stored for inspection. The zooplankton were collected using a No. 25 plankton net (200 mesh) and immediately fixed with 4% formaldehyde. For the collection of phytoplankton, a 1 L water sample (0.5 m below the surface of the water body) was collected in a plexiglass water collector, and 15 mL of Lugol's reagent was added for fixation. Sampling areas of 0.5 m × 0.5 m and 1 m × 1 m were used for vascular plants and submerged plants, respectively. In the field, the plant species were distinguished and weighed (wet weight). Fish resources were collected mainly using trawl nets, and the travel distance for each sampling point was no more than 100 m. The biomass of organic detritus was calculated based on the relationship between the primary productivity and water transparency [24] (Equation (1)).

$$
\log D = 0.954 \log PP + 0.863 \log E - 2.41 \tag{1}
$$

where *<sup>D</sup>* is the detrital biomass (g C·m−2), *PP* is the primary production volume (g C·m−2·a<sup>−</sup>1), and *<sup>E</sup>* is the average transparency (m).

#### *2.2. Ecopath Model*

2.2.1. Principle of the Ecopath Model

The Ecopath is an ecological model that can directly determine the structure of an ecological system and describe its energy flow and mass transfer through the principle of nutrition dynamics. The Ecopath model stipulates that each function group (*i*) energy input and output in the ecosystem are balanced. The model uses a set of simultaneous linear equations to define an ecosystem, and each function group is represented by a linear equation [25] (Equation (2)).

$$B\_{\dot{l}} \times \left(\frac{P}{B}\right)\_{\dot{i}} \times EE\_{\dot{i}} = \sum\_{j=1}^{n} B\_{\dot{j}} \times \left(\frac{Q}{B}\right)\_{\dot{i}} \times DC\_{\dot{i}\dot{j}} + EX\_{\dot{i}} \tag{2}$$

In the formula, *Bi* is the biomass of the *<sup>i</sup>*-th function group. *<sup>P</sup> B i* is the ratio of the annual average production of the *i*-th function group to the annual average biomass, which is the biomass turnover rate. *EEi* is the ecological nutrient conversion efficiency of the *i*-th function group and is usually obtained by the calculation of the model [26]. *<sup>Q</sup> B i* is the ratio of the consumption of the *i*-th function group to the biomass. *DCij* is the ratio of the *i*-th the prey group to the total predation of the predator group *j*. *EXi* is the output of the *i*-th function group.

#### 2.2.2. Division of Function Groups

The purpose of the establishment of function groups is to merge populations with highly overlapping niches to simplify the food web. Based on the survey results for the Lianshi Lake ecosystem, characteristics of the community, and survival habits of each species, this work identified organisms with similar ecological characteristics, which were grouped into function groups [24]. A total of 14 function groups were set up (Table 1).



#### 2.2.3. Parameter Setting

Some of the parameters refer to lakes with water environmental conditions similar to those of Lianshi Lake, such as the production/biomass (P/B) and consumption/biomass (Q/B) parameters of fish resources that were obtained by querying (Available online: http://www.fishbase.org) (accessed on 20 May 2021). The plant P/B coefficient refers to related research on Taihu Lake [27]. The P/B coefficients of zooplankton (Zoop), other

benthos (Othb) and macrocrustaceans (Macc) were estimated based on the measured data by referring to the research results for Taihu Lake [27], Qiandao Lake [28], Dianshan Lake [29] and Zhushan Bay [30]. The P/Q coefficient refers to recognized data. The values for Zoop, OthB and Macc were 0.05 [31], 0.02 [32] and 0.075 [33], respectively. The food composition (*DCi*) is shown in Table 2. In addition to ecological investigation and research, we also referred to related research results [27,34], such as those for Zhushan Lake [30] and Qiandao Lake [35,36]. According to the cited literature, the *GS* values for general OthP and herbivorous fish (HerF) were set to 0.2 and 0.41, and the *GS* values for Zoop, OthB and Macc were set as 0.65, 0.94 and 0.7, respectively [31–33,37].


**Table 2.** Food composition data entered into the model.

2.2.4. Model Balance Calculation

The Ecopath model requires the input of six basic parameters: *Bi*, *P B i* , *EEi*, *Q B i* , *DCij* and *EXi*. Model balancing was executed, which was established on the basis of conforming to objective laws. The model parameters could be slightly modified to meet the requirements of the model operation, but it was necessary to avoid changing reliable data sources [25]. The nutritional conversion efficiency (*EEi*) ranged from 0 to 1 [38], the group *EE* was close to 1 in the face of considerable predation pressure, and the underutilized function group had a lower *EE* value. The value range of *GE <sup>P</sup> Q* was generally 0.1–0.3. If the balance test had one or more *EE* > 1, it was necessary to locate which predators caused the problem for specific prey groups in the predation mortality. The level of model confidence is mainly related to the quality and reliability of the acquired data. The accuracy of the Ecopath model was measured using the Pedigree index. The higher the index, the higher the quality of the model. The input and output parameters of the balanced model are shown in Table 3.

**Table 3.** Parameters of Lianshi Lake ecosystem construction.



**Table 3.** *Cont*.

#### *2.3. Fuzzy Comprehensive Evaluation Method*

The fuzzy comprehensive evaluation method is based on fuzzy mathematics. There is a certain characteristic of *n* things to be evaluated; these *n* things comprise the object set *X* = {*x*1, *x*2,..., *xn*}, the factor set *U* = {*u*1, *u*2,..., *un*} and the evaluation set *V* = {*v*1, *v*2,... *vm*}. Suppose that the weight distribution of the factors is the fuzzy subset A on V, denoted as *A* = {*a*1, *a*2,..., *an*}. In the formula, *ai* is the weight corresponding to the *i*-th factor *ui*, and the general rule is ∑*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *ai* = 1.

The evaluation steps for the fuzzy comprehensive evaluation method are as follows:

(1) Establish a factor set.

Assuming there are a total of *i* factors of the judged object, the factor universe of the evaluated object *U* is *U* = {*u*1, *u*<sup>2</sup> ..., *ui*}.

(2) Determine the membership function.

Assuming that the evaluation level is divided into *i* levels, the set is *V* = {*v*1, *v*2, *v*3,... *vi*}.


Determine the weight vector of the judgement factor *A* = (*a*1, *a*2,..., *an*). *A* is the subordination relationship of each factor in *U* to the thing being judged; it depends on the focus of people when making fuzzy comprehensive judgments, and it is equivalent to assigning weights according to the importance of each factor in the judgment.

(5) Establish a fuzzy comprehensive evaluation matrix.

According to the calculated maximum membership value, the program selection or category evaluation is performed.

#### **3. Results**

#### *3.1. Screening of Typical Biological Chains of Fuzzy Comprehensive Evaluation Coupled with Ecopath*

The Ecopath model was used to screen the key species in the Lianshi Lake ecosystem, identify the types of algae-controlling biological chains containing the key species, and establish a biological chain evaluation index system that could simultaneously express the delivery efficiency and stability. Finally, the fuzzy comprehensive evaluation method was used to screen out the biological chain with the largest weight value.

See Flowchart 2 for specific screening ideas (Figure 2).

**Figure 2.** Screening flowchart of a typical biological chain.

(1) Analysis of key species in the ecosystem

The key function group plays an important role in the function of the Lianshi Lake ecosystem. This study used the method developed by Libralato et al. [39] to calculate the keystoneness (KS) index of each function group. The one with the largest KS value is regarded as the key function group. Compared to the calculation method for the KS value proposed by Valls et al. [40] and Power et al. [41], this calculation is simpler and comprehensively considers the effects of biomass and capacity. The equation is *KSi* = *log*[*εi*(1 − *pi*)], ε*<sup>i</sup>* = ∑*n <sup>j</sup>*=<sup>1</sup> *<sup>m</sup>*<sup>2</sup> *ij*ε,*pi* <sup>=</sup> *Bi* ∑*n <sup>i</sup> BK* , where *KSI* is the criticality index of function group *i*, *ε<sup>i</sup>* is the total impacts of function group *i* in the ecosystem (total impacts), and *mij* is the value of the mixed nutrition effect of function group *i* on function group *j*, indicating the mutual relationship between each other's strengths, and *Pi* is the ratio of the biomass *Bi* of function group *i* to the biomass of the entire ecosystem ∑*<sup>n</sup> <sup>k</sup> Bk*. Since *KSi* and *pi* are negatively correlated, the criticality index will not be too high due to the high biomass of the function groups (Table 4).

**Table 4.** Criticality data table for each function group of Lianshi Lake.


As shown in Figure 3, the first key function groups in the Lianshi Lake ecosystem are submerged macrophytes (SubM). However, other macrophytes (OthM) are redundant in the food web of the entire ecosystem. Other benthos (OthB), Phyt, other fish (OthF) and Zoop have a relative total impact second only to SubM, and their criticality index ranks in the top four.

(2) Screening of target biological chains for alga control and salt reduction

With the goal of alga control, we identified the biological chain with Phyt and SubM as primary producers, and then determined the biological chain containing key species such as Zoop, OthB and OthF. The results are shown in Figure 4.

(**a**) Alga control and saltƺcutting biological chain containing key species of OthB.

**Figure 4.** *Cont*.

(**b**) Alga control and saltícutting biological chain containing key species of OthF.ȱ ȱ

(**c**) Alga control and saltícutting biological chain containing key species of Zoop.ȱ ȱ

**Figure 4.** The biological chain of algal control in Lianshi Lake containing key species.

(3) Establishment of evaluation indicators for typical biological chains

In order to express the characteristics of the Lianshi Lake ecosystem, the ecotrophic efficiency (*EE*), production/respiration (*P*/*R*), relative total impact (*RT I*), omnivory index (*OI*) and trophic level (*TL*) are listed as typical indicators for biological chain evaluation (Table 5). The evaluation indicators have the following characteristics: (1) *EE* stands for the utilization and conversion efficiency for the energy contained in the previous trophic organism. The value of EE ranges from 0 to 1, and it is close to 1 in a group with considerable predation pressure. (2) *P*/*R* represents an important indicator of the maturity of the biological chain, which is close to 1 in a mature ecosystem. (3) *RT I* represents the relative total impact of a single function group on the entire ecosystem. (4) *OI* indicates that, when the omnivorous index is zero, the corresponding prey is oriented. (5) *TL* refers to the level of the organism in the food chain of the ecosystem. The higher the trophic level of the organism, the greater the contribution of the food chain to the stability of the ecosystem.


**Table 5.** Basic data of evaluation indicators for each function group of Lianshi Lake.

(4) Fuzzy comprehensive evaluation method for screening typical biological chains

The research did not consider the comment-level domain in the fuzzy comprehensive evaluation, and selected the typical biological chain according to the maximum weight (from Equation (3) to Equation (6)):

<sup>1</sup> The factor domain *U*, *U* = (*EE*, *P*/*R*, *RT I*,*OI*, *TL*) of the judged object was determined.

<sup>2</sup> The membership function was established.

The higher the value, the higher the membership function:

$$r\_i = \begin{cases} 0 & (x \le S\_1) \\ \frac{x - S\_1}{S\_2 - S\_1} & (S\_1 < x < S\_2) \\ 1 & (x \ge S\_2) \end{cases} \tag{3}$$

(I) The determination of the membership functions of *EE*, *P*/*R*, *RT I*, and *OI* As the variation range of *EE*, *P*/*R*, *RT I*, and *OI* was between 0 and 1, the membership function was:

$$\mathcal{U}\_A(\mathbf{x}) = \mathbf{x} \tag{4}$$

(II) The determination of the TL membership function

According to the literature, the lowest trophic level of primary producers in lake ecosystems is 1, and the trophic level of the function group that is generally at the top is 4, so the membership function is:

$$\mathcal{U}\_A(\mathbf{x}) = \frac{\mathbf{x} - 1}{\mathfrak{Z}} \tag{5}$$

<sup>3</sup> Single factor evaluation was performed, and a fuzzy relationship matrix was established.

The fuzzy evaluation matrix was established according to the following formula:

$$R = \begin{bmatrix} r\_{ij} \end{bmatrix} = \begin{bmatrix} r\_{11} & \cdots & r\_{15} \\ \vdots & \ddots & \vdots \\ r\_{51} & \cdots & r\_{55} \end{bmatrix} \tag{6}$$

In the formula, in the *i*-th row *Ri* = (*ri*1,*ri*2,...*rim*), *i* = 1, ... *m*, *i* is the degree of the membership of the *i*-th evaluation factor to the evaluation standards at all levels; for the *j*-th column *Rj* = *r*1*j*,*r*2*j*,...,*rnj* , *j* is the degree of the membership of each evaluation factor to the *j*-th evaluation standard. Taking the nine biological chains containing the Zoop group as an example, the membership degrees corresponding to different evaluation indicators of different biological chains were calculated and are shown in Table 6 below.


**Table 6.** The membership degrees corresponding to different evaluation indicators of different biological chains with the Zoop group.

The fuzzy relation matrix with the Zoop group:


<sup>4</sup> A weight vector was established.

The weight vector *A* = (*a*1, *a*2,..., *an*) of the judgement factor was determined. *A* is the subordination relationship of each factor in *U* to the thing being judged. Due to the importance of the distribution of the weight, this study used the analytic hierarchy process to determine the weight of each evaluation index; the hierarchical structure is shown in Figure 5.

**Figure 5.** Hierarchical structure model.

$$\begin{array}{rcl} \text{The jointement matrix } T & = & \begin{bmatrix} 1 & 4 & 0.33 & 0.143 & 0.25 \\ 0.25 & 1 & 0.25 & 0.143 & 0.2 \\ 3 & 4 & 1 & 0.5 & 3 \\ 7 & 7 & 2 & 1 & 8 \\ 4 & 5 & 0.33 & 0.125 & 1 \end{bmatrix} \\ \text{The calculated matrix} \end{array}$$

mum characteristic value is *λmax* = 5.42, the random consensus ratio is 0.0944 < 0.1, and the consistency test shows that the weight distribution is reasonable, Therefore, the weight of each evaluation index of a typical biological chain is *A* = (0.0712, 0.0388, 0.2224, 0.0534, 0.1336).

<sup>5</sup> The comprehensive evaluation results were analyzed.

$$B = AR = (0.2582, 0.2825, 0.2839, 0.2950, 0.2798, 0.2581, 0.2940, 0.2368, 0.2663)$$

It can be seen that the typical biological chain containing the Zoop group as the key species is <sup>4</sup> , equal to phytoplankton (Phyt)–zooplankton (Zoop)–macrocrustaceans (Macc)–other piscivorous (OthP).

According to the same calculation principle, the typical algal control biological chain containing OthB and OthF as the key species does not exist.

#### *3.2. Analysis of Typical Algal Control Biological Chain and Nutritional Structure and Function Response*

Through modeling, it was found that the production of primary producers in the Lianshi Lake ecosystem was 10,661.7 t/km2/year, of which the amount of food consumed was 3043 t/km2/year, and the remaining 71.5% was not consumed by other predators but flowed into the debris, entered recirculation, and became deposited through mineralization. Therefore, the effective removal of plant organisms from the water body plays an important role in controlling the total amount of nutrients, dissolved oxygen and transparency of the water body.

In the Lianshi Lake ecosystem, there are 16 biological chains with Phyt as the primary producer. Among them, the number of biological chains where Phyt is directly grazed by Zoop is the largest; it has nine chains. The transfer efficiency of Phyt is 0.308. The transfer efficiency of Zoop reached 0.606. It can be seen that the types of biological chains that impose significant biological density constraints on algae are mainly Zoop predation on Phyt. Combined with the screening of typical biological chains, a typical biological chain of algal control in the Lianshi Lake ecosystem is phytoplankton (Phyt)–zooplankton (Zoop)–macrocrustaceans (Macc)–other piscivorous (OthP).

This study improved the efficiency of nutrient delivery for primary producers by simulating an increase in biomass of Zoop and Macc. Under the guidance of biological manipulation theory, the key factor in the regulation of the algal control food chain is Zoop. Through the introduction of *Daphnia magna*, the food chain is opened up through insects eating algae and fish-eating insects. A symbiotic system of "*Daphnia magna*–underwater forest–aquatic animals–microbes" was constructed, and the "grass-type clean water state" self-purification system was restored.

The proliferation of ecological capacity continuously increases the biomass of the target species. By observing the changes in other function groups such as the eaten organisms in the model, when the ecological nutrition transfer efficiency of a function group in the model *EE* > 1, the model will become unbalanced; the biomass value of the released species before the imbalance of the model is the ecological capacity. The biomass of Zoop and Macc have been increased to 1.5, 2, 3 and 4 times, respectively. It can be seen from Table 7 that, within the ecological capacity of Zoop and Macc, with increases in the biomasses of both, the ecological nutrition transfer efficiency of Phyt gradually increased, from 0.308 (the current situation) to 0.458 (1.5 times), 0.607 (2 times) and 0.906 (3 times). The flow of material flowing into the second trophic level from primary producers also increased from 3043 (the current situation) to 4353 (1.5 times), 5663 (2 times) and 8283 t/km2/year (3 times).



The amount of debris flowing into the primary producers for deposition continuously decreased, from 7618 (the current situation) to 6308 (1.5 times), 4998 (two times) and

2378 t/km2/year (three times). The total primary production/total respiratory volume (TPP/TR) continued to decrease from 9.224 (the current situation) to 6.461 (1.5 times), 4.972 (two times) and 3.403 (three times). The Finn's circulation index (FCI) continued to rise from 13.59% (the current situation) to 14.33% (1.5 times), 15.23% (two times) and 17.51% (three times), whereas the Finn's average energy flow path length (FCL) rose from 2.854 (the current situation) to 2.993 (1.5 times), 3.132 (two times) and 3.410 (three times).

Based on the analysis of the nutritional structure of the Lianshi Lake ecosystem, the artificial introduction of Zoop (*Daphnia magna*) and Macc can increase the transfer efficiency and the maturity of the ecosystem. To a certain extent, it can solve the problem of excessive primary production.

#### **4. Discussion**

#### *4.1. Development Characteristics of Lianshi Lake Ecosystem*

The total primary production/total respiratory volume (TPP/TR) is an important evaluation index, which is close to 1 in mature ecosystems, far greater than 1 in developing ecosystems, and less than 1 in polluted ecosystems. The Finn's cycling index (FCI) is the ratio of the circulation flow to the total flow in the system, and the Finn's mean path length (FMPL) is the average length of each circulation through the food chain. The higher the ratio of material recycling, the longer the food chain through which the nutrient flows, and the FCI value of a mature ecosystem is close to 1. The current Lianshi Lake TPP/TR is 9.224, the FCI is 13.59%, and the FMPL is 2.854, indicating that the Lianshi Lake ecosystem is immature.

In the typical biological chain of algal control, with an increase in the biomass of Zoop (*Daphnia magna*) and Macc within the ecological capacity, the amount of phytoplankton as a primary producer flowing into the next trophic level gradually increases, the transfer efficiency gradually increases, the FCI and FMPL gradually increase, and the TPP/TR gradually approaches 1, which reduces the risk of lake eutrophication and increases ecosystem maturity to a certain extent.

Therefore, there may be two reasons for the low maturity of the Lianshi Lake ecosystem: first, the biomass of the function group that plays a key role in the Lianshi Lake ecosystem is much lower than the ecological capacity, resulting in insufficient driving force for the ecosystem to develop to a mature state. Second, the biodiversity of Lianshi Lake is low, and the flow of energy to higher levels is hindered. It is recommended to introduce Zoop (*Daphnia magna*) and indigenous herbivorous fishes to build a food chain in order to promote material and energy cycles.

#### *4.2. Prospects of Ecopath Model in Lake Ecological Restoration*

As early as 1975, Shapiro et al. [17] proposed the biological manipulation theory. Biomanipulation methods have been developed for nearly 50 years; there have been many reports in Western European countries, but this technology has not been widely promoted. This may be because biological manipulation methods involve complex biological networks, and too many factors are affected. Traditional research methods can only study the behavior of individual organisms in simple habitats or competitive environments, and there is very little research on the biological chain that significantly affects the regulation goals and ecosystem conditions. Therefore, it is generally believed that research on complex ecosystems must rely on the guidance of mathematical models or theories [42].

The Ecopath with Ecosim (EWE) model, also known as the ecological channel model, is an ecological model that can assess the true structure of the ecosystem and describe its energy flow and mass balance. The Ecopath model was originally created in 1984 [43]. After years of development, the Ecopath model has become a key tool for ecosystem research. Fuzzy comprehensive evaluation (FCE) is based on fuzzy mathematics [44], which can express the fuzzy relationship between various factors and solve the problem of ambiguity between multi-factor evaluations that cannot be solved by traditional methods, and it has been widely used in the field of policy evaluation and risk assessment. A typical

biological chain of alga control combines the advantages of the Ecopath model and FCE method, and the results show that the biomass of phytoplankton in Lianshi Lake has been effectively controlled and that the ecosystem's maturity has been improved.

Ecopath is a powerful model but is mostly used to assess the condition of the ecosystem and provide scientific management and control solutions. In the future, the Ecopath model will continue to be developed and coupled with other models, such as pollutant diffusion models and ecotoxicology models, which will have important scientific research significance for exploring the restoration of lake ecosystems.

#### **5. Conclusions**

Based on the survey results for aquatic organisms in Lianshi Lake from 2018 to 2019, this study established, for the first time, a typical biological chain screening method with fuzzy comprehensive evaluation coupled with Ecopath. Among the organisms, Phyt is the primary producer, and the typical biological chain of alga control with Zoop as the key species is phytoplankton (Phyt)–zooplankton (Zoop)–macrocrustaceans (Macc)–other piscivorous (OthP).

In a typical biological chain with significant biological density constraints on algae, when the biomass of Zoop was increased from 7.85 (the current situation) to 23.55 t/km2/year (three times) and Macc was increased from 1.58 (the current situation) to 4.74 t/km2/year (three times), the results show that the ecological nutrition efficiency of Phyt increased from 0.308 (the current situation) to 0.906 (three times), the material flow into the second trophic level from primary producers increased from 3043 (the current situation) to 8283 t/km2/year (three times), the amount of debris flowing into primary producers for sedimentation decreased from 7618 (the current situation) to 2378 t/km2/year (three times), the total primary production/total respiratory volume (TPP/TR) decreased from 9.224 (the current situation) to 3.403 (three times), the Finn's cycle index (FCI) increased from 13.59% (the current situation) to 17.51% (three times), and the Finn's average energy flow path length (FCL) increased from 2.854 (the current situation) to 3.410 (three times). In the typical biological chain of alga control, the artificial release of Zoop (*Daphnia magna*) and Macc can improve the transfer efficiency of phytoplankton as primary producers to a certain extent, reduce the harm caused by eutrophication to lake ecosystems, and improve the maturity of the lake ecosystem.

**Author Contributions:** P.Z. and X.C. designed the research content; P.Z., H.L. and W.P. processed the data and analyzed the results. All the authors participated in the writing of the manuscript. Y.G. and X.C. revised the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the Major Science and Technology Program for Water Pollution Control and Treatment (No.2018ZX07101005).

**Acknowledgments:** We thank Min Zhang for providing hydrobiological data and Su Wei for assistance with the fieldwork. We would like to thank Guanglei Qiu for his guidance on the translation of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Jianwei Wang 1, Nengzhan Zheng 1, Hong Liu 1, Xinyi Cao 1, Yanguo Teng 1,2,\* and Yuanzheng Zhai 1,\***


**Abstract:** Songnen Plain is one of the three great plains in northeast China with abundant groundwater resources. The continuous population growth and the rapid development of agriculture and economy in China has caused a series of environmental problems in the plain, such as endemic diseases caused by the accumulation of harmful substances in drinking water. This paper conducts a systematic investigation of fluorine in the groundwater of Songnen Plain. The results showed that fluorine was widespread in the groundwater of the plain in the concentration range of BDL–8.54 mg·L<sup>−</sup>1, at a mean value of 0.63 mg·L−<sup>1</sup> and detectable at a rate of 85.91%. The highest concentrations of fluorine were found in central and southwest areas of the plain. The concentration exceeded the guideline values for fluorine in drinking water and may have varying degrees of adverse effects on adults, and especially children, in the study area. The fluorine in groundwater mainly came from the dissolution of fluorite and other fluorine-containing minerals, and the concentrations and distribution of fluorine were affected by cation exchange, groundwater flow field and hydrochemical indexes (pH, TDS and HCO3 −). The study provides scientific basis for the investigation, evaluation and prevention of endemic diseases caused by fluorine.

**Keywords:** Songnen Plain; groundwater; fluorine; distribution; formation; human health risk

#### **1. Introduction**

Water is an important natural resource to guarantee normal human life and socioeconomic development. As one of the main components of water resources, groundwater is exerting a more and more significant influence on society, and groundwater is also a prerequisite for the development of other resources, especially in arid and semi-arid areas [1]. However, the accumulation of harmful elements caused by geological causes or human activities not only harms the water environment, but also seriously affects the safety of drinking water [2,3]. Fluorine is an essential element of the human body and moderate intake (0.5–1.5 mg·L<sup>−</sup>1) is beneficial to human health, according to the WHO guidelines [4–6]. However, excessive fluoride, after long-term consumption may lead to fluorosis, and it is also a serious problem for the world's geological environments [2,7]. Studies have shown that over 260 million people are at risk of fluorosis all over the world, reported in locations such as America, Argentina, China, Mexico, and India [8–13]. Therefore, the research on high fluorine groundwater have gradually become a research focus.

Drinking water is the primary route by which fluoride enters the human body [14]. The National Health Commission of the PRC [15] and the WHO [16] guidelines have set values for fluorine in drinking water at 1.0 and 1.5 mg·L−1, respectively. Fluorine is widespread in natural minerals, such as flourite, cryolite, fluorapatite, etc. [8,9,17]. Studies have shown that the main reason for the formation of high-fluorine groundwater in many regions of the world is dissolution of fluorine-bearing minerals [18–20]. However, the

**Citation:** Wang, J.; Zheng, N.; Liu, H.; Cao, X.; Teng, Y.; Zhai, Y. Distribution, Formation and Human Health Risk of Fluorine in Groundwater in Songnen Plain, NE China. *Water* **2021**, *13*, 3236. https://doi.org/10.3390/w13223236

Academic Editor: Fausto Grassa

Received: 16 September 2021 Accepted: 6 November 2021 Published: 15 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

distribution, formation and risk of high-fluorine groundwater in Songnen Plain, China have not been systematically studied.

Songnen Plain is one of the largest and most fertile plains in China, and an important agricultural production base [21,22]. The plain has a large groundwater aquifer system with multiple aquifers and abundant groundwater resources [23]. With the continuous growth of population and the rapid development of agriculture and economy, the contradiction between supply and demand of water resources is becoming more and more prominent in the area [24], and has caused a series of environmental problems, such as metal pollution [25], nitrate pollution [26,27], and endemic diseases caused by excessive content of fluorine, iodine and arsenic [28,29].

In this paper, the high content of fluorine in the groundwater of Songnen Plain is investigated systematically. Through a series of processes such as field sampling, index determination, sample preservation, pretreatment, detection and data analysis, the concentration, distribution, formation and human health risk of fluorine in groundwater in Songnen Plain are revealed. This provides scientific basis for the investigation, evaluation, and prevention of endemic diseases caused by fluorine.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Songnen Plain is in the northeastern part of China and located at longitude 121◦21 – 128◦18 and latitude 43◦36 –49◦26 . Songnen Plain is an important grain commodity production base and animal husbandry base of China. It covers an area of 103,200 km2, and is placed in the Songhua River basin. The plain evolved from the Mesozoic-Cenozoic faulted basin and has accumulated over 8 km of Cretaceous terrestrial clastic deposits. Gravel, sand and loam are the main components of strata in the study area, and cohesive soil interlayers are locally distributed [1,26]. Analysis of the strata minerals revealed that fluorine–bearing minerals are rich in the central and southwest strata of the plain, mainly including flourite, apatite, cryolite, topaz, biotite, hornblende, tourmaline [26]. The regional groundwater resources are abundant and the largest groundwater system of the entire aquifer includes Neogene fissure–pore water, Cretaceous pore–fracture water, Quaternary pore water and Paleogene fissure-pore water. Irrigation and precipitation constitute the main sources of local groundwater recharge [30]. The shallow groundwater (depth is less than 50 m) are greatly affected by anthropogenic activities, complicated and changeable chemical composition. With the rapid development of agriculture and industry in northeastern China, groundwater exploitation in this region is expanding and, coupled with decreasing precipitation, the water table is declining and the groundwater environment in the study area changed greatly.

#### *2.2. Sampling*

In this paper, a comprehensive groundwater pollution survey was carried out in Songnen Plain from 2012 to 2014. Sampling time was concentrated in May to October each year, due to the cold winters in northeast China. Groundwater sample collection in the study area relied on local mechanized wells. A total of 2683 groundwater samples were collected; their locations are shown in Figure 1. Prior to groundwater collection, the original well water was pumped more than three times to flush the well's pump [3]. The sampling bottles (500 mL) were made of polyethylene plastic and were soaked in a 10%-sodium hydroxide solution for 3 h [11], then cleaned with deionized water and distilled water in turn, and finally dried at 60 ◦C for 5 h and stored in ziplock bags. Additionally, 2 mL of concentrated nitric acid (1:1) was added to the sampling bottle for measuring heavy metals [26]. Then, 2 mL of concentrated sulphuric acid (1:1) was added to the sample bottles for the measurement of Fe and Mn [1]. The sample bottles were washed three times with the corresponding water before each sampling. Each water sample collected was refrigerated at −4 ◦C, and handled within 48 h. Unstable parameters, such as water temperature (T), pH, electrical conductivity (EC) and also water-table depth were measured in situ.

**Figure 1.** Location of sampling points in the study area.

#### *2.3. Analysis Methods*

#### 2.3.1. Data Analysis Method

In this paper, a variety of data analyses are reported; a hydrochemical analysis, simulation methods, spatial analysis and mapping software are used herein toprocess the large volume of groundwater-sample testing data collected. The comprehensive parameters that can reflect groundwater characteristics and kriging method in statistics are selected as spatial interpolation method. Based on variogram theory and structural analysis, the comprehensive parameter zoning maps of groundwater at each layer are drawn with the help of ArcGIS and Surfer professional mapping software. The SPSS software was used to carry out descriptive statistical analysis on the main components of the sampled groundwater. SUPCRTBL and PHREEQC were used to calculate the saturation index (SI) of various rocks in combination with the latest database.

#### 2.3.2. Instrumental Analysis

The instrumental analysis method used for the determinations of the water samples is referenced from a previous study [1,26]: the pH and redox potential were determined by dual-channel multi-parameter water quality analyzer (HQ40D, Field Case, cat. No:58258-00, HACH, Loveland, CO, USA); K, Ca, Na, Mg, Mnand Fe concentrations were measured by ICP-AES (IRIS Intrepid II XSP, Thermo Scientific, Waltham, MA, USA); the concentrations

of Cl− and NO3 − in the water samples were measured by ion chromatograph (Dionex2500, Dionex, Sunnyvale, CA, USA).

#### 2.3.3. Saturation Index

The saturation index (SI) is, itself, derived from the formula of a theoretical derivation and used to describe the equilibrium state of the water with respect to mineral phases therein and to determine the dissolution or precipitation of minerals in rock–water interactions [31,32]. The index was calculated by Equation (1):

$$SI = \log\left(\frac{IAP}{Ks}\right) \tag{1}$$

where *IAP* is the ion activity product of the solution and *Ks* is the solubility product of the mineral. Different *SI* values indicate different states of ion in solution: *SI* > 0 indicates oversaturation (precipitation), *SI* = 0 indicates equilibrium and *SI* < 0 indicates undersaturation (dissolution).

#### 2.3.4. Human Health Risk Assessment

A human health risk assessment consists in determining potential adverse effects of a target pollutant. The health risk assessment model (RBCA) was created to assess non-carcinogenic risks. We have mainly referred to oral exposures to pollutants in this study because the mouth is considered the primary route thereof. The calculations of non-carcinogenic risk (hazard quotient, HQ) of directly consuming water resources were extrapolated from the oral reference dose (*RfD*), hazard index (HI) represents the total non-carcinogenic risk to humans when the ADI (average daily intake) was unavailable, as shown in Equation (2). HI > 1 indicates that the exposed individual was adversely affected.

$$\text{HQ} = \frac{EDI}{RfD} = \frac{\text{CS} \times IR \times EF \times ED}{BW \times AT \times RfD} \quad \text{and} \quad \text{HI} = \sum\_{i=1}^{n} HQ\_i \tag{2}$$

where *CS* (mg·L−1) is the concentration of OPPs in the water; *IR* is the average daily water intake (1.5 and 0.7 L·d−<sup>1</sup> for adults, children, respectively); *EF* stands for exposure frequency (365 d·y<sup>−</sup>1). The *ED*s (exposure durations) for children and adults were 12 and 30 years, respectively; *BW*s (body weights) for children and adults were 10 and 60 kg, respectively [16]. *AT* represents average lifetime, and was 4380 and 1095 days for adults and children, respectively. *RfD* stands for the reference dose of the carcinogen consumed orally. The value of the *RfD* for fluoride was 0.04 mg·kg−1·day−<sup>1</sup> [33,34].

#### *2.4. Quality Assurance/Quality Control (QA/QC)*

In order to make sure the accuracy of the measurement results, blank samples and standard samples were taken from each batch during sampling. After analysis, the standard deviations ranged from 0.09% to 0.23% and the target pollutant was not detected in the blank samples, which conformed to the stated data processing standards [35,36]. In determining groundwater quality parameters, the recovery indicator was added before the water sample was processed; samples' recovery rates ranged from 84.4 to 95.7%, and the average was 89.4%; moreover, we compiled the standard curve for target objects and the results of the analysis show that the correlation coefficient of the linear equation was over 0.99; all conformed to quality assurance standards for the processing of groundwater [37,38].

#### **3. Results and Discussion**

#### *3.1. Hydrochemical Parameters and Types*

The concentrations of hydrochemical parameters are shown in Table 1. The main cations in groundwater in the study area were Ca and Na, the concentration of Ca ranged from 1.56 to 567.65 mg·L<sup>−</sup>1, with an average of 95.04 mg·L−<sup>1</sup> and the concentration of Na ranged from 4.51 to 1107.36 mg·L−1, with an average of 74.45 mg·L−1. According to the

groundwater index detection values, the Piper diagram of groundwater hydrochemistry types in the study area was drawn in Figure 2. Generally speaking, the main groundwater chemical type in the study area was HCO3–Ca type [1], accounting for 24.83% of the total water samples. Other main groundwater types include HCO3–Ca·Mg type, HCO3–Na·Ca type, HCO3–Na·Mg·Ca type and HCO3·Cl–Ca type, accounting for 19.43%, 16.78%, 13.92%, 12.11% and 9.56% of the total water samples, respectively.


**Table 1.** Concentrations of hydrochemical parameters and saturation indexes of minerals.

Unit: BDL = below detection limit.TH: total hardness; TDS: total dissolved solids.

**Figure 2.** Piper diagram of groundwater samples.

#### *3.2. Distributions of Fluorine*

The concentration range of fluorine in groundwater was BDL–8.54 mg·L−1, with a mean value of 0.63 mg·L<sup>−</sup>1, and its detection rate was 85.91%. The content distributions of fluorine are showed in Figure 3; in the study area's groundwater the highest concentrations of fluorine (over 2 mg·L<sup>−</sup>1) were found in the central and southwest areas of the Songnen Plain, such as Tongyu, Qianan, Baicheng, Lindian, Daqing, Zhaozhou, Zhaodong—and, the observed concentrations exceeded the maximum fluorine content (1.5 mg·L−1) that is beneficial to human health [6] and within the guideline values set by the National Health Commission of the PRC (1.0 mg·L<sup>−</sup>1) [15] and the WHO (1.5 mg·L<sup>−</sup>1) [16].

**Figure 3.** Distribution of fluorine in groundwater.

#### *3.3. Formation and Influencing Factors of Fluorine in Groundwater*

#### 3.3.1. Dissolution and Precipitation of Minerals

Research has shown that the dissolution of fluorine-bearing minerals and the precipitation of calcium-bearing minerals are the main influencing factors of F− enrichment in groundwater [8,20,39,40]. Analysis of the strata minerals in the study area revealed an abundance of fluorine-bearing minerals in the central and southwest strata of the plain, mainly flourite, apatite, cryolite, topaz and hornblende [1,26]. The saturation indices of SI fluorite in almost all groundwater samples in the study area were less than zero, and there was a significant positive correlation between F<sup>−</sup> concentration and SI fluorite (Figure 4a), suggesting that the dissolution of fluorite is the main source of F− in the groundwater of these areas. According to other mineral saturation indices (SI Calcite < 0, SI Halite < 0, SI Dolomite < 0), calcite, halite and dolomite had not reached the saturation state and were

easy to dissolve in reaction. SI fluorite had a logarithmic increase, with an increasing concentration of F− in groundwater (Figure 4a). Where fluorite tended to saturate, the concentration of F− reached the upper limit, indicating that the concentration of F− was restricted by the equilibrium constant of fluorite (Ksp = 10<sup>−</sup>10.059, 22 ◦C). By comparing the concentration relationship between Ca2+ and F<sup>−</sup> (Figure 4b), fluorine in the groundwater samples below the fluorite dissolution curve (dotted line in Figure 4b) mainly came from the dissolution of fluorite, and the fluorine in the groundwater samples above the dissolution curve came not only from the dissolution of fluorite, but also from other sources. The results show that the dissolution of fluorine-bearing minerals is main reason for the deposition of significant fluorine in the groundwater of Songnen Plain—similar to conditions found elsewhere in China [14,41] and the world, such as America [11], Mexico [8] and India [10].

**Figure 4.** Relationship between F<sup>−</sup> concentration and SI(Fluorite) (**a**); Ca2+ concentration (**b**).

#### 3.3.2. Cation Exchange

The formation of groundwater hydrochemistry is often closely related to cation exchange [40], and cation exchange is the significant factor affecting the formation of fluorine in this study area. Cation exchange was confirmed with chloro-alkaline CAI 1 and CAI 2, and the indices were calculated by Equations (3) and (4).

$$\text{CAI 1} = \frac{\text{CI}^- - \left(\text{Na}^+ + \text{K}^+\right)}{\text{CI}^-} \tag{3}$$

$$\text{CAI 2} = \frac{\text{Cl}^- - \left(\text{Na}^+ + \text{K}^+\right)}{\text{HCO}\_3^- + \text{SO}\_4^{2-} + \text{CO}\_3^{2-} + \text{NO}^{3-}}\tag{4}$$

where, if CAI 1 > 0 and CAI 2 > 0, it is indicated that the dissolved Na+ and K+ in the groundwater will exchange cations with the absorbed Mg2+ and Ca2+. However, when less than zero, it is indicated that the dissolved Mg2+ and Ca2+ will exchange with the absorbed Na<sup>+</sup> and K+. Moreover, the greater the absolute value, the stronger the cation exchange. Figure 5 shows the CAI 1 and CAI 2 of the groundwater samples. All values of CAI 2 do not exceed zero, and most values of CAI 1 were negative. This suggested that the cation exchange process of dissolved Mg2+ and Ca2+ exchanging cations with the absorbed Na<sup>+</sup> and K<sup>+</sup> was the driving process explaining local mineral concentrations, and it is also responsible for the decreased contents of Mg2+ and Ca2+ in the groundwater. This process may promote the hydrolysis of fluorite and other fluorine minerals (including apatite, cryolite, topaz, hornblende, tourmaline, etc.), thereby increasing the fluorine in

groundwater. This indicates that the cation exchange process can affect the fluorine content in groundwater, as is consistent with previous studies [12,19,40,42].

**Figure 5.** Relationship between CAI 1 and CAI 2.

3.3.3. Hydrochemical and Hydrological Influence Factors

Hydrochemical parameters are also one of the important factors affecting the content and distribution of fluorine [43,44]. It has been found that fluorine accumulates more easily in an alkaline environment [18], therefore, the areas with the highest pH values (pH > 8.5) in the study area (shown in Supplementary materials, Figure S1) also had the highest F− contents. Correlation analysis between fluorine and other hydrochemical parameters in groundwater showed (Figure 6) that the concentrations of F− in groundwater were positively correlated with the concentrations of TDS (total dissolved solids) and HCO3 −. The results indicated that the concentrations and distribution characteristics of F− in groundwater were closely related to the pH of the groundwater environment and the concentrations of TDS and HCO3 −, again, as is consistent with previous studies [44]. Correlation analysis of fluorine and other high concentration pollutants (I−, Mn2+, Fe) in the groundwater showed (Figure S2) that fluorine was only weakly correlated with iodine (having similar properties), indicating that the pollutants in the groundwater had little influence on each other.

Hydrological conditions can partly affect the concentration and distribution of fluorine in groundwater [1,41]. The shallow groundwater system of Songnen Plain belongs to a larger groundwater catchment basin, and groundwater gather in its central low plain [16,29]. Therefore, the groundwater in the surrounding areas, especially the high-fluorine groundwater in the southwest area, will gradually migrate to the central region through the groundwater flow field, further increasing the fluorine concentrations in the already-high-fluorine groundwater in the central plain area.

#### *3.4. Human Health Risk Assessment of Fluorine*

Health risk assessments are mainly concerned with oral exposures; to that end, the risk assessment of non-carcinogens performed was based on the concentrations of fluorine in groundwater, and the results are shown in Figure 7. The sampling points (HQ Children > 1) accounted for 92.26% of the total samples, signaling that the groundwater fluorine concentration is

high enough to have significant adverse effects on children in the study area. The sampling points (HQ Adults > 1) accounted for 32.18% of the total samples, indicating that the fluorine content was also high enough to adversely affect adults, though much less so than children. In addition, the districts where the fluorine in the groundwater showed the greatest potential influence on children and adults were roughly the same, and were concentrated in the southwest and central Songnen Plain, such as Tongyu Lindian and Daqing.

**Figure 6.** Correlation of F− and Hydrochemical Indices.

**Figure 7.** HQs of fluorine to children and Adults.

#### **4. Conclusions**

Fluorine is widespread in the groundwater of the Songnen Plain, at a concentration range of BDL–8.54 mg·L−1, with a mean value of 0.63 mg·L−<sup>1</sup> and detectable at a rate of 85.91%. The highest concentrations of fluorine (over 2 mg·L−1) were found in the central and southwest areas of the plain. The concentrations there exceeded the guideline values for fluorine in drinking water set by both the National Health Commission of the PRC (1.0 mg·L−1) and the WHO (1.5 mg·L−1), and represent varying degrees of adverse effect on adults, and especially children, in the study area. The fluorine in these groundwaters mainly came from the dissolution of fluorite and other fluorine-containing minerals in the study area; additionally, the concentrations and distribution of fluorine were shown to be affected by cation exchange, the groundwater flow field and hydrochemical indexes (pH, TDS and HCO3 −). The study provides scientific basis for the investigation, evaluation and prevention of endemic diseases caused by groundwater fluorine.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/w13223236/s1, Figure S1: Distribution of pH in groundwater, Figure S2: Correlation of F− and other pollutants.

**Author Contributions:** Data curation, X.C.; Formal analysis, H.L.; Funding acquisition, Y.T. and Y.Z.; Investigation, J.W.; Methodology, J.W.; Resources, Y.T.; Software, N.Z.; Writing—original draft, J.W.; Writing—review & editing, Y.T. and Y.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by National Natural Science Foundation of China (No. 41877355; 42077170); Major Science and Technology Program for Water Pollution Control and Treatment (2018ZX07101005-04); Beijing Advanced Innovation Program for Land Surface Science of China; the 111 Project of China (B18006).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** Authors highly appreciate the detailed and constructive comments from the editors and reviewers to help us improve the quality and English expression of the paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

