Detecting and Assessing the Spatio-Temporal Land Use Land Cover Changes of Bahrain Island During 1986–2020 Using Remote Sensing and GIS

The Kingdom of Bahrain has experienced accelerated development growth since the 1980s. These rapid land demands increased the pressure on the country area to rebuild urban centers and cities surrounding the coast. The purpose of this research is to detect and investigate changes in land use and land cover (LULC), which is one of the most critical aspects of planning and managing the use of land as a natural resource. The massive growth in land demand, particularly in small-area countries like Bahrain, forces decision-makers to re-plan the main island areas (Bahrain, Muharraq, Sitra, and Nabih Saleh). The study focuses on mapping the LULC changes detection over 1986–2020. It employs an integrated approach of remote sensing and GIS (Geographic Information System) to analyze and evaluate the changes in the LULC area in the main islands using multi-temporal and multispectral Landsat satellite imagery acquired in 1986, 1994, 2000, 2005, 2013, and 2020. In addition, high-resolution satellite images of different dates IKONOS 2000, GeoEye1 2011, 2013, Worldview3 2019, ASTER 2012, 2013, and multiresolution seamless image database-MrSID 1994, 1998 were used to enhance the LULC classification. Furthermore, different ancillary data were utilized to adjust the decision of LULC classes. The images were supervised using Maximum Likelihood Classifier (MLC) algorithms to generate the seven LULC maps. The seven-raster classification maps revealed overall accuracies exceeding 85%, and overall Kappa statistics range between 87 and 95%. The results indicate that the increment in the built-up area was dominant over the last 3 decades.


