Assessment of land ecosystem health with Monte Carlo simulation: A case study in Qiqihaer, China

Journal Pre-proof Assessment of land ecosystem health with Monte Carlo simulation: A case study in Qiqihaer, China Yijia Yang, Ge Song, Shuai Lu PII:

S0959-6526(19)34392-6

DOI:

https://doi.org/10.1016/j.jclepro.2019.119522

Reference:

JCLP 119522

To appear in:

Journal of Cleaner Production

Received Date: 25 February 2019 Revised Date:

28 November 2019

Accepted Date: 29 November 2019

Please cite this article as: Yang Y, Song G, Lu S, Assessment of land ecosystem health with Monte Carlo simulation: A case study in Qiqihaer, China, Journal of Cleaner Production (2019), doi: https:// doi.org/10.1016/j.jclepro.2019.119522. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Author Contribution Statement

Yijia Yang: Writing-Original Draft, Methodology, Conceptualization Ideas，Software, Formal analysis;

Ge Song: Writing-Review & Editing, Supervision, Funding acquisition;

Shuai Lu: Formal analysis, Data Curation.

Title : Assessment of land ecosystem health with Monte Carlo simulation: A case in Qiqihaer, China Authors: Yijia Yanga Ge Songa,b* Shuai Lua a College of Resources & Environment, Northeast Agricultural University, Harbin, 150030, China b Institute of Land Management, Northeastern University, Shenyang, 110169, China Corresponding Author：Ge Song； Xiangfang District, Changjiang Road No. 600 Harbin 150030 China Cell phone:+86-13945058282 Email: [email protected]

1

8301

2

Assessment of land ecosystem health with Monte Carlo simulation: A case study in

3

Qiqihaer, China

4

Yijia Yanga

5

a College of Resources & Environment, Northeast Agricultural University, Harbin,

6

150030, China

7

b Institute of Land Management, Northeastern University, Shenyang, 110169, China

8

Abstract： ：A better understanding of the level of land ecosystem health in a region is helpful for

9

policymakers in developing measures for eco-space management. During the evaluation process,

10

uncertainty is inevitable because of different methods of determining the related weights. In this study,

11

a Monte Carlo simulation was used to construct a sample of the indicator weights to quantify the

12

uncertainty in the process of land ecosystem health assessment to improve the accuracy of the

13

assessment result. This study used a new evaluation framework of pressure, state (vigor-organization

14

-resilience-function) and response to assess the level of land ecosystem health in Qiqihaer City from

15

2000 to 2015. Then, the spatial heterogeneity among the levels of land ecosystem health was examined

16

using hotspot spatial analyses. Our results showed the following: (1) the overall level of land

17

ecosystem health in Qiqihaer city was ordinary (Ⅲ), and the average values of the related indices were

18

0.554 (2000) and 0.563 (2015). Hotspots and coldspots identified by the land ecosystem health index

19

had obvious overlapping areas from 2000 to 2015, among them, approximately 30.6 % of the hotspot

20

area (high-value area) was located in the Zhalong Nature Reserve, and approximately 15.6 % of the

21

coldspot area (low-value area) was located in city centres and the boundary of the study area. (2) The

22

average values of the land ecosystem health pressure and response indexes were low, while the average

Ge Songa,b Shuai Lua

1

23

values of the land ecosystem health states index changed greatly. However, in terms of spatial

24

heterogeneity, the spatial differences were small for the pressure index and the states index, and the

25

spatial difference in the response index was large. The assessment results objectively reflect the health

26

status of the regional land ecosystem and have very important theoretical and practical significance for

27

safeguarding similar regional ecological security and promoting social and economic sustainable

28

development.

29

Keywords: Land ecosystem health; Pressure-state (vigor-organization-resilience-function)-response;

30

Confidence limits; Hotspot spatial analyses.

31

1 Introduction

32

With rapid urbanization and industrialization, human activities affect the structure (Li et al., 2016a)

33

and function of land ecosystems to a certain extent (Qiu et al., 2015), which lead to severe ecosystem

34

degradation and threat to the sustainable development of the social economy (Cheng et al., 2018). Land

35

ecosystem health (LEH) has become a focus issue (He et al., 2019), especially in the northwestern part

36

of the Songnen Plain in China, where there are ecological problems such as land desertification and

37

soil erosion. It is urgent to evaluate LEH to develop an effective the programme for eco-space

38

management.

39

LEH refers to a land ecosystem with a strong self-repair ability under long-term natural and

40

human disturbances (Xiao et al., 2019) that can resist external interference and maintain the stability of

41

the system state (Wu et al., 2018). The study of ecosystem health mainly focuses on the research scale

42

of ecosystem health, the framework for ecosystem health assessment and the determination of weight

43

of indicators (Mariano, 2018). (1) Many studies on the ecosystem health assessment have been

44

conducted at various scales, but the scale of the results in existing research is rarely based on grids. For 2

45

example, most of the existing research is focused on countries (He et al., 2019), provinces (Meng et al.,

46

2018), cities (Wang et al., 2018), rivers (Zhao et al., 2019), wetlands (Chi et al., 2018) and forests

47

(Ishtiaque et al., 2016). (2) The land ecosystem is a large and complex system, and the level of LEH is

48

mostly evaluated by establishing frameworks for ecosystem health assessment. At the same time,

49

existing frameworks do not comprehensively quantify the level of LEH (the natural status of the land

50

ecosystem and the interaction between the land ecosystem and human activities). For example, the

51

vigor-organization-resilience (VOR) (Yan et al., 2016) and the vigor-organization-resilience-function

52

(VORF) (Xiao et al., 2019) were mainly used to measure the status of the ecosystem, and the

53

interaction between the ecosystem and human activities was ignored. The frameworks of

54

pressure-state-response (PSR) (Sun et al., 2019) and driving force-pressure-state-impact-response

55

(DPSIR) (Flint et al., 2017) emphasize the causal relationship between human activities and changes in

56

the land ecosystem. (3) In the process of LEH assessment, determining the weight of the indicators is

57

particularly important, and it substantially affects the result of the LEH assessment. There are many

58

methods (subjective, objective and combination of subjective and objective) for determining the weight

59

of indicators, such as the entropy method, analytic hierarchy process (AHP), and coefficient of

60

variation method (Song et al., 2017). Due to the differences between the methods, inevitably,

61

uncertainty resulting from the determination of weights is introduced into the final evaluation result.

62

However, the Monte Carlo simulation can help to quantify the uncertainty in determining weights, and

63

the simulation is based primarily on the triangular probability distribution. The model was determined

64

by the Intergovernmental Panel on Climate Change (IPCC, 2006), and it minimizes parameters (Song

65

et al., 2015), and has been widely used in many fields, such as biology, sociology and ecology. For

66

example, the model can be used to quantify environmental efficiencies under uncertainty (Ewertowska 3

67

et al., 2017) and economic risk analysis (Zaroni et al., 2013).

68

In summary, LEH does not yet have a mature and operational assessment system, and the existing

69

research results needs to be supplemented and improved based on construction of the framework for

70

ecosystem health assessment and the determination of weights. Thus, this paper combines the

71

evaluation frameworks of VORF and PSR, establishes a new LEH evaluation system (P-S(VORF)-R),

72

and applies the Monte Carlo simulation model to LEH for the first time to determine the weight of

73

indicators.

74

Qiqihaer City in the northwestern part of the Songnen Plain is a typical area for LEH research due

75

to its ecological problems and land-use status. This study assessed and analyzed LEH levels over at

76

temporal and spatial scales and was conducted in Qiqihaer City using GIS, Python and detailed

77

analysis. More specifically, our goals were to study the process of LEH from 2000 to 2015 at different

78

temporal and spatial scales, explore the threshold range of LEH through Monte Carlo simulation model,

79

and improve the LEH level in Qiqihaer City to help policymakers develop measures for eco-space

80

management.

81

2 Materials and Methodology

82

This section mainly introduced the research area overview, data sources and methodology (the

83

LEH assessment).

84

2.1 Study area

85

This research uses Qiqihaer in Heilongjiang Province in the northwestern region of the Songnen

86

