Exploring effective network design for monitoring soil pH and potentially toxic elements based on geochemical surveys in economically developed area

Periodically visiting soil monitoring sites, i.e., sampling and analysis, is recognized as one of the most important ways to monitor soil quality. However, reconciling the monitoring costs with monitoring precision of the soil monitoring network (SMN) is a key technical problem to be solved. A statistically sound method, which depends on the spatial variation in monitoring indicators, was adopted to determine the number of monitoring sites and the monitoring interval as well as their ability to detect a particular change under an economically feasible scenario. The spatial variation in soil monitoring indicators was inquired from the “Multi-Purpose Regional Geochemical Survey in Zhejiang Province (MRGSZ)” project. Based on the data for soil pH and concentration of potentially toxic elements, the number of monitoring sites and the monitoring intervals that might be used for soil monitoring were determined with the administrative region as the monitoring unit. The results showed that there was great spatial variation in the MRGSZ region, which resulted in discrepancies in the minimum detectable changes (MDCs), monitoring site numbers, and temporal monitoring intervals for revisiting. Our research proposes a number of monitoring sites (nr) that could reconcile the monitoring costs, practicability and monitoring precision; thus, it was recommended for the design of SMNs. Under nr, the MDC values of each monitoring indicator were acceptable for all administrative regions, and the temporal monitoring intervals were practical with variations of 6.7–14.8 years.


Introduction
Soil is not only a type of nonrenewable resource but also the material basis for human survival and development.High-intensity industrial and mining operations, as well as the unsustainable use of chemical fertilizers and pesticides in recent decades, have resulted in widespread soil pollution all over the world which include China (Chen et al., 2014; Ministry of Environmental Protection and Ministry of Land and Resources of P. R. China, 2014, Luo et al., 2020, Shukla L et al., 2022).Based on the data from the "Multi-purpose Regional Geochemical Survey (MRGS)", which covers an area of more than 1.50 m km 2 , a geochemical investigation report of arable land in China (China geological survey, 2015) showed that the proportion of soil heavy metal contents exceeding the standard was as high as 8%, of which the sample ratios of slight to mild-and moderate-to-severe pollution were 5.7 and 2.5%, respectively.Therefore, researchers carried out extensive research on the status of soil pollution, cause, risk, etc. (Wu et al., 2020, Karimian et al., 2021, Golia et al., 2022).
Due to great significance for the early detection of quality changes, carrying out soil monitoring and periodic visits to soil monitoring sites was recognized as the most efficient way to sustainably manage and protect the availability and quality of soil (Lark et al., 2006;Desaules, 2012).Many countries have established an official national framework of soil monitoring networks, although the monitoring methods and monitoring scopes are not consistent (Table 1).Teng et al. (2014) systematically reviewed soil environmental quality monitoring in China and put forward some suggestions on the establishment of multilevel (national, provincial, and local scales) soil environmental monitoring networks.The promulgation of the action plan for soil pollution prevention and control (State Council of P.R. China, 2016) and the soil pollution control law have advanced the soil environmental monitoring network in China, but many related scientific challenges, such as the number of monitoring sites in the SMN, monitoring intervals, optimal composite subsamples and determination of methods of monitoring indicators, have not been solved effectively in the process of SMN design (Arrouays et al., 2009;Lark et al., 2019;Marchant & Arrouays, 2014).
To clarify the abovementioned scientific problems, many studies have been carried out on SMN design and improved measurements have been proposed in view of the problem existing in the monitoring network after years of trial operation (Chang et al, 2016;Jia, et al., 2021;Morvan et al., 2008;Sparling et al, 2004;Taylor et al, 2010;van Leeuwen et al., 2017;Wu et al., 2023).de Gruijter et al. (2006) described the methods that could satisfy the sampling network design systematically, including simple random sampling, stratified random sampling or systematic random sampling.The number of sampling points remarkably decreased after using an improved annealing model, which provided a theoretical basis for the optimal design of soil property monitoring (Li & Wang, 2019).Chang et al. (2016) determined the optimal sampling density and interval for monitoring soil organic carbon (Corg) changes in degraded Tibetan grasslands by using power analysis.Jia et al. (2021) proposed an optimization plan for the deployment of regional monitoring points according to the result of regional monitoring network of soil geological environment quality in Beijing.Wu et al. (2023) formed the ideas, principles and methods of provincial control basic points by combing the foundation of soil environmental monitoring and the construction of monitoring network in Hainan Province.
The most practicable method in SMN design was provided by Saby et al. (2008), who proposed the theoretical model of minimum detectable change (MDC), which was determined by the standard deviation (SD) of monitoring indicators and the number of monitoring sites.The MDC model takes economic feasibility into consideration, which has enhanced in use in many SMN designs (Chang et al., 2016;Xia et al., 2017).Based on the MDC theoretical model, our study estimate the recommended number of monitoring sites for revisiting, the monitoring interval for a given change, and the MDC under the recommended number of monitoring sites and thus recommend a sampling framework for estimating the status of and trend in soil pH and potentially toxic elements (PTE) concentration in Zhejiang Province.The spatial variation in monitoring indicators was acquired from the Multi-purpose Regional Geochemical Survey in Zhejiang Province (MRGSZ).