Introduction
The Kingdom of Bahrain's rapid development boom during the 1980s raised the demand for land, particularly in coastal areas, over the past 34 years . As a result, the pressure on the country's overall area to develop urban centers and re-design the existing cities along the shore had intensified due to this growing land demand. Therefore, different information and statistics on LULC are required for sustainable land resource management and policymaking processes to monitor continually, model, and update changes.
Remote sensing science and techniques play an essential role in producing spatial information and detailed raw data on terrestrial phenomena using space-borne platforms (Jansen and Gregorio 2004). These data on different geographical locations facilitate the characterization and evaluation of the changes (Yuan et al. 2005). Remote sensing data's spatial, temporal, and spectral characteristics are utilized to map LULC to aid in land resource management decision-making (Berlanga-Robles and Ruiz-Luna 2002). The emergence of satellite imagery has developed tools for systematic observation of land cover from space (Mollicone et al. 2003). Data required for large land areas are collected at different time intervals and used to monitor changes on the earth's surface (Jensen 2005). Landsat data help assess LULC for different years, analyze changes, understand landuse patterns and factors, and build databases for long-term monitoring tools when combined with geographic information system data and various statistical analytical tools (Jensen 2005). GIS and Remote sensings are effective techniques for obtaining precise and timely data and information on the spatial distribution of LULC changes over broad areas (Guerschman et al. 2003;Rogan and Chen 2004;Zsuzsanna et al. 2005). Remote sensing and GIS data have been used to map land uses and land cover and detect changes in many studies over the past fifty years. These data have been used to study expansion and urbanization, quantify land changes, and study the effects of land changes in faster, straightforward ways than traditional methods in survey studies (Da Costa and Cintra 1999). GIS provides a versatile framework for gathering, storing, displaying, and evaluating digital data needed for change detection (Demers 2008;Wu et al. 2006). GIS relies heavily on remote sensing imagery as a data source. Satellite imagery is employed to recognize synoptic data of the earth's surface (Ulbricht and Heckendorff 1998). Over the last 20 years, researchers concentrated on using remotely sensed image band characterization to study urban and nonurban area changes (Gadal and Ouerghemmi 2019; Vigneshwaran and Kumar 2018;Jacquin et al. 2008;Xu 2007;Zha et al. 2003). Due to their ability to provide immediate and synoptic usage of land cover, remotely sensed images are suitable tools for identifying, managing, and monitoring urban built-up areas, spatial distribution, and expansion (Hegazy and Kaloop 2015;Rawat and Kumar 2015). Many types of research used different classification techniques such as supervised, unsupervised, object-based, or deep learning classification in their studies for extracting urban areas from a multispectral satellite image (Ghosh and Siddique 2018;Bramhe et al. 2018;Forget et al. 2017;Hegazy and Kaloop 2015;Rawat and Kumar 2015;Zhang et al. 2014;Ndehedehe et al. 2013). With a spatial accuracy of 30 m, Landsat multispectral images became the backbone of large-scale land cover change studies due to the long period covered by satellite data over 45 years (Duan et al. 2020;Tewabe and Fentahun 2020;Lu et al. 2010;Wulder et al. 2008;Cohen and Goward 2004). Many satellite images are being used to monitor and study land cover changes, such as SPOT, ASTER, IR, IKONOS, QuickBird, WorldView1, 2, 3 (Allen et al. 2013;Lu and Weng 2009). The land cover information extracted from remote sensing data is an essential part of various applications, including land-use mapping and monitoring of variables. Different satellite-based LC data applications have been created in the last 30 years. The most important of which is the analysis of land cover dynamics that significantly impacts natural habitats (Yuan et al. 2005), as well as the creation of LC maps utilizing a variety of classification methods, all of which are linked to direct statistics to determine classification accuracy and error ratio. Land cover analysis and classification results are based on four main characteristics: temporal, spatial, radiometric, and spectral accuracy (Allen et al. 2013). Macleod and Congalton (1998) identified four manifestations of detection of land-use changes in the study of land resources, namely, the detection of changes that have occurred, the determination of nature of change, the measurement of the areal extent of change, and the assessment of the spatial pattern of change (Macleod and Congalton 1998). Landsat Multispectral Scanner (MSS), Thematic Mapper (TM), and Enhanced Thematic Mapper Plus (ETM+) data have been broadly employed in studies towards the determination of land cover since 1972, the starting year of the Landsat program, mainly in forest and agricultural areas (Campbell 2007). The availability of detailed data on land uses for urbanization is essential for planning processes (Jensen and Cowen 1999). Remote sensing data have become the cornerstone of all urban studies emphasizing mapping and analyzing the changes in area, extent, and patterns . Open access to earth observation satellite data provides spatially valid datasets over broad areas with great spatial detail and temporal frequency (Xiao et al. 2006). The availability of many free satellite images and GIS layers and improvements in remotesensing data gathering with increased spatial accuracy allow and improve quantitative studies of the rate and pattern of urban LULC change (Epstein et al. 2002).