Plain as the study area; specifically, this area ranges from 45° to 48° N and from 122° to 126° E (Fig.

87

1). The region is in the middle temperate zone, has a temperate continental monsoon climate, with an

88

average annual temperature of 3.6°C and an average annual rainfall of 415 mm, and its mineral 4

89

resources are abundant (non-metallic minerals are dominant). The area is dominated by plains and hills,

90

with higher areas in the northeast and west, including Qiqihaer City, Nehe City and 7 counties (Ganan

91

County, Longjiang County, Fuyu County, Yi’an County, Baiquan County, Keshan County, and Kedong

92

County), with an area of approximately 42.3 thousand km2.

93

The study area has typical ecological problems, such as land desertification and soil erosion. This

94

paper explores the health status of the land ecosystem in the study area, and the research results will

95

help to improve the regional ecological environment and include measures for eco-space management.

96

In brief, the selection of areas with typical ecological problems has certain research value.

97 98 99

Fig. 1 Location of the eastern part of Heilongjiang Province, China 2.2 Data source and processing

100

Two forms of data were used in this research, namely, spatial data and statistical data.

101

Spatial data: These data were used to represent information such as the geographical location, shape, 5

102

and attached property information. The spatial data in this study consisted of the following data: first,

103

the administrative boundary map of Qiqihaer, the GDP grid data and land-use data, which were

104

procured from the Data Centre for Resources and Environmental Sciences, Chinese Academy of

105

Sciences (RESDC); second, the normalised vegetation index (NDVI) data, which were downloaded

106

from the Geospatial Data Cloud. The time nodes of the above data were both 2000 and 2015.

107

Statistical data: The social and economic data were mainly from the Qiqihaer Statistical Yearbook,

108

Statistical Bulletin of the National Economy and Social Development of Qiqihaer (county) and

109

Environmental Statistics Bulletin. The time nodes of the above data were both 2001 and 2016. Part of

110

the evaluation indicator data was obtained from calculations of the original data, and these data were

111

obtained by kriging interpolation.

112

The spatial analyst of ArcGIS is a data analysis tool that is based on the location and morphology of

113

geographical objects and provides technical support for research. In this study, all the collected data

114

were converted to the raster data format using a raster cell of 1 km*1 km.

115

2.3 Research process

116

This paper, using Qiqihaer in the northwest of Songnen Plain as the study area, proposed the

117

framework of P-S (VORF) -R for assessing LEH. This study used five weighting methods to construct

118

weight samples and then used the Monte Carlo simulation model to repeat the calculation 500 times.

119

Finally, this study obtained the land ecosystem health index (LEHI). A flow chart of this procedure is

120

illustrated in Fig. 2.

6

121 122

Fig. 2 Flow chart for the assessment of LEH (land ecosystem health pressure (LEHP), land ecosystem

123

health state (LEHS), and land ecosystem health response (LEHR))

124

2.3.1 Index system construction

125

In addition to a broad literature review, the index framework was established based on the

126

principles of data that were quantifiable, available, objective, dynamically predictive and representative.

127

All these indicators were divided into three categories including pressure, state, and response. The

128

pressure subsystem describes the impacts that human activities have on a land ecosystem and has

129

indicators mainly based on population, environment and economic stress. The state subsystem

130

describes the natural state of the land ecosystem under multiple stresses, characterized mainly by four

131

factors: vigor, organization, resilience and function. The response subsystem describes the actions of

132

policymakers and managers under pressure and status indicators, primarily in terms of environmental

133

governance, living standards of farmers, population and food production (Table 1).

134

7

135

Table 1

136

The index system of LEH of study area Indicator

Ecosystem

Subsystems

No.

Indicators

Units

Calculation method character

P1

P2

Density of population（DP） Human activity density index

Person/km2

-

/

-

%

-

Farmland area/land area

%

-

Construction land area / Land area

kg

-

Statistical data

Population/Land area

（HAD） P3

Rate of land reclamation (LR)

LEHP P4

P5

Proportion of construction land (CP) Farmland use of pesticides and fertilizer (PF)

P6

Economic density (ED)

Yuan/km2

-

land output value of 1km2

S1(V)

NDVI

/

+

Spatial data

S2(O)

Landscape diversity index (LD)

/

+

S3(R)

Ecological resilience index (ER)

/

+

S4(F)

Ecosystem services value (ESV)

yuan

+

yuan

+

Statistical data

%

+

Statistical data

%

+

Statistical data

kg

+

Statistical data

Yuan/person

+

Statistical data

‰

+

Statistical data

Land LEHS ecosystem

E = 0.3* Resil + 0.7 * Resist

health

Total investment in the treatment R1 of environmental pollution (TIE) R2

Forest cover rate (FC) Attainment rate of the industrial

R3 LEHR

waste water discharged (IWW) R4 R5

R6

137 138 139 140

Grain Production (GP) Per capita net income of rural residents (IRR) Natural population growth rate (NPG)

(Note: For the landscape diversity (LD) index, the analysis of this indicator was performed using the Shannon index (S) (Velázquez et al., 2019); n is the number of land-use types, ai is the i-th land-use type area, pi is the human activity density parameter of the i-th land-use type, vi is the ESV per unit area of the i-th land-use type, and A is the total area of the evaluated unit.) 8

141 142

In related research on ecological elasticity, there are some differences in ecological elasticity

143

scores, but the order of the elastic scores of each land-use type are essentially the same. This paper

144

refers to the research results of Peng et al (2017), and combines the land-use and ecological

145

environment status of study area to determine the ecological elasticity scores of different land-use

146

types (Table 2).

147

Table 2

148

Ecological elasticity coefficient of land-use types in study area Cultivated Land use type

Forest land

Grassland

Construction

Unused

land

land

Water body

land Resilience coefficient

0.3

0.5

0.8

0.7

0.2

1.0

Resistance coefficient

0.6

1.0

0.7

0.8

0.3

0.2

Elasticity coefficient

0.51

0.85

0.73

0.77

0.27

0.44

149 150

This paper ignores the difference in the utilization intensities of the same land-use type, and

151

scholars (He et al., 2015) set the human activity density parameters of different land-use types in the

152

study area (Table 3).

153

Table 3

154

Human activity density parameters of different land-use types Cultivated Land use type

Unused Forest land

Grassland

Construction land

Water body

land Parameters

0.55

land 0.1

0.23

155 9

0.95

0.115

0.14

156

This paper uses ESV to indicate that a land ecosystem provides service functions for human society.

157

The calculation process of ESV per unit area is found in Section S4 of Appendix A (Table 4).

158

Table 4

159

Coefficients of the ESV in study area Unit:Yuan(RMB)/hm2

160

Ecosystem Wetland

water

Cropland

Forestland

services

161 162

Construction

Unused

land

land

Grassland

FP

285.45

420.24

792.91

261.66

340.95

0

15.86

RM

190.3

277.52

309.24

2362.88

285.45

0

31.72

GA

1910.92

404.38

570.9

3425.38

1189.37

0

47.57

CL

10743.95

1633.4

769.12

3227.15

1236.94

0

103.08

WA

10656.73

14882.95

610.54

3243.01

1205.23

0

55.5

WT

11417.92

11774.73

1102.15

1363.81

1046.64

0

206.16

SFR

1577.89

325.09

1165.58

3187.5

1776.12

0

134.79

BD

2925.84

2719.69

808.77

3576.03

1482.74

0

317.16

RCT

3718.75

3520.53

134.79

1649.26

689.83

0

190.3

Total

43427.75

35958.53

6264

22296.68

9253.27

0

1102.14

（Note: gas regulation (GA), climate regulation (CL), water regulation (WA), soil formation and retention (SFR), waste treatment (WT), recreation cultural and tourism (RCT),

biodiversity (BD), food prod- uction (FP) and raw material (RM)）

163 164

2.3.2 Standardization of indicators

165

Because it is difficult to combine all the variables to evaluate LEH, this study needed to

166

standardize the indicators. The formula for normalizing indicators is as follows (Feng et al., 2017):

10

X ij − min( Xi )

167

Positive indicator: Yij =

168

Negative indicator: Yij =