Theory and assumptions
Information on spatial variability of monitoring indicators from a single date may be a poor guide to the design of a monitoring scheme.Revisiting soil monitoring sites on a second and subsequent occasion has been considered to be the most effective method for monitoring changes in soil properties (Lark et al., Vol.: (0123456789) 2006).The theoretical basis and calculation methods in our paper conformed to those of Saby et al. (2008).
The monitoring indicator x in any particular monitoring network with n monitoring sites was measured at times t 0 and t 1 , and the mean change ( d ) in x can be estimated as follows: (1) where x i,t is the concentration of monitoring indicator x at monitoring site i at time t.The expression of the standard error of d is

√
s 2 d n , where s 2 d is an estimate of the variance of the differences ( de Gruijter et al., 2006).However, it is impossible to estimate the s 2 d directly because no revisiting was executed at the provincial scale.An alternative estimate of the standard error of d can be expressed as follows: where s 2 is an estimation of the variance of x at the first sampling occasion (baseline).The alternative estimate requires us to abide by the following assumptions: (1) the variance does not change over time; and (2) the correlation between the two sampling times is small (Saby et al., 2008).
Assuming the mean change in the concentration of monitoring indicators conforms to a normal distribution (invoking the central limit theorem), estimates of a 100(1-α)% confidence interval for the mean change can be written as follows: where Z α is the critical value of the standard normal distribution at a given significance level α (using α = 0.05 in this paper), D is the estimates of mean changed.The condition for detecting a mean change y is as follows: Hence, the MDC is as follows: where s is an estimate of the standard deviation of the monitoring indicator on the first sampling occasion.
It is assumed that the concentration of the monitoring indicator changes at an estimated rate of k during the whole time interval t so that the following is true: and the minimum time interval to detect a given rate of change k is as follows: We could estimate the number of monitoring sites required to monitor a certain level of change y by deducing Eq. ( 5) to the following:

