Impact of Rapid Urbanization on Vulnerability of Land System from Complex Networks View : A Methodological Approach

Rapid urbanization is responsible for the increased vulnerability of land systems and the loss of many crucial ecosystem services. Land systems are typical complex systems comprised of different land use types which interact with each other and respond to external environment processes (such as urbanization), resulting in dynamics in land systems. This work develops a methodology approach by integrating complex networks and disruptive scenarios and applies it to a case study area (Wuhan City in China) to explore the effects of urbanization on land system structural vulnerability. The land system network topologies of Wuhan City during five time periods from 1990 to 2015 are extracted. Our results reveal that (1) the urban land expands at a higher speed than the urban population in Wuhan City; (2) the period of 2005–2010 has witnessed more land area conversions from ecological lands to urban land than other periods; (3) the land system is more vulnerable to intentional attacks on nodes with higher integrated node centrality and larger land area, such as paddy, dryland, and lake; and (4) the network efficiency of the land system would decline sharply if the area shrinkage of paddy, dryland, and lake is larger than 30%, 50%, and 20%, respectively. The results provide some insights into building a resilient urban land system, such as increasing the efficiency of existing urban land and controlling the shrinkage rate of important land use types. This study contributes to existing literature on complex networks by expanding its application in land systems, which highlight the potential of complex networks to capture the complexity, dynamics, heterogeneity, and emergent phenomena in land systems.


Introduction
In recent decades, rapid urbanization has led to frequent conversions among different land uses and covers, especially from forests or wetlands to artificial uses.Such conversions certainly meet short-term human needs but come at the loss of biodiversity and reductions in many ecosystem services [1], resulting in increased vulnerability of land systems [2].For example, the Millennium Ecosystem Assessment reported that over 60% of the global ecosystems were degraded or used unsustainably from the 1950s to 2005 [3].Urbanization is a key contributor.As the largest developing country, China has undergone rapid urbanization since the adoption of the Reform and Opening-up Policy in 1978 and is currently in the middle of the exponential stage of urbanization [4].The overall rate of urban population climbed from 18% to 57% by 2016 [5].Meanwhile, the urban land expanded by 5.5 million hectares from the late 1980s to 2010 [6], which encroached ecological lands with crucial ecological values (e.g., farmlands, forests, and wetlands).The shrinkage in ecological lands imposes accumulative pressure on natural ecosystems and threatens regional ecological security [3,7].Therefore, knowledge of the relationships between urbanization, land use and cover change (LUCC), and land system vulnerability is imperative for sustainable land management to reduce vulnerability and may contribute to preserving the crucial ecosystem services.
Although the term of "land system" was initially defined by Christian and Stewart in 1947 [8], the study on land systems has received limited attention until 2005 when the Global Land Project (GLP) focusing on land systems was launched by global research agendas.The definition of 2 Complexity land systems was also expanded from the initial biophysical systems (consist of topography, soils, and vegetation) [8] to coupled human-natural systems (i.e., mosaics of land use and cover) [9].Integrating experiences from the Global Change and Terrestrial Ecosystem [10] and Land Use Land Cover Change programs [11], the GLP moved towards studying land systems as complex interactions between natural, social, and coupled processes at various temporal and spatial scales [12].Yet the research focus of GLP is primarily on the type, amount, distribution, and pattern of land uses and covers [13].Regarding land system vulnerability, a growing attention is attracted in the past few years when the Future Earth plan (2014-2023) was initiated by the International Council of Scientific Unions (ICSU) in 2014 [14].The new program aims to develop the knowledge for responding effectively to the risks and challenges of global environmental changes.Thus, the focus of scientific programs is shifted from the observation of changes to more integrative vulnerability and resilience analyses of land systems [15].
Since vulnerability and resilience describe the capacity of a system to cope with external perturbations, a systemlevel understanding of the dynamics in land systems is crucial.However, most previous research on land systems has focused on land use or land cover, while the integrality and systematicness of land systems are often neglected [16].Thus, there is a growing need for modeling tools that can unpack the complexity of land systems and dissect the conversion relationships and processes among land system components under perturbations or stressors, which may contribute to revealing the mechanism of land system vulnerability.
The last decade has witnessed a rapidly growing interest in adopting complex networks to analyze systems that are large, complex, and dynamic [17].Its successful applications in many real-world networks, such as social networks, brain networks, and infrastructure networks, have revealed that the complex network analysis is capable of uncovering the organizing principles that govern the formation and evolution of such complex systems [18].In particular, the network analysis converts complex systems into networks with nodes and links, which facilitate analyzing important properties of complex systems at system-level from the perspective of network topology [19].By integrating disruptive scenarios, the complex network analysis can be further extended to predict system vulnerability under a set of scenarios that represent varying extents of disturbance.For example, Ouyang et al. [20] introduced the complex network method in evaluating structural and functional vulnerabilities of interdependent infrastructures.Li and Xiao [21] extracted the topological structure of eco-industrial parks based on complex network methods and analyzed the resilience of the parks under different attack strategies.Regarding the applications of complex network analysis in land systems, existing studies are still in an early phase of extracting network topology [22], identifying key land use types or dominant land use conversions [16,23].To our knowledge, no study to date has used the complex network methods to analyze the resilience or vulnerability of land systems.Thus, we attempt to bridge the gap by integrating theories and approaches of complex network analysis to enhance the understanding of land system vulnerability, which could be potentially very relevant for building a more resilient land system.
To illustrate how complex network analysis works, we conduct a case study in Wuhan City-a megalopolis in central China that is experiencing rapid urbanization.In particular, this paper has three specific objectives to obtain: (1) to describe the topological characteristics of the case study area from 1990 to 2015 based on commonly used measures of complex network approaches; (2) to identify the important land use types and land use conversions in the land system; and (3) to explore the response of land system structural vulnerability to different urbanization scenarios.

