Spatial Pattern Evolution and Influencing Factors on Agricultural Non-Point Source Pollution in Small Town Areas under the Background of Rapid Industrialization

To promote sustainable agricultural development in small town areas during rapid industrialization, it is important to study the evolution of agricultural non-point source pollution (ANSP) and its influencing factors in small town areas in the context of rapid industrialization. The non-point source inventory method was used to study the characteristics of ANSP evolution in 14 small town areas in Gongyi City from 2002 to 2019. Using the spatial Durbin model and geographical detectors, the factors influencing ANSP in small town areas were analyzed in terms of spatial spillover effects and the spatial stratified heterogeneity. The results showed a zigzagging downward trend of ANSP equivalent emissions over time. Spatially, the equivalent emissions of ANSP showed a distribution pattern of being high in the west and low in the east. There was a significant positive global spatial autocorrelation feature and there was an inverted “U-shaped” Environmental Kuznets Curve relationship between industrialization and ANSP. Affluence, population size, and cropping structure positively contributed to the reduction of ANSP. Population size, land size, and industrialization were highly influential factors affecting the spatial variation of ANSP and the interaction of these factors was bivariate or nonlinearly enhanced. This study provides a feasible reference for policymakers and managers to develop reasonable management measures to mitigate ANSP in small town areas during rapid industrialization.


Introduction
Agricultural non-point source pollution (ANSP) is prevalent globally, affecting 30% to 50% of the earth's surface [1]; therefore, the study of ANSP has been a hot topic of academic research. Since the implementation of China's reform and opening-up policy in 1978, industrialization has accelerated technological progress in agricultural production and large-scale mechanized agricultural production has been realized, leading to a substantial increase in agricultural production efficiency [2]. China's agricultural economy has progressed remarkably, with the production of China's major agricultural products increasing exponentially or tens of times, most of which are among the highest worldwide, with grain, oilseeds, and meat all ranking first worldwide [3]. However, the great achievements of agriculture have come at the cost of environmental losses. Since the agricultural production system is still a crude production model with high inputs and outputs [4], there are many referring to the proposed IPAT model, decomposed the impact of economic growth on environmental pollution into three aspects: population size, affluence, and technical level.
There are also studies that analyze its association with ANSP in terms of natural factors [37], financial development [38], and labor quality [12] as well as those that analyze their effects on ANSP from micro perspectives such as farmers' behavior [39] and willingness [40].
The spatial spillover effect of environmental pollution refers to the first law of geography-the geographical proximity effect-which leads to the rapid expansion of pollutants with easy mobility to neighboring areas or even farther away, manifesting as a spatially contagious relationship [41,42]. ANSP is a "spillover" pollutant with a strong spatial spillover effect [42]. Existing studies on the influencing factors of ANSP are usually divided into different spatial units based on administrative boundaries, assuming that each regional unit is an independent individual, ignoring the spillover effects of ANSP between regions [43]. In contrast, the spatial econometric model assumes the existence of a spatial correlation among regions and considers the influence of surrounding areas on the study area, which can more fully reveal the effect of each influencing factor on ANSP [12].
Therefore, this study took Gongyi City, a typical representative of small towns under the background of rapid industrialization, as the research area. Based on the panel data set of 14 small towns in Gongyi City from 2002 to 2019, the non-point source inventory method that is more suitable for this study area was selected as the calculation method of ANSP load in Gongyi City. Furthermore, based on the spatial econometric model and the EKC hypothesis, the impact of various influencing factors on ANSP was empirically explored by selecting factors such as the degree of industrialization. At the same time, the driving factors behind the spatial stratification heterogeneity of ANSP, as well as the interaction between factors, were analyzed by using factor detection and interaction detection in geographical detectors. The research results can provide theoretical and practical guidance for the coordinated development of agricultural sustainable development planning in small town areas, thus enriching the existing research contents of regional development and economic geography.
The main contributions of this study are as follows: (1) The EKC model was extended to the spatial Durbin model to empirically explore the nonlinear relationship between industrialization and ANSP. Most existing studies directly assume that there is a linear relationship between industrialization and environmental pollution [44]. Given that industrialization is at a new stage of development from "high speed" to "high quality" [45], the impact of industrialization on ANSP may have a nonlinear relationship. (2) Most of the existing studies are based on traditional econometric models to analyze the factors affecting ANSP, ignoring the impact of possible spatial spillover effects and spatial heterogeneity of ANSP on the model estimation results [46,47]. In this study, the spatial spillover effect and spatial heterogeneity of ANSP were tested and estimated using spatial econometric models and geographical detectors.