Study Area
The Kingdom of Bahrain is an archipelago of 36 lowlying islands that vary in size, surrounded by shallow waters no deeper than 20 m. The country is located in the southern part of the Arabian Gulf. Its territorial water covers roughly 9200 km 2 and accounts for more than 90% of the total area. In 2015, the country's land area was estimated to be 782.4 km 2 (2020), including rocky coral atolls (Fashts), as shown in Fig. 1. The coastal length is more than 537 km, and the marine area is more than 9200 km (Information and eGovernment Authority 2020). The country has significant population and infrastructure growth, which has raised the demand for land (UNDP 2018). As a result, the country's land area has increased from 697 km 2 in 1987 to 782 km 2 in 2019 (Information and eGovernment Authority 2020). Due to the increased demand for land, the increment rate reached its maximum Published in partnership with CECCR at King Abdulaziz University point (21 km 2 /year) between 1997 and 2007 (UNDP 2018). The majority of the country's population (more than 95%) lives in densely populated coastal areas (Mesaiqer and Al-Zayani 2008). From the 1950s until 2007, the population increased by 2.7 percent per year, reaching around 435,654 in 1987 and over 544,366 in 1994. In 2013, the estimated population was 1,253,191, which was expected to grow to 1,523,084 by 2020 (CIO (1986(CIO ( -2013; Information and eGovernment Authority 2020). According to CIO's (1986CIO's ( -2013 predictions, the country's population will reach 2.2 million by 2030 (Information and eGovernment Authority 2018; UNDP 2018). Increased demand for land, notably for urbanization and creating new cities and residential facilities, has increased population. More than half of Bahrain's main island is almost entirely urbanized. Aside from considerable pressure on infrastructure, health, education, and social services, the country is vulnerable to some issues due to its arid desert climate and small land area (UNDP 2018; UNdata 2021).

Data Acquisition
In this analysis, we used 30 m dry season six Landsat cloudfree datasets: Landsat-5, Thematic Mapper (TM) from August 14, 1986, Landsat-7 from July 11, 2000, Enhanced Thematic Mapper Plus (ETM+) from July 25, 2005, and Landsat-8 Operational Land Imager (OLI) from August 24, 2013, and August 11, 2020. Table 1 shows the specifics of Landsat satellite image data downloaded from the Earth Explorer website of the United States Geological Survey (USGS) (https:// earth explo rer. usgs. gov/). Earth Explorer is an online website that allows searching, viewing, exporting information, and downloading earth science data from the USGS (EROS 2021). In addition, the images were spatially registered on the World Geodetic System (WGS84) and then re-projected the data to the Northern Hampshire zone 39N on Universal Transverse Mercator (UTM) projection. We used the Maximum Likelihood Classifier (MLC) algorithms to classify the images. All the maps were assigned to the WGS 1984 UTM Global Geodetic Coordinates System. We used the 2020 Landsat 8 OLI image as the master reference for all the images and data used in this study.

Ancillary Data
Ancillary data, also known as auxiliary data, are frequently utilized to improve image classification and accuracy. Due to the limitations of spectral information for defining change, researchers must rely on auxiliary data (Franklin 1995).
To determine LULC Changes processes, many researchers used ancillary data as input features with multi-temporal images (Hurskainen et al. 2019;Feng et al. 2018;Zhu et al. 2016;Zoungrana et al. 2015;Yu et al. 2013 1994, and ASTER 2012 to enhance the LULCclassification. We used auxiliary data, mainly Bahrain topographic scanned maps from the "Survey and Land Registration Bureau (SLRB)" scale 1: 50,000, 25,000, and 10,000, to adjust the decision of LULC-classes, Fig. 2

LULC Classes Definition
The classification of LULC is based mainly on its creation, with the levels being divided into numerous levels based on the level of accuracy and purpose necessary. We adapted the LULC level-1 classification from Anderson et al. (1976). Most of these categories were based on statistical data, aerial photos, satellite imagery data, or directly based on Anderson et al. (1976) or other land classifications. We created a classification scheme for this study based on prior knowledge of the study region and information derived from auxiliary data and information (Anderson et al. 1976). The accessible digital Bahrain Zoning maps (UPDA 2022) and other auxiliary data were the primary sources for identifying the research area's main LULC classes. As a result, seven LULC classes covering the surface study area were identified and classified, including the reclaimed area (Rec), vegetation (Veg), built-up/urban area (Bul), rock outcrop (Roc), bare ground (Bar), gypsic soil (Gyp), and wetland/sabkhas region (Wetland/Sabkhas) (Wet). Table 2 displays the seven LULC classification schemes adapted for this study.

Image Preprocessing
This study used multitemporal and multispectral Landsat satellite images acquired in 1986, 1994, 2000, 2005, 2013, and 2020. The changes in the LULC areas were detected using a Supervised Maximum Likelihood Classification algorithm for the six temporal dates. Remote sensing and GIS Software were used to handle the geometric and radiometric enhancement, mosaicking, and the sub-setting of  1994,1998 were used to enhance the LULC-classification. Imagepreprocessing operations for this study include geometric correction (rectification), ETM+ Scan Line Corrector-SLC off and gap filling, radiometric, and atmospheric calibration (Lu et al. 2004). Geometric (rectification), ETM+ Scan Line Corrector-SLC off and gap filling, radiometric, and atmospheric calibration are a few of the image-preprocessing techniques used in this study (Lu et al. 2004). All the images were co-registered and radiometrically calibrated to compare different temporal dates of the study area. Radiometric calibration was done to reduce or correct errors in the digital numbers of the images. The calibrated Digital Numbers (DNs) were first transformed to absolute units of at-sensor spectral brightness and then to Top-Of-Atmosphere (TOA) reflectance (Chander et al. 2009). The DNs stored in the original image are converted into biophysical variables of standard significance (reflectance). For Landsat 5, 7, and 8, the TOA spectral radiance is calculated using the band-specific multiplicative rescaling factor. The TOA reflectance is then adjusted for the solar angle. TOA reflectance has three advantages over at-sensor spectral radiance when comparing images from different sensors. First, the time difference between data collection reduces the cosine effect of differing solar zenith angles. Second, TOA reflectance adjusts for spectral band changes in exo-atmospheric solar irradiance, resulting in varied values of exo-atmospheric solar irradiance. Third, the TOA reflectance compensates for differences in Earth-Sun distance between data gathering dates. These variances might be significant geographically and temporally (Chander et al. 2009).

Image Classification
We employed the supervised classification methodology to classify the LULC categories (sites) into the predetermined seven classes in this study, Table 2. Combining training samples (sites) with prior knowledge and familiarity with the study area is necessary to complete the supervised classification (Jensen 2005). The training sites were chosen based on visible areas in all the image sources. We began by defining polygons around representative portions of each LULC category (Sites) in each Landsat image as identified and defined Bar Areas with exposed soils and un-vegetated land generally contain thin soil, bare soil, sand, or rocks Gypsic soil Gyp Areas covered with Gypsisols soil are a fine white powder, crystals, pebbles, and stones formed through dissolution from calcium sulfate Wetlands/sabkhas Wet In-land and near shallow coasts, marshes, mudflats, natural sabkhas, and the in-land are formed by dredging and filling 1 3 Published in partnership with CECCR at King Abdulaziz University in Table 2. Between 200 and 500 training sites were identified in each rectified satellite image. We employed the maximum likelihood classifier-MLC to produce spectral signatures to classify all the pixels in the image after the training sites had been established. Following the first classification, a majority filter was used to smooth the classification results by reducing noise from the classed raster maps, resulting in the final LULC raster maps utilized for further analysis. The 1:10,000 Bahrain topographic maps were utilized to construct ground signatures for the supervised categorization of the 1986 image. For the 1994 supervised classification, ground control signatures were created using MrSID orthophotos from 1994 and Bahrain topographic maps at a scale of 1:10,000.  Table 2. We then supervised all the images using the maximum MLC algorithm (Lillesand and Kiefer 1999).

Accuracy Assessment
The error matrix is the most used technique for determining accuracy (Congalton and Green 2019). The classification accuracy can be evaluated to produce an overall measure of the map's quality, which can then be used to compare alternative change detection systems (Foody 2002). The minimum level of accuracy in interpreting remote sensing data to identify LULC classes should be at least 85% (Anderson et al. 1976). These standards error matrices (accuracy assessment statistics) are computed based on the same data references for each image to calculate the components of overall accuracy, user's accuracy, producer's accuracy, error of commission (EC), error of omission (EO), and kappa coefficient. Overall accuracy (OA) is the total number of successes compared to the total number of samples in the categorized image. It is calculated by summing the number of correctly classified values and dividing it by the total number of values in the confusion matrix in Eq. (1). User's accuracy (UA) is the probability of classified pixel on each map representing the actual class on the ground or real-world location (Congalton 1991;Jensen 2005;Campbell 2007) and is calculated using Eq.
(2). On the other hand, the producer's accuracy (PA) measures the error of omission. It is the probability that a reference pixel is classified correctly. It is calculated by dividing the number of corrected classified samples of a specific category by the total number of reference samples using Eq.
(3). The error of commission (EC) is the proportion of a pixel that is predicted to be in a class, but it does not. It is calculated using Eq. (4). The error of omission (EO) is the proportion of observed pixels on the ground that are not classified on the map, and it corresponds to the producer's accuracy. The kappa coefficient (KC) is a measure of the difference between the actual agreement between reference data and an automated classifier and the chance agreement between the reference data and a random classifier, expressed using Eq. (5) (Lillesand et al. 2015). The percentage of correctly categorized pixels is calculated from the percentage expected by chance. It measures the difference between the actual agreement and chance (random) agreement between the map and the validation data on the ground (Congalton 2001). The higher the classification accuracy of the map, the more valuable it is for land administrators and land-use planners. The KC (Koc et al. 2012) is a discrete multivariate approach for determining the level of agreement or accuracy. Its value ranges from −1 to 1. However, it frequently falls between 0 and 1 (Zanotta 2018): where r is the rows number in the matrix, xii is the number of observations in row i and column i (the diagonal elements), x +i and x i+ are the marginal totals of row i and column j, respectively, and N is the observations' number. This study used the error matrix (confusion matrix), the most efficient accuracy assessment method for the six-time period's MLC raster classified maps: 1986, 1994, 2000, 2005(Roy and Inamdar 2019. We used this step to quantitatively assess how efficiently the classified remotely sensed data pixels were sampled into corrected LULC's category. The resulting classified images Error of Omission (EO) = ∑ of f Diagonal element of Column Column Total × 100 were verified pixel by pixel for the accuracy assessment (Keshtkar et al. 2017). The seven classes, Table 2, were considered for accuracy assessment with a minimum of 71 sample points for each defined class, as recommended by Congalton (1991). We used the "Equalized stratified random" sampling method to create randomly distributed points, where each class has the same number of points. Four hundred and ninety-seven (497) systematic point samples for each classified image were used to test the precision by comparing them with specific points from the available ancillary data for this study. We calculated the confusion error matrix considering an accuracy level equal to or above 85% to be an excellent and reliable LULC classification as recommended by previous studies (Alrababah and Alhamad 2006;Manandhar et al. 2009;Koc et al. 2012). We used the confusion matrix for obtaining descriptive and analytical statistics of the classification accuracy assessment, as explained in Sup Tables 1-6 (Foody, 2002;Congalton, 1991;Jensen, 2005).

Change Detection
Change detection aims to compare the spatial representation of two points in time while controlling all variations due to differences in the variables of interest (Green et al. 1994).
The ability to identify changes in the earth's surface features in real-time and with high accuracy lays the groundwork for a greater understanding of the linkages and interactions between human and natural events, allowing for improved resource management and use (Lu et al. 2004). The most common data type used to detect changes is geographic data, usually in digital forms such as satellite imagery, analog format (earlier aerial images), and vector format (maps).
Other data types, such as historical and economic data, can be employed (Singh 1989). Due to the availability of massive archival data sets in recent decades, several digital change detection algorithms and methodologies for evaluating and identifying LULC changes have been developed and evaluated (Dewidar 2004). These methods and procedures have been thoroughly examined, with great descriptions and summaries supplied (Haque and Basak 2017). Since digital change detection is significantly influenced by the temporal, spatial, spectral, and thematic resolutions of remotely sensed data, choosing the right change detection approach is critical for producing accurate findings (Lu et al. 2004). Image differencing, image rationing, PCA, CVA, and Post-Classification Comparison are the most common approaches for detecting changes (Xu et al. 2009;Bekalo 2009 (Talukdar et al. 2020) and comparing machine-learning algorithms (Camargo et al. 2019). In addition, some research has been conducted to determine the most appropriate and accurate algorithm for LULC mapping among the many machine-learning classifiers (Camargo et al. 2019;Jamali 2019). In this study, we used the GIS-based change detection approach, which integrates GIS and remote sensing method (Lu et al. 2004) to explore and identify the changes that have taken place to the spatial extent and pattern of our study area (Gallego 2004). We used GIS because of its ability to incorporate a different source of data with different data accuracies and formats (Petit and Lambin 2001) for long-period intervals into LULC change detection (Lu et al. 2004). This approach helps analyze the direction, rate, and spatial pattern of LULC changes (Weng 2002). The magnitude of change (MC), the percentage of change (PC), and the annual rate of change (ARC) of the classified images were calculated based on the following equations: where A i is the class area (km 2 ) at the initial time, A f is the class area (km 2 ) at the final time, and (n) is the number of years of the study period.
There are different techniques for change detection analysis, and the most common and used ones include post-classification comparison, image ratio, and manual on-screen digitization of change, principal components analysis, image regression, conventional image differentiation, and multidate image classification (Ayele et al. 2018;Lu et al. 2004).
Published in partnership with CECCR at King Abdulaziz University We used the post-classification comparison techniques to integrate the six classified images, which were classified individually to produce a land cover raster map, Fig. 3a-f, Tables 3 and 4. Then, we compared the corresponding classes to identify areas where change has occurred (El-Hattab 2016). We build a series of "from-to" matrixes by comparing them on a pixel-by-pixel basis (Jensen 2005;Aldoski et al. 2013). It involves the overlay (or "stacking") of two or more classified images. Change areas are simply those not classified the same at different times. The multidate change detection uses a binary mask applied to two dates (Jensen 1996;Dobson et al. 1995) and quantifies the change rates and magnitude (Ukor et al. 2016). The degree of change detection success depends on the accuracy of each image classification (Jensen 2005 (1986,1994,2000,2005,2013, and 2020) were compared by applying Boolean logical "AND" operation using GIS Software. Three matrices were developed to understand the changing status and magnitude rate for 1986 to 2005, 2005 to 2020, and 1986 to 2020, Table 5. In addition, six chord diagrams were created using Microsoft Power BI to visualize the classified raster maps classes ( Fig. 4a-f), and three to visualize the changes in the specified temporal dates (Fig. 5a-c). In addition, the data presented in Table 5 was used to create three-chord diagrams to visualize the changes for the specified years in this study : 1986-2005, 2005-2020, and 1986-2020 (Fig. 5a-c). These diagrams represent the changes in the seven LULC-classes in the study area: reclamation, vegetation, built-up, rock outcrop, bare soil, Gypsic soil, and wetland from the "initial year" to the "first year." The data in the three tables were reconstructed and added to the program to fit the parameters' criteria. The Chord visual offers three fields (From, To, and Values). The data from the initial year go into the first field (From), the data from the first year (year of change) go into the second field (To), and the change between classes goes into the third field (Values). Finally, the Chord diagram is displayed according to the parameters selected visually (Data Visualization | Microsoft Power BI, n.d.; Ferrari and Russo 2016). Different researchers have employed Chord-Diagram, as a graphical approach for visualizing the inter-relationships between data

Accuracy Assessment
We used six confusion matrices to assess the classification accuracy in this study. Sup Tables 1-6 indicate that the image classification achieved satisfying accuracy results. The OA of the classified images was 88.73, 89.54, 89.13, 88.93, 90.34, and 95.77% for the temporal dates 1986, 1994, 2000, 2005, 2013, and 2020, respectively, thus ascertaining that the classification was within the excellent range according to Anderson et al. (1976), as he indicated that the minimum accuracy value for reliable LC-classification is 85%. The KA reached 86.85, 87.79, 87.32, 87.09, 88.73, and 95.07% for 1986KA reached 86.85, 87.79, 87.32, 87.09, 88.73, and 95.07% for , 1994KA reached 86.85, 87.79, 87.32, 87.09, 88.73, and 95.07% for , 2000KA reached 86.85, 87.79, 87.32, 87.09, 88.73, and 95.07% for , 2005KA reached 86.85, 87.79, 87.32, 87.09, 88.73, and 95.07% for , 2013, and 2020, individually. The results of the seven classes for all the temporal times have excellent individual user's accuracy ranges between 80% and more than 94% for 1986. The user's accuracy for all the classes is outstanding except for the "Rock outcrop," about 80.28%, and "Gypsic soil," which is 81.69%. In 1994, the UA ranged from 73% to more than 97%, and the results indicate that the "Vegetated" surfaces were most accurately classified (97.18). The less accurate class is "Reclaimed" land, 73.24%. In 2000, the UA ranged between 70 and 100%. The 2000 classified image results indicate that the "Vegetated" and "Rock outcrop" surfaces were most accurately classified (95.77%). The less accurate class is Gypsic, which is 83.1%. The 2005 classified image results indicate that the "Vegetation" surface revealed that UA reached 95.55%, while all the other classes of UA ranged between 83.1% and 87.32%. The UA of 2013 and 2020 Landsat8 OLI-classified images ranged between 91.55% and 98.59% and 87.32 and 100%, respectively. The producer's accuracy (PA) has a nearly excellent result for all class-classified images. In 1986, the PA ranged between 76.47 and 100%. In addition, the results indicate that the "Reclaimed area" surfaces were most accurately classified (100%). The less accurate class is "Bare Land" land, 76.47%. In 1994, the PA ranged between 68% and more than 98.39%. Furthermore, the results indicate that the "Gypsic soils" surfaces were most accurately classified (98.39). The less accurate class is the "Bare land" area, 68%. The 2000 classified image results indicate that the "Rock outcrop." Moreover, "Gypsic soils" surfaces were most accurately classified (100%). The less accurate class is the built-up area, which is 77.11%. The 2005 classified image results indicate that the "Gypsic soils" surface revealed a PA up to 96.72%, while all the other classes of PA ranged between 77.78 and 95.24%. The UA of 2013 and 2020 Landsat8 OLI-classified images ranged between 93.15 and 100% and 92.65 and 100%, respectively.

LULC Classified Maps
The six classified raster maps were used to reach the objectives of LULC change detection in this study. Table 3, Fig. 3a-f, and chord diagram Fig. 4a-f explicate the result of the six classified raster maps of 1986, 1994, 2000, 2005, 2013, and 2020 in detail. The chord diagram represents flows (connections) between the seven classes of LULC (reclamation, vegetation, built-up, rock outcrop, bare soil, gypsic soil, and wetland/sabkhas). A fragment on the outer part of the circular layout represents each class in the chord diagram, and then arcs are drawn between each class. The size of the arc is proportional to the significance of the change. The results of the 1986 LULC classified image Fig. 3a shows that Bare soil dominates a total area (41.4%), followed by wetland and sabkhas (15.9%), as illustrated in Fig. 4a. The built-up area covered 12.6% of the total area. More than 78 km 2 of vegetation covered the ground, accounting for 11.2% of the total land surface. Bare soil class accounts for 28.4% (208.1 km 2 ) of the entire area. Nonetheless, the built-up area covers 19.6% (143.2 km 2 ) of the land. Vegetation land covers 8.8% of the total land area (64.6 km 2 ). The smallest class is the gypsic soils, which make up 5.9% of the total land area (42.9 km 2 ). In the 2013 classified image results, Figs. 3e and 4e, the land area was covered by 28.6% wetland and sabkhas (220.8 km 2 ), followed by built-up land (190.9 km 2 , 24.8% of the total area). The bare soil covered 21.2% or 163.7 km 2 . The reclaimed area covered 6.9% or 52.8 km 2 of the total area. The LULC classes for 2020, Fig. 3f, indicate that the built-up land and bare soil area almost occupied the same total area of 312.2 and 301.1 km 2 , which account for 39.8% and 38.4% of the total area, respectively, Fig. 4f. On the other hand, the gypsic soils occupied the minor area this year and reached 0.8% or 6.3 km 2 of the total area.

LULC Change Detection
The spatial analysis was conducted to assess the patterns of LC-change and overall LU-changes from 1986 to 2020. The LULC-classes of the study area were reclaimed, vegetation, built-up, rock outcrop, bare land, gypsic soils, and wetland/sabkhas. We designated three critical years: the starting year of the study, 1986, the mid-year 2005 (approximately), and 2020. The change detection analysis concentrated on three periodic times by creating three area matrices from 1986-2005 and 2005-2020, and the whole period of the temporal study dates from 1986 to 2020. The changing patterns were depicted using chord diagrams. Table 4 illustrates the LULC classes in 1986, 2005, and 2020 for each class, the total area, and the spatio-temporal changes in the study area for the last 34 years. In addition, three chords' diagrams, Fig. 5a-c presents the changes from 1986-2005, 2005-2020, and 1986-2020 1986-2005, 2005-2020, and 1986-2020. Table 4 shows that four LULC classes have registered a negative annual change from 1986 to 2020. The vegetation, rock outcrop, gypsic soil, and wetland/sabkhas classes have recoded loss annually by −0.70, −0.42, −1.51, and −2.45 km 2 , respectively. The annual rate of change (ARC) of these classes were −1.59, −0.90, −1.26, −2.26, and −2.22%, respectively, for the same period. The built-up and bare soil gain was annually 6.6 and 0.36 km 2 , respectively, and their annual rate of change was 7.49 and 0.12%, respectively. Built-up areas have registered the highest annual growth rate. The wetland/sabkhas registered the lowest annual growth rate in the last 34 years. Table 5 Table 5 Table 5 and Fig. 5c depict the LULC-classes change transitions from 1986 to 2020. Most changes in the classes' category between 1986 and 2020 were directed to built-up areas. The reclaimed areas were converted mostly to builtup and wetlands/sabkhas, and 18.2 and 18 km 2 were added to these classes. The vegetation class lost 43.6 km 2 , which converted mainly to built-up areas. The rock outcrop class area lost 14.1 km 2 , which converted mainly to bare soil in 2005. 60.8 km 2 of the bare soil areas changed to built-up by 2020. In addition, 23.9 and 34.4 km 2 from the gypsic soil wetland/sabkhas were added to built-up areas. 12.3 km 2 of the rock outcrop areas in 1986 converted to built-up areas by 2020. Around 83.4 km 2 of the wetland/sabkhas area in 1986 was lost. The change detection from 1986 to 2020 indicates the shifting of most LULC-classes to the built-up area and the reduction in vegetation and wetland/sabkhas land in the last 34 years. The change detection from 1986 to 2005 indicates the shifting to the built-up area and the reduction in vegetation land.

Conclusions
This study integrated remote sensing and GIS to quantify and analyze the LULC changes in Bahrain's main islands over 34 years from 1986 to 2020. Moreover, LULC changes between 1986, 2005, and 2020 were demonstrated and visualized using statistics and chord diagrams. Furthermore, the identified seven LULC-classes revealed substantial change patterns in the study area. LULC-changes were determined using six multitemporal Landsat satellite imagery, and the classification accuracy was measured using the confusion matrix. The overall classification accuracy was acceptable. Conferring to the quantitative evidence from our study, Bahrain Island has witnessed significant land use and land cover changes since 1986. The study results showed a significant change in the LULC during the study period. During the study period, the built-up areas showed an increasing trend of 7.5% annually, while the Gypsic, wetland/sabkhas, and vegetation showed a decreasing trend of 2.26, 2.22, and 0.9%, respectively. The changes in LULC were effectively captured by the remote sensing Landsat satellite sensors with different spectral, spatial, and temporal resolutions with GIS analysis. The change detection analysis using GIS and remote sensing delivers valuable information to understand the annual patterns of land use dynamics for planners and decision-makers; therefore, sustainable land management planning is possible. The expansion of the built-up area in the study area was mainly at the expense of vegetation and wetland/sabkhas lands. The results indicate that the builtup area covered more than 40% of Bahrain's main islands 1 3 Published in partnership with CECCR at King Abdulaziz University land. For the environment's long-term sustainability, LULC changes should be constantly monitored in the future. Moreover, Landsat TM5, ETM+7, and OLI8 images archives and the newly invented remote sensing satellite imagery should be used to create and monitor accurate maps of LULC changes in Bahrain.

Funding
The authors have not disclosed any funding.

Conflict of interest
The corresponding author states that there is no conflict of interest on behalf of all the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.