169

where Xij is the actual value of the indicator in the i-th year, min(Xi) is the minimum value,

170

max(Xi) is the maximum value, and Yij is the normalized value. After standardization, the values were

171

between 0 and 1.

172

2.3.3 Method of evaluating LEH

173

max( X i ) − min( X i ) max( X i ) − X ij max( X i ) − min( X i )

(7)

Step 1, construct a weight sample. In this paper, this study used AHP (Ameen & Mourshed, 2019),

174

maximum deviation (Li et al., 2018), entropy method (Zhao et al., 2018), mean square deviation

175

method (Pereira & Vasquez, 2017) and coefficient of variation method (Qian et al., 2014) to determine

176

the weight of the indicators and constructed a weight sample.

177 178 179

Step 2, establish a triangular probability distribution. The minimum, maximum and average of the weight index formed the triangular probability distribution of the weight (IPCC, 2006) (Fig. 3). Step 3, evaluate LEH. Based on the weighted sample and indicator normalization results, this

180

paper performed a Monte Carlo simulation (500 repetitions) with the help of the Python 3.7 platform

181

(Section S2 of Appendix A) (Song et al., 2015). The average of the results of the 500 repetitions was

182

the LEHI, the land ecosystem health stress index (LEHPI), the land ecosystem health status index

183

(LEHSI), and the land ecosystem health response index (LEHRI). This study further ranked the results

184

by five levels with the natural break method to reflect the extent of LEH (Table 5). Considering the

185

actual management needs, the breakpoint value was taken to one decimal place (He et al., 2015).

11

186 187

Fig. 3 The weight of the 16 variables used for the Monte Carlo simulations

188

Table 5

189

Grade of LEH in the study area Level Ⅰ Ⅱ Ⅲ Ⅳ Ⅴ

Health status Well Relatively well Ordinary Relatively weak Weak

Comprehensive evaluation value [0.7, 1] [0.6, 0.7) [0.5, 0.6) [0.4, 0.5) [0, 0.4)

190 191

2.3.4 Hotspot and coldspot spatial analyses of LEH

192

Hotspot spatial analyses have been widely used in ecological fields to help determine the spatial

193

locations of hotspots and coldspots in data, namely, spatial clusters of high and low values, respectively.

194