Overview of the Study Area
Gongyi City is bound by latitude 34 • 31 ~34 • 52 N and longitude 112 • 49 ~113 • 17 E. It has a total area of 1041 km 2 and is subordinate to Zhengzhou City of Henan Province. The topography is high in the southeast and low in the northwest, with the Song Mountains to the south and the Yellow River to the north. Gongyi City is in a temperate monsoon climate with high temperatures and rain in the summer and a cold and dry winter, which is suitable for growing crops.
Gongyi City is an important origin of Chinese township enterprises and a typical representative of a rapidly industrializing inland region [48]. From 1978 to 2019, the average annual growth rate of the secondary industry in Gongyi was 16.46% and the average annual share of the secondary industry was 50.37%, which was higher than the national average in the same period. The average annual growth rate in the secondary industry was also higher than that in Suzhou, a typical rapidly industrializing region in the east, in the same period. Therefore, Gongyi City can be considered a typical fast industrialization area [49]. In summary, in the context of rapid industrialization, this study of the evolution and influencing factors of ANSP in small towns, which takes Gongyi city as the study area, is of guiding significance for the sustainable development of agriculture in a large number of small town areas that are in or about to be in a rapidly industrializing area.
Gongyi City has 5 streets and 15 townships. The streets are Xinhua Road Street, Dufu Road Street, Zijing Mountain Street, Yong'an Road Street, and Xiaoyi Street. The townships are Jizuankou Town, Luzhuang Town, Huiguo Town, Zhanjie Town, Zhitian Town, Zhulin Town, Shibucun Town, Mihe Town, Xinzhong Town, Hailuo Town, Kangdian Town, Dayugou Town, Xicun Town, Xiaoguan Town, and Beishankou Town. Owing to the absence of multi-year data in Zhulin Town during the study period, as well as the lack of statistics related to agricultural production and rural life in five of the streets, the remaining 14 small towns were selected as the study area in this research ( Figure 1). average in the same period. The average annual growth rate in the secondary industry was also higher than that in Suzhou, a typical rapidly industrializing region in the east, in the same period. Therefore, Gongyi City can be considered a typical fast industrialization area [49]. In summary, in the context of rapid industrialization, this study of the evolution and influencing factors of ANSP in small towns, which takes Gongyi city as the study area, is of guiding significance for the sustainable development of agriculture in a large number of small town areas that are in or about to be in a rapidly industrializing area.
Gongyi City has 5 streets and 15 townships. The streets are Xinhua Road Street, Dufu Road Street, Zijing Mountain Street, Yong'an Road Street, and Xiaoyi Street. The townships are Jizuankou Town, Luzhuang Town, Huiguo Town, Zhanjie Town, Zhitian Town, Zhulin Town, Shibucun Town, Mihe Town, Xinzhong Town, Hailuo Town, Kangdian Town, Dayugou Town, Xicun Town, Xiaoguan Town, and Beishankou Town. Owing to the absence of multi-year data in Zhulin Town during the study period, as well as the lack of statistics related to agricultural production and rural life in five of the streets, the remaining 14 small towns were selected as the study area in this research ( Figure 1).

Accounting for ANSP
Considering both the reality of agricultural development in Gongyi City and the statistical yearbook data, this study categorized the four main sources of ANSP as agricultural fertilizer, agricultural solid waste, livestock and poultry breeding, and rural living ( Table 1). All pollutants are calculated in terms of total nitrogen (TN), total phosphorous (TP), and chemical oxygen demand (COD) for measurement and comparison on a unified scale. The non-point source inventory method was used to measure ANSP in Gongyi City [26] and the relevant coefficients were calculated by referring to the existing literature [44,50,51]. The formula is as follows: (1)

Accounting for ANSP
Considering both the reality of agricultural development in Gongyi City and the statistical yearbook data, this study categorized the four main sources of ANSP as agricultural fertilizer, agricultural solid waste, livestock and poultry breeding, and rural living ( Table 1). All pollutants are calculated in terms of total nitrogen (TN), total phosphorous (TP), and chemical oxygen demand (COD) for measurement and comparison on a unified scale. The non-point source inventory method was used to measure ANSP in Gongyi City [26] and the relevant coefficients were calculated by referring to the existing literature [44,50,51]. The formula is as follows: In Formula (1), E i is the pollutant emission of ANSP i, U ij is the indicator statistic of the j-th item in the inventory on pollutant i, P ij is the pollution production coefficient of the j-th item in the inventory on pollutant i, and C ij is the loss coefficient of the j-th item in the inventory on pollutant i. Determined by the pollution source and spatial characteristics, the loss coefficient indicates the combined effects of the regional ecological environment, precipitation, hydrology, and various management measures on ANSP.
Since the discharge standards and hazard levels of each pollutant vary, it was not possible to compare the discharge of various pollutants on the same scale; therefore, the concept of equivalent discharge was introduced to compare the discharge of different pollutants. The equivalent discharge of each pollutant was calculated uniformly according to the Class III standard in the Environmental Quality Standard for Surface Water (GB3838-2002) [51]. The formula is as follows: where Q i is the equivalent pollution load of pollutant i; E i is the pollutant discharge of pollutant i; and S i is the evaluation standard. COD was 20 mg/L, TN was 1 mg/L, and TP was 0.2 mg/L.

