Spatiotemporal Evolution of Urban Agglomeration and Its Impact on Landscape Patterns in the Pearl River Delta, China

: An urban agglomeration is the engine of regional and national economic growth, but also causes many ecological and environmental issues that emerge from massive land changes. In this study, the spatiotemporal evolution of an urban agglomeration was quantiﬁed and its impacts on the urban and regional landscape patterns were evaluated. It showed that the urbanized land area of the Pearl River Delta Urban Agglomeration (PRDUA) in China nearly quadrupled, having linearly increased from 1819.8 km 2 to 7092.2 km 2 between 1985 and 2015. The average annual growth rate presented a bimodal wave-like pattern through time, indicating that the PRDUA has witnessed two rounds of the urbanization process. The growth modes (e.g., leapfrog, edge-expansion, inﬁlling) were detected and they exhibited co-existing but alternating dominating patterns during urbanization, demonstrating that the spatiotemporal evolution of the urban development of the PRDUA follows the “spiral diffusion-coalescence” hypothesis. The morphology of the PRDUA presented an alternating dispersal-compact pattern over time. The city-level and regional-level landscape patterns changed synchronously with the spatiotemporal evolution of the PRDUA over time. The urbanization of the PRDUA increased both the complexity and aggregation of the landscape, but also resulted in an increasing fragmentation and decreasing connectivity of the natural landscape in the Pearl River Delta region. These ﬁndings are helpful for better understanding how urban agglomerations evolve and in providing insights for regional urban planning and sustainable land management.


Introduction
Cities concentrate a large amount of labor, information, and financial resources [1,2]. They contributed nearly 85% of the global gross domestic product (GDP) in 2015 and are thus the engines of economic growth and social change [3]. In 2018, 55% of the world's population inhabited urban areas, and this is predicted to rise to 68% by 2050 [4]. As the world's most populous country, China has experienced accelerating urbanization in the past few decades, leading to large built-up areas replacing the natural landscapes [5]. Such urbanization has brought tremendous stress onto the environment and natural resources [6][7][8], and has affected ecosystem functions and services from the local to global scale [9]. Therefore, patterns of the urban landscape of seven cities in the PRDUA. The landscape metrics exhibited several distinct temporal trajectories, including a monotonous increase, and single-hump shaped and double-hump shaped patterns for different cities, while urban expansion had led to the fragmentation of green spaces, farmland, and grassland in the BTHUA, China [57]. However, how the spatial and temporal evolution of the entire urban agglomeration influences the urban and regional landscape patterns remains unclear.
The PRDUA is the fastest-growing metropolitan region in China [25] and has been undergoing a rapid urbanization process since China's reform and opening-up and Deng Xiaoping's "coastal development strategy" in the late 1970s [58], which transformed large farmlands and forests into an impervious surface [59]. Thus, it is an ideal area to investigate how a Chinese urban agglomeration has evolved over the past few decades. This study aims to investigate the spatiotemporal dynamics of urban growth in the PRDUA and their impacts on the urban and regional landscapes. We specifically address the following questions: (1) What are the evolution patterns of the PRDUA over space and through time? (2) Does the PRDUA follow the spiraling diffusion-coalescence hypothesis? (3) What are the impacts of the urban agglomeration evolution on the urban and regional landscape patterns of the PRDUA?

Study Area
The PRDUA is located in the lower reaches of the Pearl River in Guangdong Province, China, and borders the Hong Kong and Macao Special Administrative Regions. According to the Outline of the Reform and Development Plan of the Pearl River Delta Region (2008-2020) [60], the PRDUA encompasses nine independent but highly interconnected cities, i.e., Guangzhou, Shenzhen, Zhuhai, Foshan, Jiangmen, Dongguan, Zhongshan, Huizhou, and Zhaoqing ( Figure 1). The PRDUA covers an area of approximately 55,300 km 2 and has a typical subtropical monsoon climate. The GDP of the PRDUA had reached CNY 6.25 trillion (nearly USD 0.983 trillion) by the end of 2015, more than seven times that in 2000 (CNY 847.13 billion). The resident population of the PRDUA increased from 42.89 million in 2000 to 58.74 million in 2015 [61], resulting in rapid urbanization in the PRDUA. At present, the PRDUA, as one of China's three major world-class urban agglomerations (the other two are the BTHUA and the YRDUA), is a pioneer region, implementing China's reform and opening-up strategies, and is an important economic hub of China [60]. As the main body of the emerging Guangdong-Hong Kong-Macao Great Bay Area, the PRDUA plays a remarkable leading role in promoting regional cooperation and integration, as well as in driving national economic development [59]. Therefore, analyzing the spatiotemporal evolution of the PRDUA is beneficial to guide sustainable urban planning and land-use management of emerging urban agglomerations.

