Spatiotemporal evolution analysis of agricultural non-point source pollution risks in 1 Chongqing, China, based on the ITO3dE model and GIS 2

7 Abstract 8 To determine the risk state distribution, risk level, and risk evolution situation of agricultural non-point source 9 pollution (AGNPS), we built an ‘Input-Translate-Output’ three-dimensional evaluation (ITO3dE) model that 10 involved 12 factors under the support of GIS and analyzed the spatiotemporal evolution characteristics of AGNPS 11 risks from 2005 to 2015 in Chongqing by using GIS space matrix, kernel density analysis, and Getis-Ord Gi* analysis. 12 Land use changes during the 10 years had a certain influence on the AGNPS risk. The risk values in 2005, 2010, and 13 2015 were in the ranges of 0.40-2.28, 0.41-2.57, and 0.41-2.28, respectively, with the main distribution regions being 14 the western regions of Chongqing (Dazu, Jiangjin, etc.) and other counties such as Dianjiang, Liangping, Kaizhou, 15 Wanzhou, and Zhongxian. The spatiotemporal transition matrix could well exhibit the risk transition situation, and 16 the risks generally showed no changes over time. The proportions of ‘no-risk no-change’, ‘low-risk no-change’, and 17 ‘medium-risk no-change’ were 10.86%, 33.42%, and 17.25%, respectively, accounting for 61.53% of the coverage 18 area of Chongqing. The proportions of risk increase, risk decline, and risk fluctuation were 13.45%, 17.66%, 7.36%, respectively. Kernel density analysis was suitable to explore high-risk gathering areas. The peak values of 20 kernel density in the three periods were around 1,110, suggesting that the maximum gathering degree of medium- 21 risk pattern spots basically showed no changes, but the spatial positions of high-risk gathering areas somehow 22 changed. Getis-Ord Gi* analysis was suitable to explore the relationships between hot and cold spots. Counties with 23 high pollution risks were Yongchuan, Shapingba, Dianjiang, Liangping, northwestern Fengdu, and Zhongxian, while 24 counties with low risks were Chengkou, Wuxi, Wushan, Pengshui, and Rongchang. High-value hot spot zones 25 gradually dominated in the northeast of Chongqing, while low-value cold spot zones gradually dominated in the 26 Midwest. Our results provide a scientific base for the development of strategies to prevent and control AGNPS in 27 Chongqing.