Spatial Autocorrelation
The premise of the spatial econometric model used was the existence of significant spatial autocorrelation of the explanatory variables; therefore, this study first used the spatial autocorrelation model to analyze the equivalent emissions of ANSP in Gongyi City. The formula is as follows: ( where I is Moran's I for the equivalent emissions of ANSP; n denotes the number of township units; y i and y j are the equivalent emissions of ANSP in neighboring areas i and j, respectively; y is the average of the observed value y i ; and w ij is the spatial weight matrix. The value range of I was −1 to 1. When I was close to 1, it indicated that the ANSP in Gongyi City was positively correlated in space; when I was close to −1, it indicated that the spatial correlation of ANSP was negative. When I approached 0, there was no spatial autocorrelation. The larger the absolute value of I, the stronger the spatial correlation of ANSP.

Spatial Econometric Model
(1) Spatial measurement model construction and selection To study the spatial effects of industrialization on ANSP, three spatial panel models-the spatial lagged model (SLM), the spatial error model (SEM), and the spatial Durbin model (SDM)-were constructed. These models are the most used because they focus on the spatiotemporal effects of the dependent variable and therefore have more accurate estimation results [52]. In this study, we selected the appropriate spatial panel model with corresponding fixed effects based on the estimation and testing framework for spatial panel econometric models proposed by Elhorst [53]. Among them, the expression of the SDM is: where y it and x it denote the dependent and independent variables of the i-th township in year t, respectively; β is the spatial lag coefficient of the dependent variable; γ is the estimated coefficient of the independent variable; α is the spatial spillover coefficient of the independent variable; W denotes the spatial matrix; u i and v t denote the spatial and temporal effects, respectively; and ε it denotes the disturbance term obeying independent distribution. When α = 0 and β = 0, Equation (4) was simplified to SLM and the model did not include the interaction effects of the explanatory variables. When α + βγ = 0, Equation (4) was simplified to SEM and the SEM considers that the spatial effects exist in the perturbation error terms; the implication of this is the error shock on these factors from the changes in the explanatory variables in the neighboring regions, which will have spatial effects on the regional explanatory variables.
(2) Weighting matrix The model in this study was a spatial model that must portray the degree of interregional dependence. The degree of interdependence between regions is usually characterized by a spatial weight matrix. To increase the robustness, this study analyzed the 0-1 adjacency spatial weight matrix as the basis and the robustness of the geographic distance spatial weight matrix was checked, which uses the inverse function of the distance between two places as the elements in the matrix. The formula is as follows: Whether the spatial locations are adjacent is the basis for judging the 0-1 adjacency space weight matrix. If two regions are adjacent, W 1 = 1; if they are not adjacent, W 1 = 0: where d is the distance between the locations of the two regional geographic centers. (

3) Selection of influencing factors
Referring to relevant studies on the analysis of factors influencing ANSP, it can be found that scholars' selection of factors influencing ANSP mostly revolves around industrialization, affluence, population size, planting scale, and planting structure [52,[54][55][56][57]. Combined with the development background of rapid industrialization in Gongyi city, this study used industrialization as the core explanatory variable. With the rapid development of industrialization, relatively high wage levels attract young and strong rural laborers from the agricultural sector to the industrial sector, and according to the theory of induced technological change, when labor is insufficient in the agricultural production process, farmers tend to choose chemical technologies that can replace labor and are relatively inexpensive, exacerbating agricultural surface source pollution [58]. Furthermore, the non-agriculturalization of arable land is an inevitable phenomenon in the process of industrialization and construction. In this process, the area of arable land decreases, the quality decreases, and the degree of intensification increases, which in turn has a certain impact on ANSP [59]. At the same time, the development of industrialization is conducive to promoting technological progress in agricultural production, and the development and promotion of some environmentally friendly technologies have reduced the generation of ANSP to a certain extent. Thus, it can be seen that the impact of industrialization on ANSP is uncertain. Therefore, based on the spatial econometric model, this study introduced the squared term of industrialization to empirically analyze the EKC of industrialization on ANSP. Considering that all towns in Gongyi City are dominated by industry and agriculture, this study expressed the degree of industrialization by dividing the total industrial output value by the total industrial and agricultural output value.
Affluence level, population size, planting size, and planting structure were selected as control variables ( Table 2). The level of per capita net income of farmers was used to characterize the affluence of small town areas. The level of affluence mainly refers to the level of economic development of an area. The level of economic development determines people's productivity and lifestyle, such as the method of land use, agricultural production and operation, and residents' awareness of environmental protection. When the level of economic development is low, driven by economic interests, people will pursue the growth of output and income and adopt rough agricultural production and management methods, resulting in an increase in ANSP emissions. At higher levels of economic development, people's higher awareness of environmental protection induces agricultural producers to improve production methods and reduce ANSP emissions. Therefore, the economic development levels at different stages will have different effects on ANSP emissions. The population size factor was expressed in terms of the population density. The impact of the population size on ANSP is mainly analyzed from two perspectives: agricultural production and operation and rural life. The rural population is the main participant in agricultural economic activities, and the larger the size of the rural population in agricultural production and operation, the greater the demand for agricultural raw materials and the greater the input of agricultural chemicals, resulting in more serious ANSP. However, some scholars argue that population growth leads to an increased demand for resources, forcing people to improve production techniques, find substitutes, and increase efficiency, making the impact of population growth on the environment neutral or even positive [60]. Rural domestic pollution emissions are one of the components of ANSP, involving the discharge of domestic sewage, domestic waste, and human waste. The larger the population living in the area, the greater the likelihood of causing ANSP. Therefore, there is uncertainty in whether the impact of population size on ANSP is positive or negative.
The sown area per capita is used to characterize the scale of cultivation. As a proxy for the supply status of demand for agricultural resources in economic development, this indicator can reflect the scarcity of agricultural resources. On the one hand, the larger the sown area per capita, the richer the agricultural resources in the region, and less pressure is placed on intensive land operations. This creates prerequisites for land transfer and is conducive to promoting large-scale land operations. According to the economies of scale theory, a moderate-scale operation can improve the utilization rate of agricultural production factors and effectively reduce ANSP emissions [61]. On the other hand, a larger sown area per capita also indicates that more agricultural production factor inputs are required, resulting in an increase in ANSP emissions. Therefore, there is uncertainty in whether the planting scale has a positive or negative impact on ANSP.
The ratio of the sown area of cash crops to the sown area of food crops was used to characterize the crop planting structure. The crops planted by farmers mainly include two types of crops: food crops and cash crops. There are large differences in the types and quantities of agricultural production materials, such as the fertilizer demanded by different types of crops; for example, vegetables and wheat in food crops require large amounts of water and fertilizer in agricultural cultivation, while fruit trees and cotton in cash crops have relatively little demand for chemical fertilizer [30]. Therefore, the crop cultivation structure is an important factor affecting ANSP and changes in the structure will lead to large changes in ANSP emissions.

Geographical Detectors
Although spatial econometric models consider spatial factors in ANSP, they lack an explanation for the spatial stratified heterogeneity in ANSP. Spatial stratified heterogeneity refers to the fact that an attribute varies between types or regions [62]. Geographical detectors can detect the spatial stratified heterogeneity of the explanatory variables to some extent [63], overcoming the limitations of the mechanisms of interaction among variables in traditional analysis methods and more clearly testing the magnitude of the effect of univariate or bivariate factors on the influence of the dependent variable [64]. Therefore, in this study, the factor and interaction detectors in the geographical detectors were adopted to quantitatively analyze the driving factors of ANSP spatial differentiation in Gongyi City and the interactions among the factors. The calculation equation is as follows: where q denotes the explanation of the driving factor, which ranged from [0, 1], where the closer the variable was to 1, the greater its explanatory power; h is the stratification of the explanatory variables Y or X; N h and N are the number of cells in stratum h and the whole region; σ h 2 and σ 2 represent the variance of the values in stratum h and the whole region Y, respectively. The interaction detector determines the characteristics of the interaction between bivariate factors by comparing the q value of a single factor and the q value of a two-factor interaction. The method is advantageous because the assumptions of the interaction are not limited to those of traditional statistical methods. The interaction of the driver is identified by the q(x i ∩x j ) value of the detection result to determine whether the joint action between the drivers increases or weakens the explanatory power of the analyzed variables. The interaction detection results can be classified into five categories: nonlinearly enhanced, independent, two-factor enhanced, single-factor nonlinearly diminished, and nonlinearly diminished.

Data Sources and Processing
The statistical data in this research were mainly from the Statistical Yearbook of Gongyi In addition, in the spatial econometric model, logarithmic treatment was applied to ANSP equivalent emissions, farmers' net income per capita, population density, sown area per capita, and the crop cultivation structure to reduce the absolute values of variables and eliminate the heteroskedasticity among variables as much as possible. Stata 16.0 software was used to establish a mixed regression model and calculate the variance inflation factor (VIF). The mean value of VIF among the five explanatory variables was 1.44 and the maximum value was 2.02; all variables were less than 10, indicating that there was no significant multicollinearity.  (Figure 3). It can be seen that ANSP in Gongyi City has significant spatial heterogeneity, with a spatial distribution pattern that is high in the west and low in the east. The areas with higher

Spatial Evolution of ANSP
The natural breakpoint method was applied in ArcGIS 10.2 to classify the ANSP equivalent emissions of each township into four categories in 2002 and 2019 (Figure 3). It can be seen that ANSP in Gongyi City has significant spatial heterogeneity, with a spatial distribution pattern that is high in the west and low in the east. The areas with higher ANSP equivalent emissions are concentrated in Huiguo and Luzhuang Towns in the west, which have larger rural populations and a better agricultural resources endowment, a larger agricultural production scale, and higher ANSP equivalent emissions. The areas with lower emissions of ANSP are concentrated in the eastern areas of Xinzhong Town and Dayugou Town, which are typical mountainous towns with sparse rural populations and relatively regressive agricultural development.

Spatial Evolution of ANSP
The natural breakpoint method was applied in ArcGIS 10.2 to classify the ANSP equivalent emissions of each township into four categories in 2002 and 2019 (Figure 3). It can be seen that ANSP in Gongyi City has significant spatial heterogeneity, with a spatial distribution pattern that is high in the west and low in the east. The areas with higher ANSP equivalent emissions are concentrated in Huiguo and Luzhuang Towns in the west, which have larger rural populations and a better agricultural resources endowment, a larger agricultural production scale, and higher ANSP equivalent emissions. The areas with lower emissions of ANSP are concentrated in the eastern areas of Xinzhong Town and Dayugou Town, which are typical mountainous towns with sparse rural populations and relatively regressive agricultural development.   (Table 3). The results show that Moran's I index of ANSP equivalent emissions is significantly positive during the study period and the change trend also indicates Moran's I index is increasing, demonstrating that ANSP in each township of Gongyi City is not randomly spatially distributed but has positive global spatial autocorrelation characteristics. Over time, the agglomeration gradually increased.  (Table 3). The results show that Moran's I index of ANSP equivalent emissions is significantly positive during the study period and the change trend also indicates Moran's I index is increasing, demonstrating that ANSP in each township of Gongyi City is not randomly spatially distributed but has positive global spatial autocorrelation characteristics. Over time, the agglomeration gradually increased.

Spatial Econometric Model Empirical Analysis
Given that spatial econometric models have different forms, this study used Stata16.0 to perform LM, Wald, LR, and Hausman [30,65] tests. The results are shown in Table 4. Note: *** and ** denote significance at the 1% and 5% levels, respectively.
As can be seen from Table 4, the LM test and R-LM test results reject the original hypothesis at the 5% significance level, indicating that the model has both spatial lag and spatial error terms. The LR test and Wald test results reject the original hypothesis at the 1% significance level, indicating that the SDM cannot be reduced to a spatial lag model and SEM, i.e., the SDM is optimal for the simulation of ANSP in Gongyi City. The SDM is optimal for the simulation of influencing factors. Meanwhile, the Hausman test results show that the original hypothesis of random effect is rejected at the 1% significance level. Therefore, the fixed effects model of the SDM was selected for the analysis and the regression results are shown in Table 5. The fixed effects model contains three types of effects: time fixed, spatially fixed, and spatiotemporally fixed. In order to compare which fixed effects model has the best fitting effect, based on the magnitude of the log-likelihood and dispersion (sigma 2 ) in Table 6, the spatiotemporally fixed effects' fitting results were optimal. Therefore, the spatiotemporally fixed SDM model was finally selected as the analytical model in this study.  Note: ***, **, and * indicate significance at the 1%, 5%, and 10% levels, respectively; numbers in parentheses are Asymptot t-stat values for the coefficients of the corresponding variables. Note: ***, **, and * indicate significance at the 1%, 5%, and 10% levels, respectively; numbers in parentheses are Asymptot t-stat values for the coefficients of the corresponding variables.

Decomposition of Spatial Spillover Effects
The estimated coefficients of the model could not explain the spatial spillover effects of ANSP. Therefore, a partial differential solution was calculated based on Lesage's study to assess the magnitude of the influence and spillover effects of the explanatory variables through direct and indirect effects [66] and the results are shown in Table 6. The direct effect is expressed as the magnitude of the effect of the explanatory variables in the region on the explanatory variables in the region; the indirect effect is expressed as the effect of the explanatory variables in neighboring regions on the explanatory variables in the region; and the total effect is expressed as the sum of the direct and indirect effects.
There is an inverted "U-shaped" EKC relationship between industrialization and ANSP. As can be seen from Table 7, the total effect of the primary term of industrialization is positive, the total effect of the secondary term is negative, and both pass the 1% significance test. Thus, as industrialization progressed, ANSP first increased and then decreased. With rapid industrialization, the relatively high wage level attracts young and strong rural laborers to transfer from the agricultural sector to the industrial sector, which leads to the emergence of part-time farming behavior. According to the theory of induced technological change, when there is a shortage of labor in the agricultural production process, farmers tend to choose chemical technologies that can replace labor and are relatively cheap, exacerbating ANSP [67]. Furthermore, the non-agriculturalization of arable land is an inevitable phenomenon in the process of industrialization construction, in which the area of arable land decreases, the quality decreases, and the degree of intensification increases, which in turn has a certain impact on ANSP [68]. At the same time, the technological progress brought on by industrialization has accelerated the evolution of traditional agriculture to modern agriculture (i.e., petrochemical agriculture), and agricultural development has entered a "chemical trap" with an increasing number of chemical products being applied to agricultural production [69]. All these effects have exacerbated ANSP. When industrialization crosses the inflection point, the positive externalities of industrialization play a dominant role and the growth of national fiscal revenue provides sufficient financial support for ANSP control. At the same time, the development of efficient and low-consumption environmental protection technologies and advanced management technologies also makes the prevention and control of ANSP more complete. The indirect effect was not significant, indicating that the increase in the industrialization level in neighboring townships had no significant effect on ANSP in this township. From the decomposition of the effects of each control variable, affluence played a positive role in promoting ANSP. The direct effect of the net per capita income of farmers on ANSP was positive and passed the test at the 5% significance level, indicating that increasing the rural economy level played a driving role in ANSP in this township. This may be because as affluence increases, increasing factors of production, such as chemical fertilizers, pesticides, and agricultural films are put into agricultural production, which in turn increases ANSP emissions. The indirect effect was positive; however, it did not pass the significance test. This indicates that the increase in the economic level of neighboring townships had a positive effect on the ANSP in the township.
Population size plays a positive role in promoting ANSP. The direct effect of the rural population density on ANSP was positive and passed the test at the 1% significance level, indicating that the increasing rural population density played a facilitating role in ANSP in this township. This is because the rural population is the main participant in agricultural economic activities and rural life; therefore, a higher population density leads to increased pressure on the agricultural environment in this township; a higher population density means that human production, living, and consumption activities in the area will all increase, and the accompanying emissions of ANSP will also increase [12]. The indirect effect of the rural population density on ANSP was negative and insignificant, indicating that, although there is some difference in the rural population density within each township, this difference did not trigger population movement in neighboring townships, which in turn has an impact on ANSP emissions in this township.
The effect of the planting scale on ANSP was not significant. The direct effect of the per capita sowing area on ANSP was positive and passed the test at the 1% significance level, indicating that increasing the per capita sowing area played a driving role in ANSP in this township. This is likely because most townships in Gongyi City are located in mountainous and hilly areas, with poor agricultural production resources and generally small land operation scales; to meet the social demand for agricultural products, they will increase the land intensification degree, increasing the input of pesticides, fertilizers, and other production materials, which to a large extent aggravates ANSP. The indirect effect was negative and passed the test at the 10% significance level, indicating that the increase in the per capita sown area in neighboring townships had a suppressive effect on the ANSP in this township. This is because the higher planting scale of neighboring townships will produce a clustering effect, attracting the inflow of agricultural production resources from this township to the neighboring area, which in turn reduces the ANSP emissions in this township. The total effect was not significant, likely because the direct and indirect effects cancelled each other out.
There is a significant positive contribution of the planting structure to the ANSP. The direct effect of the cropping structure is positive and passes the test at the 1% significance level, indicating that the increase in the proportion of the sown area of cash crops will promote the growth of ANSP emissions in this township. As the level of economic development increases, farmers pursue higher profits in agricultural production and operation, resulting in a decrease in the proportion of the sown area of grain in the plantation industry and an increase in the sown area of cash crops. While cash crops bring farmers rich economic returns, they also produce more serious ANSP. In terms of the amount of fertilizer applied, it can be seen from the 2010 Guidance on Scientific Fertilization of Major Crops that the amount of fertilizer applied to cash crops is significantly higher than that of food crops. The average utilization rate of nitrogen fertilizer on farmland planted with vegetables, fruit trees, and flowers only reached approximately 10%, which is much lower than the utilization rate of 35% for food crops [36]. The indirect effect was positive but failed the significance test, indicating that the planting structure of neighboring townships was exemplary; however, the spillover effect of this demonstration was not significant and had a very limited impact on ANSP emissions in this township.

Robustness Test
A reasonable measure of the spatial relationship between regions is a prerequisite for spatial econometric analysis and the construction of a spatial weight matrix is a key step. A subjective and single spatial weight matrix may not truly reflect the complex relationships among different regions [70]. Therefore, the accuracy and robustness of the spatial econometric modeling results are highly dependent on the spatial weight matrix. To ensure the reliability of the above empirical results, the model was tested by replacing the weight matrix. The direct, indirect, and total effects under the SDM model of spatiotemporally fixed effects based on the geographic distance weight matrix are shown in Table 7. Compared with the 0-1 adjacency weight matrix, the regression coefficients and significance changes of each effect based on the geographic distance weight matrix were smaller, indicating that the results of this study were robust.
The model may have endogeneity and omitted variables that lead to biased estimation results since ANSP is a dynamic and continuous process; consequently, the current increase in ANSP will be influenced by both current-and prior-period factors. Therefore, a dynamic spatial Durbin model (DSDM) with a first-order lag containing the explanatory variable ANSP was simultaneously constructed. The model results are shown in Table 8, and the results indicate that the direction and significance of the regression coefficients are basically consistent with those of previous studies after considering endogeneity. Furthermore, from the total effect of short-term and long-term effects, each influencing factor exhibits a stronger temporal effect of short-term effects than long-term effects. Note: ***, **, and * indicate significance at the 1%, 5%, and 10% levels, respectively; numbers in parentheses are Asymptot t-stat values for the coefficients of the corresponding variables.

Analysis of ANSP Drivers Based on Geodetectors
To further explore the driving factors of spatially stratified heterogeneity of ANSP in Gongyi City, this study employed geodetector software for driving factor detection. In this study, the five detection factors of industrialization (X1), per capita net income of farmers (X2), population density (X3), per capita sown area (X4), and grain crop planting structure (X5) were selected from five dimensions of industrialization, affluence, population size, land size, and planting structure to analyze the spatial heterogeneity of ANSP in Gongyi City. RStudio software and the R language "GD" package were used to detect the geographic probe factors (Table 9) and interactions ( Figure 4).  Table 9 shows that the spatial heterogeneity of ANSP is mainly the result of the combined effect of multiple factors and the driving factors are in the order: X3 > X4 > X1 > X5 > X2. This indicates that the population size, land size, and industrialization degree were the most influential factors affecting ANSP in Gongyi City. In the process of ANSP management in small town areas in the future, it is more important to strengthen the population and land planning to reduce the pressure of the environmental carrying capacity. The design and promotion of agricultural green production technology should be accelerated in the meantime, crossing the industrialization EKC inflection point as early as possible, releasing the positive externalities of industrialization, reducing ANSP, and achieving green and sustainable agricultural development. farmers (X2), population density (X3), per capita sown area (X4), and grain crop planting structure (X5) were selected from five dimensions of industrialization, affluence, population size, land size, and planting structure to analyze the spatial heterogeneity of ANSP in Gongyi City. RStudio software and the R language "GD" package were used to detect the geographic probe factors (Table 9) and interactions ( Figure 4).   Table 9 shows that the spatial heterogeneity of ANSP is mainly the result of the combined effect of multiple factors and the driving factors are in the order: X3 > X4 > X1 > X5 > X2. This indicates that the population size, land size, and industrialization degree were the most influential factors affecting ANSP in Gongyi City. In the process of ANSP management in small town areas in the future, it is more important to strengthen the population and land planning to reduce the pressure of the environmental carrying capacity. The design and promotion of agricultural green production technology should be accelerated in the meantime, crossing the industrialization EKC inflection point as early as possible, releasing the positive externalities of industrialization, reducing ANSP, and achieving green and sustainable agricultural development.
ANSP is the result of the joint action of multiple factors; therefore, the interaction factor exploration module was used for in-depth analysis ( Figure 4). The interaction of the factors was bifactor-enhanced or nonlinearly enhanced, i.e., the influence of factors on the interaction was greater than the sum of the independent forces of two factors, reflecting the complex characteristics between the compound influencing factors and ANSP. The explanatory power of the interaction between the population density and the sown area per capita was significantly stronger than that of the interaction among other factors, which was the controlling factor for the significant spatial variation of ANSP.

Discussion and Recommendations
This study examined the evolution of ANSP and the factors influencing it in small town areas in the context of rapid industrialization from a geographic perspective. Although some scholars have already conducted related studies, this study still had some innovative features. First, in comparison to Mahmood et al. [59] who suggested that industrialization would lead to increased environmental pollution, this study showed an ANSP is the result of the joint action of multiple factors; therefore, the interaction factor exploration module was used for in-depth analysis ( Figure 4). The interaction of the factors was bifactor-enhanced or nonlinearly enhanced, i.e., the influence of factors on the interaction was greater than the sum of the independent forces of two factors, reflecting the complex characteristics between the compound influencing factors and ANSP. The explanatory power of the interaction between the population density and the sown area per capita was significantly stronger than that of the interaction among other factors, which was the controlling factor for the significant spatial variation of ANSP.

Discussion and Recommendations
This study examined the evolution of ANSP and the factors influencing it in small town areas in the context of rapid industrialization from a geographic perspective. Although some scholars have already conducted related studies, this study still had some innovative features. First, in comparison to Mahmood et al. [59] who suggested that industrialization would lead to increased environmental pollution, this study showed an inverted "U" shaped EKC relationship between industrialization and ANSP in small towns in a typical area of ANSP in China. This study is helpful to clarify the relationship between industrialization and ANSP in small towns with contradictory agricultural development and to extend the boundary of environmental Kuznets theory to a certain extent. Furthermore, compared with the traditional regression model used by Zhang et al. [26] to analyze the influencing factors of ANSP, this study analyzed the influencing factors of ANSP in terms of both spatial spillover effects and spatial stratified heterogeneity, which was conducive to improving the credibility of the results. Meanwhile, it was found that affluence, population size, and cropping structure positively contributed to ANSP; these findings are consistent with those of existing studies [35,36,38]. The spatial spillover effect of industrialization on ANSP was not significant, which is likely due to the small scale of the study area in this study and the low exchange of resource factors between regions. The results of this study can provide a scientific basis for the prevention and the control of ANSP in similar rapidly industrializing small town areas and play a supporting role in sustainable agricultural development.
Based on the results of the empirical analysis, combined with the basic ideas of agroecological economics and sustainable agricultural development, the following suggestions for the prevention and control policies of ANSP are proposed.
(1) Control at the source to prevent the risk of ANSP is necessary. Temporally, ANSP in Gongyi City from 2002 to 2019 showed a zigzagging downward trend, indicating that the agricultural development in Gongyi City was positive; however, there are still some inefficiencies and environmental losses and there is substantial room for agricultural pollution reduction. From the perspective of pollution sources, ANSP emissions mainly come from farmland solid waste and rural life, livestock and poultry breeding, and agricultural fertilizer; pollutant emissions were mainly in the order of TN > TP > COD. Thus, real-time prediction and assessment of production and discharge sources of pollution is crucial.
This will also improve the effective use of agricultural solid waste storage and processing equipment, promote the updating and upgrading of relevant technologies, and enhance the standardization of rural living environments via targeted pollution risk mitigation.
(2) Strengthening cooperation for mutual benefits and creating "eco-efficient" situations is key. According to the spatial econometric model results, ANSP shows a significant spatial correlation and clustering distribution pattern. Therefore, policymakers should fully consider the spillover effect of ANSP when formulating policies, eliminate local protectionism, conduct comprehensive planning from the perspective of regional integration, establish a complete agricultural policy linkage mechanism, and strengthen inter-regional agricultural production resource sharing, exchange, and cooperation.
(3) Increasing the research, development, and promotion of green agricultural production technologies will continue to be important. From analyzing the impact of industrialization on ANSP, we can see that the improvement of agricultural green production technology can effectively promote an early crossing of the inflection point of the EKC. Technological innovation is the only path for China's ongoing agricultural development. In terms of China's agricultural growth, although the contribution of science and technology to agricultural development has reached 42%, it is still far behind that of other developed countries (which tend to be closer to 70%). Therefore, we should increase investment in agricultural science and technology innovation and vigorously develop and promote new resource-saving and environmentally friendly agricultural technologies. For example, by promoting soil testing and fertilization techniques, encouraging farmers to apply more organic fertilizers, developing new pesticides such as bionic pesticides, and promoting biodegradable agricultural films. Furthermore, as China is limited by farmers' education levels, learning abilities, and production input levels, there is a greater need to develop simple, easy-to-operate, and low-cost agricultural green production technologies.
In general, this study has some limitations. First, the non-point source inventory method was chosen for the measurement of ANSP, which means the research results may deviate from the actual situation. Second, the spatial Durbin model was used for the influence factor analysis, which lacked consideration of local spatial regression models. Therefore, in future studies on ANSP, field ANSP measurements should be increased as much as possible and local spatial regression models should be tried in the research method, so as to improve the accuracy of the research conclusions.

Conclusions
In this study, panel data from 14 small towns in Gongyi City from 2002-2019 were used to study the evolution of ANSP and influencing factors in small town areas in the context of rapid industrialization using the non-point source inventory method, spatial econometric model, and geographical detectors. The main conclusions are as follows. (2) Spatially, the equal standard emissions of ANSP showed a distribution pattern of being high in the west and low in the east, with higher areas concentrated in the western Huiguo and Luzhuang Towns and lower areas concentrated in the eastern Xinzhong and Dayugou Towns. Meanwhile, there was a positive global spatial autocorrelation feature. (3) There was an inverted "U-shaped" EKC curve between industrialization and ANSP in Gongyi City, i.e., as industrialization progressed, the equivalent emissions of ANSP first increased and then decreased; affluence, population size, and planting structure contribute positively to the reduction in ANSP. Furthermore, all influencing factors showed stronger short-term temporal effects. (4) Population size, land size, and industrialization were the high-impact factors on ANSP in Gongyi City and the interactions between these factors showed a two-way or nonlinear enhancement. Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The associated dataset of the study is available upon request to the corresponding author.