Using multisource data and the V-I-S model in assessing the urban expansion of Riyadh city, Saudi Arabia

ABSTRACT This paper examines the application of remote sensing, based on the Vegetation-Impervious surface-Soil (V-I-S) model and spatial metrics, in an urban analysis for promoting sustainability and understanding urban growth theory. In order to improve the accuracy of land-cover classification, spectral angle mapping (SAM), spectral mixture analysis (SMA) and band ratioing were applied on satellite images for land-cover classification and comparison of the discrimination efficiency of these techniques. For the SMA, subsets of the Landsat (2, 3, 4, 5, and 7) and ASTER (1, 2, and 3N) images were selected. After endmember extraction and purification, the Bayesian probability of each component was computed and used for spectral unmixing. The classified images of different years were compared to analyze the changes in land-use and spatial pattern using V-I-S, a form of percentage of landscape (PLAND) and annualized urban sprawl index (AUSI). The result indicates that the performance of band ratioing (69% accuracy) is not as good as that of SAM (75%) and SMA (86%) in discriminating between vegetation and agricultural land. It is concluded from the land-use analysis that the growth dynamics substantiate the urban theory of diffusion and coalescence and urban growth management strategies have not been completely successful.


Introduction
The city of Riyadh is one of the fastest-growing cities in the Middle East. The population of the city is currently about 6.5 million and the annual population growth rate of Riyadh, although already declining, stands at around 4% (ArRiyadh Development Authority (ADA), 2016). However, the size of growth is a major influencing factor that threatens the sustainability of the city. The resultant spatial footprint pattern and expansion have been extensive and the city's population density is declining in a typical sprawl format. The spatial extent of Riyadh is currently about 1000 km 2 and the population density has decreased from 38000 people per square km in 1950 to 6825 in 1999(Al-Sahhaf, 2000 and 4300 in 2013 (http://www.atlasofurbanexpansion.org/cities/view/ Riyadh). Urban growth management strategies (urban growth boundaries -UGB) have been adopted by ArRiyadh Development Authority to reduce spatial growth and foster sustainable urban development (ArRiyadh Development Authority (ADA), 2003, ArRiyadh Development Authority (ADA), 2009, UN-HABITAT, 2017). This is in congruence with the Saudi national spatial strategy of promoting a decentralized and spatially balanced urban system (Aldalbahi & Walker, 2015;Alkhedheiri, 2002). However, there are indications that the strategies have not been totally successful (Aldalbahi & Walker, 2015;Al-Hathloul & Mughal, 2004;Mubarak, 2004). Urban land has extended beyond the urban growth boundaries and natural habitats and vegetation are being changed and encroached upon. The prevailing urban growth needs to be monitored to ensure sustainable urban development.
The Vegetation-Impervious surface-Soil (V-I-S) model is being applied in deriving the biophysical composition of the urban environment and thereby in analyzing urban environmental change (Kuang, 2019;Li et al., 2019). The application of the V-I-S model might be driven by data derived from remote sensing since different approaches have been developed to extract the biophysical composition of urban areas from satellite images (Ridd, 1995). The practical application of the V-I-S model, especially in urban planning, is still restricted due to unsatisfactory levels of accuracy in discriminating and quantifying V-I-S components in urban areas and the mixed pixel problem (Weng, 2012;Guo, Huang, Li, Sun, & Zhang, 2014). Though the level of accuracy has been improving, the mixed pixel problem is still a concern in using medium spatial resolution images (Weng, 2012). Thus, there is a need to verify the discriminating capability of classification methods before adopting them for the V-I-S model. The model can be applied at different scales from the sub-pixel level to an entire city. The focus of this study is on the city scale.
Band ratioing is one of the existing change detection and land-use mapping techniques (Borra, Thanki, & Dey, 2019;Lu, Mausel, Brondizio, & Moran, 2004a). Band ratioing can reduce the effect of illumination variation resulting from environmental factors such as topography (Borra et al., 2019;Mather & Koch, 2011). Al-Saud (2005), in a study of Al-Hasa in Saudi Arabia, applied a set of band ratios to discriminate land features in an arid environment. In the same vein,  used a set of indices based on ratioing to classify land covers in some Saudi cities. The set of band ratios by Al-Saud (2005), which has not been rigourously tested, could be further examined in a similar environment (Al-Saud, 2005).
The spectral angle mapping (SAM) algorithm computes the similarity between a known (reference) spectrum and an unknown (target) spectrum and a group of similar spectra is classified into a specific class (Jensen, 2015). The spectra (both known and unknown) are considered as vectors and the angles between the vectors are computed to evaluate their similarity (Renza, Martinez, & Molina, 2017). The angles usually vary from 0°to 90°and 0°indicates high similarity while 90°indicates dissimilarity (Mather & Koch, 2011). Similarly, spectral mixture analysis (SMA) is based on spectral analysis. It is used to extract features from an image by computing the fraction of each feature in a pixel (sub-pixel analysis) (Mather & Koch, 2011). Spectral mixture analysis can be used to generate urban impervious surfaces and land covers without elaborate training data if the endmembers are adequately defined (Van de Voorde, Jacquet, &Canters, 2011). MacLachlan, Roberts, Biggs, andBoruff (2017) asserted that SMA is increasingly used for representing urban land cover using the V-I-S framework due to the capability of extracting fractions of urban materials at sub-pixel level. However, the endmembers might not adequately represent the complex spectral heterogeneity of urban areas (MacLachlan et al., 2017). Thus, the selection of appropriate endmember for each feature is critical in spectral mixture analysis because endmembers are used as a reference to determine the spectral components of the pixels. This study uses image subset and typicality threshold to select endmembers for spectral mixture analysis and compares the discrimination efficiency of the SMA with that of the other techniques (SAM and Band ratioing) to ascertain the best method to use with the V-I-S model. The SAM and the SMA methods were selected because the SMA is widely used for extracting impervious surfaces and suitable for the V-I-S model (Kuang, 2019) while the SAM is based on the same concept of analyzing spectra.
Urban and remote sensing researchers have asserted that the analysis of land-use change and urban spatial expansion needs to be linked with urban form and processes for it to lead to effective urban intervention and planning (Dahal, Benner, & Lindquist, 2016;Dietzel, Herold, Hemphill, & Clarke, 2005;Longley & Mesev, 2001). Dietzel et al. (2005) and Dahal et al. (2016) applied spatial metrics to show, by presenting empirical results for some North American cities, that the spatial evolution of urban areas oscillates between two phasesdiffusion and coalescence (theory of urban growth). The question still remains whether other cities, especially in the developing countries, will follow the same pattern.
Recent studies have focused on urban growth in Riyadh due to the rapid urbanization and spatial growth occurring in the city. Mashee and Malthus (2009) examined the urban expansion in Riyadh from 1994 to 2007 by using an object-oriented technique. Al-Ahmadi, See, Heppenstall, and Hogg (2009a, 2009b, 2009c and Alqurashi, Kumar, and Sinha (2016) focused on developing, calibrating and validating predictive urban growth models of Riyadh city. Although the studies presented some of the changes in land use and spatial structure, the focus was on modeling and predicting urban spatial growth by 2020. This study explores the application of band ratioing and spectral mapping techniques for land-use change and urban growth analysis and how the results of urban growth analysis are linked with urban growth theory. The objectives of the study are to (1) examine the spatial and temporal change in land use and urban growth in Riyadh using the V-I-S model at the city scale; (2) compare the efficiency of band ratioing and spectral mapping methods in discriminating and mapping urban land cover in an arid environment; (3) link land-use change and urban growth analysis in Riyadh with urban growth theory using a spatial metric; and (4) assess the effectiveness of the urban growth boundaries implemented in Riyadh.

V-I-S model
There is a growing interest in adopting the V-I-S model as a simple means of analyzing changes in the biophysical composition of urban areas. Ward, Phinn, and Murray (2000) applied the V-I-S model in southeast Queensland, Australia to examine the changes in land-cover types and urban growth between 1988 and 1995 and noted the need for further work on developing methods for evaluating urban composition and comparing V-I-S model results from different cities. Setiawan, Mathieu, and Thompson-Fawcett (2006) demonstrated the use of the model in developing countries and asserted that the model can be applied in comparing land-cover changes in an urban area over a period of time or comparing different urban areas. They noted that the model is meant for characterizing urban land cover but can also be used for land-use study because there is a relationship between V-I-S components and land use. Lin, Guo, Yan, and Heng (2018) used the V-I-S model in classifying the land-use components of an urban area in China to landscape pattern changes.

Improving the accuracy of image classification
In order to improve the accuracy of extracting urban composition from satellite images and also to resolve the mixed pixel problem, some studies (Fan, Fan, & Qihao, 2015;Jacobson, 2014;Kwarteng & Small, 2010;Weng, 2010) had used spectral mixture analysis for image classification. The mixed pixel problem, which is resolution dependent, arises from the basic assumption of pixel-based classification techniques that the pixel is homogenous. In reality, most pixels in moderate-and low-resolution images comprise two or more components that are found in spatially varying fractions of the pixels in an image (Sliuzas, Kuffer, & Masser, 2010). So, a pixel might be classified as a feature while it comprises two or more features having equal fractions. The mixture analysis method is effective in handling mixed pixel problems and providing spatially continuous biophysical variables (Fan et al., 2015;Feizizadeh & Blaschke, 2013;Weng, 2010;Wetherley, Roberts, & McFadden, 2017) if appropriate endmembers are selected. Mather and Koch (2011) asserted that spectral angle mapping could be considered as an alternative to mixture analysis to avoid the critical step of selecting endmembers in the mixture analysis method since they are both based on the analysis of spectra. Dennison, Halligan, and Roberts (2004), in an earlier study, showed by using a spectral library (compilation of spectral signatures of surface materials) that multiple endmember SMA and SAM produced similar results. Though the results were similar, the outcome might have been different if the endmembers were taken from an image and not from a spectral library. Yonezawa (2007) and Pan, Hu, and Wang (2019) demonstrated the capability of SAM by using a SAM algorithm to complement other methods in order to improve the classification accuracy. They asserted that the overall accuracy of the land-cover classification was improved by the SAM algorithm.
Band ratioing is another technique that had been used to provide urban composition data for the V-I-S model. Setiawan et al. (2006) applied the hierarchical classification scheme (H-VIS) in mapping urban land use in Indonesia. The classification scheme was conducted by applying thresholds to the Normalized Difference Vegetation Index (NDVI). The study demonstrated the applicability of band ratioing for land-use classification. Shrestha and Conway (2011) also successfully delineated ex-urban development by using band ratioing (NDVI) and binary recoding.

Analysis of urban landscape
The combination of spatial metrics with urban mapping and analysis methods, such as remote sensing, can provide new information about urban dynamics and spatial structure (Herold, Couclelis, & Clarke, 2005;Fichera et al. 2012). One of the spatial metrics applied by Dietzel et al. (2005); Weng, Liu, and Lu (2007); and Aina et al. (2008) in analyzing urban spatial growth is the percentage of landscape (PLAND). PLAND is a metric that shows the percentage of the landscape occupied by each land-cover class per unit area. Though PLAND is a simple metric, as a measure of landscape composition and dominance, it can capture the characteristics of urban landscape as shown in the studies by Kong, Yin, Nakagoshi, and James (2012); Chen, Yao, Sun, and Chen (2014); and Sertel, Topaloğlu, Şallı, Yay Algan, and Aksu (2018). A multidate computation of the metric depicts the dynamics of landscape composition (Shrestha et al. 2012). A drawback of the metric is the lack of spatial explicitness (Liu et al., 2010). That is, the metric does not depict urban spatial pattern or changes in spatial structure in a location-specific manner. Liu et al. (2010) suggested a new spatially explicit index called the land expansion index (LEI). The index can explicitly indicate if spatial expansion is due to infilling, edge-expansion or outlying growth. It is computed by creating a buffer around the newly grown urban patch and determining if the patch is infilling, edge-expansion or outlying growth based on the percentage of intersection between old and new patches. The value ranges between 0 and 100. The computation of the index is more complex than that of PLAND and it is also sensitive to buffer distance. A different method of calculating PLAND, based on a changing circular urban extent (Urban Field), is explored in this study to make the metric more amenable to depicting changes in urban spatial structure. An annualized urban sprawl index (Hepinstall-Cymerman, Coe, & Hutyra, 2013) is also explored in studying landuse change. It is based on land consumption by the population. Moreover, it is relevant to the United Nations Sustainable Development Goals (SDGs) as indicator 11.3.1 focuses on land consumption (https:// unstats.un.org/sdgs/metadata/?Text=&Goal=&Target= 11.3).

Study area
The study area covers 1488 km 2 of the Riyadh city area. Riyadh is the capital city of Saudi Arabia (see Figure 1) and plays a major role in the Saudi urban system as an administrative, legislative, and commercial centre. It is bounded in the north and east by the Ad-Dhana sand dune landscape, in the west by the Arabian shield and in the south by the Rub Al-Khali desert (empty quarters). The city is located in a desert landscape about 400 km from the east coast (Gulf) and 1000 km from the west coast (Red Sea) of Saudi Arabia. Riyadh has experienced rapid urbanization especially during the oil boom period in the 1970s and 1980s. Rapid spatial expansion of the city has been driven by the oil economy and automobile dependency.). The ecological footprint coefficient of the city, which is 4.5, is higher than the coefficients of most cities in the neighbouring countries (Aina, 2017;Choguill, 2006). The sustainability of the city is threatened by the increasing demand for water and infrastructure, traffic pollution and encroachment on natural vegetation (ArRiyadh Development Authority (ADA), 2003, UN-HABITAT, 2017). Currently, urban infrastructures such as electricity, telephone, water, and sewage are being overstretched and the percentage of population provided with public utilities is reducing (Abubakar & Aina, 2016).
In order to manage rapid urban growth, two phases of urban growth boundary had been implemented by ArRiyadh Development Authority -ADA (2003,2009). That is, the boundary beyond which new urban development is not officially supported with infrastructure and facilities. The urban growth boundary (UGB) for Phase I (632 km 2 ) was delineated for 1995 while the boundary for Phase II (1,149 km 2 ) was delineated for 2005 and Phase III (2,130 km 2 ) was delineated for 2015 (see Figure 2). There is an additional boundary regarded as the development environs. It is about 3,120 km 2 in extent and it can be allocated to cater for future growth. The boundary might become Phase IV of the UGB series. The main  objectives of the growth boundaries are to make the footprint of the city more compact, protect agricultural land and natural areas and reduce pressures on the provision of additional infrastructure to new haphazard settlements (Alkhedheiri, 2002).

Data and image preprocessing
The data sources for the analysis include multidate and multisource satellite images from 1972 to 2014. Vector data (shapefile format) for the two phases of UGB and land use were also acquired for the study from the ArRiyadh Development Authority. The land-use shapefile contained the land-use data of Riyadh and it was used for classification accuracy assessment (complemented by Google Earth TM images of the area). The satellite images include Landsat MSS (1972), Landsat TM (1985 and1995), ASTER (2005) and Landsat OLI (2014) images (see Figure 3). ASTER image was used because of the SLC-off problem of Landsat data acquired after May 2003 that have resulted in data gaps in Landsat images. Anniversary date images were chosen to reduce errors due to seasonal changes (see Table 1). The satellite images were all georeferenced with a root-mean-square error of less than 0.25 pixel and projected to UTM (Universal Transverse Mercator) Zone 38 North. The 1972 image has a coarser resolution of 60 m. The images, except the 1972 image, were resampled to 30 m using the nearest neighbour method. The satellite image coverages for 1985 and 1995 comprised two scenes each, since the study area was not fully covered by single scenes. The images were atmospherically corrected by using the dark object subtraction method (Mather & Koch, 2011) implemented in IDRISI TM Selva 17.2 software (Eastman, 2012). Image mosaic and subset operations were carried out on images with two scenes to extract the development area of Riyadh for the study. The location of the testing samples (275) for accuracy assessment is shown in Figure 4. The training samples, which are the endmembers, are separate from the testing samples. The training samples were selected from the image, five locations for each class and 10 pixels from each location.

Image classification and analysis
A post-classification technique for land-use change detection and analysis was adopted in this study. Three different classification techniquesspectral angle mapping, spectral mixture analysis, and band ratioingwere applied for the classification. The ratios presented by Al-Saud (2005) were used to classify the images into four classes. The ratios were based on Landsat TM imagery but they can be used to classify images with similar bands. However, in this study, the ratios were assessed on classification of Landsat imagery only. The ratios are Green Band/ Red Band, Near-Infrared Band/Shortwave Infrared I Band, Shortwave Infrared I Band/Thermal Band, and Shortwave Infrared I Band/Shortwave Infrared II Band. The ratio Green Band/Red Band is very good for discriminating natural vegetation while ratio Near-Infrared Band/Shortwave Infrared I Band is good for extracting barren land. Band ratios Shortwave Infrared I Band/Thermal Band and Shortwave Infrared I Band/ Shortwave Infrared II Band are good for discriminating urban area and natural vegetation, respectively (Al-Saud, 2005). After the ratios had been computed, different thresholds were applied to the results to get the most suitable threshold (based on trial and error and visual assessment) for each feature.
The IDRISI TM Selva 17.2 software package was used in performing the sub-pixel classification (Eastman, 2012). The most important aspect of spectral mixture analysis is the selection of endmembers. Authors such as Small (2006); Lu, Batistella, Moran, and Mausel (2004b); and Weng (2010) mentioned that applying principal component analysis (PCA) or minimum noise fraction (MNF) or using a subset of Landsat images can reduce the correlation between images. This study used subsets (2, 3, 4, 5, and 7) of the Landsat images and subset (1, 2, and 3N) of the ASTER image. Endmembers were selected from the images and the parametric option of the PURIFY module (Eastman, 2012) of IDRISI TM was applied to the selection to improve the purity of the endmembers. The module calculates the typicality of each training pixel by computing the distance between the pixel and the mean of its class. Pixels with low typicality are dropped from the sample and pixels with high typicality are included in the training sample based on a threshold chosen by the user (0.7 was chosen for this study). The pixels with lower typicality might not be pure enough for the study and choosing higher typicality would result in rejecting most of the pixels.
After the purification of the endmembers, the UNMIX module (Eastman, 2012;Warner & Campagna, 2013) was used to generate the fractions of land-cover types in the study area. The module uses the Linear Spectral Unmixing (LSU) method to compute the fractions. The method assumes a linear combination of the pure training samples (endmembers). The mathematical equation of LSU can be written as where R is the spectral reflectance of a mixed pixel in the spectral bands, N is the number of endmembers, f i is the fraction of ith endmember, r i is the reflectance of ith endmember, e is the error term. The approach is most effective if the endmembers are few and distinct (Eastman, 2012). This study applied the probability guided unmixing method. In this method, the Bayesian probability of each component is computed prior to mixture analysis and the user can also select the number of components to consider in each pixel (Eastman, 2012;Song, 2005). The Bayesian probability method uses the endmember signature probability distribution rather than the constant signature used by a typical LSU analysis (Song, 2005). Then, the HARDEN module (Eastman, 2012;Warner & Campagna, 2013) was used to assign the fractions to classes. The HARDEN module assigns a pixel to the class with the maximum fraction. After the classification, a majority filter (3x3) was applied on the classified images to remove isolated pixels. The flowchart of the classification process is presented below (see Figure 5). The results of the image classifications of different years were compared to depict the changes in land use and spatial pattern of urban growth. The changes were analyzed by computing the PLAND, the sprawl index and percentages of land-cover types. PLAND is the percentage of landscape or the sum of areas of a land-cover type divided by the total area. Classification accuracies were also computed for the different techniques to ascertain the most appropriate method for land-cover classification in the study area. The SMA method was chosen for the analysis due to the level of accuracy and being suitable for extracting the urban materials for the V-I-S model. The urban growth boundary shapefiles were used to assess, in spatial context, the effectiveness of the urban growth boundary policy. A metric (PLAND) (Dietzel et al., 2005) was computed for urban extent in each classified image in order to examine the linkage between urban growth and spatial processes (changes in spatial structure of a city).
Traditionally, the total area is computed from the administrative boundary of the area and is thereby constant. This study explored a varying total area by computing the urban field area for each year. It is the area of a circle that bounds the extent of the study area (Longley & Mesev, 2001). The "Urban Field" PLAND (UFPLAND) could depict the process of urban spatial growth as diffusing or coalescing just like the PLAND (Dietzel et al., 2005). If the growth is diffusing, it is expected that there will be scattered "seeds" of centers of growth over the landscape and UFPLAND value will be low (Sertel et al., 2018). On the other hand, if the process is in a state of coalescence, urban extent will be more compact and UFPLAND value will be high (Sertel et al., 2018). The UFPLAND might not be sensitive to spatial pattern and processes within the urban field especially if the UFPLAND is low. Thus, in a typical urban setting, the PLAND will increase with time showing mere urban expansion (Dietzel et al., 2005). On the other hand, UFPLAND has a varying total landscape area. If the total landscape area is expanding over a period of time, it means more area outside the spatial extent has been developeddiffusion. If the total land area is fixed over a period of time, it means development is mainly taking place within the spatial extentcoalescence. Though the metric will not indicate the differences in pattern of change within the spatial extent, it will indicate the overall pattern.
An annualized urban sprawl index (AUSI) was also computed for the study area to further depict the changes in urban land consumption throughout the study period. AUSI was computed by dividing an annualized change in urban extent (in km 2 ) by an annualized change in population size for the same period (Hepinstall-Cymerman et al., 2013). Changes in urban land cover within and outside the urban growth boundaries were analyzed by calculating the amount in km 2 of new urban land cover within and outside the urban growth boundaries.

Results and discussion
Accuracy assessment The classification result shows that the band ratioing method can be applied to discriminate landcover types especially if an appropriate threshold is used. But the performance is not as effective as the SAM and SMA methods in discriminating between vegetation and agricultural land. The performance of all three methods is actually poorer in discriminating between agricultural land and vegetation than mapping other features (see Tables 2-4). The overall accuracy for band ratioing is 69% while the overall accuracy for SAM is 75%. The overall accuracy for spectral mixture analysis ranges from 83% to 86% for the different years. The highest kappa statistics for SMA is 81%. It should be noted that the urban extent tends to be slightly overestimated by the SAM algorithm in 2005 while the band ratio method slightly underestimates the urban extent (upon visual assessment) (see Figure 6). The urban extent in 2005 according to the SAM technique is 638 km 2 while the SMA urban extent for the same year is 622 km 2. A limitation of the study    Figure 6. Image classification (1985)band ratio, SAM, and SMA methods.
is that resampling the images might introduce some errors that might affect the accuracy of the classifications. Nevertheless, the nearest neighbor used in this study preserves the pixel values of the images (Borra-Serrano, Peña, Torres-Sánchez, Mesas-Carrascosa, & López-Granados, 2015). The band ratio method outperforms the SAM in depicting the gaps (undeveloped lands) within the urban area. It is also more suitable for extracting the road network if the appropriate threshold is applied. Despite having lower accuracy than the SMA, the SAM technique depicted similar results upon visual inspection (see Figure 6). In their study of five Saudi cities,  obtained an accuracy of 87.5% (for 1985 image) to 92.8% (for 2014 image) by using object-oriented image analysis (OBIA) method. They classified vegetation, water, urban, and bare soil . Since 85% is widely used as sufficient accuracy (Foody, 2006) and the use of urban fractions is more appropriate to the V-I-S model, LSU analysis might be preferred. The study by Lu and Weng (2006) yielded 83.78% as the overall accuracy. A recent study by Guo et al. (2014) achieved the overall accuracy of 96.69% using a support vector machine (SVM). However, the classes are different from the ones used by this study. The classes are high albedo, low albedo, vegetation, and soil.

Change analysis
As mentioned above, the Spectral Mixture Analysis technique was used in the change analysis. The analysis shows that the urban extent in Riyadh increased from about 43 km 2 in 1972 to 780.6 km 2 in 2014 (see Table 5). The period between 1972 and 1985 had the highest spatial expansion from 43 km 2 to 396 km 2 . The growth rate of urban extent which had started tapering as shown by the urban land expansion from 1995 to 2005 has picked up again from 2005 to 2014. The lowest rate of expansion was recorded in 2005. The urban spatial extent compiled from the satellite images is actually less than the 1000 km 2 reported by ArRiyadh Development Authority (ArRiyadh Development Authority (ADA), 2003, UN-HABITAT, 2017) but it is consistent with findings of Mashee and Malthus (2009) (see Figure 7). This might be due to leap-frog development outside the city fringe. Some land areas of the city are just subdivided without being developed immediately and they remain vacant for many years.
The vegetation cover had a different trend from the urban area. Vegetation cover increased from 13.4 km 2 in 1972 to 39.6 km 2 in 1995 but later contracted to 29.8 km 2 in 2005. The decrease in vegetation cover indicates unprecedented built encroachment upon vegetation and agricultural land by urban areas. It was only in 2005 that the land area gained by urban land use (105.6 km 2 ) exceeded the land area lost by barren open land (−95.8 km 2 ). The findings corroborate the earlier works by Al-Sahhaf (2000), Al-Ahmadi et al. (2009a, 2009b, 2009c and Mashee and Malthus (2009) that had highlighted the rapid urbanization in Riyadh. The changes in vegetation cover have implications on the environment in the Riyadh area. The increase in vegetation cover might have implications on water consumption and underground water usage while the decrease might reduce the mitigating effect of vegetation on the microclimate as suggested by Stabler, Martin, and Brazel (2005). The trend of reduction in vegetation cover has changed in 2014 with an increase in the area covered by vegetation from 29.8 km 2 in 2005 to 34.5 km 2 . The fractions of impervious surface and vegetation also highlight the changes discussed above. The mean fractions changed with an increase or decrease in the spatial extent of the urban components (see Table 6). The annualized sprawl index indicates a decreasing trend in land consumption from 1972 to 2014 (see Table 7). There is an indication that the urban growth boundary implemented by the Riyadh Authority has not been completely successful. Urban development had occurred outside the 1995 and 2005 boundaries in 1995 already (see Figure 8 and Table 8). Even the urban field area had reached 2300 km 2 since 1985 but had remained constant till 2014. Despite the shortcomings, the policy has helped to curb haphazard and outlying (leapfrog) spatial growth in the city by achieving a constant urban field since 1985 and limiting urban development to this area. However, with the new trend (2014) (see Table 8) of an increase in land development outside the growth boundaries (even the 2015 boundary), the Riyadh Authority might need to increase efforts in curbing sprawl. The new spate of land development might be attributed to the commencement of the third phase of UGB and increasing population pressure.
The result of the spatial metric computation indicates that the urban growth process in Riyadh followed the urban theory of diffusion and coalescence highlighted by Dietzel et al. in 2005and Dahal et al. in 2016. Liu et al. (2010 and Dahal, Benner, and Lindquist (2017) mentioned three types of growth during the diffusion-coalescence phases of urban growth; in-filling, edge-expansion, and outlying growth. Outlying growth occurs mainly at the diffusion phase while in-filling and edge-expansion occur mainly at the coalescence phase. The studies (Dietzel et al., 2005;Dahal et al., 2016Dahal et al., , 2017Liu et al., 2010) did not use LSU in their classifications but adopted spatial metrics in monitoring urban expansion. The growth of Riyadh started within the city walls in the 1950s (ArRiyadh Development Authority (ADA), 2003). Thereafter, the city expanded outward in a gradual manner (ArRiyadh Development Authority (ADA), 2003). In the 1970s, due to an unprecedented increase in government spending the spate of development increased, which led to the leap-frog growth witnessed later (Garba, 2004). The growth between 1972 and 1985 constituted the diffusion phase whereby expansion was mainly through outlying growth. After 1985, the spatial expansion of the city started tapering as new development was largely due to edge-expansion and "in-filling" rather than isolated and outlying growth at the edge of the city. This was the coalescence phase.
The lowest PLAND figure is 3% for 1972 and thereafter, the figure increases to 42% in 2005 and 52% in 2014. This indicates the rapid change of land cover from barren open land to impervious surface due to urbanization. The change in PLAND value for vegetation depicts a low percentage of coverage by the green area in Riyadh. Unlike the PLAND values, the UFPLAND depicts the phases of the changes occurring in the land use and land cover of the study area. The highest UFPLAND figure is 32% for 1972 which indicates coalescence. Thereafter, the figure declines to 17% in 1985 (diffusion) and picks up again in 1995 and 2014 to reach 34% (coalescence) (see Figure 9). The increases in UFPLAND values from 1985 to 2014 are due to the fact that the urban field area is constant (2300 km 2 ) while urban land-cover area is increasing. It should be noted that the PLAND and UFPLAND of 1972 might have been affected by the differences in sensor resolution. Application of the V-I-S model to Riyadh also depicts the same process. The land-cover composition in 1972 was desert according to V-I-S ternary because desert had the largest share of the landscape. It changed to low-density residential in 1985 and medium to high density in 1995, 2005 and 2014 as the impervious surface increased. This also corroborated by the fractions of the impervious surface. The mean increased from 0.35 in 1972 to 0.63 in 2014 indicating increasing imperviousness. The mode of the fraction of the impervious surface also increased from 0.05 to 1 indicating an increase in the size of the fractions thereby depicting coalescence. The percentages of barren open land and urban area were 63% and 35% respectively in 1995 and 56% and 42% respectively in 2005. This indicates that Riyadh became a medium-density residential area in 1995 and 2005 and it is progressing towards high density.

Conclusion
This paper presents a form of a spatial metric (percentage of landscape) and the V-I-S model being applied to the city of Riyadh. The metric highlights the changes in urban spatial pattern in a better way than the traditional version. It depicts the changes in the diffusion and coalescence phases of urban spatial growth. The V-I-S model and AUSI were applied in highlighting the changes in land cover and the effectiveness of UGBs in Riyadh. In regard to Riyadh, this study has shown the spatial impact of urban expansion, especially on vegetation loss. There had been an increase in vegetation cover in Riyadh since the 1950s until around 1995 when encroachment on vegetation cover became apparent. Though the vegetation loss has been reversed in 2014, the task at hand is to develop a means of adequately forestalling increases in vegetation loss. On the other hand, the metrics of urban expansion are improving as indicated by the reduction in the rate of spatial expansion and increase in the urban area's percentage of landscape.
The analysis has, therefore, shown that despite a few failures, the urban growth boundary strategy in Riyadh has been working to some extent. The hitherto rapid spatial expansion has slowed down and the "gaps" and brownfields within the city are being developed. There is a caveat. How can it be confirmed that the changes are due to planning intervention and not the normal processes of urbanization, suburbanization, and re-urbanization or oscillatory diffusion and coalescence? Having used these techniques to understand the spatial processes of the city, how do planners know when and where to intervene? Arriyadh Authority had made some efforts through the initiation of the Urban Growth Strategy for Riyadh and this could be a starting point towards better urban growth management.
This work could still be extended further by looking into other spatial metrics (Taubenböck, Wurm, Geiß, Dech, & Siedentop, 2019) that were used in testing the hypothesis of diffusion and coalescence. Such a study will help in confirming the universality of the theory.

Acknowledgments
The authors are grateful to the anonymous reviewers for their useful comments. The authors would like to acknowledge the assistance of the following institutions in the acquisition of the necessary data for this work -King Abdul Aziz City for Science and Technology (KACST), ArRiyadh Development Authority (ADA), Japan Aerospace Exploration Agency (JAXA), United States Geological Survey (USGS) and Global Land Cover Facility of University of Maryland. The

Disclosure statement
No potential conflict of interest was reported by the authors.