Complexity Analysis of Land Systems
Complex systems consist of interconnected components that interact in a network or a geographic space following simple rules and lead to all kinds of interesting dynamics [24,25].By definition, the complexity of complex systems originates from three key sources: heterogeneous components, dynamic interactions, and emergent system functions.Land systems have most of the key sources of complexity described above and thus are widely deemed as typical complex systems [12,26].
First, land systems are huge systems that consist of hierarchical and heterogeneous components.From a macroscopic perspective, a land system is composed of cropland, water, forest, grassland, and urban subsystems.Each of the subsystems is constituted by subsystems of the lower level; for example, the terrestrial water system is made up of river, lake, wetland, reservoir, and aquifer subsystems.From a microscopic point of view, a land system is an assembly of land units (e.g., a site or satellite image pixel) which are geographically related [26].In the paper, we are mainly concerned about different land types, focusing on the ecological lands (cropland, grassland, forest, water, and so on) to analyze the roles of different land types in LUCC processes.
Second, different land types that interact with each other across both spatial and temporal scales form complex conversions between different land types.The subsystems exchange fluxes and flows, such as nutrients, energy, carbon, and information [26].These interactions occur at different spatial locations and take place within the internal environment of land systems.In addition, as land systems are human-influenced, the interactions among subsystems also drive and respond to external environment changes, such as climate, urbanization, and macroeconomy [9].In this paper, we extract the conversion relationships between different land types incurred by rapid urbanization area as edges.It is noted that the conversions in land systems have directed and weight features.
Last, land systems generate emergent phenomena and nonlinear behaviors due to the dynamic interactions among the heterogeneous components.Emergent phenomena are described as the aggregate outcomes that cannot be predicted by inspecting the components of the system in isolation [27].An example of emergent phenomena includes the famous spatial segregation model by Schelling.The model specifies that a spatial setup comprises households from different races that form racially mixed neighborhoods.The households prefer residing with households of similar race and can relocate until the ratio of neighborhoods of similar race is above a satisfactory threshold [28].The model finally exhibits high levels of segregation, illustrating how local interactions can lead to surprising aggregate spatial patterns.Another example is the aggregate distribution of commercial and residential areas, which can be identified as emergent properties of land markets.To conclude, land systems can be denoted as directed weighted complex systems composed of different land types (nodes) and conversions between different land types (edges), which is interconnected in LUCC processes.Complex network theory provides a powerful tool to gauge the structural vulnerability within land systems and contributes to preserving the stable ecological functions in order to build resilient cities in rapid urbanization area.

An Introduction to Complex
Networks.Complex network analysis (both theories and methods), which has its origin in graph theory and statistical physics, is considered as an important tool to study complex evolving systems [29].It is capable of modeling real-world networks that comprise a large number of components interacting with each other in a complicated manner [30] and dynamically evolves in time [31].Specifically, the rapidly growing interest in complex networks was triggered by the emergence of the smallworld model (1998) [32] and the scale-free networks (1999) [29].Since then, complex network analysis has been widely applied to the investigation of real-world networks such as, for example, the World Wide Web [33], electric power grids [34], transportation systems [35], ecological networks [21], genetic regulatory networks [36], and brain networks [17].
The fundamental units of a network include a collection of nodes that identify the components of a system and edges that denote relationships between pairs of components.Edges in complex networks can have directions and weights (e.g., the length, thickness, capacity, load, or strength associated with an edge).In general, networks can be divided into four categories according to whether directions and weights are considered, namely, undirected unweighted, directed unweighted, undirected weighted, and directed weighted networks.
Based on complex properties of land systems, we extract the topology structure of land systems and established directed weighted networks for land systems.Therefore, an introduction of basic concepts and important network measures of complex network methods is given focus on weighted complex networks based on several review papers [31,37,38] and a book on complex networks [24].
Node Degree.The number of edges incident with the node is defined in terms of the adjacency matrix A as follows: In directed networks, it is possible to distinguish two different degrees: the in-degree  in,푖 = ∑ 푛 푗=1  푗푖 , referred to as the incoming edges linked to node , and the out-degree  out,푖 = ∑ 푛 푗=1  푗푖 , referred to as the outgoing edges from node .The total degree and the average total degree of a network are then defined as ( Node Strength.It is a key metric that measures the importance and the connectivity of a node in a network by taking into account all the weights linked to the node.Similarly, for a directed weighted network, the node strength consists of two parts: out-strength and in-strength.The out-strength ( 푤 out,푖 ), in-strength ( 푤 in,푖 ), total node strength ( 푤 푖 ), and average node strength of a network ( 푤 ) can be defined as Node Betweenness.Two nonadjacent nodes (e.g.,  and ℎ) can be linked via a path starting from node , passing through several pairs of nodes and eventually connecting to node ℎ.Thus, nodes with more edges passing through to connect other nodes play more important roles in a network to maintain its connectedness.Here, we introduced the node betweenness ( 푖 ) to identify key nodes in complex systems: where  푗ℎ is the total number of geodesics connecting  and ℎ in the network and  푗ℎ () is the number of these geodesics that pass through node i.Specifically,  ̸ = ℎ ̸ = .example, for electricity network and railway network, the shortest path provides an optimal pathway to achieve a fast transfer of loads and passengers.The shortest path length ( 푖 푗 ) between node  and node  is defined as the number of edges along the shortest path connecting  and .Accordingly, the average shortest path length () is the average value of  푖 푗 for all the possible pairs of nodes in the network: Clustering Coefficient.It measures the extent to which nodes in a network tend to cluster together.The weighted clustering coefficient C푖 푤 is the ratio of the number of edges among the neighbors of node  to the maximum number of possible edges.C푖 푤 measures the centralization degree of a system.
Therefore, the formula for C푖 푤 is defined as where  [1/퐾] = { 푖 푗 1/푘 }; t퐷 푖 is the number of directed triangles actually formed with its weighted counterpart;  푖 퐷 is the maximum number of possible edges among the neighbors of node .As defined before,  푖 tot is the node degree.
Mathematically, node  can be possibly linked to a maximum of  푖 tot ( 푖 tot − 1)/2 pairs of edges.Node  and each pair of edges can form up to two triangles as the edge between them can be oriented in two ways; thus, the maximum number of possible triangles is  푖 tot ( 푖 tot − 1).However, the "false" triangles formed by  and by a pair of directed edges pointing to the same node are included, for example,  →  and  → .Therefore, we estimate  푖 퐷 by subtracting 2 푖 ↔ .The average clustering coefficient for the whole network ( 푤 ) is then defined as follows:

A Complex Network Method to Analyze Land System
Vulnerability.To investigate system vulnerability that arises from land system structural dynamics, we proposed a hybrid method with three steps.First, we extracted the topology structure of the land system and established a set of networks to represent the land system at different time-steps.Then the topological properties of land system networks were studied based on the important measures of complex network analysis as introduced above.In particular, an integrated node centrality indicator was adopted to identify key land use types and conversions.Finally, we evaluated the structural vulnerability of the land system using a network efficiency indicator by integrating the scenario analysis.The process of how we analyzed land system dynamics and vulnerability is displayed in Figure 1.

Establishment of Land System
Networks.The land system network is identified as a directed weighted network, in which land use types are identified as nodes, and edges Complexity 5 represent the transitions between land use types.Edges in land systems have directions; that is,  푖 푗 is a connection going from node  to node , representing the transition from land use type  to type .Specifically,  푖 푗 = 1 if a land system undergoes the conversion from type  to , and  푖 푗 = 0 otherwise.The edges form an adjacency matrix, which exhibits all the transitions that take place during a time-step.
In addition, the weights ( 푖 푗 ), associated with each edge, are measured by the area converted from land use type  to type .If the transition does not exist,  푖 푗 is set to zero.A dynamic complex network with  time-steps can be represented as a sequence of networks  = { 1 ,  2 , . . .,  푇 }.

Identification of Key Land Use Types and Transitions.
The identification of critical components (e.g., nodes and edges) in a network can be assessed by degree, strength, and betweenness [39,40].
Here, we constructed an integrated node centrality (INC 푖 ) indicator [21] by combining total node degree ( tot 푖 ) and node betweenness ( 푖 ) to detect the important nodes in land system networks: where  푖 is the node degree of ;  푖 and  푖 are the normalized values of  tot 푖 and  푖 , respectively;  1 and  2 represent the weight of node degree and node betweenness.In the paper, we assumed that  tot 푖 and  푖 contribute equally to the domain of topological structure.Therefore, both  1 and  2 take the value of 0.5.
The dominant land use transitions in different stages of urbanization process were identified based on total node strength ( 푤 푖 ), in-strength ( 푤 in,푖 ), and out-strength ( 푤 out,푖 ).Specifically, the value of node strength measures all conversions associated with a land use type.The division of the outstrength and in-strength measures whether a land use type has more area gains or losses during a reference period.If the out-strength of a node is larger than the in-strength, then the node is denoted as an output dominated land use type.Otherwise, it is defined as an input dominated one.
The statistic measures used in analyzing network topological characteristics, including the definitions, equations, and their meanings in land systems, are summarized in Table 1.

Analysis of Land System Structural Vulnerability.
The structural vulnerability is the degree to which a system is likely to experience "damage" due to exposure to a disturbance that affects the topological structure of the system [41].In this study, we are particularly interested in how urbanization affects structural vulnerability in land systems; thus, urbanization is considered as a source of disturbances.It is obvious that urbanization leads to area changes of different land use types in land systems.However, the effects of urbanization on structural vulnerability remain unclear.Thus, we evaluated the structural vulnerability of land systems by attacking network structure and then evaluate how much the performance of the network is affected [42].This method involves three steps: (1) define the load capacity each node can take in a system; (2) define the attack schemes.As each node contributes to the connectedness of a complex system, attack scenarios can be defined by partly or completely removing a node from the network; and (3) define the indicator to measure the performance degradation due to the attacks.
Node Load and Load Capacity.In general, nodes in complex systems carry an initial quantity of some type of information or resource, known as load ( 푝,V ) in complex network theory, such as traffic load in a transport network; power load in an electric network; and information load in a virtual network.Node is one type of weight.The load can flow in networks via links (i.e., edges).Additionally, each node is also associated with a maximum load capacity.If the load of each node equals to or is smaller than its maximum load capacity, a system reaches a dynamic balance state and thus can function normally.However, when the network is exposed to an attack, for example, one or more nodes are affected, the initial load balance of nodes in the system will be destroyed, which leads to the load redistribution over the adjacent nodes.If the load reassigned to a node exceeds its load capacity, then the exceeding load will be further passed down to its linked nodes, which may trigger another round of load redistribution and eventually lead to dramatic changes in the performance of the network.
In land system networks, the load ( 푝,V ) associated with each node is represented by land use area.The maximum load capacity ( V ) is defined as the upper limit of land use area for each land use type in land systems.To maintain the diversity of landscape as well as the important functions provided by different land uses, one land use type cannot undergo unlimited area increase.In reality, the upper limit for land use types is often regulated by land use planning, city planning, and natural resources protection plans.In complex networks, the load capacity of each node is positively proportional to its initial load  푝,V (0) and constrained by a tolerance parameter () with its value varying from 10% to 100%.Load capacity can be defined as follows: Attack Scheme.In land systems, we design attack schemes to simulate the effects of different urbanization scenarios on land system structural vulnerability.Two types of attacks are widely considered in complex networks: random or intentional attacks.Due to rapid urbanization, land use transitions in land systems have certain directions, for example, from ecological land use types to built-up land.Here, we only take into account intentional attacks.In general, the loss of ecological lands has a positive relationship with the pace of urbanization; that is, a higher speed of urbanization would lead to a greater loss in ecological lands.Hence, we design a series of attack scenarios by reducing the area of each attacked node by 10% to 100% with 10% increment to simulate the effects of different urbanization speeds.Then the decreased land use area (load) of the attacked node (i.e., land use type) will be assigned to its adjacent nodes (via the outgoing edges from the node).There are many rules to assign the load.The rule we adopted in this study is to assign the area of the Average number of conversions it takes to link one land use type to another in the land system.

Node betweenness 𝐵𝑃 푖
The ratio of number of edges passing through node  to connect two nonadjacent nodes (e.g.,  and ℎ) to the total number of geodesics connecting  and ℎ in the network.An indicator to measure the importance of a node in the network.
The ratio of links passing through land use type  to connect two land use types (e.g.,  and ℎ) that have no direct conversions, to the total number of links connecting  and ℎ in the network.An indicator to measure the importance of a land use type in the land system.

Clustering coefficient 𝐶 푤
In a weighted network, clustering coefficient is the ratio of the number of edges among the neighbors of node  to the maximum number of possible edges.It measures the centralization degree of a system, that is, the degree to which nodes in a network tend to cluster together.
It measures whether the network structure of the land system is more clustered.
Load is the initial quantity of some type of information or resource carried by a node in complex systems, such as traffic load in transport network; power load in electric network; and information load in virtual network.It is one type of weight.
Load associated with each land use type is represented by land area.

Maximum load capacity 𝐴 V
Each node is also associated with a maximum load that the node can handle.
To maintain the diversity of landscape as well as the important functions provided by different land uses, one land use type cannot undergo unlimited area increase.In reality, the upper limit for land use types is often regulated by land use planning, city planning, and natural resources protection plans.

Topological efficiency 𝐸(𝐺)
The topological efficiency is defined as the inverse of geodesic distance; it measures the efficiency of information transmission among nodes in a network.The network has larger efficiency if it has smaller geodesics between pairs of nodes.
Land systems with high topological efficiency will be more resilient to larger area shrinkage in key land use types since there are more conversions among land use types in the network.Thus, the attacked land use types with area shrinkages will be more likely to be compensated by the conversions from other land use types, which lowers the negative effects of intentional attacks.
affected node to its adjacent nodes according to the original area share.We assumed that, after a node is attacked, it has no in-flows from other land use types.The impact of the attack will continue until the updated load for each remaining node equals the load capacity as specified in (10).
Topological Efficiency.An indicator of the performance of networks can be evaluated by topological efficiency, which measures the efficiency of information transmission among components in a network.The topological efficiency is defined as follows: where () denotes the topological efficiency of the network and 0 ≤ () ≤ 1, and  푖 푗 is the geodesic distance between node  and node .From the formula, it is obvious that the network has larger efficiency if it has smaller geodesics between pairs of nodes.Similar to other real-world networks, it is also crucial to maintain a high topological efficiency in the land system network.In general, land systems with high topological efficiency will be more resilient to larger area shrinkage in key land use types since there are more conversions among land use types in the network.Thus, the attacked land use types with area shrinkage will be more likely to be compensated by the conversions from other land use types, which lowers the negative effect of intentional attacks.

Software Platforms.
In this study, several software tools are used to implement the proposed framework as illustrated in Figure 1, including ArcGIS Desktop 10.4,MATLAB R2014a, NetDraw 2.160, and NetLogo 6.0.2.Specifically, the ArcGIS geoprocessing tool "Intersect" was used to identify changes within the land system between two time points of interest (e.g., 2010 and 2015) and generate the adjacency matrix during the period (e.g., from 2010 to 2015).The calculation of important measures for complex networks is processed in MATLAB.The visualization of land system networks is processed in NetDraw (Figures 3 and 8-11).
Additionally, the agent-based modeling (ABM) was applied to simulate structural vulnerability in land systems under different urbanization scenarios in NetLogo platform.ABM is a microlevel simulation approach to model complex systems comprised of "agents" which both interact with each other and environment and can change their attributes as a result of interactions [43].Here, in land system networks, different land use types are identified as agents and land area is the attribute of the agents.Directed links are established between two land use types that have conversions during a reference period.The initial land area for each land use type is entered into the model to initialize the agent's attribute.The scenario settings with varying parameters (attack rate and tolerance parameter) characterize the external environment that affects conversions among land use types in the land system.Thus, the attribute (i.e., land area) of each agent (i.e., land use type) changes over time through the interactions with other agents and environment.The number of conversions and the network efficiency are coded into the model in NetLogo platform as indicators of land system structural vulnerability.The "Behavior Space" tool of NetLogo is used to experiment with the model with parameters of interests being systematically varied.Figure 5 is a series of screenshots of NetLogo interface showing the response of the land system to different scenario settings, including which node is attacked, the values of attack rate, and tolerance parameter.

Study Area and Data.
To illustrate how complex networks can be applied in evaluating land system structural vulnerability, we conduct a case study based on the methodology described above.Wuhan City is a metropolis situated in the central part of China (latitudes 29 ∘ 58 耠 N to 31 ∘ 22 耠 N, longitudes 113 ∘ 41 耠 E to 115 ∘ 05 耠 E) (Figure 2).The region covers ∼8450 km 2 with an average annual temperature of 15.8-17.5 ∘ C and precipitation of 1150-1450 mm.The terrain is dominated by flat plains and nearly 25.6% of the city area is covered by rivers and lakes.In this paper, the urbanization rate was measured by the ratio of urban population to total population.The population data was collected from historical Statistical Bureau of Wuhan City.Although China implemented the economic reform and opening-up policy   in 1978, most resources and opportunities were provided to the eastern coastal area.In the recent decade after the implementation of the strategy of central China rising, the development of Wuhan City has entered a new era with rapid economic growth and urbanization process [44].We selected a relatively long time period (from 1990 to 2015) to capture the historical trends of urbanization in Wuhan City and project its impact on land system vulnerability.During 1990-2000, the urbanization rate increased from 55.9% to 58.9% with an annual growth rate of 0.3%.After entering the new century, the urbanization process of Wuhan City is accelerating and the urban population reached 70.6% by 2015 with an annual growth rate of 0.8% during 2000-2015.The land use data of Wuhan City in 1990City in , 1995City in , 2000City in , 2005City in , 2010, and 2015 was obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC).This dataset adopted the classification system developed by Chinese Academy of Sciences [45], which classified the land system into 25 subtypes.However, six of them, including permanently snow/ice-covered land, sea beach, sandy land, Gobi, saline-alkali land, and other unused land such as desert and tundra do not exist in the land system of Wuhan City, which resulted in a total of 19 land use types for the analysis in this study (Table 2).The trends of urban land area expansion were consistent with that of the urban population, however, with higher annual growth rate.For example, the urban land area grew from 263.2 ha (1990) to 295.8 ha (2000) with an annual growth rate of 1.0%.The trend of urban sprawl accelerated during 2000-2015 as the annual growth rate reached 6.0%.Based on the preliminary analysis and the land use change map (Figure 2), we can see that the landscape in Wuhan City has been largely modified by the rapid urbanization process and the land system structure is also highly dynamic.The complex network analysis can provide more information about the key land change processes that take place within the land system along with the rapid urbanization.

4.2.
Land System Topological Structure.Figure 3 shows the directed weighted land system network for the period of 2010-2015.The visualization of land system networks for other four stages is displayed in Figures 8-11.All nodes are connected in the networks.Since the networks are complex and highly dynamic, we used the complex network measures to further reveal the information hidden in the networks.We first interpret the value of average node degree, which represents the number of conversions among land use types that take place within the land system.As shown in Table 2, the average node degree in the land system of Wuhan City had experienced a downward trend in the first four periods from 1990 to 2010 and then increased during 2010-2015.This indicates the land system tends to be more stable before Reservoir and ponds (43) High cover grass (31) Sparse woods ( 23) Urban ( 51) Woods ( 21 3, the 19 land use types are divided into three major categories, that is, ecological lands, artificial lands, and other lands.Generally, the artificial land types (i.e., urban land, rural residential land, and other built-up land) are input dominated (i.e., outstrength/in-strength < 1), suggesting that these lands have more area gains than losses.Among them, other built-up land (53) has the largest node strength, followed by urban land (51).These two land use types have much lower out-strength/instrength values than rural residential land (52).It is noticeable that the out-strength/in-strength of rural residential land is close to 1.0 in most periods, which suggest that the area of rural residential land is in dynamic balance.The increase in the rural residential land area comes from household splits when children grow up and set up new families.The decrease is mainly contributed by the "increasing versus decreasing balance" land use policy in China [46].Recent decades have witnessed an increasing out-migration of rural labors [47], especially in rapid urbanization areas, which have led to large abandoned rural residential land.Some of the abandoned land was reclaimed for agricultural use to compensate for cropland loss as required by the policy.In contrast, the other two artificial land uses are growing rapidly and have primarily encroached cropland.Both paddy and dryland have the largest node strength values and are output dominated.More importantly, the divisions of the out-strength and in-strength of both paddy and dryland increase dramatically after 2000, indicating the loss of cropland because rapid urbanization has accelerated over time.Although rivers, lakes, and wetlands have experienced fluctuations, they are roughly output dominated land use types with more losses than gains.This is consistent with the observations in previous studies that the area of lakes in Wuhan City has shrunk drastically since the 1980s [48,49].In contrast, reservoirs and ponds belong to the input dominated type in the former four periods but change to the output type during 2010-2015.The transformation may attribute to the decline in the importance of aquaculture industry in Wuhan City.Additionally, the dynamics in the area of river beach is caused by the river level fluctuations, instead of urban expansion.Other land use types, such as woods and grasslands, account for a small share of the total area in Wuhan City.They also present smaller changes and are mostly output dominated.Overall, the main transitions among land use types are dominated by conversions from ecological lands to artificial lands, especially from croplands and water bodies to urban land and other built-up land.The trend of ecological land losses is becoming more evident after 2000, particularly for the period of 2005-2010.

Structural Vulnerability Assessment.
The attacks of more important nodes have larger impacts on network structure than attacking less important ones [50].Taking the land system network in the latest period (i.e., 2010-2015) as an example, we explored how intentional attacks lead to changes in the structure of land system networks.As discussed above, land use changes in rapid urbanization areas are dominated by transitions from ecological land use types to urban land.In general, a faster urbanization speed results in a greater loss in the area of ecological lands.Thus, we used the area declining at varying levels (from 10% to 100%) in ecological lands to simulate the impacts of different paces of urbanization on land system structural vulnerability.In addition, the network efficiency is an indicator of system structural vulnerability.
The simulation was processed in NetLogo with the "Behavior Space" tool with varying parameter settings, that is, attack rate and tolerance parameter vary from 10% to 100% with 10% increment.Figure 5 is an example showing the responses of the land system network when each of the five nodes with the highest INC values is attacked.Specifically, the attack rate (the percentage of area reduction) and tolerance parameter () for this visualization are set at 30%.We can see that before attacking, there are 289 links in total in the land system network, and the nodes are tightly connected with a high network efficiency of 0.93 (Figure 5(a)).However, when the important nodes are attacked, the connectivity and efficiency of the land system network decrease drastically (Figures 5(b)-5(f)).For example, the removal of paddy land area by 30% contributes to a decline in network connectivity with the remaining number of Sparse woods (23) Other woods (24) High coverage grassland (31) Moderate coverage grassland (32) Low coverage grassland (33) River (41) Lake (42) Reservoir and pond (43) River beach (46) Urban land (51) Rural residential land (52) Other built-up land (53) Wetland (64) Bare soil (65) Rocky land (66) links decreasing from 289 to 178 and the network efficiency decreasing from 0.93 to 0.53 (Figure 5(b)).Figure 6 exhibits the changing trends of the network efficiency when the area of each node decreases from 10% to 100%.First, the shrinkage of paddy area has the largest negative impact on the network efficiency of land system networks, which is followed by dryland and lake.In contrast, attacks on other land use types have the smallest impacts on the network efficiency of land system topological structure, even when their land area reduced by 40% or larger rates.This confirms that the land system is more vulnerable to intentional attacks on nodes with higher INC and larger load capacity.Second, the Sparse woods (23) Other woods (24) High coverage grassland (31) Moderate coverage grassland (32) Low coverage grassland (33) River (41) Lake (42) Reservoir and pond (43) River beach (46) Urban land (51) Rural residential land (52) Other built-up land (53) Wetland (64) Bare soil (65) Rocky land (66) shrinkage in paddy area leads to the decrease of network efficiency across different attack rates, whereas attacking dryland or lake brings about fluctuations (both rises and drops) in network efficiency.Upon attack, the decreased load of the attacked node will be redistributed to its adjacent nodes via the outgoing links.Thus, the declines in network efficiency of dryland or lake can attribute to the conversion to nodes with lower connectivity, such as urban land and other built-up land, while the increase trend is contributed by the conversion to nodes with higher connectivity, such as paddy.Finally, given that paddy, dryland, and lake have the largest impacts on the network efficiency and are most likely to experience land area loss due to urban sprawl, we further identified a critical point for each of them based on the result of disruptive scenarios as shown in Figure 6.Here, the critical point refers to an attack level over which an accelerating downward trend in the network efficiency occurs.It is obvious that 30% is a critical point for paddy, because once the critical point is passed, the network efficiency would experience sharp shrinkages and drop to zero when the paddy area declines by 60%.Since the relationship between attack rates and network efficiency of dryland and lake exhibits fluctuations, we select the first point witnessing a sudden loss of network efficiency as the critical point.Accordingly, 50% and 20% are identified as the critical points for dryland and lake, respectively.Hence, land management interventions are needed to protect paddy, dryland, and lake by controlling their area shrinkage lower than their critical point to avoid the accelerated loss of resilience in the land system.Figure 7 illustrates the relationship between network efficiency and the tolerance parameter () when the attack rate is set to 30%.Generally, the network efficiency has a positive association with .It is noteworthy that paddy, dryland, and lake are most sensitive to variations in , indicating that the land system network is more susceptible to intentional attacks on nodes with larger load capacity.Specifically, 50% is a critical point since the increases in  from 0 to 50% lead to remarkable gains in network efficiency.

Conclusions
This work develops a methodology approach by integrating complex networks and disruptive scenarios to study the responses of land systems to different urbanization levels.To illustrate the method, we conduct a case study and choose Wuhan City a case study area, because it is a typical representative of the rapid urbanization areas in China that involves both economic opportunities and environmental challenges.On the one hand, owing to its advantageous geographic  growth engine of China by the "Development Plan for City Clusters along the Middle Reaches" (referred to as the "new plan") [51] designed by the Chinese Central Government in 2015.Thus, Wuhan City has a great development potential and is currently in an accelerating stage of urbanization.On the other hand, the rapid urbanization in the past two decades have led to serious environmental problems, such as the overconsumption of ecological lands, the irreversible damage in the city's water systems (both qualitatively and quantitatively) [52], and the frequent episodes of heavy haze [53].Thus, knowledge of the historical dynamics in the land system and its vulnerability responses to different urbanization scenarios is crucial for developing a resilient city for Wuhan City, as well as for other areas that are undergoing rapid urbanization processes all over the world.
In this study, we first extract the topology structure of the land system in Wuhan City during five time periods, that is, 1990-1995, 1995-2000, 2000-2005, 2005-  results for node strength reveal that this period has the largest losses of paddy and dryland but the largest gains in other built-up land and urban land, indicating that rapid urban expansion during 2005-2010 consumed a large amount of cropland.These trends slow down during the last period (2010-2015) but are still much greater than before 2005.Second, we adopt two indicators, that is, the integrated node centrality and node strength to identify key land use types and major land use conversions.Paddy emerges as the most important node in land system networks as it has the highest number of conversions and largest area involved in these transitions.Time trends for the indicator of out-strength/instrength reveal that paddy, dryland, lake, river, wetland, woods, and low coverage grassland are primary output dominated land use types, whereas urban land and other built-up land are input dominated.This can be confirmed by a previous study by Li et al. [54] that the urbanization in Wuhan City has consumed large amount and many types of ecological lands, which largely affect the structural stability in land systems.Finally, disruptive scenarios are used to explore how land systems structural vulnerability respond to attacking different nodes with different attack rate that Complexity represents varying urbanization levels.Our results suggest that the land system is more vulnerable to intentional attacks on nodes with higher integrated node centrality and larger load capacity (land area), such as paddy, dryland, and lake.The area reductions of an important land use type will lead to the area redistribution among its linked nodes since the total area in the study area remains unchanged.Thus, the attack on important nodes will exert cascading impacts on other nodes in the network, finally affecting the connectivity and network efficiency of the land system.This is similar to the mechanism of structural vulnerability in other complex networks [21,50].In addition, we find that the network efficiency and connectivity of the land system decline sharply if the area shrinkage of an important node exceeds its critical point.Specifically, the critical points for paddy, dryland, and lake are 30%, 50%, and 20%, respectively.Meanwhile, an efficient way to enhance network efficiency of the land system is to increase the area of important nodes by 0-50%.
The results of this study yield three insights to build a resilient urban land system.First, in the land system of Wuhan City, paddy, dryland, and lake are the most vulnerable nodes from a network topology perspective; thus, greater efforts are needed to secure these nodes as they may have the largest cascading impacts on the structural stability of land system networks.Second, it is crucial to control the shrinkage rate of each important land use type by keeping it less than the critical point.This information can be integrated into land use planning and regulations to guide the setting of upper and lower limits for each land use type towards a more sustainable land system.Additionally, it is noteworthy that the urban land expands at a much higher speed than of urban population in Wuhan City; that is, the annual growth rate of urban land is 5 times larger than of urban population during 2000-2015.Thus, a plausible way to meet the land requirement for urban development lies in increasing the efficiency of existing urban land by redeveloping inefficient urban land uses, such as shanty towns, "villages inside city," and discarded factories as suggested by Liu et al. [55].
This study contributes to existing literature on complex networks by expanding its application in land systems, which highlight the potential of complex networks to capture the complexity, dynamics, heterogeneity, and emergent phenomena in land systems.Furthermore, we explore the use of Agent-Based Modeling in NetLogo to simulate how the land system responds to different urbanization scenarios.Specifically, we find that the "Behavior Space" tool in NetLogo is especially powerful and have broad potential applicability to reveal the mechanism behind the vulnerability formation in other complex systems by coupling with complex networks.However, some limitations of this study should be acknowledged.First, we simplify the design of maximum load capacity (the upper limit of land use area) by assuming it is proportional to its initial load and constrained by a tolerance parameter ().However, in reality, this value should be drawn from related plans, such as land use planning, city planning, and natural resources protection plans, which may provide more practical and useful results for policy-makers.Second, since structure affects function [56], the structural durability of the land system is crucial for ensuring its functionality.However, only structural vulnerability is addressed in this study.The quantitative details on how structural vulnerability leads to functional vulnerability need to be established.Our future study will investigate the impacts of urbanization on the functional vulnerability of land systems by integrating ecological indicators.Third, land systems are geographic systems with rich spatial information.However, this study focuses on the topological properties of land systems but neglects the spatial aspect, which will be another direction of our future study.

Figure 2 :
Figure 2: Study area: geographic location, a Landsat image, and land use changes during 1990-2015 of Wuhan City.

Figure 3 :
Figure 3: The topology of land system networks showing the conversions among land cover types in 2010-2015.

Figure 5 :
Figure 5: Land system networks after attacking the nodes and connected links.Note.In this example, tolerance parameter and area shrinkage are set at 0.3.The original network and five nodes with largest INC are shown here.NetLogo interface shows the network efficiency and remaining total links in each network graph after the corresponding node is attacked.

Figure 6 :
Figure 6: Network efficiency versus different shrinkage in area of each land use type ( = 0.3).

Figure 7 :
Figure 7: Network efficiency versus different tolerance parameter (area shrinkage due to attack is set as 30%).

Figure 9 :
Figure 9: The topology of land system networks showing the conversions among land cover types in 1995-2000.

Figure 11 :
Figure 11: The topology of land system networks showing the conversions among land cover types in 2005-2010.

Table 1 :
Mathematical definitions of complex network measures and their representation in land systems.
푖 Number of links connected to a node .Number of land use conversions linked to land use type .In-degree  in,푖 Number of ingoing links to node .Number of land use conversions from other land use types to land use type .Out-degree  out,푖 Number of outgoing links from node .Number of land use conversions from land use type  to other land use types.Shortest path length  푖 푗 Number of edges along the shortest path connecting nodes  and .Minimum sum of conversions that are needed to connect land use type  and type ; reflect the internal structural stability of the land system.A smaller  푖 푗 value indicates the transitions between land use type  and type  are easier.

Table 2 :
Topological indicators of the directed weighted land system network.
Figure 4: The identification of key land use types based on INC (ranked by average INC value).2010as the number of conversions decreases and becomes more active in the last period.Then we move to the average node strength, which further provides the average land area that is converted during the conversions.Interestingly, the period of 2005-2010 has the lowest average node degree but the largest average node strength, indicating that land use conversions during 2005-2010 are more concentrated than in other periods.In other words, more land area was converted via fewer pairs of land use types.Regarding average shortest path, a larger value indicates the network has a lower level of connectivity; that is, more links are needed to join two nodes that are not adjacent.We find that nodes in the land system are less connected during 2000-2005 and 2005-2010.In addition, the network efficiency in the period 2005-2010 is lower, whereas the clustering coefficient is much larger than other periods.All these indicators reveal that the land system network in the period 2005-2010 is more clustered with fewer numbers of conversions among land use types but larger area being converted.Thus, it is important to figure out what are the key nodes that have most in-flows and out-flows and what are the most important transitions that dominate the dynamics of the land system of Wuhan City.
(41)their INC values decline dramatically in later periods.The reason is that the areas of these two land use types are relatively small and have shrunk largely due to urban expansion before 2005.The remaining high coverage grass or sparse woods are mostly distributed at fringes of water bodies or in hilly areas, which lower the likelihoods to be converted into built-up uses.In contrast, urban land (51) and river(41)are becoming more important in the land system in later periods, since their INC values increase significantly from 1990 to 2010 and peak in the period of 2005 to 2010.This result further explains why the land system network is highly clustered during 2005 to 2010, as there is a drastic growth in the number of conversions associated with urban land and river in this period.4.4.The Dominant Land Use Transitions.After identifying the key land use types, we further investigate what are the transitions associated with these key nodes that dominate the dynamics in land systems.As shown in Table

Table 3 :
Node strength † and the division of the out-strength and in-strength in land system networks in different periods.
Note.†The unit of node strength is km 2 ; * rocky land has no incoming weights (no other land transfer to rocky land during 2005-2010).
location, Wuhan City is considered as a new economic The topology of land system networks showing the conversions among land cover types in 1990-1995.
2010, and 2010-2015.We find that the land system of Wuhan City in the period of 2005-2010 is most noteworthy, as the network structure is highly clustered with most of the land use transitions taking place between less land use types.The The topology of land system networks showing the conversions among land cover types in 2000-2005.