Data Sources and Re-Processing
Urban growth is defined here as the expansion of urbanized land, i.e., the conversion of non-urban areas to urban areas [62]. Urbanized land refers to the impervious surface within urban administrative boundaries, mainly including residential, industrial, transportation, public facilities, and other land-use types [18,63]. The proportion of urbanized land has been widely used as another kind of index for measuring the urbanization degree in a city or a country [30,[64][65][66]. In this study, the yearly urbanized land data covering the entire PRDUA between 1985 and 2015 were derived from the global annual urban dynamics (GAUD) dataset with a 30 m spatial resolution, which was developed by Liu et al. [67]. The original dataset includes yearly urbanized land, non-urbanized land, and green recovery land, and has an overall accuracy of 76% (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and 82% (2000-2015) for the humid regions (accounting for over 90% of global urban lands) [67]. The dataset has been widely utilized to study urban growth [68][69][70][71] and can be freely downloaded [67]. To investigate the spatiotemporal evolution of the landscape of the PRDUA, the original dataset was processed using the following procedures.

Data Sources and Re-Processing
Urban growth is defined here as the expansion of urbanized land, i.e., the conversion of non-urban areas to urban areas [62]. Urbanized land refers to the impervious surface within urban administrative boundaries, mainly including residential, industrial, transportation, public facilities, and other land-use types [18,63]. The proportion of urbanized land has been widely used as another kind of index for measuring the urbanization degree in a city or a country [30,[64][65][66]. In this study, the yearly urbanized land data covering the entire PRDUA between 1985 and 2015 were derived from the global annual urban dynamics (GAUD) dataset with a 30 m spatial resolution, which was developed by Liu et al. [67]. The original dataset includes yearly urbanized land, non-urbanized land, and green recovery land, and has an overall accuracy of 76% (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and 82% (2000-2015) for the humid regions (accounting for over 90% of global urban lands) [67]. The dataset has been widely utilized to study urban growth [68][69][70][71] and can be freely downloaded [67]. To investigate the spatiotemporal evolution of the landscape of the PRDUA, the original dataset was processed using the following procedures.
Firstly, the time-series data of urban and non-urban land use during 1985-2015 was divided into six equidistant time periods, including 2015-2010, 2010-2005, 2005-2000, 2000-1995, 1995-1990, and 1990-1985, to match the China national "5-Year Plan" [21]. The land-use maps of urban and non-urban areas were then extracted from the original raster images in the years of 1985,1990,1995,2000,2005,2010, and 2015, respectively. Secondly, we used administrative boundaries to clip the land-use map to generate land-use types of urban and non-urban for the study area. The administrative boundaries and divisions of the cities in the PRD region were downloaded from the National Geographic Information Resources Catalogue Service System of the Ministry of Natural Resources of China [72], and the Ministry of Civil Affairs of the PRC [73], respectively. Thirdly, the original landuse map contains numerous small patches, some of them may be noise pixels and must be excluded. We used the method proposed by Wu et al. [21] and adopted the 3 × 3 majority filter in ArcMap 10.7 to reprocess the maps of 2015, 2010, 2005, 2000, 1995, 1990, and 1985 to exclude the noise pixels. These seven reproduced maps of the PRDUA were used for calculating various landscape metrics and urban growth rates. The Figure 2 showed the steps of the data-processing in this study. Firstly, the time-series data of urban and non-urban land use during 1985-2015 was divided into six equidistant time periods, including 2015-2010, 2010-2005, 2005-2000, 2000-1995, 1995-1990, and 1990-1985, to match the China national "5-Year Plan" [21]. The land-use maps of urban and non-urban areas were then extracted from the original raster images in the years of 1985,1990,1995,2000,2005,2010, and 2015, respectively. Secondly, we used administrative boundaries to clip the land-use map to generate land-use types of urban and non-urban for the study area. The administrative boundaries and divisions of the cities in the PRD region were downloaded from the National Geographic Information Resources Catalogue Service System of the Ministry of Natural Resources of China [72], and the Ministry of Civil Affairs of the PRC [73], respectively. Thirdly, the original land-use map contains numerous small patches, some of them may be noise pixels and must be excluded. We used the method proposed by Wu et al. [21] and adopted the 3 × 3 majority filter in ArcMap 10.7 to reprocess the maps of 2015, 2010, 2005, 2000, 1995, 1990, and 1985 to exclude the noise pixels. These seven reproduced maps of the PRDUA were used for calculating various landscape metrics and urban growth rates. The Figure 2 showed the steps of the data-processing in this study.

Quantifying Urban Growth Rate of the PRDUA
Equation (1) was used to calculate the annual growth rate of urbanized land (AGRUL) and the urban growth rate over the past 30 years [16]: where UL t+n and UL t are the areas of urbanized land in year t + n and t, respectively.

Measuring Growth Modes of the PRDUA
Urban growth generally exhibits three major types of growth modes, i.e., leapfrog, edge-expansion, and infilling [35,36,38]. Leapfrog refers to newly increased urbanized land occurring in areas beyond the existing urbanized land. Edge-expansion means that the newly grown urbanized land extends at the edge of the existing urban land. Infilling denotes that the new urban patches emerge in areas surrounding the existing urbanized land [16,36]. The Landscape Expansion Index (LEI) was used to model the urban growth modes in the PRDUA, and was calculated using Equation (2) [16,36]: where A 0 is the intersection between a predefined buffer zone of the new urban patch and an existing urban patch. A v is the intersection between the existing non-urban area and the buffer zone of the newly grown urban patch [16,21]. We set 1 m as a buffer distance according to the previous study [36]. The LEI value ranges between 0 and 100. When a patch has an LEI value of 0, it belongs to the leapfrog mode; with an LEI value between 0 to 50, it belongs to the edge-expansion mode; otherwise, it belongs to the infilling mode when the LEI is within 50 and 100 [16,36].

Quantifying Urban Growth Rate of the PRDUA
Equation (1) was used to calculate the annual growth rate of urbanized land (AGRUL) and the urban growth rate over the past 30 years [16]: where ULt+n and ULt are the areas of urbanized land in year t + n and t, respectively. The mean expansion index (MEI), reflecting the aggregate properties of the patch mosaic, is the integration of the LEI of all the patches over the full extent of the data and is calculated using Equation (3) [36]:

of 19
To describe the characteristic of morphological variation, the area-weighted mean expansion index (AWMEI) was also computed using the following Equation (4) [16]: where LEI i refers to the value of the LEI for a newly grown urban patch i, and a i is the area of the patch. N is the count of the newly grown urban patches. A refers to the total area of all the newly grown urban patches. A larger MEI or AWMEI value signals a more substantial compacting trend in addition to the landscape expansion [16,36].

Quantifying the Impact of Urban Agglomeration Evolution on Landscape Pattern Changes
To further quantify the spatial and temporal changes in urban morphology and the impact of urban agglomeration evolution on the city and regional landscape pattern, the landscape metrics were calculated at the class and landscape levels for the whole PRDUA using the Fragstats 4.2 software. Eight frequently used landscape metrics were selected and calculated, including Edge Density (ED), Largest patch index (LPI), Area-weighted Mean Patch Fractal Dimension (FRAC_AM), Area-weighted Mean Contiguity Index Distribution (CONTIG_AM), Contagion index (CONTAG), Aggregation Index (AI), Patch density (PD), and Landscape shape index (LSI). These landscape metrics can accurately capture and quantify urban morphological variations in shape complexity, density, dominance, and aggregation, and can be classified into three types, as shown in Table 1 [50,74]. We also profiled the landscape metrics over time and fitted the change trend of each metric to detect the landscape pattern dynamic along with the evolution of the PRDUA landscape. Describes the complexity of landscape shape by the relationships between patch perimeter and area [52].

Class, Landscape None
Edge density (ED) Describes the length of edges between patches scaled to the area of the landscape [52].

Class, Landscape Meter per ha
Landscape shape index (LSI) Describes the deviation of patch structure from regular shape (square) [52].

Class, Landscape None
Aggregation Contagion index (CONTAG) Describes the juxtaposition and dispersion of landscape elements [52].

Urban Growth and Spatial Expansion of the Urban Agglomeration in the PRDUA
The urbanized land in the PRDUA increased from 1819.8 km 2 to 7092.2 km 2 between 1985-2015, nearly quadrupling over the 30-year period. The top three cities with the highest increase in urbanized land were Guangzhou, Foshan, and Dongguan, with areas of net increase up to 1036.7 km 2 , 985.2 km 2 , and 958.1 km 2 , respectively. Meanwhile, the top three cities with the fastest increase rates are Dongguan, Huizhou, and Zhongshan, whose increase rates were 5.4, 5.3, and 4.8 times, respectively. The urbanized land area increased linearly at the urban agglomeration scale ( Figure 3a) and either linearly or non-linearly at the city scale ( Figure 3b).

Aggregation
Contagion index (CON-TAG) Describes the juxtaposition and dispersion of landscape elements [52].
Class, Landscape none Aggregation index (AI) Describes aggregation levels of spatial patterns [76].

Urban Growth and Spatial Expansion of the Urban Agglomeration in the PRDUA
The urbanized land in the PRDUA increased from 1819.8 km 2 to 7092.2 km 2 between 1985-2015, nearly quadrupling over the 30-year period. The top three cities with the highest increase in urbanized land were Guangzhou, Foshan, and Dongguan, with areas of net increase up to 1036.7 km 2 , 985.2 km 2 , and 958.1 km 2 , respectively. Meanwhile, the top three cities with the fastest increase rates are Dongguan, Huizhou, and Zhongshan, whose increase rates were 5.4, 5.3, and 4.8 times, respectively. The urbanized land area increased linearly at the urban agglomeration scale (Figure 3a) and either linearly or non-linearly at the city scale (Figure 3b).  The spatial expansion of the urbanized land exhibited a continuous trend between 1985 and 2015, and the urban expansion was mainly concentrated in the central and southeastern parts of the PRDUA, where Guangzhou, Foshan, Dongguan, and Shenzhen are located (Figure 4). In the earlier period of the urbanization process, between 1985 and 2000, the urban growth was characterized by a rapid expansion, with a large number of scattered newly urbanized patches and the emerging spatial framework of the urban agglomeration. While in the later period of the urbanization process, between 2001 and 2015, the urban growth exhibited a pattern of spatial coalitions, which was characterized by the edge-expansion and infilling modes from the old urban areas and newly grown patches, eventually forming the current morphology of the urban agglomeration. located ( Figure 4). In the earlier period of the urbanization process, between 1985 and 2000, the urban growth was characterized by a rapid expansion, with a large number of scattered newly urbanized patches and the emerging spatial framework of the urban agglomeration. While in the later period of the urbanization process, between 2001 and 2015, the urban growth exhibited a pattern of spatial coalitions, which was characterized by the edge-expansion and infilling modes from the old urban areas and newly grown patches, eventually forming the current morphology of the urban agglomeration. The average annual growth rate of urbanized land (AGRUL) in the PRDUA was 4.6% between 1985 and 2015. The AGRUL exhibited a bimodal change pattern, indicating that the PRDUA experienced two rounds of urbanization processes. The first round started in 1985 and covered three periods. During the earlier period of 1985-1990, the AGRUL of the whole PRDUA was 4.0%, and in 1990-1995, it quickly increased to the first peak of 10.4%, which is also the maximum value throughout the whole study period. The value then sharply declined to the minimum value of 2.5% in the period 1995-2000. The second round exhibited a relatively slower growth rate and came to the second peak value of 5.  The average annual growth rate of urbanized land (AGRUL) in the PRDUA was 4.6% between 1985 and 2015. The AGRUL exhibited a bimodal change pattern, indicating that the PRDUA experienced two rounds of urbanization processes. The first round started in 1985 and covered three periods. During the earlier period of 1985-1990, the AGRUL of the whole PRDUA was 4.0%, and in 1990-1995, it quickly increased to the first peak of 10.4%, which is also the maximum value throughout the whole study period. The value then sharply declined to the minimum value of 2.5% in the period 1995-2000. The second round exhibited a relatively slower growth rate and came to the second peak value of 5.

Spatiotemporal Evolution of the Urban Agglomeration in the PRDUA
The spatiotemporal evolutions of the PRDUA were well quantified by the spatial and temporal distribution of three urban growth modes over time, respectively. During the first round of urbanization, in the earlier period of 1985-1990, the spatial evolution of the

Spatiotemporal Evolution of the Urban Agglomeration in the PRDUA
The spatiotemporal evolutions of the PRDUA were well quantified by the spatial and temporal distribution of three urban growth modes over time, respectively. During the first round of urbanization, in the earlier period of 1985-1990, the spatial evolution of the PRDUA was characterized by the co-occurrence of numerous newly grown small urban patches in the leapfrog and edge-expansion growth modes (Figure 6a); the edgeexpansion and infilling modes dominated in the middle period of 1990-1995 (Figure 6b The temporal trajectories of the urban agglomeration evolution in the PRDUA during the past three decades can be explicitly depicted by the area proportion and the number of newly increased urban patches in the growth modes of infilling, edge-expansion, and leapfrog ( Figure 7). Firstly, the three growth modes coexisted in all periods, but presented a pattern of alternating relative dominance, which was indicated by their proportion of area and patch number during the six time periods. Secondly, each urban growth mode  The temporal trajectories of the urban agglomeration evolution in the PRDUA during the past three decades can be explicitly depicted by the area proportion and the number of newly increased urban patches in the growth modes of infilling, edge-expansion, and leapfrog ( Figure 7). Firstly, the three growth modes coexisted in all periods, but presented a pattern of alternating relative dominance, which was indicated by their proportion of area and patch number during the six time periods. Secondly, each urban growth mode The evolution of the PRDUA was also reflected in the other two indices of landscape expansion-the MEI and the AWMEI-both of which presented a wave-like pattern over time. They both started from a small value in the early stage of urbanization and increased gradually to a peak value, and then fell back to a small value and started a new round with a similar increase pattern. This meant that the urban morphological structure presented a dispersal form, indicated by small MEI and AWMEI values in the early urbanization stage, and gradually came to a compact form, indicated by a large and peak value in the late urbanization stage. This followed a process of diffusion and coalescence in the first round of urbanization, and then came to the second round and followed a similar diffusion-coalescence pattern (Figure 8). The evolution of the PRDUA was also reflected in the other two indices of landscape expansion-the MEI and the AWMEI-both of which presented a wave-like pattern over time. They both started from a small value in the early stage of urbanization and increased gradually to a peak value, and then fell back to a small value and started a new round with a similar increase pattern. This meant that the urban morphological structure presented a dispersal form, indicated by small MEI and AWMEI values in the early urbanization stage, and gradually came to a compact form, indicated by a large and peak value in the late urbanization stage. This followed a process of diffusion and coalescence in the first round of urbanization, and then came to the second round and followed a similar diffusion-coalescence pattern (Figure 8).

Landscape Pattern Changes with the Evolution of the Urban Agglomeration in the PRDUA
The variations in the different landscape metrics effectively characterize the landscape pattern changes with urban agglomeration evolution over time in the PRDUA. The class-level landscape metrics reflected the landscape pattern of the urban agglomeration and presented three distinct characteristics. Firstly, the metrics that reflected the shape complexity of the cities, such as the edge density (Figure 9a), the area-weighted mean fractal dimension (Figure 9c), the patch density (Figure 9e), and the landscape shape index (Figure 9g), exhibited an M-shape change pattern, which nearly remained in sync with the wave-like patterns of the urban growth rate (Figure 5), the growth mode (Figure 7), and the landscape expansion indices, the MEI and AWMEI (Figure 8). Secondly, the largest patch index (Figure 9b), the aggregation index of the contiguity (Figure 9d), and the aggregation index ( Figure 9f) present a stepwise increase over time. Thirdly, all the landscape metrics generally displayed an increasing trend over time, with the exceptions of the patch density and landscape shape index, which presented a decreasing trend over time.

Landscape Pattern Changes with the Evolution of the Urban Agglomeration in the PRDUA
The variations in the different landscape metrics effectively characterize the landscape pattern changes with urban agglomeration evolution over time in the PRDUA. The class-level landscape metrics reflected the landscape pattern of the urban agglomeration and presented three distinct characteristics. Firstly, the metrics that reflected the shape complexity of the cities, such as the edge density (Figure 9a), the area-weighted mean fractal dimension (Figure 9c), the patch density (Figure 9e), and the landscape shape index (Figure 9g), exhibited an M-shape change pattern, which nearly remained in sync with the wave-like patterns of the urban growth rate (Figure 5), the growth mode (Figure 7), and the landscape expansion indices, the MEI and AWMEI (Figure 8). Secondly, the largest patch index (Figure 9b), the aggregation index of the contiguity (Figure 9d), and the aggregation index (Figure 9f) present a stepwise increase over time. Thirdly, all the landscape metrics generally displayed an increasing trend over time, with the exceptions of the patch density and landscape shape index, which presented a decreasing trend over time. The spatiotemporal evolution of the urban agglomeration also influenced the regional landscape pattern, which can be quantified by the changes in the landscape-level metrics over time. Firstly, the shape complexity indices of the edge density (Figure 10a), areaweighted mean fractal dimension (Figure 10c), landscape shape index (Figure 10g), along with the dominance index of patch density (Figure 10e), all presented an M-shaped pattern over time and remained synchronous with the wave-like pattern of urban growth, which was very similar to that of the class-level metrics. Secondly, the temporal pattern of the area-weighted mean contiguity index and the aggregation index showed a W-shaped pattern (Figure 10f,h), while the largest patch index (Figure 10b) and contagion index ( Figure 10d) demonstrated a monotonic decline over time. Thirdly, the landscape became more complex, while the landscape dominance and aggregation presented an opposite trend, i.e., a declining trend over time.  The spatiotemporal evolution of the urban agglomeration also influenced the regional landscape pattern, which can be quantified by the changes in the landscape-level metrics over time. Firstly, the shape complexity indices of the edge density (Figure 10a), area-weighted mean fractal dimension (Figure 10c), landscape shape index (Figure 10g), along with the dominance index of patch density (Figure 10e), all presented an M-shaped pattern over time and remained synchronous with the wave-like pattern of urban growth, which was very similar to that of the class-level metrics. Secondly, the temporal pattern of the area-weighted mean contiguity index and the aggregation index showed a W-shaped pattern (Figure 10f,h), while the largest patch index (Figure 10b) and contagion index (Figure 10d) demonstrated a monotonic decline over time. Thirdly, the landscape became more complex, while the landscape dominance and aggregation presented an opposite trend, i.e., a declining trend over time.

Quantifying Spatiotemporal Evolution of Urban Agglomerations
An urban agglomeration is the outcome of urban growth and spatial aggregation in cities within a region over time. The evolution of an urban agglomeration is normally driven by multiple forces and follows a spatiotemporal pathway from a city cluster to a

Quantifying Spatiotemporal Evolution of Urban Agglomerations
An urban agglomeration is the outcome of urban growth and spatial aggregation in cities within a region over time. The evolution of an urban agglomeration is normally driven by multiple forces and follows a spatiotemporal pathway from a city cluster to a megalopolis [2]. However, quantifying such a complex spatiotemporal process is a challenge. Our study established a framework consisting of a series of metrics of growth rate, growth modes, and landscape pattern metrics to quantify the evolution of the urban agglomeration in the PRD region over a three-decade period.
Firstly, our temporal change detection, using the AGRUL and growth curve fitting of the urbanized area, explicitly reveals the wave-like urban growth pattern and the two urbanization processes of the PRDUA between 1985-2015 ( Figure 5). This method has been applied to detect the urban growth pattern changes in four major cities in the central YRD region [16], as well as individual cities such as Shanghai, Nanjing, Hangzhou, etc. [21,77]. Our method is a more straightforward way to demonstrate the wave-like urban growth pattern than that of previous studies, which utilized landscape metrics along the urbanization gradient to test the wave-like urban growth pattern [42,43]. Our results showed diverse fitting patterns of the urbanized land in the PRD region; the entire urban agglomeration exhibited a linear growth trajectory, while the nine major cities under the PRDUA presented linear, exponential, logarithmic, and logistic growth trajectories over time (Figure 4), which have never been reported for urban agglomerations in other regions. Previous studies showed that three American metropolises outperformed three Chinese urban agglomerations in terms of the degree of compact development or the quality of urban expansion measured by population density [29]. Due to the following reasons, the diverse growth trajectories revealed in this study imply that the Chinese government should adopt different approaches to better manage land-use efficiency and improve the urbanization quality at the regional scale. The AGRUL of the entire PRDUA reached 4.6% in the past 30 years, which is nevertheless lower than that of 7.4% in the YRDUA between 2000-2015 [78].
Secondly, the spatial evolution of an urban agglomeration over time can be quantified and the exact spatial extents can be explicitly delineated using growth modes and landscape expansion indices. Based on the urban growth modes identified using the landscape expansion index, we found that the leapfrog and edge-expansion modes dominated in the early stage of urbanization, which facilitated fast urban expansion, while the infilling mode dominated in the later stages and promoted the coalescence of urban patches. Each mode presented its alternating relative dominance, measured using the proportion of patch area and the patch number in different time periods during the urbanization processes in the PRD region ( Figure 6). Other studies that only used limited time period data could not find such an urban growth change pattern of an urban agglomeration. For instance, He et al. [29] and Zhou et al. [78] both identified that edge-expansion was the only primary growth mode in the PRDUA and the YRDUA, respectively, because they only used two and four time-period data, respectively. Yu and Zhou [79] revealed diverse patterns of growth modes in three Chinese urban agglomerations, such as edge-expansion dominating the urban growth in the BTHUA, edge-expansion and infilling dominating in the YRDUA, and infilling dominating in the PRDUA when using the land-use data of 2000, 2005, and 2010. Therefore, changing patterns of long-term urban growth modes can be helpful to explicitly detect the spatial and temporal evolutions of an urban agglomeration over time.

Morphological Dynamics with Urbanization Processes of Urban Agglomeration
Urban morphology can display highly dynamic and dramatic changes in spatial patterns during rapid urbanization processes. The spatial dynamic of urban morphology with urbanization processes has been widely detected in many countries and regions at the city scale using landscape metrics and/or urban expansion indices [16,21,30,38,42,49,55,70]. A wave-like change pattern of the urban fringe, which was proposed firstly by Blumenfeld [39] and tested by Boyce [40] to demonstrate metropolitan expansion, has been observed. The diffusion-coalescence theory has also been widely tested in individual cities to demonstrate the urban morphological change processes [16,21,38,42]. For urban agglomerations, previous studies successfully detected the morphological changes of urban agglomerations worldwide, but they only partially revealed the urban morphology, such as the landscape pattern [26,78] and the diffusion-coalescence processes of urban growth [25,30]. The earlier study on urban morphological change dynamics in the central YRD detected morphological changes of the urban agglomeration through multiple perspectives, simultaneously, in the landscape pattern, diffusion-coalescence processes, and dispersal-compact processes [16]. Our study revealed that the morphological change of the PRDUA exhibited a wave-like expansion pattern, co-existing diffusion-coalescence processes of the urban form, and an alternating switch between the compact urban form and dispersal urban form that corresponded to the coalescence and diffusion processes, respectively, during the urban form evolution of the past 30 years. Our study confirmed the spiral diffusion-coalescence hypothesis proposed by Li et al. [16], which demonstrated that the three growth modes of leapfrog, edge-expansion, and infilling alternated in relative dominance patterns in a spiral manner during the urbanization processes. Similarly, a study conducted in three Swiss urban agglomerations also supported the spiral diffusion-coalescence hypothesis [28].
Combining the LEI and AWMEI illustrates how the urban morphology changed during rapid urbanization in the PRDUA. The temporal change of the MEI and the AWMEI exhibited a wave-like pattern with two peaks, which implies that the PRDUA experienced two rounds of alternating between dispersal and compact urban forms, which correspond to the two rounds of diffusion and coalescence processes. The temporal variations of the MEI and AWMEI (Figure 7) exhibited a synchronous match with that of the proportions of the patch area and the patch number of infilling and edge-expansion urban patches (Figure 8), but were inverse to that of leapfrog. The diffusion process corresponds to the proliferation of the leapfrog in the early urbanization processes, which results in the dispersal of the urban morphology. The coalescence corresponds to the edge-expansion and infilling in the later urbanization processes, which make the urban form compact. This pattern of diffusion and compaction was also shown in Xinjiang, where the edge-expansion and infilling modes led to the compact growth of urbanized land [34].

Dual Reflectance of Landscape Metrics to Urban Form and Landscape Pattern Dynamics of Urban Agglomeration
Landscape metrics have been widely used to characterize the urban form and its changes during urbanization [49,53,80,81]. Most of the previous studies were from the perspective of individual cities to identify the compactness or sprawl of the urban form and classified the cities into different groups or types. For instance, Huang et al. [53] used seven landscape metrics and successfully identified the 77 selected cities from Asia, Europe, the USA, Latin America, and Australia as "sprawl" and "compact", and classified these cities into different groups in terms of the spatial metrics. They also revealed that the urban form of the cities in the developing world was more compact and denser than their counterparts in Europe and North America. Similarly, Schwarz [49] utilized landscape metrics to quantify the urban form of 231 European cities and classified them into eight clusters. Another study used landscape metrics to identify the urban form of 194 cities of the world into four types, namely "compact-grey, transitional, ragged-small, and fragmented-complex", and they found that the urban forms of these cities tended to be more homogeneous and most of them were in transitional states as a result of their fragmentation and compactness [80]. However, these studies only used one or two time series of land-use data to derive the landscape metrics and could not capture the temporal trajectories of the urban form evolutions over time. Using the landscape metrics derived from the long-term time series land-use data, our study demonstrated that, on one hand, the urban form of the entire PRDUA generally exhibited an increasingly complex evolution trend, reflected by the increasing edge density (Figure 9a), area-weighted mean fractal di-mension index (Figure 9b), and regularity indicating by landscape shape index (Figure 9g), and the compact, demonstrated in the declining patch density (Figure 9e), increasing largest patch index (Figure 9b), contiguity (Figure 9d), and aggregation index (Figure 9f). They are similar to the increasing trends in the fragmentation, irregularity, and complexity of the urban form of the YRDUA when measured using landscape metrics [16,78]. On the other hand, the changing trajectory of each landscape metric showed a synchronous pattern with the diffusion and coalescence that resulted in the compactness or dispersal of the urban form of the entire urban agglomeration, which was also observed for the urban agglomeration in the central YRD region [16].
The evolution of the PRDUA influenced the regional landscape pattern, which exhibited a similar change trajectory to that of the urban form itself. Our results demonstrated that the regional landscape of the PRD region also presented temporal change patterns associated with the spatial and temporal evolution of urban agglomerations. For instance, the regional landscape patterns also exhibited wave-like trajectories (Figure 10a,c,e,g), which corresponded to the temporal change patterns of the urban expansion (Figures 5 and 8) and urban form (Figure 9a,c,e,g). On one hand, the shape complexity (Figure 10a,c) and irregularity (Figure 10g) of the regional landscape increased over time, and that of the connectivity of the regional landscape decreased (Figure 10f,h), while the fragmentation increased (Figure 10b,d), which were similar to those observed in the YRDUA [16,78] and those in the BTHUA [57].

Limitations
There are some limitations in our research that should be addressed in future studies. Firstly, the more detailed urban expansion or the temporal change of urbanized land may not be fully captured using the five-year interval due to the high frequency of land surface changes in the PRDUA, especially for Shenzhen and Guangzhou [82,83]. Shortening the inter-annual interval for the analysis could render a more detailed view of the temporal features of urbanization, which is helpful in understanding the evolution of an urban agglomeration.
Secondly, the details of the spatial changes in the landscape patterns may not be well captured by the 30 m resolution urbanized land derived from the Landsat images. In general, the grain size and extent have been determined by default when the study area was determined; in particular, the data used were derived from remote sensing images. The higher the spatial resolution, the more accurate and finer the scale information of the landscape patterns presented [84,85]. Therefore, fine resolution is preferred when computation efficiency is not a limitation. However, it is hard to acquire satellite data with high resolution at periodical intervals [86], particularly at the regional scale. Therefore, a reliable spatiotemporal data fusion method, such as STARFM or FSDAF [87,88], to improve the spatial resolution and frequency of the images can be further explored.
Thirdly, the dataset we used allows the classification of urban, non-urban, and green recovery, which could not fully describe the complex features of landscapes, such as residential landscape and commercial landscape, and their dynamics. Thus, more detailed classification schemes for time-series remote sensing data for the analysis of landscape changes may provide more accurate characterization. Recently, increasing land-use and land cover datasets at national and global scales have been made freely available and accessible [89], providing us with the possibility to conduct more detailed urban expansion monitoring in the future.
Lastly, most of the landscape expansion indices that were proposed in the previous studies [33,34] only described the temporal change of urban expansion, but overlooked the spatial changes of the evolution process of urbanization. For example, how the cities and in-between areas spatially integrate into an agglomeration and how the interconnected transportation develops [37] still require further studies. Thus, it is suggested to conduct more spatial analysis, and even multiple spatial and temporal scales analyses, to provide more valuable insights for understanding the evolution of an urban agglomeration.

Conclusions
Based on the spatiotemporal trajectories of the urban agglomeration and its impacts on the urban and regional landscape patterns in the Pearl River Delta region between 1985 and 2015, we found that the urbanized land area of the entire urban agglomeration increased linearly, with an annual growth rate of 4.6%, while the nine major cities within the PRDUA presented diverse growth trajectories, varying between linear and non-linear. The PRDUA had experienced two cycles of urbanization processes, which were well reflected by a bimodal wave-like pattern of the average annual growth rate of the PRDUA over time. Three urban growth modes-infilling, edge-expansion, and leapfrog-existed simultaneously, but only one mode dominated at a particular period, and the dominant mode alternated during the urbanization processes, demonstrating that the spatiotemporal evolution of the urban growth of the PRDUA follows the "spiral diffusion-coalescence hypothesis". The morphology of the PRDUA correspondingly showed an alternating dispersal-compacting pattern over time. The city-and regional-level landscape patterns changed synchronously with the spatiotemporal evolution of the PRDUA over time. The response of the landscape patterns agreed well with the change in the growth rate and growth modes. The urbanization of the PRDUA had increased the landscape complexity at the city and regional scales, but also resulted in the increasing fragmentation and decreasing connectivity of the natural landscape in the PRD region. These findings help to better understand how an urban agglomeration evolves and provide insights for regional urban planning and the sustainable development of urban agglomerations.