In this paper, ArcGIS 10.2 software (hotspot analysis tool (Getis-Ord Gi*) was used to explore the 12

195

spatial distribution of LEH, LEHP, LEHS and LEHR (2000a and 2015a) (Li et al., 2016b).

196

3 Results and Analysis

197

To gain a deeper understanding of the LEH in the study area, this paper analyzed the levels of

198

LEH, LEHP, LEHS and LEHR, and their uncertainty.

199

3.1 Comprehensive LEH analysis

200 201 202

This part mainly analysed the level of LEH in the study area at temporally and spatially. 3.1.1 LEH The average value of the LEH in the study area was 0.554 (2000) and 0.563 (2015), indicating that

203

the status of the LEH in the past 15 years was generally stable and slightly increased. Among the

204

various levels of LEH, the area that was considered at the relatively weak ( ) level decreased the most,

205

from 22.20 % (2000) to 11.30 % (2015); the increase in the area that was considered at the ordinary ( )

206

level increased the most obvious, from 50.33 % (2000) to 61.39 % (2015). The area of that was

207

considered well ( ), relatively well ( ) and ordinary ( ) in 2000 and 2015 was more than 70 %,

208

indicating that the overall level of the LEH in the study area was relatively good (Table 6 and Fig. 4).

209

The phenomenon of LEH with simultaneously improved and degraded conditions existed between

210

2000 and 2015 (Fig.4). The area where the level of LEH was well (LEH level turning well) accounted

211

for 22.61 % of the study area and was concentrated in the eastern part (Baiquan, Keshan and Kedong)

212

and the western part (Longjiang and Qiqihaer) of the study area. The area where the level of LEH was

213

weak (LEH level turning weak) accounted for 11.09 % of the study area and was concentrated in the

214

southwestern part of the study area (Tai Lai and Longjiang) and the northern part (Gannan and Nehe).

215

Table 6

216

Overview of the area and proportion of the level of LEH in the study area from 2000 to 2015 13

Ⅰ

Ⅱ

Ⅲ

Ⅳ

Ⅴ

Year Area

Porportion

2

(km ）

(%)

Area

Porportion

2

(km ）

(%)

Area

Porportion

2

(km ）

(%)

Area 2

(km ）

Porportion (%)

Area

Porportion

2

(km ）

(%)

2000

728.62

1.72

10764.09

25.49

21250.05

50.32

9374.66

22.4

106.85

0.25

2015

964.16

2.28

10456.22

24.76

25920.59

61.39

4773.33

11.4

109.97

0.26

217 218

Fig. 4 Spatial distribution pattern of LEH level and their level changes in the study area from 2000 to

219

2015

220

3.1.2 LEH spatial heterogeneity

221

Fig. 5 shows the spatial distribution of hotspots and coldspots based on the LEHI. For 2000 to

222

2015, the spatial distributions of hotspots and coldspots were roughly the same. Hotspots identified by

223

the LEHI had obvious overlapping areas from 2000 to 2015, and approximately 30.6 % of the

224

identified hotspots were located in Zhalong Nature Reserve. The spatial distribution of the coldspot

225

areas (2000) was relatively dispersed, and approximately 15.6 % of the identified coldspots were

226

located in city centres and the boundary of the study area.

14

227 228

Fig. 5 Spatial heterogeneity of the LEHI in the study area in 2000 and 2015

229

3.2 Comprehensive LEHP, LEHS and LEHR analysis

230

To more fully explore the level of LEH in the study area, we needed to analyse not only the

231

current status of the entire system but also the current status of the subsystem. Thus, this part of the

232

study mainly analysed the levels of LEHP, LEHS and LEHR in the study area at temporally and

233

spatially.

234

3.2.1 LEHP, LEHS and LEHR

235

From the perspective of spatial patterns, there were significant differences in LEHP, LEHS and

236

LEHR in the study area (Fig.10). The mean LEHP values were 0.229 (2000) and 0.227 (2015), and the

237

average LEHS values were 0.125 (2000) and 0.135 (2015). The average LEHR values were 0.199

238

(2000) and 0.200 (2015) (Fig. 6).

15

239 240 241

Fig. 6 Spatial distribution pattern of the values of LEHP, LEHS and LEHR in 2000 and 2015 From 2000 to 2015, the overall change in the LEHP in the study area was small. The average of

242

DP (P1), HAD (P2), LR (P3), CP (P4) and PF (P5) showed a decreasing trend, while ED (P6) showed

243

an increasing trend. For example, the mean value of HAD (P2) decreased from 0.600 to 0.510, the

244

mean value of LR (P3) decreased from 0.312 to 0.282, the mean value of PF (P5) decreased from

245

0.506 to 0.467, and the mean value of ED (P6) increased from 0.870 to 0.983. The low-value area of

246

LEHP was the economic core area of the study area. Compared with other indicators of pressure in the

247

study area, the gradual development of human activities (farming), urban construction, population and

248

economy had had greater pressure on the LEH in the study area (Fig. 7 and Fig. 8).

249

From 2000 to 2015, the LEHS in the study area varied greatly. The average values of NDVI (S1)

250

and ER (S3) showed an increasing trend, while the mean values of LD (S2) and ESV (S4) showed a

251

decreasing trend. For example, the mean value of ER (S3) increased from 0.468 to 0.563, and the 16

252

average value of ESV (S4) decreased from 0.243 to 0.203. The high-value areas of LEHS were mainly

253

distributed in the middle, west, and northeast of the study area; for example, the Zhalong Nature

254

Reserve, which was located in a high-value area, mainly contained wetlands and woodlands. These

255

land-use types were conducive to increasing the values of ESV (S4), ER (S3) and NDVI (S1) (Fig. 7

256

and Fig. 8).

257

From 2000 to 2015, the spatial heterogeneity of LEHR in the study area was significant, mainly

258

due to the combined effects of six indicators. Among them, the average values of IRR (R5) showed a

259

decreasing trend, while, the average values of TIE (R1), FC (R2), IWW (R3), GP (R4) and NPG (R6)

260

showed an increasing trend; for example, the mean value of TIE (R1) decreased from 0.308 to 0.422,

261

the mean value of GP (R4) decreased from 0.471 to 0.545, and the mean value of NPG (R6) increased

262

from 0.431 to 0.506. The high-value area of LEHR was distributed in the southeastern part of the study

263

area. Because the southeastern part (Baiquan, Keshan and Kedong) is a typical demonstration site for

264

soil and water conservation, the investment in environmental pollution control was high, and the

265

artificial afforestation area increased (Fig. 7 and Fig. 8).

17

266 267

Fig.7 The spatial distribution of 16 variables (standardized value) in 2000

18

268 269

Fig.8 The spatial distribution of 16 variables (standardized value) in 2015

270

3.2.2 LEHP, LEHS and LEHR spatial heterogeneity

271 272 273 274

Fig. 9 shows the extensive spatial distribution of hotspots and coldspots based on the LEHP, LEHS and LEHR. For the LEHP, the coverage area of hotspots decreased from 16.7 % to 16.2 %, and the coverage area of coldspots decreased from 9.5 % to 6.4 %. Hotspots identified by the LEHPI had obvious 19

275

overlapping areas from 2000 to 2015, and approximately 74.2 % of the identified hotspots were located

276

in Qiqihaer (Zhalong Nature Reserve) and marginal areas of the study area. The coldspots identified by

277

the LEHPI had obvious overlapping areas from 2000 to 2015, and approximately 57.3 % of the

278

identified coldspots were located in Longjiang, Nehe and the city centre of the study area.

279

For the LEHS, the coverage area of hotspots decreased from 15.9 % to 15.2 %, and the coverage

280

area of cold spots increased from 1.6 % to 1.9 %. Hotspots identified by the LEHSI had obvious

281

overlapping areas from 2000 to 2015, and approximately 72.6 % of the identified hotspots were located

282

in Qiqihaer and the edge area of the study area. Compared with that identified by the LEHP, the spatial

283

distribution of hotspots identified by LEHSI was relatively scattered, and the coverage area was

284

smaller. The spatial distribution of the coldspot area was relatively dispersed, and the coldspots

285

identified by the LEHSI had obvious overlapping areas from 2000 to 2015, and approximately 44.1 %

286

of the identified coldspots were located in the city centre of the study area.

287

For the LEHR, the coverage area of hotspots increased from 11.9 % to 12.6 %, and the coverage

288

area of coldspots increased from 9.3 % to 13.8 %. Hotspots identified by the LEHRI had obvious

289

overlapping areas from 2000 to 2015, and approximately 22.3 % of the identified hotspots were located

290

in Qiqihaer, Longjiang and Nehe. The overall distribution of the hotspots identified by the LEHRI was

291

relatively concentrated. The spatial distribution of the coldspot areas was relatively dispersed, and the

292

coldspots identified by the LEHRI had obvious overlapping areas from 2000 to 2015, and

293

approximately 11.3 % of the identified coldspots were located at the edge area of the study area.

20

294 295

Fig. 9 Spatial heterogeneity identified by the LEHPI, LEHSI and LEHRI in the study area in 2000 and

296

2015

297

3.3 Uncertainty

298

In the process of Monte Carlo simulation, the uncertainty of the weight was gradually transferred

299

to the evaluation results with the calculation process, so there was uncertainty in this evaluation result.

300

This part used points A, B, C and D as examples to illustrate the uncertainty in the evaluation of the

301

LEH, LEHP, LEHS and LEHR.

302

3.3.1 LEH uncertainty

303

For point A, the values of LEH ranged from 0.71 to 0.78 (2000) and 0.74 to 0.82 (2015); the mean

304

values of LEH and the 95 % confidence limits were 0.747±0.016 (2000) and 0.793±0.010 (2015); the

305

mean values of LEH and the 99 % confidence limits were 0.747±0.021 (2000) and 0.793±0.013 (2015);

306

and the adjusted R2 was 0.825 (2000) and 0.925 (2015) (Fig. 10, Fig. 11 and Table 7). 21

307

For point B, the values of LEH ranged from 0.55 to 0.60 (2000) and 0.59 to 0.64 (2015); the mean

308

values of LEH and the 95 % confidence limits were 0.582±0.014 (2000) and 0.615±0.017 (2015); the

309

mean values of LEH and the 99 % confidence limits were 0.582±0.018 (2000) and 0.615±0.022 (2015);

310

and the adjusted R2 was 0.867 (2000) and 0.806 (2015) (Fig. 10, Fig. 11 and Table 7).

311

For point C, the values of LEH ranged from 0.47 to 0.55 (2000) and 0.46 to 0.53 (2015); the mean

312

values of LEH and the 95 % confidence limits were 0.509±0.021 (2000) and 0.501±0.022 (2015); the

313

mean values of LEH and the 99 % confidence limits were 0.509±0.027 (2000) and 0.501±0.029 (2015);

314

and the adjusted R2 was 0.866 (2000) and 0.784 (2015) (Fig. 10, Fig. 11 and Table 7).

315

For point D, the values of LEH ranged from 0.73 to 0.79 (2000) and 0.69 to 0.75 (2015); the mean

316

values of LEH and the 95 % confidence limits were 0.761±0.013 (2000) and 0.721±0.013 (2015); the

317

mean values of LEH and the 99 % confidence limits were 0.761±0.017 (2000) and 0.721±0.017 (2015);

318

and the adjusted R2 was 0.889 (2000) and 0.903 (2015) (Fig. 10, Fig. 11 and Table 7).

22

319 320

Fig. 10 Spatial pattern of the levels of LEH in 2000, with four typical sites (A, B, C and D) to illustrate

321

the embedded uncertainty of the Monte Carlo simulation

23

322 323

Fig. 11 Spatial pattern of the levels of LEH in 2015, with four typical sites (A, B, C and D) to illustrate

324

the embedded uncertainty of the Monte Carlo simulation

325

24

326

Table 7

327

Uncertainty analysis of LEHI (A, B, C, and D) 2000year

2015year

the mean of LEHI the mean of LEHI Adjust R2 and the 95% and the 99% confidence limits confidence limits

the mean of LEHI the mean of LEHI Adjust R2 and the 95% and the 99% confidence limits confidence limits

A

0.747±0.016

0.747±0.021

0.825

0.793±0.010

0.793±0.013

0.952

B

0.582±0.014

0.582±0.018

0.867

0.615±0.017

0.615±0.022

0.806

C

0.509±0.021

0.509±0.027

0.866

0.501±0.022

0.501±0.029

0.784

D

0.761±0.013

0.761±0.017

0.889

0.721±0.013

0.721±0.017

0.903

328 329

3.3.2 LEHP, LEHS and LEHR uncertainty

330

LEH was calculated by the three indexes of LEHP, LEHS and LEHR. The uncertainty of LEH

331

was directly determined by them; thus, it was necessary to understand the uncertainty of LEHP, LEHS

332

and LEHR. This study took point A (LEHI’s high-value) and point C (LEHI’s low-value) as examples

333

(Fig.12 and Table 8).

334

For point A, the mean values of LEHP were 0.280 (2000) (99 % confidence limits: 0.254-0.306;

335

95 % confidence limits: 0.260-0.300) and 0.293 (2015) (99 % confidence limits: 0.266-0.320; 95 %

336

confidence limits: 0.273-0.313). Its adjusted R2 was 0.804 (2000) and 0.759 (2015). The mean values

337

of LEHS were 0.227 (2000) (99 % confidence limits: 0.192-0.262; 95 % confidence limits: 0.201-

338

0.253) and 0.235 (2015) (99 % confidence limits: 0.202-0.268; 95 % confidence limits: 0.210-0.260).

339

Its adjusted R2 was 0.814 (2000) and 0.818 (2015). The mean values of LEHR were 0.243 (2000)

340

(99 % confidence limits: 0.225-0.261; 95 % confidence limits: 0.229-0.257) and 0.268 (2015) (99 %

341

confidence limits: 0.246-0.290; 95 % confidence limits: 0.251-0.285). Its adjusted R2 is 0.833 (2000)

342

and 0.774 (2015). 25

343

For point C, the mean values of LEHP were 0.221 (2000) (99 % confidence limits: 0.198-0.244;

344

95 % confidence limits: 0.204-0.238) and 0.213 (2015) (99 % confidence limits: 0.189-0.237; 95 %

345

confidence limits: 0.195-0.231). Its adjusted R2 was 0.853 (2000) and 0.804 (2015); the mean values of

346

LEHS were 0.097 (2000) (99 % confidence limits: 0.075-0.119; 95 % confidence limits: 0.081-0.113)

347

and 0.108 (2015) (99 % confidence limits: 0.083-0.133; 95 % confidence limits: 0.089-0.127). Its

348

adjusted R2 was 0.888 (2000) and 0.802 (2015); the mean values of LEHR were 0.191 (2000) (99 %

349

confidence limits: 0.176-0.206; 95 % confidence limits: 0.180-0.202) and 0.177 (2015) (99 %

350

confidence limits: 0.165-0.189; 95 % confidence limits: 168-0.186). Its adjusted R2 was 0.833 (2000)

351

and 0.912 (2015).

352 353 354

Fig. 12 The LEHP, LEHS and LEHR in 2000 and 2015, with four typical sites (A and C) to illustrate the embedded uncertainty of the Monte Carlo simulation

355

Table 8

356

Uncertainty analysis of LEHP, LEHS and LEHR (A and C) 26

2000year 95%

Subsystem

2015year 99%

95%

99%

confidence

confidence

limits

limits

Adjust R2

Mean

0.266-0.320

0.759

0.293

0.195-0.231

0.189-0.237

0.804

0.213

0.227

0.210-0.260

0.202-0.268

0.818

0.235

0.888

0.097

0.089-0.127

0.083-0.133

0.802

0.108

0.225-0.261

0.833

0.243

0.251-0.285

0.246-0.290

0.774

0.268

0.176-0.206

0.833

0.191

0.168-0.186

0.165-0.189

0.912

0.117

Adjust R

2

confidence

confidence

Mean

limits

limits

A

0.260-0.300

0.254-0.306

0.804

0.280

0.273-0.313

C

0.204-0.238

0.198-0.244

0.853

0.221

A

0.201-0.253

0.192-0.262

0.814

C

0.081-0.113

0.075-0.119

A

0.229-0.257

C

0.180-0.202

LEHP

LEHS

LEHR

357 358

4 Discussion

359

This study mainly highlights the following five aspects: the rationality of the evaluation results, the

360

scientific nature of the selection model, the main factors affecting the evaluation results, the formulation

361

of the eco-space management programme, and limitations and prospects.

362

4.1 Rationality of the assessment of LEH in the study area

363

For LEH, the increase in the LEHI in study area was small (from 0.554 (2000) to 0.563 (2015)),

364

and its grade did not changed significantly (2000 (III) and 2015 (III)), mainly due to the existence of

365

weak areas (Tai Lai, Longjiang and Gannan) in the study area, where there was land desertification,

366

low investment in environmental protection and other issues. However, the area where LEH was

367

considered well (22.61 %) was larger than the area where LEH was considered weak (11.09 %), which

368

was also the reason for the growth of the LEHI. Hotspot areas (high-value LEH) were mainly

369

distributed in the Zhalong Nature Reserve. This area is rich in wetland resources, with a high LEHS

370

value, a low degree of human disturbance, and a high level of social protection. Coldspot areas

371

(low-value LEH) were mainly distributed in the central area of the city, which had a high economic 27

372

level, a high degree of human interference, and a low LEHS value.

373

For subsystem of the LEH, this study observed that the average of LEHS in the study area

374

changed more than the average of LEHP and LEHR, so LEHS could determine the level of LEH. For

375

example, on one hand, in 2000, the mean values of LEHS were in the order of D (0.246) >A

376

(0.227) >B (0.128) >C (0.097); the mean values of LEH were in the order of D (0.761) >A (0.747) >B

377

(0.582) >C (0.509); their values were in the same order, which further demonstrates that LEHS had a

378

high contribution rate to the level of LEH. On the other hand, the spatial distributions of the LEHP and

379

LEHS were similar to that of LEH, but the spatial distributions of LEHR and LEH were substantially

380

different mainly because the indicators of subsystems of LEHR were different between 2000 and 2015;

381

therefore, the subsystem of LEHR did not form a relatively stable spatial distribution.

382

4.2 Monte Carlo simulation uncertainty to improve the accuracy of evaluating LEH

383

This study selected typical points (A, B, C and D) to illustrate the uncertainty of the assessment of

384

the LEH based on the Monte Carlo simulation. The principle of the Monte Carlo simulation model was

385

that the uncertainty of results could be well quantified by the uncertainty of parameters, and the range

386

of the uncertainty of results and the statistical distribution were obtained; that is, the results with

387

uncertainties were more reasonable and reliable with the Monte Carlo simulation model than with other

388

models (Song et al., 2015).

389

Monte Carlo simulations were more reasonable than other methods of determining weight (Jiang

390

et al., 2013), because they fully considered the uncertainty and provided the scope and probability

391

distribution of the assessment to demonstrate its advantages. Other methods (for example, fuzzy

392

analytic hierarchy process (Koulinas et al., 2019) and membership functions (Lu et al., 2012) also

393

involved uncertainty analysis) were not as comprehensive in terms of uncertainty analysis (Schader et 28

394

al., 2019), while Monte Carlo simulations model established weighted sample data (subjective

395

weighting methods and objective weighting methods), associated the data with spatial information

396

using the Python environment, considered all the uncertainties and provided both the ranges and

397

probability distributions of the assessments, obtaining a fitted normal distribution curve of LEHI,

398

LEHP, LEHS and LEHR. According to the research results, the adjusted R2 of fitted normal

399

distribution curves were greater than 0.75, indicating that the normal distribution curve was well fitted

400

(Shi, 2014).

401

4.3 Mechanism analysis of factors affecting LEH in the study area

402

According to the ordinary least squares (OLS) model and the geographically weighted regression

403

(GWR) model (Section S1 of Appendix A), this study found that the influencing factors affecting the

404

LEH of the study area were ESV (S4), HAD (P2) and TIE (R1) (Fig. 13).

405

The area of HAD (P2) as the main negative driving influence factor was 18375.72 km2 (2000)

406

and 14902.81 km2 (2015), accounting for 43.52 % (2000) and 35.29 % (2015) of the study area. The

407

area with this indicator as the negative driving factor was mainly distributed in the central and

408

northeastern part of the study area. These areas were mainly cultivated land, construction land and

409

unused land, which also leaded to a higher value of HAD (P2). The factor was a negative effect factor,

410

after standardization, its value was low, which was not conducive to the improvement of the level of

411

LEH in the study area.

412

The area of ESV (S4) as the main positive driving influence factor was 10680.55 km2 (2000) and

413

11725.13 km2 (2015), accounting for 43.52 % (2000) and 35.29 % (2015) of the study area. The area

414

with this indicator as the positive driving factor was mainly distributed in the eastern, northeastern and

415

western parts of the study area. These regions are located in areas with relatively slow economic 29

416

development, where the stress index and response index values were relatively low; that is, the ESV

417

(S4) factor (the indicator of state) had a significant impact on them. The higher the ESV (S4) index, the

418

greater the improvement in the LEH was in the study area.

419

The area of TIE (R1) as the main positive driving influence factor was 8259.22 km2 (2000) and

420

8303.36 km2 (2015), accounting for 19.56 % (2000) and 19.66 % (2015) of the study area. The regions

421

with this indicator as the positive driving factor were mainly distributed in the southwestern, central

422

and eastern regions of the study area, and the ecological conditions of this coverage area were poor (for

423

example, in the western part of the study area, land desertification was prominent). The higher the

424

capital investment in environmental management is, the greater the improvements LEH (Fig.13).

425 426

Fig. 13 Spatial distribution of the major influencing factors of LEH in the study area in 2000 and 2015

427

4.4 Management recommendations based on LEH evaluation results

428

The LEH in the study area is of great significance to the development of the main grain-producing

429

areas. Therefore, the level of LEH should be improved by implementing eco-space management in

430

conjunction with cleaner production techniques. According to the results of the assessment of the LEH

431

in the study area, the key areas of LEH management were reasonably determined (Shi &Yang, 2014)

432

(Fig. 14).

433

Region III was the Zhalong Nature Reserve, and the management of the region should focus on 30

434

maintaining the current state of ecosystem health. The region had the best current LEH status in the

435

study area, mainly due to increased government funding for wetland resources. For future management,

436

the region should strictly reduce the planning of human productive activities through legislative means

437

in accordance with the requirements of the overall plan of the nature reserve, and the region should

438

also increase the construction of ecological corridors to avoid the “island effect” of nature reserves.

439

Region II was the area where the LEH of the study area deteriorated and was mainly located west

440

of Tailai and south of Longjiang. The management goal of this area was to reverse the current

441

deterioration of LEH mainly caused by the indicators related to pressure and state; therefore,

442

management measures should also focus on the load on the land ecosystem and establish a clean

443

agricultural production technology ("One clear" and "two reductions"), for example, cleaning field

444

production waste, reducing the amount of agricultural fertilizer input, reducing agricultural non-point

445

source pollution, and improving land resource utilization efficiency. At the same time, management

446

measures for the region should also increase funding for environmental management on land

447

desertification and vegetation ecological restoration projects. Ultimately, an overall strategy of

448

eco-space management from source to process should be formed to improve LEH levels in the region.

449

The health status of Region

deteriorated significantly between 2000 and 2015, mainly due to

450

an increase in the human activity density index and a decline in ecosystem services. This area should

451

avoid implementing human activities on these ecological lands, optimize the spatial layout of land-use,

452

and maintain ecosystem service functions such as production capacity.

31

453 454 455

Fig. 14 Four key areas for LEH management 4.5 Research limitations and prospects

456

In this study, there are some limitations in the assessment of LEH. First, this study only

457

constructed the LEH evaluation system from the pressure-state (VORF)-response, but ignores some

458

equally important factors (filed surveys can be conducted to obtain specific data, such as soil erosion,

459

waste treatment capacity and environmental quality). Therefore, the accuracy of the evaluation

460

framework has certain limitation. Second, an issue with the grid approach was how to more accurately

461

convert socio-economic data within an administrative unit to a spatial grid. Finally, based on the trend

462

of LEH, this paper divided typical area of management and combined the needs of the cleaner

463

production to proposed management measures (only from the aspects of pollution reduction and

464

efficiency improvement), therefore, management measures still needed to be improved.

465

Focusing on the existing limitations, we should actively determine the future research direction.

466

LEH is a complex concept, involves many factors and is the interaction between people and the land

467

system. In the future, this study will further improve the evaluation framework of LEH to better reflect

468

the interaction between natural factors and socioeconomic factors, develop more comprehensive 32

469

management measures based on the perspective of cleaner production and link socio-economic

470

statistics with land-use types (grid units) to ensure the evaluation results closer to the actual situation.

471

5 Conclusion

472

In order to comprehensively and scientifically evaluate the level of LEH, this study constructed a

473

new evaluation framework (P-S(VORF)-R). At the same time, the Monte Carlo simulation

474

(quantitative uncertainty in the process of determining weights) was introduced for the first time, and

475

the level of LEH in the study area was analyzed (temporal-space scale) with 1km*1km as the

476

evaluation unit. The results obtained in this paper were as follows:

477

(1) The status of integrated LEH can be reflected by the LEHI value in the target layer. At a

478

temporal scale, the LEHI values of study area were 0.554 (2000) and 0.563 (2015), and their grades

479

were all ordinary (III). At a spatial scale, the areas with high LEH levels were distributed in the middle

480

of the study area, the western part of Fuyu, the northwest and northeast of Longjiang, and the central

481

part of Nehe. The areas with low LEH levels were distributed in the southwestern part of Tailai and

482

Longjiang.

483

(2) The detailed impact mechanism can be reflected by the subsystems. In the three subsystems,

484

their magnitude of change was LEHS (0.01) > LEHP (-0.002) > LEHR (0.001), and the LEHS

485

subsystem was an important condition for determining the quality of LEH. At the same time, given the

486

results of the GWR model analysis, this study should focus on improving ESV (S4), HAD (P2) and

487

TIE (R1).

488

(3) The Monte Carlo simulation quantified the uncertainty that was passed to the results due to

489

differences in indicator weights. This simulation also provided the range of uncertainty for the LEHI,

490

LEHP, LEHS and LEHR so that the uncertain evaluation results were reliable. 33

491

The research on LEH can help policymakers to determine the contribution of different factors, and

492

based on these contributions, policymakers can draw on the concept of the cleaner production to

493

formulate an eco-space management measures (controlling the use of pollutants from the source and

494

sustainably using land resources) to achieve regional sustainable development；Second, after adjusting

495

the weights and specific coefficients, the established method system of LEH evaluation can also be

496

applied to other areas. Overall, this assessment of LEH is practical, and in the future, we will consider

497

more potential factors and more scientific methods to improve the theoretical system of LEH.

498

Acknowledgements

499

This research was financially supported by the National Natural Science Foundation of China (No.

500

41571165 and No. 41071346). We all grateful to all the data providers ,including the Data Centre for

501

Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn),

502

the Geospatial Data Cloud (http://www.gscloud.cn/), the Chinese Meteorological Science Data Sharing

503

Service, the Harmonised World Soil Database (V1.2) and the Cold and Arid Regions Sciences Data

504

Center at Lanzhou (http://westdc.westgis.ac.cn).

505

Appendix A, Supplementary data

506

Supplementary data associated with this article can be found ,in the online version.

507

References

508

Ameen, R.F. M., Mourshed, M., 2019. Urban sustainability assessment framework development: The

509

ranking and weighting of sustainability indicators using analytic hierarchy process. Sustain. Cities.

510

Soc., 44, 356-366. https://doi.org/10.1016/j.scs. 2018. 10.020

511

Cheng, X., Chen, L.D., Sun, R.H., Kong, P. R., 2018. Land use changes and socio-economic

512

development strongly deteriorate river ecosystem health in one of the largest basins in China. Sci.

513

Total. Environ., 616-617, 376-385. https://doi.org/10.1016/j.scitotenv. 2017. 10.316

514

Chi, Y., Zheng, W., Shi, H.H., Sun, J.K., Fu, Z.Y., 2018. Spatial heterogeneity of estuarine wetland 34

515

ecosystem health influenced by complex natural and anthropogenic factors. Sci. Total. Environ.,

516

634, 1445-1462. https://doi.org/10.1016/j.scitotenv.2018.04.085

517

Ewertowska, A., Pozo, C., Gavalda, J., Jimenez, L., Guillen-Gosalbez, G., 2017. Combined use of life

518

cycle assessment, data envelopment analysis and Monte Carlo simulation for quantifying

519

environmental efficiencies under uncertainty. J. Clean. Prod., 166, 771-783.

520

https://doi.org/10.1016/j.jclepro.2017. 07.215

521

Feng, Z., Liu, X.P., Zhang, J.Q., Wu, R., Ma, Q.Y., Chen, Y.N., 2017. Ecological vulnerability

522

assessment based on multi-sources data and SD model in Yinma River Basin, China. Ecol. Model.,

523

349, 41-50. https://doi.org/10.1016/j.ecolmodel.2017.01.016

524

Flint, N., Rolfe, J, Jones, C.E., Ssellens, C., Johnston, N.D., Ukkola, L., 2017. An ecosystem health

525

Index for a large and variable river basin: methodology, challenges and continuous improvement

526

in Queensland's Fitzroy Basin. Ecol. Indic., 73, 626-636. https://doi.org/10.1016/j.ecolind.2016.

527

10.007

528

He, X., Jiang, G.H., Zhang, R.J., Ma, W.Q., Zhou, T., 2015. Temporal and spatial variation of land

529

ecosystem health based on the pressure-state-response model: A case study of Pinggu District,

530

Beijing. J. Nat. Resour., 30(12), 2057-2068. DOI:10.11849/zrzyxb.2015.12.008 (in Chinese with

531

English abstract).

532

He, J.H., Pan, Z.Z., Liu, D.F., Guo, X.N., 2019. Exploring the regional differences of ecosystem health

533

and its driving factors in China. Sci. Total. Environ., 673, 553-564. https://doi.org/10.1016/ j.scito

534

tenv.2019.03.465

535

Intergovernmental Panel on Climate Change (IPCC). 2006. IPCC Guidelines for national greenhouse

536

gas inventories. In: Eggleston, H.S., Buendía, L., Miwa, K., Ngara, T., Tanabe, K. (Eds.), Prepared

537

by the National Greenhouse Gas Inventories Programme. IGES. Japan.

538

Ishtiaque, A., Myint, S.W., Wang, C.Y., 2016. Examining the ecosystem health and sustainability of the

539

world's largest mangrove forest using multi-temporal MODIS products. Sci. Total. Environ.,

540

569-570,1241-1254. https://doi.org/10.1016/ j.scitotenv. 2016.06.200

541

Koulinas, G.K., Marhavilas, P.K., Demesouka, O.E., Vavatsilkos, A.P., Koulouriotis, D.E., 2019. Risk

542

analysis and assessment in the worksites using the fuzzy-analytical hierarchy process and a

543

quantitative technique- A case study for the Greek construction sector. Safety. Sci., 112, 96-104.

544

https://doi.org/10.1016/j.ssci.2018.10.017 35

545

Li, B.J., Chen, D.X., Wu, S.H., Zhou, S.L., Wang, T., Chen, H., 2016a. Spatio-temporal assessment of

546

urbanization impacts on ecosystem services: A case study of Nanjing City, China. Ecol. Indic., 71,

547

416-427. https://doi.org/10.1016/j.ecolind.2016.07.017

548

Li, G.D., Fang, C.L., Wang, S.J., 2016b. Exploring spatiotemporal changes in ecosystem-service values

549

and hotspots in China, Sci. Total. Environ., 545, 609-620. https://doi.org/10.1016/ j.scitotenv.2015.

550

12.067

551

Li, L.X., Li, Y.N., Ye, F., Li, Z., 2018. Carbon dioxide emissions quotas allocation in the Pearl River

552

Delta region: Evidence from the maximum deviation method. J. Clean. Prod., 177, 207-217.

553

https://doi.org/10.1016/j.jclepro.2017.12.205

554

Lu, C.Y., Gu, W., Dai, A.H., Wei, H.Y., 2012. Assessing habitat suitability based on geographic

555

information system (GIS) and fuzzy: A case study of Schisandra sphenanthera Rehd. et Wils. in

556

Qinling Mountains. China. Ecol. Model., 242, 105–115, https://doi.org/10.1016/j.ecolmodel. 2012.

557

06.002

558

Mariano, D.A., dos Santos, C.A.C., Wardlow, B.D., Anderson, M.C., Schiltmeyer, A.V., Tadesse, T.,

559

Svoboda, M.D., 2018.Use of remote sensing indicators to assess effects of drought and

560

human-induced land degradation on ecosystem health in Northeastern Brazil. Remote. Sens.

561

Environ. 213, 129-143. https://doi.org/10.1016/j.rse.2018.04.048

562

Meng, L.G., Huang, J., Dong, J.H., 2018. Assessment of rural ecosystem health and type classification

563

in Jiangsu province, China. Sci. Total. Environ., 615, 1218–1228. https://doi.org/10.1016/j.

564

scitotenv.2017.09.312

565

Peng, J., Liu, Y.X., Li, T.Y., Wu, J.S., 2017. Regional ecosystem health response to rural land use

566

change: A case study in Lijiang City, China. Ecol. Indic., 72, 399-410. https://doi.org/10.1016/j.

567

ecolind. 2016.08.024

568 569

Pereira, J., Vasquez, O.C., 2017. The single machine weighted mean squared deviation problem. Eur. J. Oper. Res., 261, 515-529. https://doi.org/10.1016/j.ejor.2017.03.001

570

Qian, C.J., Zhang, M.G., Chen, Y.T., Wang, R., 2014. A Quantitative judgement method for safety

571

admittance of facilities in Chemical Industrial Parks based on G1-Variation coefficient method.

572

Procedia Eng., 84, 223-232. https://doi.org/10.1016/j.proeng.2014.10.429

573

Qiu, B.K., Li, H.L., Zhou, M., Zhang, L., 2015. Vulnerability of ecosystem services provisioning to

574

urbanization: A case of China. Ecol. Indic., 57, 505–513. https://doi.org/10.1016/j.ecolind.2015.04. 36

575

025

576

Schader, C., Curran, M., Heidenreich, A., Landert, J., Blockeel, J., Baumgart, L., Ssebunya, B.,

577

Moakes, S., Marton, S., Lazzarini, G., Niggli, U., Stolze, M., 2019. Accounting for uncertainty in

578

multi-criteria sustainability assessments at the farm level: Improving the robustness of the

579

SMART-Farm tool. Ecol. Indic., 106, 105503. https://doi.org/10.1016/j.ecolind.2019.105503

580

Shi, X.Q., Yang, J.X., 2014. A material flow-based approach for diagnosing urban ecosystem health. J.

581

Clean. Prod., 64, 437-446. https://doi.org/10.1016/j.jclepro.2013.09.051

582

Shi, G.G., 2014. Data mining and knowledge discovery for geoscientists, Amsterdam, Netherlands.

583

Reference to a chapter in an edited book: Chapter 2-Probability and Statistics. pp. 22-53.

584

https://doi.org/10.1016/B978-0-12-410437-2.00002-3

585

Song, G.B., Li, Z., Yang, Y.G., Semakula, M., Zhang, S.S., 2015. Assessment of ecological

586

vulnerability and decision-making application for prioritizing roadside ecological restoration: A

587

method combining geographic information system, Delphi survey and Monte Carlo simulation.

588

Ecol. Indic., 52, 57-65. https://doi.org/10.1016/j.ecolind.2014.11.032

589

Song, D.B., Gao, Z.Q., Zhang, H., Xu F.X., Zheng, X.Y., Ai, J.Q., Hu, X.K., Huang, G.P., Zhang, H.B.,

590

2017. GIS-based health assessment of the marine ecosystem in Laizhou Bay, China. Mar. Pollut.

591

Bull., 125(1-2), 242-249. https://doi.org/10. 1016/j.marpolbul.2017.08.027

592

Sun, B.D., Tang, J.C., Yu, D.H., Song, Z.W., Wang, P.G., 2019. Ecosystem health assessment: A PSR

593

analysis combining AHP and FCE methods for Jiaozhou Bay, China. Ocean. Coast. Manage., 168,

594

41-50. https://doi.org/10.1016/j.ocecoaman.2018.10.026

595

Velázquez, J., Gutiérrez, J., García-Abril, A., Hernando, A., Aparicio, M., Sánchez, B., 2019. Structural

596

connectivity as an indicator of species richness and landscape diversity in Castilla y León (Spain).

597

Forest Ecol. Manag., 432, 286-297. https://doi.org/10.1016/ j.foreco.2018.09.035

598

Wang, P., Deng, X.Z., Zhou, H.M., Qi, W., 2018. Responses of urban ecosystem health to precipitation

599

extreme: A case study in Beijing and Tianjin. J. Clean. Prod., 177, 124-133. https://doi.org/10.

600

1016/j.jclepro.2017.12.125

601

Wu, L.Y., You, W.B., Ji, Z.R., Xiao, D.J., 2018. Ecosystem health assessment of Dongshan Island

602

based on its ability to provide ecological services that regulate heavy rainfall. Ecol. Indic., 84,

603

393-403. https://doi.org/10.1016/j.ecolind.2017.09.006

604

Xiao, R., Yu, X.Y., Shi, R.X., Zhang, Z.H., Yu, W.X., Li, Y.S., Chen, G., Gao, J., 2019. Ecosystem 37

605

health monitoring in the Shanghai-Hangzhou Bay Metropolitan Area: A hidden Markov modeling

606

approach. Environ. Int., 133, 105170. https://doi.org/10.1016/j.envint.2019.105170

607

Yan, Y., Zhao, C.L., Wang, C.X., Shan, P., Zhang, Y.J., Wu, G.., 2016. Ecosystem health assessment of

608

the Liao River Basin upstream region based on ecosystem services. Acta. Ecol. Sin., 36(4),

609

294-300. https://doi.org/10.1016/j.chnaes.2016.06.005

610

Zaroni, H., Maciel, L.B., Carvalho, D.B., Pamplona, E.D., 2019. Monte Carlo Simulation approach for

611

economic risk analysis of an emergency energy generation system. Energy, 172,498-508.

612

https://doi.org/10.1016/j.energy.2019.01.145

613

Zhao, J.C., Ji, G.X., Tian, Y., Chen, Y.L., Wang, Z., 2018. Environmental vulnerability assessment for

614

mainland China based on entropy method. Ecol. Indic., 91, 410-422. https://doi.org/10.1016/j.

615

ecolind.2018.04.016

616

Zhao, Y.W., Zhou, L.Q., Dong, B.Q., Dai, C., 2019. Health assessment for urban rivers based on the

617

pressure, state and response framework-A case study of the Shiwuli River. Ecol. Indic., 99, 324-

618

331. https://doi.org/10.1016/j.ecolind.2018.12.023

38

Table 1 The index system of LEH of study area Indicator Ecosystem

Subsystems

No.

Indicators

Units

Calculation method character

P1

P2

Density of population（DP） Human activity density index

Person/km2

-

/

-

%

-

Farmland area/land area

%

-

Construction land area / Land area

kg

-

Statistical data

Population/Land area

（HAD） P3

Rate of land reclamation (LR)

LEHP P4

P5

Proportion of construction land (CP) Farmland use of pesticides and fertilizer (PF)

P6

Economic density (ED)

Yuan/km2

-

land output value of 1km2

S1(V)

NDVI

/

+

Spatial data

S2(O)

Landscape diversity index (LD)

/

+

S3(R)

Ecological resilience index (ER)

/

+

S4(F)

Ecosystem services value (ESV)

yuan

+

yuan

+

Statistical data

%

+

Statistical data

%

+

Statistical data

kg

+

Statistical data

Yuan/person

+

Statistical data

‰

+

Statistical data

Land LEHS ecosystem

E = 0.3* Resil + 0.7 * Resist

health

Total investment in the treatment R1 of environmental pollution (TIE) R2

Forest cover rate (FC) Attainment rate of the industrial

R3 LEHR

waste water discharged (IWW) R4 R5

R6

Grain Production (GP) Per capita net income of rural residents (IRR) Natural population growth rate (NPG)

(Note: For the landscape diversity (LD) index, the analysis of this indicator was performed using the Shannon index (S) (Velázquez et al., 2019); n is the number of land-use types, ai is the i-th land-use type area, pi is the human activity density parameter of the i-th land-use type, vi is the ESV per unit area of the i-th land-use type, and A is the total area of the evaluated unit.)

Table 2 Ecological elasticity coefficient of land-use types in study area Cultivated Land use type

Forest land

Grassland

Construction

Unused

land

land

Water body

land Resilience coefficient

0.3

0.5

0.8

0.7

0.2

1.0

Resistance coefficient

0.6

1.0

0.7

0.8

0.3

0.2

Elasticity coefficient

0.51

0.85

0.73

0.77

0.27

0.44

Table 3 Human activity density parameters of different land-use types Cultivated Land use type

Unused Forest land

Grassland

Construction land

Water body

land Parameters

0.55

land 0.1

0.23

0.95

0.115

0.14

Table 4 Coefficients of the ESV in study area Unit:Yuan(RMB)/hm2 Ecosystem Wetland

water

Cropland

Forestland

Construction

Unused

land

land

Grassland

services FP

285.45

420.24

792.91

261.66

340.95

0

15.86

RM

190.3

277.52

309.24

2362.88

285.45

0

31.72

GA

1910.92

404.38

570.9

3425.38

1189.37

0

47.57

CL

10743.95

1633.4

769.12

3227.15

1236.94

0

103.08

WA

10656.73

14882.95

610.54

3243.01

1205.23

0

55.5

WT

11417.92

11774.73

1102.15

1363.81

1046.64

0

206.16

SFR

1577.89

325.09

1165.58

3187.5

1776.12

0

134.79

BD

2925.84

2719.69

808.77

3576.03

1482.74

0

317.16

RCT

3718.75

3520.53

134.79

1649.26

689.83

0

190.3

Total

43427.75

35958.53

6264

22296.68

9253.27

0

1102.14

（Note: gas regulation (GA), climate regulation (CL), water regulation (WA), soil formation and retention (SFR), waste treatment (WT), recreation cultural and tourism (RCT),

biodiversity (BD), food prod- uction (FP) and raw material (RM)）

Table 5 Grade of LEH in the study area Level Ⅰ Ⅱ Ⅲ Ⅳ Ⅴ

Health status Well Relatively well Ordinary Relatively weak Weak

Comprehensive evaluation value [0.7, 1] [0.6, 0.7) [0.5, 0.6) [0.4, 0.5) [0, 0.4)

Table 6 Overview of the area and proportion of the level of LEH in the study area from 2000 to 2015

Ⅰ

Ⅱ

Ⅲ

Ⅳ

Ⅴ

Year Area 2

(km ）

Porportion (%)

Area 2

(km ）

Porportion (%)

Area 2

(km ）

Porportion (%)

Area 2

(km ）

Porportion (%)

Area 2

(km ）

Porportion (%)

2000

728.62

1.72

10764.09

25.49

21250.05

50.32

9374.66

22.4

106.85

0.25

2015

964.16

2.28

10456.22

24.76

25920.59

61.39

4773.33

11.4

109.97

0.26

Table 7 Uncertainty analysis of LEHI (A, B, C, and D) 2000year

2015year

the mean of LEHI the mean of LEHI Adjust R2 and the 95% and the 99% confidence limits confidence limits

the mean of LEHI the mean of LEHI Adjust R2 and the 95% and the 99% confidence limits confidence limits

A

0.747±0.016

0.747±0.021

0.825

0.793±0.010

0.793±0.013

0.952

B

0.582±0.014

0.582±0.018

0.867

0.615±0.017

0.615±0.022

0.806

C

0.509±0.021

0.509±0.027

0.866

0.501±0.022

0.501±0.029

0.784

D

0.761±0.013

0.761±0.017

0.889

0.721±0.013

0.721±0.017

0.903

Table 8 Uncertainty analysis of LEHP, LEHS and LEHR (A and C) 2000year 95%

Subsystem

2015year 99%

95%

99%

confidence

confidence

limits

limits

Adjust R2

Mean

0.266-0.320

0.759

0.293

0.195-0.231

0.189-0.237

0.804

0.213

0.227

0.210-0.260

0.202-0.268

0.818

0.235

0.888

0.097

0.089-0.127

0.083-0.133

0.802

0.108

0.225-0.261

0.833

0.243

0.251-0.285

0.246-0.290

0.774

0.268

0.176-0.206

0.833

0.191

0.168-0.186

0.165-0.189

0.912

0.117

Adjust R

2

confidence

confidence

Mean

limits

limits

A

0.260-0.300

0.254-0.306

0.804

0.280

0.273-0.313

C

0.204-0.238

0.198-0.244

0.853

0.221

A

0.201-0.253

0.192-0.262

0.814

C

0.081-0.113

0.075-0.119

A

0.229-0.257

C

0.180-0.202

LEHP

LEHS

LEHR

1

Highlights

2

1, Monte Carlo simulation are applied to the evaluation of land ecosystem

3

health;

4

2. Building a land ecosystem health evaluation system based on

5

P-S(VORF)-R;

6

3. Identify the healthy and unhealthy areas of land ecosystem in 2000 and

7

2015.

8

4. The level of land ecosystem health in Zhalong Nature Reserve is

9

relatively well.

10

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: