Integrating a Three-Level GIS Framework and a Graph Model to Track, Represent, and Analyze the Dynamic Activities of Tidal Flats

: Tidal ﬂats (non-vegetated area) are soft-sediment habitats that are alternately submerged and exposed to the air by changeable tidal levels. The tidal ﬂat dynamics research mainly utilizes the cell-level comparisons between the consecutive snapshots, but the in-depth study requires more detailed information of the dynamic activities. To better track, represent, and analyze tidal ﬂats’ dynamic activities, this study proposes an integrated approach of a three-level Geographic Information Science (GIS) framework and a graph model. In the three-level GIS framework, the adjacent cells are assembled as the objects, and the objects on different time steps are linked as lifecycles by tracking the predecessor–successor relationships. Furthermore, eleven events are deﬁned to describe the dynamic activities throughout the lifecycles. The graph model provides a better way to represent the lifecycles, and graph operators are utilized to facilitate the event analysis. The integrated approach is applied to tidal ﬂats’ dynamic activities in the southwest tip of Florida Peninsula from 1984 to 2018. The results suggest that the integrated approach provides an effective way to track, represent, and analyze the dynamic activities of tidal ﬂats, and it offers a novel perspective to examine other dynamic geographic phenomena with large spatiotemporal


Introduction
The study of dynamic geographic phenomena is a challenging topic in Geographic Information Science (GIS). In this topic, the core issue is the spatiotemporal characteristics of the study subjects, and the prerequisite is to identify and delineate the study subjects from the source data. Compared with the delineation, a more critical problem is how to track better, represent, and analyze the dynamic activities from the delineated results, which can help obtain the relevant information to fully describe the spatiotemporal characteristics of dynamic geographic phenomena and explain the reasons behind them.
In the GIS community, modeling the dynamic activities of geographic phenomena has been discussed, developed, and refined during the past several decades. As summarized by Worboys (2005) [1], this procedure can be divided into three stages: (1) the temporal snapshots, (2) the object changes, and (3) the events and action. As the delineated results can represent the static state of the dynamic geographic phenomena at one single moment, the first stage makes it possible to find the temporal sequences of objects, their attributes, and relationships by associating the states at different moments along the time sequence. For instance, Armstrong (1988) [2] proposes the snapshot view, which allows the query for the temporal information of a single object and all items associated with a given time. In addition, Worboys (2005) [1] points out that the time domain should allow interpolation between measurements in case of the continuous movements. However, the researchers often feel frustrated when considering the changes, because the first stage focuses on the temporal sequences of snapshots only [3]. In response to the new challenges, the second stage shifts the focus to the changes retrieved from the series of temporal snapshots. Hornsby and Egenhofer (2000) [4] emphasize that the change identification should be based on the comparison between the subsequent snapshots. As one of the earliest attempts, Peuquet and Duan (1995) [5] propose a data model that explicitly represents the changes based on the temporal comparisons in sequence. On the other hand, Claramunt and Thériault (1995) [6] propose a temporal framework focused on the geographic entities' topological changes, which provides a novel perspective to examine the complex evolution of dynamic geographic phenomena. However, the second stage's changes can only contribute to the preliminary representation and analysis for the dynamic geographic phenomena, and therefore, the third stage is developed in response to the new challenges. As a breakthrough, Yuan (2001) [7] systematically organizes a framework in which a temporal sequence of states forms a process, and a temporal aggregation of the processes makes an event.
In the third stage, Worboys (2005) [1] highlights "the continuants that endure through time" and "the occurrents that happen or occur and are then gone". This perception further deepens and expands the concept of objects in the second stage, and it yields the eventbased perspective, which provides a possible way to better model the dynamic geographic phenomena. The objects can be tracked as lifecycles through which the dynamic activities can be observed [8]. McIntosh and Yuan (2005) [9] propose a four-level framework (zone, sequence, process, and event), which facilities the spatiotemporal query and analysis for the distributed dynamic geographic phenomena, and a case study of rainstorms is implemented. In addition, the following studies of marine dunes [10], ocean eddies [11], convective storms [12], urban heat islands [13], and Land Use and Land Cover (LULC) [14] have also endorsed the good performance of the event-based model. However, due to the data availability and computing performance, the previous studies often face the challenges of domain selection and spatiotemporal scale limitation. The in-depth study of event-based model is expected to eliminate these shortcomings with the aid of high-performance computing.
As a following issue, a good representation strategy could not only contribute to the visualization of the lifecycles but also facilitate the further analysis of the dynamic activities. More recently, the graph theory has been applied to GIS studies, which provides a possible solution to better represent the lifecycles. Thibaud et al. (2013) [10] utilizes a spatiotemporal graph model, in which "the graph edges are associated with the concepts of spatial relations and filiation relations". Cheung et al. (2015) [15] suggests that a sequence of snapshots at particular time steps can be threaded together into a single spatiotemporal graph, which provides an enhanced way to visualize the evolutionary trajectories over time. Furthermore, Wu et al. (2016) [16] proves that the graph theory can demonstrate a concise and intuitionistic way to capture and represent the storm-induced coastal changes. As the lifecycles are represented as graphs, the dynamic activities have great potential to get further analyzed with graph operators and algorithms. Previous studies mainly concentrate on the representation and visualization of dynamic geographic phenomena via the graph model, so it needs further study on how to adapt the graph algorithms to fully explore the dynamic activities in different domains.
In this study, we select the tidal flat as the subject of research. It is the coastal sediment that is alternately and periodically submerged in water and exposed to the air by the changeable tidal levels. As the land-sea interactions take place, the temperature, salinity, and acidity in the tidal flat region become variable [17]. A unique type of wetland ecosystem is the product under the above physical and chemical conditions, which is the homeland of shorebird [18], fungus [19], plankton [20], coastal fish [21], and so on. Aside from its contribution to biodiversity, the tidal flat also plays an important role in carbon preservation and global warming prevention [22][23][24][25][26]. In addition, the tidal flat has economic benefits. The unvegetated tidal flat all over the world contributes a total of $2.44 trillion USD per year as of the value in the year 2011 [27]. Hence, the dynamic activities of tidal flats have attracted the attention of researchers from different backgrounds.
It is essential to clarify the scope of the tidal flat region in this study because it can be defined in two ways. In a broad sense, the tidal flat region includes the bare intertidal flats, It is essential to clarify the scope of the tidal flat region in this study because it can be defined in two ways. In a broad sense, the tidal flat region includes the bare intertidal flats, as well as the vegetated flats covered by salt marshes, mangroves, and seagrasses [28,29]. In a narrow sense, the tidal flat region only includes the unvegetated intertidal flats [30,31]. This study explicitly aims at the narrow sense of tidal flats ( Figure 1). Following the existed theoretical basis, it is expected to propose a novel approach that facilitates the acquisition and summarization for the spatiotemporal characteristics of tidal flat dynamic activities. More specifically, the goal of this study is to integrate a three-level GIS framework and graph model. The three levels are as follows: (1) the cell level, which is the satellite image pixel of the tidal flat; (2) the object level, which is the set of contiguous tidal flat cells with a specific state at a specific moment; and (3) the lifecycle level, which is the set of tidal flat objects at different time steps linked as a chain based on the relative position relationships. In addition, different events to depict the dynamic activities can be derived from the lifecycle level. Therefore, this three-level GIS framework should be capable of tidal flat lifecycle tracking, which is the foundation of events capture. On the other hand, the graph model can convert the spatiotemporal information in the lifecycle level to the directed graphs, and its derivative spatiotemporal database offers the storage and query features for the morphologic attributes of the tidal flats throughout each individual lifecycle. Consequently, the graph can provide an enhanced way to represent the lifecycles and analyze the dynamic activities of tidal flats.
A case study is applied to the tidal flats in the southwest tip of the Florida peninsula from 1984 to 2018. Owing to the development of high-performance cloud computing in recent years, it has become feasible to access, process, and analyze geospatial big data, and therefore, our integrated approach can be applied to the larger spatiotemporal scales. With the available big datasets, the machine learning algorithms are widely used in cloud computing platforms. The random forest (RF) machine learning algorithm is used to delineate the cell level tidal flat information in this study.

Tidal Flat Cell Delineation
To study the spatiotemporal characteristics, it is essential to delineate the cells of the dynamic phenomena from the satellite data [11]. As different types of LULC have different spectral characteristics, the remote sensing-derived spectral indices are commonly used for water body identification, which is the prerequisite for tidal flat cell delineation. These indices include the Normalized Difference Water Index (NDWI) [32], Land Surface Water Index (LSWI) [33], Modified Normalized Difference Water Index (MNDWI) [34], Automated Water Extraction Index (AWEI) [35], and so on.
Based on water body identification results, many efforts have been spent on delineating tidal flats from the satellite images. The conventional methods include linear re- Following the existed theoretical basis, it is expected to propose a novel approach that facilitates the acquisition and summarization for the spatiotemporal characteristics of tidal flat dynamic activities. More specifically, the goal of this study is to integrate a three-level GIS framework and graph model. The three levels are as follows: (1) the cell level, which is the satellite image pixel of the tidal flat; (2) the object level, which is the set of contiguous tidal flat cells with a specific state at a specific moment; and (3) the lifecycle level, which is the set of tidal flat objects at different time steps linked as a chain based on the relative position relationships. In addition, different events to depict the dynamic activities can be derived from the lifecycle level. Therefore, this three-level GIS framework should be capable of tidal flat lifecycle tracking, which is the foundation of events capture. On the other hand, the graph model can convert the spatiotemporal information in the lifecycle level to the directed graphs, and its derivative spatiotemporal database offers the storage and query features for the morphologic attributes of the tidal flats throughout each individual lifecycle. Consequently, the graph can provide an enhanced way to represent the lifecycles and analyze the dynamic activities of tidal flats.
A case study is applied to the tidal flats in the southwest tip of the Florida peninsula from 1984 to 2018. Owing to the development of high-performance cloud computing in recent years, it has become feasible to access, process, and analyze geospatial big data, and therefore, our integrated approach can be applied to the larger spatiotemporal scales. With the available big datasets, the machine learning algorithms are widely used in cloud computing platforms. The random forest (RF) machine learning algorithm is used to delineate the cell level tidal flat information in this study.

Tidal Flat Cell Delineation
To study the spatiotemporal characteristics, it is essential to delineate the cells of the dynamic phenomena from the satellite data [11]. As different types of LULC have different spectral characteristics, the remote sensing-derived spectral indices are commonly used for water body identification, which is the prerequisite for tidal flat cell delineation. These indices include the Normalized Difference Water Index (NDWI) [32], Land Surface Water Index (LSWI) [33], Modified Normalized Difference Water Index (MNDWI) [34], Automated Water Extraction Index (AWEI) [35], and so on.
Based on water body identification results, many efforts have been spent on delineating tidal flats from the satellite images. The conventional methods include linear regression [36], Gaussian function [37], Otsu thresholding [38,39], and frequency thresh-olding [29,40,41]. Compared with the above conventional methods, the RF could classify the cells with respect to multiple indices and their statistical distributions and provide higher robustness [42]. To date, the RF is widely used in LULC classification issues, such as cropland [43,44], woody vegetation [45], human settlement [46,47], coastal lands [48], and so on. For tidal flat cell delineation, several methods have been developed based on the RF. Zhang et al. (2019)'s RF-based method [28] used MNDWI, Normalized Difference Vegetation Index (NDVI) [49], LSWI, Enhanced Vegetation Index (EVI) [50], Modified Soil-Adjusted Vegetation Index (MSAVI) [51], and Soil Brightness [52]. Their methods extract both the intertidal flats and supratidal vegetated flats, but we focus on the unvegetated intertidal zone. Following Murray et al. (2019) [31], the indices used in this study include AWEI, NDWI, MNDWI, and NDVI. The four indices are expressed as Equations (1)-(4): (1) where ρ blue , ρ green , ρ red , ρ nir , ρ swir1 , and ρ swir2 are the surface reflectance values of the corresponding bands of the Landsat images. In addition to the above indices, the Landsat near-infrared band (NIR) and Landsat short-wave infrared (SWIR) band are also considered.
The Landsat images are categorized as stacks consisting of all images acquired within each year. For any cell with 30 meters' spatial resolution, we can get the results by computing AWEI, NDWI, and MNDWI from all cells at the same location on different images within the same stack. To get the spectral change pattern in a stack, we also compute the maximum, minimum, median, standard deviation, the 10th, 25th, 50th, 75th, and 90th percentiles, and the means of all inputs in the specified percentile ranges (0-10, 10-25, 10-90, 25-50, 25-75, 50-75, 75-90, and 90-100). However, for NDVI, NIR, and SWIR, we only calculate the means of interval between the 10th and 90th percentiles. At this point, there are 54 variables for each cell.
In addition, another two variables are used. (1) The first is the ETOPO1 Global Relief Model, which is developed by the National Oceanic and Atmospheric Administration (NOAA) of the US. It is a global relief model of the Earth's surface, which integrates land topography and ocean bathymetry [53]. (2) The second is the global surface water occurrence data. This dataset has generated the location and temporal distribution of surface water, as well as the statistics on the extent and change of these water surfaces with 30 meters' spatial resolution based on a total of 3,066,102 tiles of Landsat 5, 7, and 8 images, which were acquired between 16 March 1984 and 10 October 2015 [54]. These two variables combining with the 54 statistical results for each cell form 56 predictor variables. A pre-computed and pre-classified (tidal flat, permanent water, and other) sample point dataset is used as a training dataset for the RF machine learning algorithm.
After the classification and post-processing, we can get the binary images, in which the cells with the value of 1 correspond to the tidal flats and the cells with the value of 0 correspond to the non-tidal flats. To validate the result, we use the validation dataset via random stratified sampling [31]. As the tidal flat cell delineation is accomplished, the foundation for tidal flat object level has been laid.

Tidal Flat Object Delineation
In the two-dimensional space, one cell may have up to eight neighboring cells (up, down, left, right, upper-left, upper-right, lower-left, and lower-right) [55]. For the tidal flat ISPRS Int. J. Geo-Inf. 2021, 10, 61 5 of 23 cells satisfying this eight-neighboring definition, the component-labeling algorithm [56] can assemble them as a contiguous patch, which is a tidal flat object in the context of this study. In addition, we can define a threshold to discard those too small tidal flat objects that may potentially add unnecessary complexity to the further analysis [12]. The threshold selection should be with respect to the spatial resolution and the size of the study area.
The objects generated from different domains, such as ocean eddies [11], convective storms [12], and urban heat islands [13], can be tracked and linked as time-series chains based on the relative positions between two objects on different time steps. These previous studies prove that the object modeling is the foundation of the lifecycle tracking. Specific to this study, the object delineation plays as the foundation of capturing the morphological dynamics of tidal flats, which is invisible in the conventional (cell-level) assessment.

Tidal Flat Lifecycle Tracking
For the tidal flat lifecycle tracking, the strategy should always be determined by the characteristics of the specific subject of research [11]. For the tidal flat lifecycle tracking, the overlapping method is realized, which finds the overlapping region between the two objects at the two adjacent time steps, considering the area ratios of the overlapping region to the two objects; then, it determines whether the two objects can be associated or not. In our study, for a pair of tidal flat objects, α is the object at the previous time step, β is the object at the current time step, and γ is the overlapping region between α and β. The overlapping ratio, R, is formularized as: where A α , A β , and A γ are the areas of α, β, and γ, respectively. The value of R should not be less than 0, which corresponds to no overlapping case, and it should not be greater than 2, which corresponds to the exactly overlapping case. A threshold is defined to determine whether the two objects can get linked or not. As α and β are linked, a predecessorsuccessor relationship is established, in which α is the predecessor and β is the successor. The predecessor-successor relationships links the tidal flat objects between consecutive images, which refers to the lifecycle in this study. For the objects within the derived lifecycle, the numbers of their predecessors and successors may vary, and the areas of each pair of predecessor and successor may be different. To categorize the dynamic activities of tidal flats, we firstly define six events according to the numbers of the predecessors and successors ( Figure 2). In addition, we further define five events based on the area comparison between the predecessor and successor ( Figure 3).
The object(s) on the first, second, and third time step are represented as A(s) (in red), B(s) (in blue), and C (in yellow), respectively. The solid line means the actual-existed object, and the dashed line means the pseudo-object ( Figure 2). The six events according to the numbers of predecessors and successors include:  To demonstrate the tidal flat lifecycle tracking procedure, three consecutive snapshots taken at different time steps (t1, t2, and t3) are provided as an example (Figure 4), which form a three-dimensional space (X, Y, T). There are a total of three tidal flat objects  To demonstrate the tidal flat lifecycle tracking procedure, three consecutive snapshots taken at different time steps (t1, t2, and t3) are provided as an example (Figure 4), which form a three-dimensional space (X, Y, T). There are a total of three tidal flat objects To demonstrate the tidal flat lifecycle tracking procedure, three consecutive snapshots taken at different time steps (t 1 , t 2 , and t 3 ) are provided as an example (Figure 4), which form a three-dimensional space (X, Y, T). There are a total of three tidal flat objects (A 1 , C 1 , and E 1 ) at t 1 , six tidal flat objects (A 2 , B 1 , C 2 , D 1 , E 2 , and E 3 ) at t 2 , and four tidal flat objects (B 2 , B 3 , C 3 , and E 4 ) at t 3 . For better reference, the locations of some tidal flat objects are projected on the adjacent snapshots, which are marked with undertint colors and dashed lines such as A 1 ' of A 1 .  As Figure 4 shows, A1' is the projection of A1 at t2, which shares sufficient overlapping area with A2; therefore, the predecessor-successor relationship is established, where A1 is the predecessor and A2 is the successor, and a continuation happens to them. The other predecessor-successor relationships in this figure are all established in the same way, which contribute a total of four tidal flat lifecycles. All other types of events based on the numbers of predecessors and successors are also visible in this figure: A2' is the projection of A2 at t3, and we can find nothing overlaps with it, so A2 has no successors and a disappearance occurs. Similarly, B0 is the projection of B1 at t1, and we can find nothing overlaps with it, so B1 has no predecessors and an appearance occurs. On the other hand, B1' is the projection of B1 at t3. It overlaps with B2 and B3, which indicates that B1 has two successors, and a splitting comes out. C2' and D1' are the projections of C2 and D1 at t3, and both overlap with C3. It means C3 has two predecessors and a merger happens. In the same way, we could find E2 and E3 are the successors of E1, and they are also the predecessors of E4. These four objects form a recovery.

Graph Representation and Analysis for Tidal Flat Lifecycles
In graph theory, the graph is defined as "a representation of a set of points (the nodes) and of how they are joined up (the edges), and any metrical properties are irrelevant" [57]. Since the graph is a concise way to deliver the information, it has been widely used in GIS to represent the evolutionary procedures of dynamic geographic phenomena. Moreover, the graph operators provide an effective way to capture the events automatically and exhaustively from a large amount of tidal flat lifecycles, which also high- As Figure 4 shows, A 1 ' is the projection of A 1 at t 2 , which shares sufficient overlapping area with A 2 ; therefore, the predecessor-successor relationship is established, where A 1 is the predecessor and A 2 is the successor, and a continuation happens to them. The other predecessor-successor relationships in this figure are all established in the same way, which contribute a total of four tidal flat lifecycles. All other types of events based on the numbers of predecessors and successors are also visible in this figure: A 2 ' is the projection of A 2 at t 3 , and we can find nothing overlaps with it, so A 2 has no successors and a disappearance occurs. Similarly, B 0 is the projection of B 1 at t 1 , and we can find nothing overlaps with it, so B 1 has no predecessors and an appearance occurs. On the other hand, B 1 ' is the projection of B 1 at t 3 . It overlaps with B 2 and B 3 , which indicates that B 1 has two successors, and a splitting comes out. C 2 ' and D 1 ' are the projections of C 2 and D 1 at t 3 , and both overlap with C 3 . It means C 3 has two predecessors and a merger happens. In the same way, we could find E 2 and E 3 are the successors of E 1 , and they are also the predecessors of E 4 . These four objects form a recovery.

Graph Representation and Analysis for Tidal Flat Lifecycles
In graph theory, the graph is defined as "a representation of a set of points (the nodes) and of how they are joined up (the edges), and any metrical properties are irrelevant" [57]. Since the graph is a concise way to deliver the information, it has been widely used in GIS to represent the evolutionary procedures of dynamic geographic phenomena. Moreover, the graph operators provide an effective way to capture the events automatically and exhaustively from a large amount of tidal flat lifecycles, which also highlights that the object and lifecycle levels could provide more diversified dynamic information than the cell level.
In this study, the tidal flat objects are represented by their centroids, and the predecessor-successor relationships are represented by the directed edges from the predecessors to the successors. The centroids form the node set V(G) and the directed edges form the set E(G), which constitute the directed graph G in two-dimensional space (X, Y). To demonstrate how the graph works in this study, the information of Figure 4 is derived as graphs and drawn in Figure 5. In this figure, it is easy to read the predecessor-successor relationships and the corresponding six events based on the numbers of predecessors and successors including appearance, disappearance, splitting, merger, continuation, and recovery.
As the graphs are derived, the information of each individual tidal flat lifecycle is collected, organized, and stored in a spatiotemporal database simultaneously ( Figure 6). This database can help us better understand and query the morphological changes between the adjacent tidal flat objects, which are invisible in Figure 5, because tidal flat objects are represented by their centroids as the nodes. Especially, the five events defined in Section 2.3 are determined by the area comparison between the predecessor and successor, and this spatiotemporal database can help us capture these events. Moreover, for the six events based on the number of predecessors and successors, we introduce two graph operators to assist the analysis for tidal flat dynamic activities. lights that the object and lifecycle levels could provide more diversified dynamic information than the cell level.
In this study, the tidal flat objects are represented by their centroids, and the predecessor-successor relationships are represented by the directed edges from the predecessors to the successors. The centroids form the node set V(G) and the directed edges form the set E(G), which constitute the directed graph G in two-dimensional space (X, Y). To demonstrate how the graph works in this study, the information of Figure 4 is derived as graphs and drawn in Figure 5. In this figure, it is easy to read the predecessor-successor relationships and the corresponding six events based on the numbers of predecessors and successors including appearance, disappearance, splitting, merger, continuation, and recovery.
As the graphs are derived, the information of each individual tidal flat lifecycle is collected, organized, and stored in a spatiotemporal database simultaneously ( Figure 6). This database can help us better understand and query the morphological changes between the adjacent tidal flat objects, which are invisible in Figure 5, because tidal flat objects are represented by their centroids as the nodes. Especially, the five events defined in Section 2.3 are determined by the area comparison between the predecessor and successor, and this spatiotemporal database can help us capture these events. Moreover, for the six events based on the number of predecessors and successors, we introduce two graph operators to assist the analysis for tidal flat dynamic activities.    lights that the object and lifecycle levels could provide more diversified dynamic information than the cell level.
In this study, the tidal flat objects are represented by their centroids, and the predecessor-successor relationships are represented by the directed edges from the predecessors to the successors. The centroids form the node set V(G) and the directed edges form the set E(G), which constitute the directed graph G in two-dimensional space (X, Y). To demonstrate how the graph works in this study, the information of Figure 4 is derived as graphs and drawn in Figure 5. In this figure, it is easy to read the predecessor-successor relationships and the corresponding six events based on the numbers of predecessors and successors including appearance, disappearance, splitting, merger, continuation, and recovery.
As the graphs are derived, the information of each individual tidal flat lifecycle is collected, organized, and stored in a spatiotemporal database simultaneously ( Figure 6). This database can help us better understand and query the morphological changes between the adjacent tidal flat objects, which are invisible in Figure 5, because tidal flat objects are represented by their centroids as the nodes. Especially, the five events defined in Section 2.3 are determined by the area comparison between the predecessor and successor, and this spatiotemporal database can help us capture these events. Moreover, for the six events based on the number of predecessors and successors, we introduce two graph operators to assist the analysis for tidal flat dynamic activities.    For any node in a graph, its degree is defined as the number of edges connected to this node. For a directed graph, the degree of the node is further divided into indegree and outdegree. The indegree is defined as the number of arrows ending at this node, and the outdegree is defined as the number of arrows starting from this node [57]. In the context of this study, for any object, the indegree of its corresponding node equals the number of its predecessors, and the outdegree of its corresponding node equals the number of its successors.
Based on this perception, the five events based on the numbers of predecessors and successors correspond to the following cases: • Appearance: when the node's indegree is 0; • Disappearance: when the node's outdegree is 0; • Splitting: when the node's outdegree is greater than or equal to 2; • Merger: when the node's indegree is greater than or equal to 2; and • Continuation: when the successor node's indegree is 1 and the predecessor node's outdegree is also 1.

Find Non-Cut Vertex to Identify the Event of Recovery
A node v of a graph G is a cut vertex if the removal of v increases the number of components of G [57]. Apparently, the starting and ending nodes, which correspond to appearance and disappearance, are not the cut vertices. For the tidal flat lifecycles, which do not include the event of recovery, the rest nodes must be cut vertices. On the other hand, according to the definition of recovery, there are at least two parallel paths that exist between (but do not include) the splitting and merger nodes, and the nodes on the parallel paths are not the cut vertices. Therefore, we can justify whether a tidal flat lifecycle includes the recovery event or not by calculating the number of non-cut vertices other than the starting and ending nodes. If the number equals to zero, the tidal flat lifecycle has no recovery event; if the number is positive, the tidal flat lifecycle must include at least one recovery event.

Data and Study Area
The following datasets are used in this case study: (1) Landsat TM (Thematic Mapper), ETM+ (Enhanced Thematic Mapper Plus), OLI (Operational Land Imager), and TIRS (Thermal Infrared Sensor) acquired by Landsat 4, 5, 7, and 8 during 1984 to 2018; (2) the ETOPO1 Global Relief Model; (3) the global surface water occurrence data; (4) the training point dataset; and (5) the validation point dataset. We use Google Earth Engine (GEE), a high-performance cloud computing platform, to access and process the above datasets and implement the RF machine learning algorithm. More than thirty years of geospatial datasets have been stored in the GEE virtual archives, which occupy over twenty petabytes of space and are updated and expanded daily. These large datasets can be instantly invoked and processed by the high-performance computing resources [58].
The case study area is the southwest tip of Florida peninsula, which belongs to the counties of Collier, Monroe, and Miami-Dade. This area is located within the coverage of the Landsat tiles of Path=15 Row=42, Path=15 Row=43, and Path=16 Row=42 (Figure 7). From 1984 to 2018, a total of 2906 tiles of Landsat 4, 5, 7, and 8 images in this area have been pre-processed to surface reflectance Tier 1, and all of them are used to delineate the tidal flat cells in this study. The annual distribution of available Landsat images is shown in Figure 8  A major part of this study area belongs to Everglades National Park. Geologically, the study area is divided into two subprovinces: the west coastal area belongs to the Ten Thousand Islands Subprovince, while the south coastal area belongs to the Southern Atlantic Coastal Strip Subprovince. The Ten Thousand Islands Subprovince is the transitional area between the Everglades and Gulf of Mexico, which consists of tidal flats, lagoons, uninhabited islands, mangrove forests, and highly productive estuaries. As the climate changes, the sea level fluctuates, while the sediment grows and recedes over time. Tidal flat is one of the products during this procedure, which plays an important role in protecting the inland areas from the hurricane-caused floods and breaking waves in the Gulf of Mexico. The Southern Atlantic Coastal Strip Subprovince consists of marine limestone, which is covered by thin sheets of quartz sand. The limestone stops water in the Everglades from flowing into the Atlantic Ocean. However, there is a marshy channel system across the limestone ridge (i.e., Taylor Slough), which drains the water  A major part of this study area belongs to Everglades National Park. Geologically, the study area is divided into two subprovinces: the west coastal area belongs to the Ten Thousand Islands Subprovince, while the south coastal area belongs to the Southern Atlantic Coastal Strip Subprovince. The Ten Thousand Islands Subprovince is the transitional area between the Everglades and Gulf of Mexico, which consists of tidal flats, lagoons, uninhabited islands, mangrove forests, and highly productive estuaries. As the climate changes, the sea level fluctuates, while the sediment grows and recedes over time. Tidal flat is one of the products during this procedure, which plays an important role in protecting the inland areas from the hurricane-caused floods and breaking waves in the Gulf of Mexico. The Southern Atlantic Coastal Strip Subprovince consists of marine limestone, which is covered by thin sheets of quartz sand. The limestone stops water in the Everglades from flowing into the Atlantic Ocean. However, there is a marshy channel system across the limestone ridge (i.e., Taylor Slough), which drains the water A major part of this study area belongs to Everglades National Park. Geologically, the study area is divided into two subprovinces: the west coastal area belongs to the Ten Thousand Islands Subprovince, while the south coastal area belongs to the Southern Atlantic Coastal Strip Subprovince. The Ten Thousand Islands Subprovince is the transitional area between the Everglades and Gulf of Mexico, which consists of tidal flats, lagoons, uninhabited islands, mangrove forests, and highly productive estuaries. As the climate changes, the sea level fluctuates, while the sediment grows and recedes over time. Tidal flat is one of the products during this procedure, which plays an important role in protecting the inland areas from the hurricane-caused floods and breaking waves in the Gulf of Mexico. The Southern Atlantic Coastal Strip Subprovince consists of marine limestone, which is covered by thin sheets of quartz sand. The limestone stops water in the Everglades from flowing into the Atlantic Ocean. However, there is a marshy channel system across the limestone ridge (i.e., Taylor Slough), which drains the water in the Everglades into Florida Bay. As a result, the mud is deposited along the coast, and the tidal flat region is produced [59][60][61][62]. Hence, the study area has a sufficient tidal flat region, and it is regarded as an ideal place to implement the case study.

Implementation
As all available Landsat images are found and sorted by year, the case study can be implemented following Figure 9. The Landsat images acquired at the same time are mosaicked and clipped according to the study area's boundary. The clouds, cloud shadows, cirrus, and scan-line corrector (SLC)-off gaps are classified as bad-quality observations. To identify and eliminate them, the Function of mask (FMask) Algorithm is used in this study [63,64]. in the Everglades into Florida Bay. As a result, the mud is deposited along the coast, and the tidal flat region is produced [59][60][61][62]. Hence, the study area has a sufficient tidal flat region, and it is regarded as an ideal place to implement the case study.

Implementation
As all available Landsat images are found and sorted by year, the case study can be implemented following Figure 9. The Landsat images acquired at the same time are mosaicked and clipped according to the study area's boundary. The clouds, cloud shadows, cirrus, and scan-line corrector (SLC)-off gaps are classified as bad-quality observations. To identify and eliminate them, the Function of mask (FMask) Algorithm is used in this study [63,64]. Next, we use the cells of the same location from all mosaicked-denoised images under the same annual stack to calculate the 56 predictor variables as described in Section 2.1. As the predictor variable calculation and the RF classifier training are finished, the RF classification could delineate the tidal flat cells and have an overall accuracy of 95.23%. In the end, every annual stack contributes one binary image. The following processes will be operated based on this tidal flat cell delineation result.
As shown in Figure 9, the following processes are tidal flat object assemblage, lifecycle tracking, the graph-based representation, and the spatiotemporal analysis of events. Proceeding with the work, a prototype system is developed via MATLAB. A threshold value of 3 is used to discard the tidal flat objects, which are smaller than three cells' coverage area (2700 m 2 ), and the value of 0.6 is used as the overlapping threshold (Equation (5)) to determine the predecessor-successor relationship.
Based on the implemental procedures as illustrated above, we get the graphs of tidal flat lifecycles, and the spatial, temporal, and morphologic attributes of tidal flats are stored in the database. The graph operators are used to analyze the generated graphs and get the spatiotemporal characteristics for each type of event. The results of this case study are shown and discussed in Section 3.3.

Results and Discussion
3.3.1. Characteristics of Tidal Flat Cells, Objects, and Lifecycles Next, we use the cells of the same location from all mosaicked-denoised images under the same annual stack to calculate the 56 predictor variables as described in Section 2.1. As the predictor variable calculation and the RF classifier training are finished, the RF classification could delineate the tidal flat cells and have an overall accuracy of 95.23%. In the end, every annual stack contributes one binary image. The following processes will be operated based on this tidal flat cell delineation result.
As shown in Figure 9, the following processes are tidal flat object assemblage, lifecycle tracking, the graph-based representation, and the spatiotemporal analysis of events. Proceeding with the work, a prototype system is developed via MATLAB. A threshold value of 3 is used to discard the tidal flat objects, which are smaller than three cells' coverage area (2700 m 2 ), and the value of 0.6 is used as the overlapping threshold (Equation (5)) to determine the predecessor-successor relationship.
Based on the implemental procedures as illustrated above, we get the graphs of tidal flat lifecycles, and the spatial, temporal, and morphologic attributes of tidal flats are stored in the database. The graph operators are used to analyze the generated graphs and get the spatiotemporal characteristics for each type of event. The results of this case study are shown and discussed in Section 3.3.

Characteristics of Tidal Flat Cells, Objects, and Lifecycles
As described in Section 3.2, we get a total of 35 binary images, which are the annual records of the tidal flat cells from 1984 to 2018. To demonstrate the delineation result, Figure 10 shows the annual tidal flat area from 1984 to 2018. As shown in the figure, the maximum record (262.5 km 2 ) happened in 1992, which is significantly greater than the mean annual area (129.39 km 2 ) and the second maximum record (173.1 km 2 in 2008). Hurricane Andrew (Category 5) may be used to explain the unusual record in 1992. According to Risi et al. (1995) [65], this hurricane created a substantial amount of sediment deposits along the shallow subtidal coastline in our study area. As described in Section 3.2, we get a total of 35 binary images, which are the annual records of the tidal flat cells from 1984 to 2018. To demonstrate the delineation result, Figure 10 shows the annual tidal flat area from 1984 to 2018. As shown in the figure, the maximum record (262.5 km 2 ) happened in 1992, which is significantly greater than the mean annual area (129.39 km 2 ) and the second maximum record (173.1 km 2 in 2008). Hurricane Andrew (Category 5) may be used to explain the unusual record in 1992. According to Risi et al. (1995) [65], this hurricane created a substantial amount of sediment deposits along the shallow subtidal coastline in our study area. To better understand the spatial distribution of the delineated tidal flat cells, we aggregate the occurrences of tidal flat cells based on all 35 binary images ( Figure 11). As shown, the tidal flat cells in different years are heavily overlapped, and they are mostly distributed in four regions: (A) the northwestern corner of the study area; (B) the middle portion of the west coast; (C) the turning point between the west coast and south coast; and (D) the eastern portion of the south coast. Region A belongs to Rookery Bay National Estuarine Research Reserve, which is one of the few remaining undisturbed mangrove estuary system in North America. This system plays an important role for the resuspension and transportation of sediments and contributes to an advantageous environment for tidal flat dynamic activities [66]. Similarly, region C is located around Cape Sable, in which a man-made canal project was carried out in 1922. This project has largely changed the hydrologic environment of the surrounding area. As a result, an active tidal flat region was produced around Sandy Key Basin, which is located off the coast of Cape Sable [67]. On the other hand, region B is located around the estuary of Shark River Slough, while region D is located around the estuary of Taylor Slough. These two rivers are the principal natural drainages that keep carrying the overland water and mud from the Everglades to the sea. Consequently, the tidal flat regions are intensive in the two places. To better understand the spatial distribution of the delineated tidal flat cells, we aggregate the occurrences of tidal flat cells based on all 35 binary images ( Figure 11). As shown, the tidal flat cells in different years are heavily overlapped, and they are mostly distributed in four regions: (A) the northwestern corner of the study area; (B) the middle portion of the west coast; (C) the turning point between the west coast and south coast; and (D) the eastern portion of the south coast. Region A belongs to Rookery Bay National Estuarine Research Reserve, which is one of the few remaining undisturbed mangrove estuary system in North America. This system plays an important role for the resuspension and transportation of sediments and contributes to an advantageous environment for tidal flat dynamic activities [66]. Similarly, region C is located around Cape Sable, in which a man-made canal project was carried out in 1922. This project has largely changed the hydrologic environment of the surrounding area. As a result, an active tidal flat region was produced around Sandy Key Basin, which is located off the coast of Cape Sable [67]. On the other hand, region B is located around the estuary of Shark River Slough, while region D is located around the estuary of Taylor Slough. These two rivers are the principal natural drainages that keep carrying the overland water and mud from the Everglades to the sea. Consequently, the tidal flat regions are intensive in the two places. By comparing the 35 cell level maps, we can summarize the area changes of tidal flats in the above four regions and the rest of the region, as shown in Figure 12. Region A has the smallest annual average area among the four regions (8.22 km 2 ). The tidal flat area in this region periodically changes, as it reaches the peak every eleven years (1986, 1997, and 2008). The maximum area happened in 1997, which is 19.62 km 2 , while the minimum area happened in 2016, which is only 1.87 km 2 . Region B has the second smallest annual average area among the four regions (18.04 km 2 ). Its tidal flat area change also has regularity, as it reaches the peak every five to seven years (1986, 1991, 1997, 2002, 2008, and 2015). Similar to region A, the maximum tidal flat area of region B happened in 1997 (32.75 km 2 ), and the minimum tidal flat area (7.78 km 2 ) happened in 2017, which is only one year after that of region A. Region C has the second largest annual average area among the four regions (42.28 km 2 ). Although the regularity of the tidal flat area change is not significant, it is obvious that two peaks exist in 1989 and 2009, while there is a valley between these two peaks. The maximum area (72.61 km 2 ) happened in 2009, while the minimum area (8.84 km 2 ) happened in 1995. Region D has the largest annual average area among the four regions (51.43 km 2 ). Similar to region C, its tidal flat area change does not have significant regularity. The maximum area (100.3 km 2 ) happened in 1992, while the minimum area (14.96 km 2 ) happened in 1989. The annual average area of the tidal flat in the rest of the region is 9.42 km 2 , which is only 7.3% of the annual average area of the tidal flat in the whole study area. It further proves that the tidal flat cells in the study area are mainly distributed in the four intensive regions. By comparing the 35 cell level maps, we can summarize the area changes of tidal flats in the above four regions and the rest of the region, as shown in Figure 12. Region A has the smallest annual average area among the four regions (8.22 km 2 ). The tidal flat area in this region periodically changes, as it reaches the peak every eleven years (1986, 1997, and 2008). The maximum area happened in 1997, which is 19.62 km 2 , while the minimum area happened in 2016, which is only 1.87 km 2 . Region B has the second smallest annual average area among the four regions (18.04 km 2 ). Its tidal flat area change also has regularity, as it reaches the peak every five to seven years (1986, 1991, 1997, 2002, 2008, and 2015). Similar to region A, the maximum tidal flat area of region B happened in 1997 (32.75 km 2 ), and the minimum tidal flat area (7.78 km 2 ) happened in 2017, which is only one year after that of region A. Region C has the second largest annual average area among the four regions (42.28 km 2 ). Although the regularity of the tidal flat area change is not significant, it is obvious that two peaks exist in 1989 and 2009, while there is a valley between these two peaks. The maximum area (72. 61  Based on the delineated tidal flat cells, a total of 3672 objects are assembled, and 962 tidal flat lifecycles are tracked from the study area. Every tidal flat lifecycle has a duration of at least one year. Figure 13 shows the duration distribution of these 962 tidal flat lifecycles and gives the mean duration of 2.  To demonstrate the graph representation and analysis, we select the 178 th tidal flat lifecycle as an example, and it is illustrated in Figure 14. This tidal flat lifecycle lasts six years (1984)(1985)(1986)(1987)(1988)(1989)(1990) and includes a variety types of events, so it is considered as an ideal case to be analyzed. The location of this tidal flat lifecycle within the whole study area is shown in Figure 14a, which is in region D. As described in Figure 6, the starting and ending coordinates of each tidal flat lifecycle are stored in the spatiotemporal database, so we can easily query the location of any specified tidal flat lifecycle. Moreover, the locations of the tidal flat objects (represented by their centroids as the nodes) and the movement track (represented as the edges) within this lifecycle are shown in Figure 14b, which are also queried from the spatiotemporal database. The derived graph of Figure Based on the delineated tidal flat cells, a total of 3672 objects are assembled, and 962 tidal flat lifecycles are tracked from the study area. Every tidal flat lifecycle has a duration of at least one year. Figure 13 shows the duration distribution of these 962 tidal flat lifecycles and gives the mean duration of 2.  To demonstrate the graph representation and analysis, we select the 178 th tidal flat lifecycle as an example, and it is illustrated in Figure 14. This tidal flat lifecycle lasts six years (1984)(1985)(1986)(1987)(1988)(1989)(1990) and includes a variety types of events, so it is considered as an ideal case to be analyzed. The location of this tidal flat lifecycle within the whole study area is shown in Figure 14a, which is in region D. As described in Figure 6, the starting and ending coordinates of each tidal flat lifecycle are stored in the spatiotemporal database, so we can easily query the location of any specified tidal flat lifecycle. Moreover, the locations of the tidal flat objects (represented by their centroids as the nodes) and the movement track (represented as the edges) within this lifecycle are shown in Figure 14b, which are also queried from the spatiotemporal database. The derived graph of Figure   Figure 13. The duration distribution of the 962 tidal flat lifecycles.
To demonstrate the graph representation and analysis, we select the 178th tidal flat lifecycle as an example, and it is illustrated in Figure 14. This tidal flat lifecycle lasts six years (1984)(1985)(1986)(1987)(1988)(1989)(1990) and includes a variety types of events, so it is considered as an ideal case to be analyzed. The location of this tidal flat lifecycle within the whole study area is shown in Figure 14a, which is in region D. As described in Figure 6, the starting and ending coordinates of each tidal flat lifecycle are stored in the spatiotemporal database, so we can easily query the location of any specified tidal flat lifecycle. Moreover, the locations of the tidal flat objects (represented by their centroids as the nodes) and the movement track (represented as the edges) within this lifecycle are shown in Figure 14b, which are also queried from the spatiotemporal database. The derived graph of Figure 14b is drawn in Figure 14c. Compared with Figure 14b, the derived graph in Figure 14c is a succinct version to illustrate the predecessor-successor relationships within the 178th tidal flat lifecycle. By calculating the indegree and outdegree of each object and querying the areas of each tidal flat object from the spatiotemporal database, we can easily summarize the event types existing in this lifecycle, which include appearance, continuation (stability), splitting (separation), merger (annexation), recovery, and disappearance. The corresponding event information is also stored in the database. This graph representation and analysis provide a straightforward way to examine the dynamic activities within an individual tidal flat lifecycle. 14b is drawn in Figure 14c. Compared with Figure 14b, the derived graph in Figure 14c is a succinct version to illustrate the predecessor-successor relationships within the 178 th tidal flat lifecycle. By calculating the indegree and outdegree of each object and querying the areas of each tidal flat object from the spatiotemporal database, we can easily summarize the event types existing in this lifecycle, which include appearance, continuation (stability), splitting (separation), merger (annexation), recovery, and disappearance. The corresponding event information is also stored in the database. This graph representation and analysis provide a straightforward way to examine the dynamic activities within an individual tidal flat lifecycle. As the first effort of the multi-level GIS framework on the issue of tidal flats, we examine the dynamic activities of tidal flats from a novel perspective. The temporal and spatial characteristics of tidal flat events are presented in the below sections. Since the cases of the continuation event are fully covered by three types of events based on the area comparison between the predecessor and successor (contraction, stability, and expansion), the below sections do not discuss the continuation event. As the first effort of the multi-level GIS framework on the issue of tidal flats, we examine the dynamic activities of tidal flats from a novel perspective. The temporal and spatial characteristics of tidal flat events are presented in the below sections. Since the cases of the continuation event are fully covered by three types of events based on the area comparison between the predecessor and successor (contraction, stability, and expansion), the below sections do not discuss the continuation event.

Temporal Characteristics of Tidal Flat Events
Annual numbers of appearance, disappearance, splitting, merger, and recovery are drawn in Figures 15 and 16. For better reference, the descriptive statistics derived from these two figures is shown in Table 1. As shown in Figure 15, the initial year (1984) in this case study has the historical maximum record of appearance (104). Apparently, it could be overestimated, because we do not have available data before 1984. Similarly, the tidal flat lifecycles that end in the final year (2018) may exist longer and disappear in the future years. To eliminate the interference, these two years' records of appearance and disappearance are both ignored in the statistics in Table 1.
From Figure 15, we can see that the annual occurrence of appearance is more stable than that of disappearance, and this preliminary observation can be verified in Table 1. As the smaller value of standard deviation implies higher stability, the event of appearance has the smaller standard deviation (σ = 7.76) than that of disappearance (σ = 9.93). In particular, there were 60 disappearance events in 2004 and 46 disappearance events in 2011, while the mean value of annual occurrence of disappearance is 27.51. We can conclude that the years of 2004 and 2011 have unusually higher occurrences of disappearance. Annual numbers of appearance, disappearance, splitting, merger, and recovery are drawn in Figures 15 and 16. For better reference, the descriptive statistics derived from these two figures is shown in Table 1. As shown in Figure 15, the initial year (1984) in this case study has the historical maximum record of appearance (104). Apparently, it could be overestimated, because we do not have available data before 1984. Similarly, the tidal flat lifecycles that end in the final year (2018) may exist longer and disappear in the future years. To eliminate the interference, these two years' records of appearance and disappearance are both ignored in the statistics in Table 1.
From Figure 15, we can see that the annual occurrence of appearance is more stable than that of disappearance, and this preliminary observation can be verified in Table 1. As the smaller value of standard deviation implies higher stability, the event of appearance has the smaller standard deviation (σ = 7.76) than that of disappearance (σ = 9.93). In particular, there were 60 disappearance events in 2004 and 46 disappearance events in 2011, while the mean value of annual occurrence of disappearance is 27.51. We can conclude that the years of 2004 and 2011 have unusually higher occurrences of disappearance.  In Table 1, we can also find the annual mean occurrences of splitting, merger, and recovery are 16.09, 15.12, and 2.50. Based on these mean values, we can find a total of five peaks with significantly higher event numbers than the corresponding mean values in Figure 16, which are in the years of 1988, 1996, 2004, 2011, and 2015. Aside from the peak Annual numbers of appearance, disappearance, splitting, merger, and recovery are drawn in Figures 15 and 16. For better reference, the descriptive statistics derived from these two figures is shown in Table 1. As shown in Figure 15, the initial year (1984) in this case study has the historical maximum record of appearance (104). Apparently, it could be overestimated, because we do not have available data before 1984. Similarly, the tidal flat lifecycles that end in the final year (2018) may exist longer and disappear in the future years. To eliminate the interference, these two years' records of appearance and disappearance are both ignored in the statistics in Table 1.
From Figure 15, we can see that the annual occurrence of appearance is more stable than that of disappearance, and this preliminary observation can be verified in Table 1. As the smaller value of standard deviation implies higher stability, the event of appearance has the smaller standard deviation (σ = 7.76) than that of disappearance (σ = 9.93). In particular, there were 60 disappearance events in 2004 and 46 disappearance events in 2011, while the mean value of annual occurrence of disappearance is 27.51. We can conclude that the years of 2004 and 2011 have unusually higher occurrences of disappearance.  In Table 1, we can also find the annual mean occurrences of splitting, merger, and recovery are 16.09, 15.12, and 2.50. Based on these mean values, we can find a total of five peaks with significantly higher event numbers than the corresponding mean values in Figure 16, which are in the years of 1988, 1996, 2004, 2011, and 2015. Aside from the peak In Table 1, we can also find the annual mean occurrences of splitting, merger, and recovery are 16.09, 15.12, and 2.50. Based on these mean values, we can find a total of five peaks with significantly higher event numbers than the corresponding mean values in Figure 16, which are in the years of 1988, 1996, 2004, 2011, and 2015. Aside from the peak in 2015, the other peaks imply a persistent and stable periodicity of about seven to eight years. Recovery 10 (1988, 2014, 2015) 0 (1985-1986, 1990-1992, 1995, 1997, 1999, 2001-2002, 2006, 2008-2010, 2014 The annual number and descriptive statistics for the events based on the area comparison between the predecessor and successor are shown in Figures 17 and 18, and Table 1. Figure 17 and Table 1 10 (1988, 2014, 2015) 0 (1985-1986, 1990-1992, 1995, 1997, 1999, 2001-2002, 2006, 2008-2010, 2014 The annual number and descriptive statistics for the events based on the area comparison between the predecessor and successor are shown in Figures 17 and 18, and Table 1. Figure 17 and Table 1   of annexation is 4.41, which is significantly larger than the overall annual mean of the event occurrence of annexation (x̄ = 3.91).
In addition to the above discoveries, another interesting finding is that the year of 2004 is a unique year. As suggested in Section 3.3.1, the hurricane activities may be the reason. There have been three hurricanes that affected the study area in that year, which are Hurricane Charley (Category 4), Hurricane Frances (Category 4), and Hurricane Ivan (Category 5) [68].

Spatial Characteristics of Tidal Flat Events
In addition to the temporal characteristics, the spatial characteristics of each type of event can be revealed by plotting the event occurrence places as points on the map. These event occurrence places are determined by the definitions of the specified types of events. In detail, the disappearance, splitting, and separation events are only defined by the number of successors. Therefore, each event must have exactly one predecessor node, which is regarded as the event occurrence place. However, for the rest event types, the terminal node is unique (i.e., the successor node of appearance, continuation, contraction, stability, expansion, merger, annexation, and the successor node of the merger is part of recovery). So, the occurrence places of these events are determined as the terminal node. For better visualization, a fishnet is created via ArcGIS to cover the whole study area, in which the size of each grid is 4 km 2 . The number of event occurrences is counted for each grid, and the grid-based summary maps for all types of events are drawn in Figure 19. The preliminary observation of Figure 19 suggests that these event types can be classified as four categories according to their occurrences, in which the events of appearance (Figure 19a (Figure 19i) belong to the third highest occurrence category, and the event of recovery (Figure 19j) belongs to the lowest occurrence category. This finding is consistent with the comparison of the annual means (x) in Table 1 and it is easy to understand, because the event of annexation is the special case of merger, event of separation is the special case of splitting, and the event of recovery is the special case of both merger and splitting.
According to Section 3.3.1, the tidal flat cells are intensively distributed in four regions within the study area, which should potentially affect the spatial distribution of each event. This is proved with Figure 19, as all types of events are mainly distributed within or around the four regions. Accordingly, the occurrence places of these events are In addition to the above discoveries, another interesting finding is that the year of 2004 is a unique year. As suggested in Section 3.3.1, the hurricane activities may be the reason. There have been three hurricanes that affected the study area in that year, which are Hurricane Charley (Category 4), Hurricane Frances (Category 4), and Hurricane Ivan (Category 5) [68].

Spatial Characteristics of Tidal Flat Events
In addition to the temporal characteristics, the spatial characteristics of each type of event can be revealed by plotting the event occurrence places as points on the map. These event occurrence places are determined by the definitions of the specified types of events. In detail, the disappearance, splitting, and separation events are only defined by the number of successors. Therefore, each event must have exactly one predecessor node, which is regarded as the event occurrence place. However, for the rest event types, the terminal node is unique (i.e., the successor node of appearance, continuation, contraction, stability, expansion, merger, annexation, and the successor node of the merger is part of recovery). So, the occurrence places of these events are determined as the terminal node. For better visualization, a fishnet is created via ArcGIS to cover the whole study area, in which the size of each grid is 4 km 2 . The number of event occurrences is counted for each grid, and the grid-based summary maps for all types of events are drawn in Figure 19. The preliminary observation of Figure 19 suggests that these event types can be classified as four categories according to their occurrences, in which the events of appearance (Figure 19a This finding is consistent with the comparison of the annual means (x) in Table 1 and it is easy to understand, because the event of annexation is the special case of merger, event of separation is the special case of splitting, and the event of recovery is the special case of both merger and splitting. est records among all types of events and verify the above finding. Another interesting finding is that even though the tidal flat area in region C is not the largest among the four regions, a total of eight types of events (appearance, disappearance, contraction, expansion, merger, annexation, splitting, and separation) are more in region C than elsewhere. Similarly, although the tidal flat area in region A is the smallest among the four regions, the number of event occurrences in region A is ranked as the third most among the four regions in cases of appearance, disappearance, and stability events and the second most among the four regions in cases of the rest event types. Therefore, we can conclude that the tidal flat dynamic activities are more frequent in regions A and C than regions B and D. This may be explained by regarding the fact that the tidal flat areas in regions B and D are mostly the products of the principal natural drainages, which are different from the cases of regions A and C.  According to Section 3.3.1, the tidal flat cells are intensively distributed in four regions within the study area, which should potentially affect the spatial distribution of each event. This is proved with Figure 19, as all types of events are mainly distributed within or around the four regions. Accordingly, the occurrence places of these events are counted by five classes (the four intensive regions and the rest region), and the summary is drawn in Figure 20. This figure shows that 89.47% of annexation events and 75.44% of separation events are distributed within the four regions, which are the highest and lowest records among all types of events and verify the above finding. Another interesting finding is that even though the tidal flat area in region C is not the largest among the four regions, a total of eight types of events (appearance, disappearance, contraction, expansion, merger, annexation, splitting, and separation) are more in region C than elsewhere. Similarly, although the tidal flat area in region A is the smallest among the four regions, the number of event occurrences in region A is ranked as the third most among the four regions in cases of appearance, disappearance, and stability events and the second most among the four regions in cases of the rest event types. Therefore, we can conclude that the tidal flat dynamic activities are more frequent in regions A and C than regions B and D. This may be explained by regarding the fact that the tidal flat areas in regions B and D are mostly the products of the principal natural drainages, which are different from the cases of regions A and C. The above spatiotemporal characteristics for the events prove that the graph representation in our case study can visualize the tracked tidal flat lifecycles and make it possible to capture the events from them. In addition, the graph operators in our case study provide an effective way to automatically and exhaustively capture the events from a large amount of tidal flat lifecycles, which also highlights that the object and lifecycle levels could provide more diversified dynamic information than the cell level.

Conclusions
An integrated approach of three-level GIS framework and graph model is proposed in this study, which provides an effective way to track, represent, and analyze the dynamic activities of tidal flats from a novel perspective. The three-level GIS framework consists of the cell level, the object level, and the lifecycle level. Furthermore, a variety types of events are defined based on (1) the numbers of predecessors or successors of each tidal flat object; and (2) the area comparison between the predecessor and successor. These events provide an enhanced way to describe the dynamic activities of tidal flats throughout their lifecycles.
On the other hand, the graph model conceptualizes the tidal flat objects as the nodes, and the predecessor-successor relationships as the edges. As the derivative of the graph model, a spatiotemporal database is derived simultaneously, which offers the storage and query features for a variety types of morphologic attributes of tidal flats. In addition, graph operators are introduced in this study, which assist the proposed graph model to capture the events automatically and exhaustively.
In the case study, all GEE-archived Landsat imageries from 1984 to 2018 in the study area are invoked, and the RF machine learning algorithm is used to extract tidal flat cells from them. The tidal flat cells are assembled as tidal flat objects, which are further linked as tidal flat lifecycles with respect to the predecessor-successor relationships. The MATLAB prototype system is implemented to track, represent, and analyze the dynamic activities of tidal flats. The results imply that the geology, hydrology, and meteorology factors may affect the dynamic activities of tidal flats in the southwest tip of Florida Peninsula.
Our research is the first effort that systematically organizes the multi-level GIS framework and graph theory to track, represent, and analyze the dynamic activities of tidal flats, which also provides a novel perspective to examine the other dynamic geo- The above spatiotemporal characteristics for the events prove that the graph representation in our case study can visualize the tracked tidal flat lifecycles and make it possible to capture the events from them. In addition, the graph operators in our case study provide an effective way to automatically and exhaustively capture the events from a large amount of tidal flat lifecycles, which also highlights that the object and lifecycle levels could provide more diversified dynamic information than the cell level.

Conclusions
An integrated approach of three-level GIS framework and graph model is proposed in this study, which provides an effective way to track, represent, and analyze the dynamic activities of tidal flats from a novel perspective. The three-level GIS framework consists of the cell level, the object level, and the lifecycle level. Furthermore, a variety types of events are defined based on (1) the numbers of predecessors or successors of each tidal flat object; and (2) the area comparison between the predecessor and successor. These events provide an enhanced way to describe the dynamic activities of tidal flats throughout their lifecycles.
On the other hand, the graph model conceptualizes the tidal flat objects as the nodes, and the predecessor-successor relationships as the edges. As the derivative of the graph model, a spatiotemporal database is derived simultaneously, which offers the storage and query features for a variety types of morphologic attributes of tidal flats. In addition, graph operators are introduced in this study, which assist the proposed graph model to capture the events automatically and exhaustively.
In the case study, all GEE-archived Landsat imageries from 1984 to 2018 in the study area are invoked, and the RF machine learning algorithm is used to extract tidal flat cells from them. The tidal flat cells are assembled as tidal flat objects, which are further linked as tidal flat lifecycles with respect to the predecessor-successor relationships. The MATLAB prototype system is implemented to track, represent, and analyze the dynamic activities of tidal flats. The results imply that the geology, hydrology, and meteorology factors may affect the dynamic activities of tidal flats in the southwest tip of Florida Peninsula.
Our research is the first effort that systematically organizes the multi-level GIS framework and graph theory to track, represent, and analyze the dynamic activities of tidal flats, which also provides a novel perspective to examine the other dynamic geographic phenomena with large spatiotemporal scales. The future work will be conducted in three aspects. First, we plan to involve the whole coastal area of the conterminous US to implement this integrated approach, to not only verify its universality and robustness, but also find more diversified spatiotemporal regularities of tidal flat dynamic activities. Second, we are also interested in the in-depth explorations for the natural and anthropogenic factors, which may potentially affect the spatiotemporal characteristics of tidal flat dynamic activities such as urbanization and sea level rise. Third, we noticed that the NOAA tide stations (https://tidesandcurrents.noaa.gov/) are producing high temporal resolution data of the water level along the US coast, which enables and encourages us to contribute to the interdisciplinary research collaborated with the colleagues of coastal hydrology and oceanography backgrounds. Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://developers.google.com/earth-engine/datasets/.