30
Agricultural non-point source pollution (AGNPS) refers to water pollution caused by nitrogen 31 and phosphorus, pesticides, and other contaminants through farmland water runoff and percolation 32 (Smith and Siciliano, 2015). In China, in the prevention and control of water pollution, point sources 33 have mainly been taken into consideration, while non-point sources are largely being ignored (Liu 34 topography, a high proportion of rural areas in the urban and rural dual structure, a hot and rainy 48 season, and a high and concentrated precipitation, which result in a high potential threat, a wide 49 coverage area, and a high driving energy of AGNPS in Chongqing (Zhou et al., 2019;Li et al., 50 2018). The application amount of chemical fertilizer, pesticide and agricultural film in Chongqing 51 is increasing constantly. In 2014, the application level of chemical fertilizer, pesticide and 52 agricultural film reached 411 kg ha, 9.5 kg ha and 81.9 kg ha respectively, which was far higher 53 than the international standard upper limit. In addition, due to the influence of topography and 54 climate on soil erosion, the total amount of soil erosion is 146 million tons per year. In addition, 55 there are other problems such as high multiple cropping index of agricultural land in Chongqing 56 (Hu, 2017). Chongqing is located in the center of the Three Gorges Reservoir area and represents 57 the connection point of the "Belt and Road" and the Yangtze River economic belt. Hence, Chongqing 58 has an important status in the national regional development pattern and is an important ecological 59 barrier of the upstream of the Yangtze River; in this sense, strict water management strategies are 60 indispensable. To ensure the ecological security of the Yangtze River basin, the issue of AGNPS 61 needs to be resolved. In view of this, it is of great importance to understand the spatiotemporal 62 evolution situation of AGNPS risks in districts and counties of Chongqing. To sum up, there are 63 7 rate is only about 35%. Pesticide use intensity is 9.5 kg ha, with a use rate of only 30%. Livestock 130 and poultry stocks are currently increasing significantly. In addition, the hilly and mountainous areas 131 in Chongqing account for 94% of the total area, while the cultivated land area with a slope greater 132 than 25 degrees accounts for 16% of the total cultivated land area, which is 11 percentage points 133 higher than the national average level. Rainfall in Chongqing is large and concentrated, and the area 134 stable at 25%. The artificial surface area has increased from 1.42% to 3.18%, while the proportion 144 of forested land showed a steady increase (Fig. 1). Overall, the farmland area of Chongqing showed 145 a significantly decreasing trend, with a reduction from 39.35% to 33.03%. Farmland is the main 146 source of AGNPS. To sum up, the actual regional conditions of land use types, topography, climate, 147 fertilizer and pesticide application lead to the aggravation of AGNPS in Chongqing.

Research methods 151
The risk assessment of AGNPS in Chongqing was carried out by constructing the ITO3dE 152 model, and the specific content was analyzed by using GIS technology. The overall framework of 153 the study was as follows (Fig. 2).

ITO3dE model building 157
The ITO3dE model was built from three dimensions of "Input-Translate-Output" (Table 1). 158 Here, the I dimension (A1) includes fertilizer use intensity (I1) (Liu, 2014), pesticide use intensity  The computational formula of positive indices is as follows: 167 The computational formula of negative indices is as follows: 169 In the above equations, Ii denotes the calculation result of a certain index of i, Ci is the true value of 171 the i index, and Ei is the reference value. The calculation results for each index can be classified into 9 five grades of no risk, low risk, medium risk, high risk, and extremely high risk, corresponding to 173 the values of ≤ 0.7, 0.7-1.0, 1.0-3.0, 3.0-5.0, and ≥ 5.0. This grading standard refers to the relevant 174 technical planning of China's agricultural sector (Chen and Zhang, 2010), and the grading standard 175 of this study was written based on Chongqing local standard "Specification for evaluation risk of 176 AGNPS in Chongqing". 177

Delphi method 178
The weights of the evaluation indices of AGNPS risks were determined by the Delphi method 179 (Strand et al., 2017). We distributed 120 questionnaires to the leaders of related departments 180 engaged in agricultural and rural work, environmental work, and water conservancy programs, as 181 well as to experts at colleges, universities, and research institutes, and to technical personnel at the 182 grass-roots level; in total, we received 115 valid questionnaires. After the first round of weights 183 assignment, we returned the calculation results to the experts for adjusting the index weights,

Spatiotemporal transition matrix of risk index 196
Transition matrices are widely used in analyzing the spatiotemporal variations and can clearly 197 exhibit the variations of risk indices at different periods (Zhu et al., 2015). In this study, we assigned 198 the five risk grades to the values of 1, 2, 3, 4, and 5, respectively, and subsequently employed the 199 following formula to analyze the risk status:

Kernel density analysis of risk index 205
Kernel density is mainly used to analyze the spatial concentration of an event and is widely 206 certain elements is to assume that there always exists an element intensity at an arbitrarily position 211 within a particular region. The density intensity of geographical elements in a specific location can 212 then be estimated by measuring the number of elements per unit area. By determining the element 213 density at different locations and the spatial differences, the relative concentration degree of the 214 spatial distribution of elements can be depicted, and the hotspot distribution regions can be identified. 215 11 grids at risk grades of high risk and extremely high risk and explore the extreme point position of 217 AGNPS in our study area, Chongqing. 218

Getis-Ord Gi* analysis 219
The Getis-Ord Gi* analysis is widely used in crime analysis, epidemiology and economic 220 geography, and is used to identify spatial gathering of high values (hot spots) and low values (cold

Statement 231
We confirm that all experimental protocols were approved by College of Resources and 232 Environment of Southwest University, all methods were carried out in accordance with relevant 233 guidelines and regulations of questionnaires, and confirming that informed consent was obtained 234 from all subjects or, if subjects are under 18, from a parent and/or legal guardian. 235

Results of risk assessment by the ITO3dE model 237
The results in the I dimension show that, overall, the distribution was high in the west and low 238 12 in the northeast and the southeast in all three periods (Fig. 3, Ⅰ, Ⅱ, Ⅲ; Table 2 Considering there are too many single-factor graphs, we omitted these graphs, but provide the 247 following description: Among the three single factors, I1 has the highest value, and I1 and I2 both 248 risk grade of I2 was relatively lower than that of I1, but overall, the spatiotemporal variation trend 254 was consistent with that of I1, except for the increasing trend of the risk level of Qianjiang. Basically, 255 the risk grade of I3 was zero; only the risk level of Bishan was in the medium risk status, while those 256 of Hechuan and Fengdu were low. 257 Spatially, the results in the T dimension presented, overall, an opposite distribution pattern 258 when compared to the I dimension, that is, with low levels in the western regions and high levels in 259 the northeastern and southeastern regions (Fig. 3, Ⅳ,Ⅴ,Ⅵ; Table 2). The annual differences in the were 1.42-5.78, 0.84-6.12, and 0.14-6.93, respectively, while those of I7 were 0, 0-5.38, and 0-5.06, 263 respectively. Because Chongqing is a typical mountainous city with purple soil (Wang et al., 2016), 264 high-risk and extremely high-risk regions, I5 and I6, are widely distributed across the city. In addition, 265 due to the introduction of the factor I8, the water areas had a higher risk level, which is consistent 266 with the actual situation of AGNPS. 267 The results in the O dimension showed a smaller interannual variation, with a low overall risk 268 level (Fig. 3, Ⅶ, Ⅷ, Ⅸ; Table 2). The O dimension levels were mainly affected by the spatial 269 changes in the paddy field area. As mentioned above, during the 10 years, the area of paddy fields 270 in Chongqing was nearly reduced by half, which led to the decrease in the spatial distribution of I12

Spatiotemporal change results of risk by transition matrix analysis 300
By assigning no risk, low risk, and medium risk levels with 1, 2, and ,3 respectively, in GIS, 301 we can obtain the spatiotemporal transition matrix according to the formula of the transition matrix. 302 Figure 5 shows the spatiotemporal transition situation of the AGNPS risk evaluation in Chongqing. 303 Basically, high levels show no changes, and the proportions of 'no-risk no-change', 'low-risk no-304 15 change', and 'medium-risk no-change' situations were 10.86%, 33.42%, and 17.25%, respectively, 305 accounting for 61.53% of the total area of Chongqing. Among these, the 'no-risk no-change' 306 situation was mainly distributed in Rongchang, the east of Nanchuan, Shizhu, Pengshui, and 307 Qianjiang; the 'low-risk no-change' situation was widely distributed in Wulong, the southeast of 308 Fengdu, the south of Nanchuan, and the northeastern counties of Chongqing, while the 'medium-309 risk no-change' situation was mainly distributed in Shapingba, Yongchuan, Dianjiang, the north of 310 Nanchuan, and Kaizhou. 311  Figure 6 shows the kernel density analysis results of the medium-risk regions. As seen in the 322 figures, the peak values of the kernel density at these three periods were all around 1,110, suggesting 323 that the maximum gathering degree of medium-risk pattern spots basically showed no changes. The 324 spatial distribution of kernel density at these three periods showed a consistent trend, but the 325  To further explore the distribution of regions with the high-risk gathering zones (Table. 3), we 339 conducted a separate analysis on the regions with kernel density values higher than 1,000 (the kernel 340 density values of these regions were divided into 10 grades with equal intervals, and the 10th grade 341 had values from 1,000 to 1,110). 342 Table 3 Distribution of regions with high-risk gathering zones 343 <Please insert Table 3 about here> 344

Results of hot and cold spots by Getis-Ord Gi* analysis 345
Applying Getis-Ord Gi* analysis is helpful to clearly identify high-value hot spots (Hot Spot-346 99% Confidence) and low-value cold spots (Cold Spot-99% Confidence). Figure 7 shows the Getis-347 Ord Gi* analysis results; the overall variation trends of high-value hot spots and low-value cold 348 spots were consistent in all periods, with significant distribution differences. The regions located in 349 the high-value hot spot zones in all three periods were Yongchuan, Shapingba, Dianjiang, Liangping, 350 northwestern Fengdu, and Zhongxian, while those located in the low-value cold spot zones were

The I, T, and O dimensions could better comprehensively analyze the risk situation of AGNPS. 369
There were many achievements in AGNPS risk research using different methods. For instance, 370 18 Huang has analyzed the risk of NPS in Taihu Lake from multiple perspectives, and the method was 371 only applicable to small-scale studies (Huang et al., 2013). In addition, it focused on the distribution 372 of pollution sources, and on waste disposal, but did not consider the process from pollution sources 373 to water bodies. Blankenberg has analyzed the relationship between wetland and pesticide 374 concentration in streams from the perspectives of pesticide spraying and changes in pesticide 375 concentration in water bodies (Blankenberg et al., 2006). and low values. This method was usually applied in studies that require the analysis of high-value 407 and low-value changes simultaneously (Stopka et al., 2017). In our research, we need to consider 408 both high-risk and low-risk changes in AGNPS. 409

Land variations influence AGNPS risks 410
The analysis of land use changes indicates that the land area in Chongqing has changed greatly 411 during the 10-year-period from 2005 to 2015. In particular, the paddy field area was reduced by 412 nearly 50%. Our observations are in agreement with previous findings (Liu et al., 2017). found that land use changes had an impact on water quality (Mehdi et al., 2015). The water quality 424 factors in this study also reflected that the water quality in the southeast and northeast areas with 425 high forest coverage was better, indicating that land use change had a great impact on water quality, 426 especially on the cultivated land, which was closely related to the use of fertilizers and pesticides. 427 In this study, we analyzed the spatiotemporal variation in AGNPS in three dimensions of 428 "Input-Transition-Output", with the conclusion that the risk levels in the Input and Output 429 dimensions were lower, while that in the Translate dimension was higher. Our conclusion accords 430 with the characteristics of Chongqing as a mountainous city; the widely distributed purple soil, the 431 larger differences in slope gradient and slope length, and the widely distributed sloping fields are 432

Focus on high-risk agglomeration area in nuclear density results 455
The high values of kernel density changed only slightly over time, indicating that the maximum 456 scope of gathering zones showed no significant changes. The spatial position variation of kernel 457 density at different periods is closely related to the changes of the factors in the three dimensions.

Focus on high-value hot spots area in Getis-Ord Gi* analysis 476
The results of the Getis-Ord Gi* analysis clearly reflect the spatiotemporal evolution situation 477 of high-value hot spot zones and low-value cold spot zones in the three periods in Chongqing. Over . For example, if the I dimension value of an area is high, but the T dimension is low, that is, 497 the measured data value is high, this does not automatically mean a high risk. The results of risk 498 assessment are based on the reliability of the data source, factor weight and grading reference value 499 of risk factors. It is based on the historical data to judge the trend of regional pollution risk, so as to 500 prevent and control the risk, and this is the advantage of risk assessment in pollution prevention. closely related with the widely distributed purple soil, the large differences in slope gradient and 512 slope length, and the widely distributed sloping fields. The risk evaluation results obtained via 513 spatiotemporal transition matrix, kernel density analysis, and Getis-Ord Gi* analysis specifically 514 quantified the regions that require close attention in the prevention and control of AGNPS from the 515 perspectives of risk level variation, the distribution of risk highly-concentrated zones, and the shift 516 of high-value hot spot zones and low-value cold spot zones; we therefore provide data to support 517 land use planning strategies as well as measures to prevent and control AGNPS. 518 Taken together, there would be significant influences on the Input and Output dimensions of 519 AGNPS that promote the "one regulation, two reductions, and three basics" policies promulgated 520 by the Chinese government, indicating areas suitable for livestock breeding in Chongqing and 521 implementing regulatory requirements in 2018. The implementation of related policies would 522 effectively reduce the risk levels of these two dimensions. In this study, the results suggest that there 523 is a higher risk level in the Translate dimension. Among its influencing factors, field slope is 524 25 important and can be improved by land use planning or artificial handling. Hence, related 525 departments should focus on the farmland protection in the high-risk gathering regions, avoid 526 farmland reclamation as far as possible, and reduce the proportion of sloping fields. At the same 527 time, the migration of pollutants into water bodies is mainly affected by the degree of soil erodibility, 528 so there would a positive impact that strengthening the control measures of regional soil erosion on 529 reducing AGNPS. 530

531
This work was funded by Chongqing Science and Technology Commission (cstc2018jszx-532 zdyfxmX0021, cstc2018jxjl20012 and cstc2019jscx-gksbX0103). We are thankful to all experts 533 involved in the Delphi survey. Thanks are expressed to the Chongqing Agricultural and Rural 534 Committee for providing information on fertilizers, pesticides, and livestock and to the Chongqing 535 Eco-environment Bureau for providing information on water-related issues. We also thank the and 536