Study area
Zhejiang Province is in Southeast China, covering 105,500 square kilometers with a population of over 60 million.It has a marked subtropical monsoon climate, with average temperatures ranging from 3 to 9 ℃ in January to 26-28.8 ℃ in July and a mean annual precipitation of 1600 mm.The areas of hills, plains and water account for 74.6, 20.3 and 5.10% of the total area, respectively (http:// tjj.zj.gov.cn/).Yellow soil and red soil are the main soil types in the hilly mountainous area.Paddy soils are mostly found in the plains and river valleys, while saline soils and desalinated soils are distributed in the coastal areas.
With rapid industrialization and urbanization, large amounts of hazardous substances have been introduced into soil by industrial and agricultural activities, causing serious environmental pollution in Zhejiang Province (Ren et al., 2019).Therefore, it is vital to monitor the soil environmental status to protect soil and food security.

MRGSZ data as the first occasion
The fundamental soil data were derived from 4 (four) MRGSZ projects carried out during 2002-2005, 2003-2008, 2014-2016 and 2015-2016.In view of the unified sampling method and sample preparation (details can be found in Li et al., 2014), analysis method (CGS, 2005) and quality control regulation (Ye, 2005;Zhang, 2005), the data from the 4 MRGSZ projects are comparable and are regarded as the first occasion data for the soil monitoring scheme design.A total of 15,305 samples were obtained from the dominant arable land in Zhejiang Province, which adopted a regular grid of the 2 × 2 km sampling method and covered an area of ~ 6.6 × 10 4 km 2 (Fig. 1).Fifty-two elements (Ag, As, Au, B, Ba, Be, Bi, Br, Cd, Ce, Cl, Cr, Co, Cu, F, Ga, Ge, Hg, I, La, Li, Mn, Mo, N, Ni, Nb, P, Pb, Rb, S, Sb, Sc, Se, Sn, Sr, Th, Ti, Tl, U, V, W, Y, Zn, Zr, SiO 2 , Al 2 O 3 , TFe 2 O 3 , MgO, CaO, Na 2 O, K 2 O, TC), pH, and Corg were analyzed in the composite samples, which provided the environmental quality, nutrient status, and physical and chemical properties of the soil (CGS, 2005;Teng et al., 2014).

Resampled data in 2017 as the subsequent occasion
In total, 1281 soil samples were subsequently collected from randomly selected regular grids of 2 × 2 km within the MRGSZ in 2017 (Fig. 1).Except for the lower sampling density in the revisit in 2017, the griding, sampling, and analytical methods were the same as those in the MRGSZ.The data for the 1281 soil samples were used in this paper to validate the hypothesis made in the methodology.

Validation of the hypothesis proposed in the methodology
According to the description in the "Theory and assumptions", to obtain the standard error of d , we must know the variance of difference ( s 2 d ).However, it is impossible to estimate s 2 d directly because of the shortage of resampled data.Considering that only the resampled data in 2017 were available, the hypothesis that the correlation between two batches of monitoring indicator concentration data is small could be validated with the data (Saby et al., 2008).Based on the SD of the original MRGSZ data and the resampled data in 2017 (Table 2), it can be assumed that the assumption is reasonable for pH, Cr, Ni and Pb, with overestimations of 16.3, 21.3, 18.7 and 8.9%, respectively, but the SD is underestimated by 2.0, 19.1, 4.8 and 14.7% for As, Cu, Zn, and Hg, respectively.As there was only a small amplitude, the estimate of MDC would be comparable to the same methodology that was applied across the provincial scale.

Spatial pattern of pH and heavy metals concentration
There is a high degree of spatial heterogeneity of pH and PTEs in the soil, manifested in the SD varied sharply in each region in the MRGSZ (Fig. 2).For pH, Ningbo had the highest SD value, while Lishui had the lowest SD value.The values varied from 0.44 to 1.35 with a difference of ~ 4 times.For Pb, Taizhou has the highest SD value of 36.78, and Jiaxing had the lowest SD value of 6.70.For Hg, Wenzhou had the highest SD value of 0.60, and Lishui had the lowest SD value of 0.06.There was a tenfold difference between the highest and lowest SD values.For Cr, Cu and Zn, Wenzhou had the highest SD values of 38.58, 30.05 and 104.68, respectively, and Lishui, Lishui, Jiaxing had the lowest SD values of 10.87, 7.22 and 18.12, respectively.The SD value of Cd varied from 0.06 in Lishui to 0.70 in Hangzhou, with a difference of ~ 12 times.The SD values of As varied from 1.47 in Jiaxing to 16.55 in Hangzhou, with a difference of ~ 12 times.Compare the SD of the soil monitoring indicators of different river basins (Fig. 3) also shows a significant difference.For pH, the values varied from 0.64 in Jiaojiang river basin to 1.41 in Yongjiang river basin with a difference of ~ 2 times.The SD values of Cd varied from 0.05 in Aojiang river basins to 0.48 in Qiantangjiang river basins, with a difference of ~ 10 times.While SD of other indicators shows a difference from 2 to 8 times.This phenomenon also appears in the SD of the soil monitoring indicators of different landforms or soil types (Figs. 4  and 5).
Because SD plays a decisive role in the determination of the MDC, the number of monitoring sites and the monitoring intervals in the revisited sampling, the significant difference in SD values for monitoring indicators in the study regions is bound to have a significant impact on the practicability of monitoring changes by using the revisiting method.
MDC under sample numbers of MRGSZ Equation ( 5) was adopted to calculate the MDC of monitoring indicators, in which the number of n sites was the number of samples in each MRGSZ (n 0 ) (Table 3).Because a high-density sampling scheme was used in the MRGSZ projects, the MDC values under the scenario of n 0 samples adopted in most administrative regions were low enough to be accepted for monitoring changes in As, Hg, Cd, Cr, Cu, Ni, Pb and Zn contents compared to the risk screening values in the "Soil environmental quality risk control standard for soil contamination of agricultural land (GB 15618-2018)" in China (CMEE,

Estimation of recommended number of monitoring sites (n r ) for revisiting
The number of monitoring sites are needed in the revisits is a significant factor that can strongly affect accuracy of soil monitoring (Khaledian & Miller, 2020;Long et al., 2018).Molla et al. (2022) proposed that adding 350 new monitoring points to the existing sampling design by spatial simulated annealing algorithm increased the prediction accuracy of PTE in urban green space soil by 64.35%.Based on Eq. ( 8), the number of monitoring sites required is controlled by (1) the SD (s) of monitoring indicator concentrations and (2) a given MDC (y).In this paper, a relative change of 10% in the provincial mean of monitoring indicator concentration was defined as a primary value of y.
We can see from Table 4 that there was a wide range of monitoring sites required for revisiting in each region.The number of As varied from 51 in Lishui to 1870 in Hangzhou.For Ni, Cu, Zn, Cr, Cd, Pb and Hg, n r varied from 47, 43, 54, 36, 48, 66 and 164 in Jiaxing, Lishui, Huzhou, Jailing, Jiaxing, Obviously, it may not be feasible to use the minimum number of monitoring sites for some administrative regions under the given primary value of MDC, as it does not consider the economic efficiency in sampling and analysis.Therefore, a higher secondary value was assigned to y during the calculation to reach a compromise between monitoring precision and economic feasibility.The higher secondary value of y was defined as 0.03 mg kg −1 and 0.03 mg kg −1 for Cd and Hg, respectively.The n r was recalculated, and the results are shown in Table 4. Finally, the largest accepted value of each monitoring indicator, which compromised the practicability and monitoring precision, was considered as nr.The difference between the n r and n 0 was show in Fig. 6, which shown that the n r were reduced in almost all regions except Hangzhou and Wenzhou.We also compared with the interpolation maps of As in Jinhua and Cd in Ningbo under scenario of n 0 and n r (random sampling from n 0 ) (Table 5).The two interpolation maps of As and Cd have kappa correlation coefficients are 0.72 (p < 0.01) and 0.78 (p < 0.01), respectively, display similar hot spots, although there are differences in the extents of the high-value regions, which reflect that the n r is able to capture the spatial patterns of PTEs content of the first occasion(Fig.7).

MDC under recommended number of monitoring sites
If the recommended number and locations for detecting concentration of soil pH and heavy metals were revisited for resampling, the MDC of soil pH and heavy metals in different administrative region is shown in Table 6.We can see that the MDC values of heavy metals in each administrative region are acceptable compared to the risk screening values in the "Soil environmental quality risk control standard for soil contamination of agricultural land (GB 15618-2018)" in China (CMEE, 2018), as the n r values adopted in this study were optimized by the SD.

Monitoring temporal interval for a given change
Besides number of monitoring sites, sampling design in time, i.e., when to make observations is another key issue for the implementation of monitoring programs (Brus & de Gruijter, 2013).According to Eq. ( 7), the calculation of the monitoring temporal interval (t) was controlled by (1) an estimated rate of change (k) of monitoring indicators; (2) the SDs of monitoring indicators in the monitoring region; and (3) the number of monitoring sites (n) adopted in the revisit.
In this study, k was obtained from Wang et al. (2010) for Ni, from Zhou et al. (2004) for As, Cu, Zn, Cr and from Xia et al. (2017) for Cd, Pb and Hg.Regarding the number of monitoring sites in revisits, some scholars (Xia et al., 2017) have suggested a ratio of the sample number in the MRGSZ.The minimum monitoring interval ranges from several years to decades depending on the monitoring indicator in different administrative regions (Fig. 8).It is not practical for some regions, such as Hangzhou, Huzhou and Wenzhou, to use the scheme because the temporal interval is too long.Therefore, the number of monitoring sites adopted in revisiting must be optimized, especially for the regions that have high SD values.There should be more monitoring sites in those regions to reduce the monitoring interval to a practical level.
According to Eq. ( 8), the n r in Table 4 is the one that has been optimized by the SD, as the higher the spatial heterogeneity of monitoring indicators is, the larger the n r is.Therefore, the temporal monitoring interval was recalculated under the scenario of n being equal to n r , and the results are shown in Table 5 and Fig. 9.The monitoring temporal interval ranged from 1.3 to 14.8 years for the monitoring indicators.For each administrative region, the maximum value of monitoring indicators was recommended as it met the requirements of all the monitoring indicators.The most recommended monitoring temporal intervals are longer than 10 years, varying from 10.0 to 14.8 years, which is approximately equal to "carrying out at least one national soil pollution survey every 10 years" stated in the soil pollution control law of China.

Limitations and prospect
There is no standardized approach for evaluating the optimum monitoring site number and monitoring interval due to the high spatial heterogeneity and small temporal variability relative to PETs content in soil.Currently, the determination of reasonable sampling numbers for monitoring PTEs in soil is mainly include empirical et al., 2018; Zhu et al.,  The unit of the pH is dimensionless, the others is mg•kg −1 .Values with underlines are unacceptable ones with reference of "Soil environmental quality Risk control standard for soil contamination of agricultural land (GB 15618-2018) 2008), statistical (Lin et al., 2018;Szatmári et al., 2019) and empirical statistics (Tenenbein, 2017;Xie et al., 2016).Among them, traditional statistical methods, such as power analysis (Liu et al., 2021) and Cochran method (Liu et al., 2021), have been well applied in many fields due to its ability to solve the problem of replacing the spatial mean and spatial variation of PTEs in soil with prior knowledge, thus meeting the optimal sampling quantity design of soil monitoring networks.Although the MDC method, an alternative expression of statistical power, has several benefits over other techniques, the prerequisite for its successful application depends on the following assumptions being satisfied by the changes of PTEs concentration in soil: (1) a normal distribution of the mean change of PETs concentration, (2) repeated sampling is carried out by the same method, (3) that the variance of PETs concentration remains the same on successive sampling occasions, (4) the variation of PETs within a unit can estimated using the MRGSZ.These assumptions have been investigated as far as possible using the available data in this paper, and it has been shown that a reasonable estimate of the MDC for soils, in each administrative region, can be made.
There are also has limitation of scale issue in this study, which could result in different baseline and increase rate of soil pH and PTEs.Shao et al. (2016) found that the increase rate of agricultural topsoil Cd in the Zhejiang province have similar amplitude to that in the whole country, but it is significantly lower the one in Yangtze River Delta Region.On a broad scale (> 10,000 km 2 ), sampling density tended to decrease with increasing spatial extent (Chen et al., 2022).Chinese government has made great efforts to cut off the sources of pollution since the passing of the Soil Pollution Prevention and Control Law of China in 2018, and relatively effective restoration strategies have been adopted and promoted on the contaminated land.Following this scenario, a smaller variance of soil pH and PTEs could occur at resampling, and sampling schemes based on baseline variance alone may underestimate the true ability to detect statistically significant change in pH and concentration of PTEs.On the other hand, there may be high variability in the implementation of management practices by farmers at field scale, which may increase the variance of soil pH and PTEs at resampling.So, it is harmful to assume identical variance of the baseline and resampled data.
In addition, the temporal change rate of PTE content (k) were quoted from relevant literatures, which calculated using temporal trend analysis based on arithmetic mean, being neither median nor geometric mean.However, to date, no reliable method for the calculation of provincial mean values of soil pH and PTE is known (Li et al., 2020), which could lead to the calculation results of the change rate inaccurate.Huang et al (2019) proposed that the content data of PTEs collected in each administrative unit should be weighted according to many factors, such as the size of the study area, but these factors are quite difficult to quantify.Further research on variability of pH value and PTEs concentration changes over time at the different spatial scales is required.
In view of these, it is recommended to conduct research on the method for determining the n r in the  future from the following aspects: (1) Before determining the n r , it is not only necessary to pay attention to the spatial variation of PETS in soil but also to clarify the soil sampling scale or sampling unit.(2) Due to the lack of a reasonable sampling method for determining the n r of PTEs in soil under lognormal distribution, there is a significant difference in the number of samples for PTEs under the same variation and accuracy requirements, and there is a lack of unified and effective methods for determining the n r as guidance.Therefore, starting from statistics and based on comparative analysis, exploring methods for determining the reasonable number of samples under skewed distribution, especially lognormal distribution and applying them to the monitoring process of soil PTEs will become a new trend in soil monitoring.

Conclusions
Using the prior knowledge of the spatial pattern of soil pH and PTE content based on the data obtained from the MRGSZ, Our research highlights the feasibility of using preliminary surveys to estimate the necessary monitoring sites to obtain a desired level of precision and to estimate the monitoring interval required to detect an expected change rate when implementing a paired sampling on the provincial scale.Further work is needed to determine the variability of pH value and PTE concentration changes over time at the different spatial scales.The conclusions are summarized as follows.
1.There was a significant difference in the monitoring indicator concentrations in the MRGSZ regions in Zhejiang Province, China.The highest SD values were ~ 3. 08, ~ 11.26, ~ 10.94, ~ 3.55 , ~ 11.66, ~ 4.16, ~ 6.07, ~ 5.49 and ~ 5.78 times the lowest SD values for pH, As, Hg, Cr, Cd, Cu, Ni, Pb and Zn, respectively.This result was bound to result in the high spatial heterogeneity of the MDC, monitoring site assignment, and monitoring temporal intervals in different administrative regions.2. Under the scenario of the same number of samples as MRGSZ (n 0 ), the MDC values of Cu, Zn and Pb were not acceptable in some regions, as their precision was lower than the standard issued in the "Soil environmental quality risk control standard for soil contamination of agricultural land (GB 15618-2018)" in China.In addition, the minimum temporal intervals were not practical for Hangzhou, Huzhou and Wenzhou, as they were unbearably long at a ratio of 1/16 MRGSZ (n 0 ). 3.There was a wide variation in the number of monitoring sites (n r ) for subsequent visits under a simulated relative change of 10% in the provincial mean concentration of monitoring indicators.
Under n r , the MDC values were acceptable for all administrative regions, and the temporal monitoring intervals were practical, with a range of 6.7 to 14.8 years.Because of the well-realized unification of economic feasibility, practicability and monitoring precision, the n r was recommended for the design of SMNs.
Data processingA total of 15,305 sampling points with coordinates were overlaid with a polygon file of the administrative region of Zhejiang Province by using GIS software ArcGIS 10.2, and the sampling points in each administrative region were obtained.The basic statistical parameters of the monitoring indicator concentrations of the sample points in each administrative region, including standard deviation (SD) and variance, were calculated and input to Microsoft Excel 2007 for the following analysis.The 1281 soil samples in 2017 were linked to those in the MRGSZ using the point-to-point "spatial join" function in ArcGIS 10.2 to form paired samples for statistical analysis of changes in monitoring indicator concentrations.

Fig. 1
Fig. 1 Sketch map for the soil sampling locations in the two dates

Fig. 2
Fig. 2 SD for monitoring indicators concentration in each administrative region calculated by the indicators concentration from MRGSZ (The unit of the pH is dimensionless, the others is mg/kg)

Fig. 3
Fig. 3 SD for monitoring indicators concentration in different river basin calculated by the indicators concentration from MRGSZ

Fig. 4
Fig. 4 SD for monitoring indicators concentration in different landforms calculated by the indicators concentration from MRGSZ

Fig. 5
Fig. 5 SD for monitoring indicators concentration in different soil types calculated by the indicators concentration from MRGSZ

Fig. 7
Fig. 7 Comparison of the interpolation maps of As in Jinhua (a, b) and Cd in Ningbo (c, d) under scenario of n 0 and n r

Fig. 8
Fig.8Temporal intervals (years) calculated using 1/16 MRGSZ sample numbers for each region (The unit of the pH is dimensionless, the others is mg/kg/a.)

Fig. 9
Fig. 9 Minimum temporal intervals (years) under n r for revisit

Table 1
Comparison of main monitoring parameters in different SMNs

Table 2
Estimates of SD of the monitoring indicator concentrations using two visits data in Zhejiang province

Table 4
The number of samples (n r ) needed for revisits