Estuarine Mapping and Eco-Geomorphological Characterization for Potential Application in Conservation and Management: Three Study Cases along the Iberian Coast

Geomorphological changes in recent decades in three estuaries along the Iberian coast were analysed using aerial orthophotographs. A hierarchical classification scheme, based on a literature review representing 26 estuarine eco-geomorphological features relevant to estuarine dynamics and functioning, is described. The estuaries selected were San Vicente de la Barquera (N Spain), Guadiana River (SW border between Spain and Portugal) and the Ebro River Delta mouth (NE Spain). For these systems, a 60-year time series of high-resolution maps was developed, analysing the changes in feature surfaces. The main subsystems analysed were beach, dunes, saltmarshes and the drainage network. The results of the cartographies showed general behaviour common to all transitional systems, relationships among main subsystems and processes inherent to each one. This work illustrates how beaches and dunes serve as a protective barrier for the tidal flats, acting as a sediment buffer for the entire system. The subsystems are connected by the drainage network responsible for the exchange of matter and energy between them. Furthermore, an accuracy assessment was performed in one of the study zones to identify the limitations of mapping with aerial photographs. The results explain the changes with time but also the processes and relationships between the estuarine features at a long-term scale. This work adds an important perspective towards a general understanding of their dependence on intrinsic and boundary conditions.


Introduction
The abundance of fertile land and natural resources has always been a claim for the settlement of human populations in estuarine areas [1][2][3]. Besides providing physical support for the location of human settlements, coastal ecosystems in general, and estuaries in particular, provide important ecosystem services, such as resource supply (extractive fishing, aquaculture or salt and aggregate extraction), culture (tourism and leisure activities) and regulation services (energy and morphosedimentary regulation or biological control; [4]). Estuaries are amongst the most productive ecosystems worldwide [5]. Nonetheless, for decades, estuarine management has been based mainly on the economic exploitation of the resources offered by these ecosystems [6]. This had led to the loss of a large part of the surface occupied by intertidal areas [7], directly affecting the general functioning of the system and, therefore, associated subsystems and the processes taking place on them. These changes generate problems such as eutrophication, alteration of hydrological regimes by the construction of dams, and the loss of areas occupied by the main systems within the estuaries [8,9]. Because of the above, in combination with short-to medium-term projections, the nature and dynamics of the substrate will be affected and, therefore, the vegetation cover. These changes may have consequences for the extension and distribution of the habitats, and situations of natural risk will increase [10].
These pressures and the imminent changes have led administrations to develop environmental management policies to evaluate the state of conservation of coastal habitat. They largely respond to the European Union (EU) mandate through the "Habitat Directive", which states that EU member states must monitor the conservation status of the habitats on their territories. The interest in coastal habitats has led the EU to include specific aspects of these environments in various directives, in addition to the previous ones (e.g., Water Framework Directive, Marine Strategy Framework Directive). To maximize the usability of these directives, the EU has also developed the INSPIRE Directive, which aims for the standardization of the environmental databases and for the development of methodologies that allow each institution to carry out standard gathering of environmental information. Despite various initiatives in recent years [11], the development of holistic methods for monitoring complex and extensive systems, such as estuaries, is still insufficient in Spain. Only a few works have included a multidisciplinary approach on the study of estuarine dynamics, but multidisciplinarity is needed to apply an integral perspective to estuarine management [12,13]. As a consequence, the lack of integral management methods, together with the scarcity of time series, often results in inefficient management measures.
In recent years, numerous studies worldwide have focused on the morphological and ecological changes that estuaries have undergone [5,[14][15][16][17][18]. However, these studies have rarely merged both perspectives [19]. The term eco-geomorphology encompasses the biogeographic approach to the study of interactions between the biosphere and physical systems [20,21]. However, most of these studies focus on plant-substrate relations [22] and rarely address a landscape scale. Although it is an increasingly common approach, it is difficult to find studies that address the study of an estuary as a complex and interconnected system. As described by [23], it is necessary to study the relationships of the estuaries with other adjacent systems, describing the importance of active processes that take place in these areas. The relative balance of waves, tides and river processes exerts considerable control over morphology but also over the ecological dynamics of deltas and estuaries [24]. These conditions, together with the broad climatic and geographical setting, are boundary conditions controlling the formation and evolution of the estuarine habitats. Thus, the essence of estuarine eco-geomorphological functioning is the connectivity between subsystems [5] and their boundary conditions. A priori, this connectivity can be considered obvious and well-described, although the reality is not so simple [25]. The difficulty of combining these perspectives, together with the inherent difficulties of studying complex systems at a whole level, requires significant effort. However, in the context of coastal conservation, it is necessary to evaluate this connectivity to develop multiple zoning and functioning schemes.
The use of aerial imagery and Geographical Information Systems (GIS) for detecting landscape changes is widely used [17,26,27] and, in the case of coastal systems, it can be challenging because of the difficulties derived from imagery quality and atmospheric and oceanographic conditions. Moreover, the older the image, the more difficult it is to analyse due to the storage and scanning qualities. However, analysis of databases of old images is indispensable for long-term studies.
This paper aims to assess changes in three estuaries under different boundary conditions in order to understand the mechanisms involved in estuarine processes at a timescale of decades. To cover this objective, the present work looks for the existence of a general pattern in estuarine systems through high-resolution evolutionary mapping and the local patterns inherent to each system. This work not only explains the estuarine changes but also the processes and relationships between the main estuarine features, adding the contribution of dependency to the relevance of boundary conditions.

Study Sites
This work focuses on three study sites ( Figure 1): (1) the San Vicente de la Barquera estuary (SVB), (2) the Guadiana River estuary (GUA) and (3) the Ebro River Delta mouth (EBR). These zones are located on different biogeographical regions along the Iberian coast: (1) the North Atlantic region, (2) the South Atlantic region, and (3) the Mediterranean region, which covers most of the environmental range and oceanographic conditions that occur on the Iberian Peninsula coasts.

Materials and Methods
For every study site, the historical evolution for the last 60 years was estimated from a variety of cartographic data managed using a Geographic Information System (GIS).

Estuarine Mapping and Area Calculation
Remote sensing imagery is a common resource to create thematic maps [18] and references therein. However, the use of satellite images often requires the definition of specific classifications for single images due to the differences in intrinsic characteristics (e.g., specific external conditions at the moment of the shot and angle and type of sensor, among others). Besides, publicly available remote sensing imagery have lower spatial and temporal resolution than aerial orthophotos (30 m/pixel compared with 0.25 m/pixel in some cases). For these reasons, this study developed on-screen digitization using visual image interpretation of different elements, i.e., colour, texture and plant association [41].
The three study areas ( Figure 1) were mapped using historical aerial ortophotographs obtained from the Spanish National Plan of Aerial Orthophotography (PNOA), Andalusian Environmental Information Network (REDIAM) and Geological and Cartographic Institute of Catalonia (ICGC). The geodetic reference system used was the European Terrestrial Reference System 1989 (ETRS89-H30), according to the guidelines of the Geographic Superior Council, and following the INSPIRE European Directive (2007/2/CE). The orthophotographs of the first years of the study period (i.e., before the 2000s) were scanned and georeferenced in ArcGIS 10.2 through a second-order polynomial function, using at least 15 well-distributed control points on every image. The coordinates of these control points were taken from correctly georeferenced recent orthophotographs. The accuracy of the georeference process was determined by calculating the total Root Mean Square (RMS), less than 0.5 m in all cases. This type of analysis requires ortophotographs taken in low tide-conditions on clear days to easily identify and record the submerged features. Unfortunately, the number of images fulfilling these requirements is limited and it was necessary to soften the requirements to get enough data for this work. Therefore, some maps have a larger uncertainty than others and, as a consequence, The first study site, San Vicente de la Barquera estuary (SVB; Figure 1 (1)), is located on the Cantabrian coast (N Spain, Gulf of Biscay). It is made up of two shallow estuarine subsystems and linked to a short river with a low hydrological regime [28]. Like most of the Cantabrian estuaries, it forms an incised valley and presents a sandy confining barrier, artificially channelled by two lateral jetties which change the morphodynamics of the entire estuary [28]. Regarding hydrodynamic conditions, SVB is a tide-dominated estuary, according to the classification proposed by [23]. It shows a semi-diurnal mesotidal regime with a mean spring tidal range (MSTR) of 3.94 m and a Highest Astronomical Tide (HAT) that reaches 4.85 m (Santander tide gauge, [29]). Tide is the most important dynamic factor in this zone, leading to well-developed tidal flats, extending upstream about 6 km from the river mouth. The predominant waves come from NW, with the most frequent significant wave height ranging between 0.5 and 1 m, peak periods greater than 9 s in most cases (frequency of 65%), and an eastward-dominant longshore transport. In this particular area, the coast orientation favours the refraction of the incident waves and forms a local coastal drift current moving westward. This region is considered a high-energy area, with strong and persistent winds from various directions but mainly blowing from W-NW [28]. The regional climate is oceanic, with no dry season.
The second study site is the Guadiana River estuary (GUA, Figure 1 (2)). It is located on the southern border between Spain and Portugal, consisting of a single channel (50-700 m wide; [30]). It is also incised in a bedrock valley [31]. The marine area of this estuary, i.e., the lowest part of the estuary influenced by seawater, extends 10 km upstream from the mouth [32]. According to [23], it is a wave-dominated system. It has a mesotidal regime, with an MSTR of 3.06 m and a HAT of 3.87 m (Isla Canela tide gauge, [29]). Predominant waves come mainly from W and SW with an average significant height between 0.5 and 1 m [33], which produces a main littoral drift moving eastward. The winds in this region come from NW and SW, with an average wind speed of 3 to 4.5 m/s [29]. The regional climate is semi-arid, with the exception of the summer and winter seasons, which are considered as arid and temperate-humid, respectively [31]. During the mid-20th century, the Portuguese margin was transformed by the construction of jetties and revetments, leading to a severe rigidization of the shoreline; for this reason, in this study, only the more natural, Spanish margin will be evaluated.
The Ebro River Delta mouth (EBR, Figure 1 (3)) is located on the Mediterranean coast of the Iberian Peninsula and is one of the most important deltas in the western Mediterranean [34]. The extension of the emerged area is 325 km 2 , which represents only 15% of the total area of the delta [35]. It is a low-lying area, with elevations not exceeding 5-6 m, and around 50% of the total surface lies below 0.5 m above MSL [36]. The Ebro Delta has a microtidal regime, with a MSTR of about 0.20 m [37], which gives a greater importance to the meteorological tide component. It is a wave-dominated system [35], and the waves come mainly from three directions: E, S and NW. The first is the dominant direction due to its frequency and energetic components, and in fact is responsible for the main transport of sediment in the mouth of the Ebro River Delta. The wind is also an important component in this area, with strong winds during spring, matching with E wave events [38]. The last component is associated with the Mestral winds, corresponding to wave calm periods [39], but energetic wind events. Mestral winds are the ones that blow most intensely in this area, although less frequently, channelled through the Ebro Valley, sometimes reaching more than 100 km/h. They mainly affect the zone during autumn and winter seasons [40]. The regional climate is littoral Mediterranean, with moderate temperatures during winter and sub-arid weather during summer, due to the high temperature and low precipitation. It is also important to take into account the average subsidence of the delta in the medium-long term, which was estimated as 1 to 3.2 mm/year in recent decades by [34].

Materials and Methods
For every study site, the historical evolution for the last 60 years was estimated from a variety of cartographic data managed using a Geographic Information System (GIS).

Estuarine Mapping and Area Calculation
Remote sensing imagery is a common resource to create thematic maps [18] and references therein. However, the use of satellite images often requires the definition of specific classifications for single images due to the differences in intrinsic characteristics (e.g., specific external conditions at the moment of the shot and angle and type of sensor, among others). Besides, publicly available remote sensing imagery have lower spatial and temporal resolution than aerial orthophotos (30 m/pixel compared with 0.25 m/pixel in some cases). For these reasons, this study developed on-screen digitization using visual image interpretation of different elements, i.e., colour, texture and plant association [41].
The three study areas ( Figure 1) were mapped using historical aerial ortophotographs obtained from the Spanish National Plan of Aerial Orthophotography (PNOA), Andalusian Environmental Information Network (REDIAM) and Geological and Cartographic Institute of Catalonia (ICGC). The geodetic reference system used was the European Terrestrial Reference System 1989 (ETRS89-H30), according to the guidelines of the Geographic Superior Council, and following the INSPIRE European Directive (2007/2/CE). The orthophotographs of the first years of the study period (i.e., before the 2000s) were scanned and georeferenced in ArcGIS 10.2 through a second-order polynomial function, using at least 15 well-distributed control points on every image. The coordinates of these control points were taken from correctly georeferenced recent orthophotographs. The accuracy of the georeference process was determined by calculating the total Root Mean Square (RMS), less than 0.5 m in all cases. This type of analysis requires ortophotographs taken in low tide-conditions on clear days to easily identify and record the submerged features. Unfortunately, the number of images fulfilling these requirements is limited and it was necessary to soften the requirements to get enough data for this work. Therefore, some maps have a larger uncertainty than others and, as a consequence, the results must be interpreted with caution.
Every feature was mapped through the study period with the Habitat Digitizer [42] extension to ArcGIS (ESRI © ). Each feature was recorded with a unique ID, according to a hierarchical classification scheme (see below, Table 1), which was based on an extensive literature review. The use of unique IDs allows the habitat digitizer tool to easily define and assign attributes to polygons [43]. The maps were developed with a minimum scale of 1:2500. This delineates features with a limited surface but with enough relevance to have a role in the estuarine system. In addition, this tool allows setting a Minimum Mapping Unit (MMU, [44]), i.e., the minimum object identifiable in the image, and areas lower than 100 m 2 were not mapped. The hierarchical classification scheme was designed following three categories or levels of organization, from general to specific. (a) Category 1 represents the main subsystems within the estuary, common to every mid-latitude estuarine zone: Beaches, Dunes, Tidal Flats/Wetland, Channel and Rocky units. The first two subsystems encompass sandy environments, the third includes muddy environments (named with two different terms in order to differentiate between micro-and mesotidal environments), the fourth corresponds to the tidal drainage subsystem (which is a mix of sandy and muddy environments) and, lastly, the rocky units are only present in the San Vicente de la Barquera estuary. (b) Category 2 covers all of the features present in each subsystem of Category 1, defined mainly by Habitat Directive 92/43/CEE, EUNIS habitat classification and CORINE Land Cover classification. (c) Category 3 encompasses those features which need to be more specific than level 2. Thus, the hierarchical scheme is intended to be a standardized procedure to be applied in every estuarine zone of mid-latitude (Table 1).
To reduce spatial uncertainties, the features defined in Table 1 were first verified in the field to improve their definition and secondly validated with field surveys after mapping. The main active processes for each study area were mapped based on reviewed literature. In the next step, the main ebb and flow directions through the channel networks were identified. These directions are responsible for the development of sandy and muddy environments in the estuarine systems. Table 1. Hierarchical classification applied to the study sites with types of features and brief definitions. Numbers between brackets after the name of the unit indicate the level of the category (Category 1: subsystem, Category 2: features, Category 3: more specific features).

Feature Description, Eco-Geomorphological Role
Shoreline Sandy Environments (1): gently sloping sand-covered shorelines affected by waves and current actions, just above the normal tide limit (EUNIS classification).
Beach (2) Accumulation of shore material formed into distinctive shapes by waves and currents. The beach form is generally a seaward-sloping boundary between a water body and mobile sediment, and a flat or landward-sloping surface at the upper limit of the beach [45].
Estuarine beach (2) Unvegetated or partially vegetated, sand, gravel, or shell intertidal beach in partially enclosed areas (estuaries, bays, lagoons and similar features) connected to the ocean dynamics [46].
Overwash deposit (2) The continuation of the uprush over the crest of the most landward (storm) berm. The resulting deposit is not subject to reworking on the active beach by normal wave and tidal action [48].

Sand bar (2)
Intertidal ridge-shaped accumulations of sand with an associated horn landward. It can be bare, without associated vegetation, or present seagrasses and algae.
Fixed dunes (2) Aeolian deposits whose mobility is impeded by a consistent plant cover (shrub or tree nature). Typical species of this band of the dune system are Crucianella maritima, Helicrysum stoechas, Sporobolus arenarius and Armeria spp., among others.
Estuarine dunes (2) Accumulations of sand, gravel, or shell on the back part of estuarine beaches in partially enclosed areas (estuaries, bays, lagoons and similar features) connected to the ocean dynamics. They present the same species as external dunes and, exceptionally, halophytes typical of tidal environments (e.g., Sarcocornia spp. and Arthrocnemun spp.) Relict dune ridges (2) Non-active dune ridges formed by aeolian processes, usually from fine-to-medium sand, by accretion of embryonic backshore/berm dunes behind the high tide zone [47]. The associated vegetation is typical of mature and established soils, i.e., tree and shrub vegetation.
Wet dune slacks (2) Damp or wet hollows left between dunes where the groundwater reaches or approaches the surface of the sand [51]. Rushes and wet grasses are the main plant communities present in these features.
Tidal flats1 (1): Sedimentary plains developed in areas affected by tides, with a predominance of fine sediment transported by water and stabilized by vegetation [52].
Wetlands2 (1): Shallow aquatic environments, from brackish to hypersaline, isolated or partially connected to the sea in environments with low tidal ranges.
1 Saltmarshes/2 Vegetated wetlands (2) 1 Coastal ecosystem in the upper coastal intertidal zone between land and open saltwater or brackish water, which are intermittently flooded by the tides and colonized by halophytes.
2 Wetlands are ecosystems that depend on constant or recurrent shallow inundation or saturation at or near the surface of the substrate. Common diagnostic features of wetlands are hydric soils and hydrophytic vegetation [53].
1 Mudflats/2 Submerged meadows (2) 1 Muddy platforms common to the intertidal zone of most estuaries (zone below neap high-tide level). The vegetation of these features is poor, mainly seagrasses and pioneer saltmarsh plants [54].
Coastal lagoon (2) Inland water body, usually oriented parallel to the coast, separated from the ocean by a barrier, connected to the ocean by one or more restricted inlets, and having depths which seldom exceed a couple of meters [56].

Feature Description, Eco-Geomorphological Role
Drainage Network (1): Estuarine channel system, including minor forms related to their bottom and margins.
Channel 1 (2) Main river channel. It is the estuary in the strict sense, which extends upstream to the edge of the salt wedge [31]. Channel 2 (2) Distributary channels directly derived from the main channel. Channel 3 (2) Distributary channels derived from the secondary distributary channels. Artificial channel (2) Non-natural channels within the fluvial network.
Sand flat (2) Sandy environment located preferentially in the lower intertidal zone. It is an unstable area characterized by the constant resuspension of sediment by tidal flood and ebb currents [57].
Flood tidal delta (2) Accumulation of sand on the shoreward side of an inlet, initially formed during storm surges and maintained by flood currents [58].
Ebb tidal delta (2) Accumulation of sand on the seaward side of an inlet formed by the ebb tidal current [58].
Cliff (2) Very steep, vertical, or overhanging rock slopes Rocky platform (2) Flat (planar) platforms carved into rocks formed by weathering and wave erosive action.
Once every unit was defined, the changes in the occupied surface over the entire study period were calculated with respect to a reference surface. In addition, the evolution of each subsystem was quantified by setting the surface occupied on the first year as a reference (RS). In this way, all subsequent changes were related to the same area. In features with large changes, they were analysed in more detail (Categories 2 and 3).
For the SVB and GUA estuaries, the reference surface, i.e., 100%, corresponded to 1956. However, the EBR site suffered great changes from 1956 to 1984, with large differences in perimeter and, therefore, a different reference surface. As a consequence, although the initial state is showed in order to better understand the deltas' evolution, the 1956 image was not included in the quantitative calculations and the reference surface was decided to be that of 1984.
Thus, for each feature, the percentage of occupied surface with respect to the reference one (RS) was calculated following (1): In addition, to quantify the rates of change of the features (RC, Category 3) in the total period, the calculation of the rates was made using (2):

Accuracy Assessment
Field campaigns were carried out on dates close to the last available orthophoto to validate the last cartographies. The collection of topographic profiles on the major geomorphological subsystems (beach, dunes and marshes) helped to delineate the most problematic features.
The system with the largest data availability, the SVB estuary, was used as a test site to carry out an accuracy assessment. This procedure allows for an estimation of errors associated with the results. Although these errors increase with the age of the orthophoto, the method is useful to identify the areas with a large probability of errors. The accuracy assessment compares the values on test pixels recorded in the field (reference points), with the corresponding pixels in the map. Besides recording reference points in the field, a flight with an Unmanned Aerial System (UAS) was also performed.
According to the Spanish regulations for UAS operations, the flying height was 120 m above ground level, obtaining an orthomosaic with a 0.03 m pixel size.
The accuracy of the maps was assessed with a confusion matrix. Confusion matrices are a common tool used for the assessment of cartography accuracy. They are widely accepted, although the bias of the test pixels also biases the accuracy of the confusion matrix [27] and references therein. Therefore, to minimize the bias of the text pixels, at least ten test pixels were considered for each feature.
According to [43], the proportion of area (p ij ) of a mapped feature i and reference feature j is estimated according to (3): where n ij is the number of reference points classified as mapped feature i and reference feature j, n i+ and N i+ are the number of pixels per sample and in the entire map in stratum i, and N is the number of pixels in the entire map. The overall proportion of pixels correctly classified (P c ) is the sum of the corrected classified pixels, i.e., the main diagonal of the confusion matrix (4), and represents the probability that a random pixel is correctly classified in the map: The user's accuracy (P Ui ) for feature i is the conditional probability that a random point classified as feature i by the map is classified as feature i by the reference data (5): The producer's accuracy (P Pi ) for feature i is the conditional probability that a random point classified as feature i by the reference data is classified as feature i by the map (6): The kappa coefficient (k) indicates, for a stratified random sample of test pixels, the proportionate reduction in classification error, compared with the error of a random assignment of features (7): where p k+ and p +k are the map and reference sums for feature k, respectively. It ranges between −1 to +1, where 1 indicates perfect agreement between reality and the classified image and 0 indicates complete randomness.

Estuarine Mapping and Area Calculation
The geomorphological features identified in the study sites encompass all the subsystems of Category 1: shoreline sandy environments, dunes, tidal flats/wetlands, channel and rocky units (Table 1) over almost 60 years. Although only 5 years of evolution maps are shown (Figures 2-4), the calculation of the areas were made on every site for all the years indicated in Figure 1.     The evolution of the SVB estuary during the last 60 years (Figure 2) shows a reduction in the area occupied by saltmarshes and, consequently, an increase in the area of mudflats. In the early 1990s, the anthropic features increased markedly at both sides of the estuary, in parallel with the general expansion of urbanization on the Spanish coast. On the eastern side of the estuary, the development of dune fields was directly limited by the increase in anthropic features. In this area, the littoral drift goes to the east. However, as indicated above, there is a local counter-drift conditioned by the morphology of the coast and subsequent wave refraction, making the formation of an external beach and dune fields possible. The beach remained almost the same thorough the entire study period, with slight variations due to differences in tidal level at the time the orthophotos were taken. The ebb tidal delta had little variation in shape, with no relevant effects on area. Regarding the sand flats and the sand bars, especially the ones in the channels of the eastern tidal flat, both increased in The evolution of the SVB estuary during the last 60 years (Figure 2) shows a reduction in the area occupied by saltmarshes and, consequently, an increase in the area of mudflats. In the early 1990s, the anthropic features increased markedly at both sides of the estuary, in parallel with the general expansion of urbanization on the Spanish coast. On the eastern side of the estuary, the development of dune fields was directly limited by the increase in anthropic features. In this area, the littoral drift goes to the east. However, as indicated above, there is a local counter-drift conditioned by the morphology of the coast and subsequent wave refraction, making the formation of an external beach and dune fields possible. The beach remained almost the same thorough the entire study period, with slight variations due to differences in tidal level at the time the orthophotos were taken. The ebb tidal delta had little variation in shape, with no relevant effects on area. Regarding the sand flats and the sand bars, especially the ones in the channels of the eastern tidal flat, both increased in area with time, particularly near the main bridge. These processes were, very probably, facilitated by the main flow directions of the ebb and flood currents (see arrows in Figure 2).
From the three estuaries, the GUA estuary (Figure 3) is the most stable, without notorious natural changes during the study period, but with a remarkable variation in the area of anthropic features around the saltmarsh due to touristic pressures. The urbanization of Isla Canela led to a loss of saltmarsh surface, which was transformed into urban area, crops or grasslands (transformed saltmarsh in Figure 3). Moreover, the area of sandy environments increased, with the formation of new dunes and the growth of the beach. The evolution of the sandy hooks upstream reveals the large availability of sediment in this area. With time, the initial hook was consolidated as a relict dune ridge, and a new active dune ridge was established with the development of a recent saltmarsh between both dune ridges. Additionally, a new hook developed at the end of the external beach in recent years (early 2000s to present), with an associated sheltered area (mudflats* in Figure 3). Very relevant for this system is the presence of a huge sediment bank formed in the outer part of the estuary and mapped as a group of sand bars to unify terms. This sand bank fluctuates over time, making it difficult to map. The drainage network became more developed and denser over the years, partly due to an increase in the quality of the orthophotos, but also to the natural evolution of the saltmarsh. Regarding the main active processes, the littoral drift in this zone moves eastward, although waves and tides have a crucial role in the development of morphologies in the lower part of the estuary, as indicated by the evolution of the main sandy structures (see black arrows in Figure 3).
The EBR estuary suffered great changes between 1956 and 1984 (Figures 4 and 5). The first available aerial photo of the zone (1947, top right image in Figure 5) shows the old mouth of the river, presently located in the southern part of San Antonio island (abandoned channel in Figure 4), and the current active mouth derives from the breakage of the old channel on its left side as a consequence of a river flooding in 1937 [59]. In 1956 ( Figure 5), the original mouth was already closed, and the active mouth opened very widely, initiating the shift towards the present shape. Huge sand sheets covered San Antonio island in 1956, which could be considered as a complex sandy system, i.e., shoreline sandy environments and dunes, with a virtual absence of muddy features, except for the proper deltaic plain. This area did not suffer remarkable changes in anthropic features due to its qualification as a Natural Park in 1989. However, it has undergone important eco-geomorphological changes during the last 30 years (Figure 4). The most evident change was the development of the beach and dune ridges in El Garxal (a submerged meadow in the left side of the main channel). These sandy features protect and make possible the formation of a wide wetland, which opens to the main channel. During the study period, an important input of sediment arrived to El Garxal, supplied by the main littoral drift transport and enabling the formation of new sandy hooks (see left black dotted arrow in Figure 4). These hooks slowly connected to land and enabled the development of a new coastal lagoon (dark blue feature in Figure 4).
Besides, on the western side of the wetland, the parabolic dunes increased in area, even burying part of the vegetated wetland. On the eastern side of the channel (San Antonio island), an erosive trend was maintained in the last years. Despite the main littoral drift flowing to the south, other remarkable currents affect the zone. Prevailing winds and waves come from the east, being responsible for the main sediment transport and shoreline morphology in this area (black dashed arrows in Figure 4). Besides, on the western side of the wetland, the parabolic dunes increased in area, even burying part of the vegetated wetland. On the eastern side of the channel (San Antonio island), an erosive trend was maintained in the last years. Despite the main littoral drift flowing to the south, other remarkable currents affect the zone. Prevailing winds and waves come from the east, being responsible for the main sediment transport and shoreline morphology in this area (black dashed arrows in Figure 4).
The total areas on the GUA ( Figure 3) and EBR ( Figure 4) estuaries are dependent on the state of the sublittoral sand bars attached to the coast and the tide level at the time of the photo (only relevant for GUA as EBR is a microtidal environment). Seasonality is also relevant for the degree of the feature development; however, it is difficult to include this variable in the general calculation of areas.

Changes in Occupied Surface
Changes on the perimeter of occupied surfaces can be used as a proxy for net size fluctuations with time ( Figure 6). This reveals that the SVB estuary increased its total occupied area with respect to 1957 ( Figure 6). This was mainly due to an increase in anthropic features, which grew until the 2000s and stabilized afterwards. As expected for its geographical constraints (see the study sites section), the rest of the SVB subsystems mostly maintained their surfaces, as the possibilities to change their perimeters are limited. The GUA estuary, although it seems to be a fluctuating system (Figure 6), was the most stable case, showing a similar surface at the initial and the final years. The changes were related to variations in the sand bars in the intertidal zone, which are very dynamic features, the emergence of which depends on tidal and seasonal conditions. Lastly, the Ebro River Delta seems to be the most vulnerable, showing a decrease in the total occupied surface larger than 10%. This trend was supported by strong changes in the erosion process on San Antonio Island, which was counterbalanced with growth on the left side of the estuary. The total areas on the GUA ( Figure 3) and EBR (Figure 4) estuaries are dependent on the state of the sublittoral sand bars attached to the coast and the tide level at the time of the photo (only relevant for GUA as EBR is a microtidal environment). Seasonality is also relevant for the degree of the feature development; however, it is difficult to include this variable in the general calculation of areas.

Changes in Occupied Surface
Changes on the perimeter of occupied surfaces can be used as a proxy for net size fluctuations with time ( Figure 6). This reveals that the SVB estuary increased its total occupied area with respect to 1957 ( Figure 6). This was mainly due to an increase in anthropic features, which grew until the 2000s and stabilized afterwards. As expected for its geographical constraints (see the study sites section), the rest of the SVB subsystems mostly maintained their surfaces, as the possibilities to change their perimeters are limited. The GUA estuary, although it seems to be a fluctuating system (Figure 6), was the most stable case, showing a similar surface at the initial and the final years. The changes were related to variations in the sand bars in the intertidal zone, which are very dynamic features, the emergence of which depends on tidal and seasonal conditions. Lastly, the Ebro River Delta seems to be the most vulnerable, showing a decrease in the total occupied surface larger than 10%. This trend was supported by strong changes in the erosion process on San Antonio Island, which was counterbalanced with growth on the left side of the estuary.
When analysing more specifically by Category 1 (Table 1), the percentage of changes allows for the identification of the most dynamic subsystems (Figure 7). To facilitate this identification, the tidal flat/wetland category was divided into Category 2 units (saltmarshes/vegetated wetland and mudflats/submerged meadows). The analysis revealed that, in the SVB estuary, the main changes affected the muddy environments, whereas the saltmarshes' extension decreased >20% (corresponding to approximately 85 ha), which gradually became mudflats (Figure 7). For the GUA and EBR estuaries, the changes mostly affected the sandy environments. In particular, the GUA estuary showed a stable saltmarsh, with a slight increase with time due to the protective role played by the sandy hooks upstream, and the mudflats almost maintained their initial surface. In the natural part of the estuary, new sandy features, especially dunes, increased. However, this effect was masked by the loss associated with anthropogenic constructions on the Isla Canela side (Figure 7). When analysing more specifically by Category 1 (Table 1), the percentage of changes allows for the identification of the most dynamic subsystems (Figure 7). To facilitate this identification, the tidal flat/wetland category was divided into Category 2 units (saltmarshes/vegetated wetland and mudflats/submerged meadows). The analysis revealed that, in the SVB estuary, the main changes affected the muddy environments, whereas the saltmarshes' extension decreased >20% (corresponding to approximately 85 ha), which gradually became mudflats (Figure 7). For the GUA and EBR estuaries, the changes mostly affected the sandy environments. In particular, the GUA estuary showed a stable saltmarsh, with a slight increase with time due to the protective role played by the sandy hooks upstream, and the mudflats almost maintained their initial surface. In the natural part of the estuary, new sandy features, especially dunes, increased. However, this effect was masked by the loss associated with anthropogenic constructions on the Isla Canela side (Figure 7).   When analysing more specifically by Category 1 (Table 1), the percentage of changes allows for the identification of the most dynamic subsystems (Figure 7). To facilitate this identification, the tidal flat/wetland category was divided into Category 2 units (saltmarshes/vegetated wetland and mudflats/submerged meadows). The analysis revealed that, in the SVB estuary, the main changes affected the muddy environments, whereas the saltmarshes' extension decreased >20% (corresponding to approximately 85 ha), which gradually became mudflats (Figure 7). For the GUA and EBR estuaries, the changes mostly affected the sandy environments. In particular, the GUA estuary showed a stable saltmarsh, with a slight increase with time due to the protective role played by the sandy hooks upstream, and the mudflats almost maintained their initial surface. In the natural part of the estuary, new sandy features, especially dunes, increased. However, this effect was masked by the loss associated with anthropogenic constructions on the Isla Canela side (Figure 7).  The shoreline sandy environments in the GUA estuary include large areas of beach and sand bars (Category 2). These features experienced huge variations, what explains the corresponding temporal curve in the general trend of the system (Figure 6). In the case of the EBR estuary, the sandy subsystems (shoreline sandy environments and dunes) increased in the last few years by almost 3% (Figure 7). As described for the GUA estuary, this gain was partially masked and counterbalanced by the strong erosion suffered on the eastern side of the estuarine system. Despite this, the sandy environments increased by between 10 and 20 ha during the last 5 years (Figure 8). Doing a similar analysis for the GUA estuary (Figure 8), it can be concluded that the saltmarsh has an equilibrated balance between sediment supplies and hydrodynamic processes. The interaction of external and internal agents favours sediment remobilization [32], resulting in an increment in the surface of around 0.3 ha/year (Figure 8). An equilibrated balance of sediment inputs allows the tidal subsystem to prograde (horizontally) and aggrade (vertically) following a natural behaviour [13]. The large availability of sandy sediments in this comes from the river itself, which erodes deeply weathered granites in its basin, and by the erosion of westward sandy cliffs, transported to the estuary by the prevailing eastward littoral drift [32,69]. Thus, the sediment supply, the coastal circulation system, the areas of wave refraction and the resulting shore currents are crucial to achieve a coastal equilibrium and to maintain the balance of the Guadiana River estuary [17,32].
The relationship between sandy and muddy environments is clearly illustrated in the SVB and the GUA estuaries. In the last one, the progradation of the sandy spits reduces the input of wave energy and leads to the generation of new sheltered muddy environments (indicated as the mudflat* in Figure 3). This also occurred in El Garxal in the EBR estuary ( Figure 4). The oldest information available for the EBR estuary ( Figure 5) allows dating the formation of El Garxal between 1956 and 1984 (no orthophotographs are available between those years), sheltered by sandy environments, namely sand bars and San Antonio Island. Together with the action of waves and winds, this island/hemi-delta, mainly formed by sandy features, facilitated the formation of littoral/dune ridges on the other side of the channel, where they were controlled by wave and wind action [59]. Nevertheless, this island has been gradually eroded through decades by wave action and also because of the lack of river sediment input, as the discharges of this river have been intensively regulated since the beginning of the 20th century [35]. As a consequence, San Antonio Island is presently recording erosion rates of almost 30 m/year in the river mouth area [35], causing the strong negative rates in overwash deposits shown in Figure 8 (more than -40 ha/year), which were recorded To maintain the same x-scale for the three study sites, the total rates of change in the relict littoral ridge, relict dune ridge and overwash deposit in the EBR, which showed rates of change greater than ± 2 ha, are indicated numerically on each bar.
The channels are also very dynamic features in the three estuarine systems (Figure 7), whose differences due to the tidal regimes (meso-and microtidal) seem to affect the magnitude of the oscillating pattern. For the mesotidal sites, these variations were related to the effect of the tidal level on the periodical emergence of certain features (especially in the GUA estuary) at the moment the orthophoto was taken.

Accuracy Assessment
The accuracy assessment was performed on the SVB estuary, as an example, and compares field data with test pixels on the defined features ( Table 1). The most recent cartography (2017) and field and UAS data (2018) made it possible to estimate an overall accuracy (P c ) of 0.87 and a k coefficient of 0.85. Incorrectly classified limits were more abundant between beach and embryo dune features, but also between saltmarshes and mudflats ( Table 2). The differences in elevation and the presence of vegetation were both difficult to detect in these environments. Table 2. Confusion matrix resulting from classifying training set pixels in the San Vicente de la Barquera estuary. Shaded cells indicate areas with larger confusion coefficients between defined features. Main diagonal, in bold, indicates features confirmed by field inspection (reference points). Final column and row, in bold, show P Ui and P Pi for the corresponding feature. Column numbers match with row numbers.

Discussion
The hierarchical classification of the estuarine features (Table 1) has allowed the criteria for drawing up eco-geomorphological maps to be unified, allowing the comparison of mapped features, which is necessary to assess changes in occupied surface in accordance with the European Union guidelines The present work reveals evolutionary trends in spatio-temporal variations and eco-geomorphological processes for three Iberian estuaries under contrasting hydrodynamic conditions. By quantifying the total areas and mapping main subsystems in each study area (Figures 6 and 7), the most vulnerable subsystems can be identified, as well as possible drivers and potential consequences for the entire system.
At the level of the entire system ( Figure 6), results illustrate that estuaries with strong physical constrictions, such as the SVB, may modify the area occupied, mainly by changing some of their features to the detriment of others, suggesting that in this type of environment the connection between subsystems is the most important driver of long-term dynamics and that an integral management of the different subsystems is fundamental.
On other estuaries lacking such physical constrains, such as GUA and EBR, temporal changes are more evident. However, changes associated with intertidal sandy subsystems have to be analysed with caution, as their presence is strongly affected by tidal and seasonal conditions. In these systems, the boundary conditions may generate contrasting processes (erosion vs. sedimentation) that may mask the magnitude of the dynamic processes. Therefore, in any of these cases, a smaller-spatial-scale analysis will help to evaluate the processes involved in the dynamics of the whole system.
Identification of the most vulnerable features sometimes requires the analysis of the features in Categories 2 and 3 (Figure 8), supporting the need for a cartography with more detail (Figures 2-4) to identify the most vulnerable features in the study areas and their associated habitats, according to European guidelines (year-by-year changes in every feature are detailed in Appendix A; Tables A1-A3). A detailed cartography reveals that main changes in San Vicente de la Barquera estuary ( Figure 2) were related to changes in muddy environments (Figure 8). The smooth change in total occupied area was due to the growth of the mudflats to the detriment of the saltmarshes. In fact, the saltmarshes in the SVB estuary suffered a strong negative trend, an around 1.4 ha/year loss (Figure 8), which implies a surface loss of around 20% since 1956. The persistence of this negative trend may lead to the disappearance of the saltmarsh and its ecosystem services, with the risk of transforming the system into a vast muddy intertidal plain partially covered by seagrasses and macroalgae in some points. The loss of the saltmarshes, and their associated communities, will reduce the capacity of the system to retain sediment [60]. One possible cause of this pattern may be the lack of sediment input into the estuary due to the small fluvial supply. Previous studies have shown that the sediment balance of these systems, both internally and externally, may enable the adjustment of the bed elevation to the sea level rise (SLR) [61]. However, if the source of sediment is reduced and the system is unable to make this adjustment, it will cause the migration of the vegetation to higher elevation areas when available [62]. In particular, for this area, [63] described SLR rates of 2.08 ± 0.33 mm yr −1 from 1943 to 2004 for Santander (nearest tidal gauge), with an accelerated rate at the end of the analysed period [14]. This, together with the impossibility of landward migration due to the orography of this area, promotes the squeeze of the saltmarsh, reducing the probability of plant survival [64]. Hence, the persistence over time of the saltmarsh depends on a balance between the supply of sediment and nutrients and external energetic forcing, including wave energy, storm surges and tidal flooding [65][66][67][68].
Doing a similar analysis for the GUA estuary (Figure 8), it can be concluded that the saltmarsh has an equilibrated balance between sediment supplies and hydrodynamic processes. The interaction of external and internal agents favours sediment remobilization [32], resulting in an increment in the surface of around 0.3 ha/year ( Figure 8). An equilibrated balance of sediment inputs allows the tidal subsystem to prograde (horizontally) and aggrade (vertically) following a natural behaviour [13]. The large availability of sandy sediments in this comes from the river itself, which erodes deeply weathered granites in its basin, and by the erosion of westward sandy cliffs, transported to the estuary by the prevailing eastward littoral drift [32,69]. Thus, the sediment supply, the coastal circulation system, the areas of wave refraction and the resulting shore currents are crucial to achieve a coastal equilibrium and to maintain the balance of the Guadiana River estuary [17,32].
The relationship between sandy and muddy environments is clearly illustrated in the SVB and the GUA estuaries. In the last one, the progradation of the sandy spits reduces the input of wave energy and leads to the generation of new sheltered muddy environments (indicated as the mudflat* in Figure 3). This also occurred in El Garxal in the EBR estuary ( Figure 4). The oldest information available for the EBR estuary ( Figure 5) allows dating the formation of El Garxal between 1956 and 1984 (no orthophotographs are available between those years), sheltered by sandy environments, namely sand bars and San Antonio Island. Together with the action of waves and winds, this island/hemi-delta, mainly formed by sandy features, facilitated the formation of littoral/dune ridges on the other side of the channel, where they were controlled by wave and wind action [59]. Nevertheless, this island has been gradually eroded through decades by wave action and also because of the lack of river sediment input, as the discharges of this river have been intensively regulated since the beginning of the 20th century [35]. As a consequence, San Antonio Island is presently recording erosion rates of almost 30 m/year in the river mouth area [35], causing the strong negative rates in overwash deposits shown in Figure 8 (more than −40 ha/year), which were recorded only for San Antonio Island. The possible disappearance of San Antonio Island in future years could expose El Garxal to eastern, energetic waves and very probably would cause its disappearance and that of its associated features and, as consequence, the elimination of the current river mouth complex, i.e., the two hemi-deltas.
Sandy/muddy environment relationships have been also described in other estuarine systems [5], and it is widely described that the weakening of a protective dune system after external impacts, natural or artificial, has consequences for local saltmarshes' functioning [70,71]. Building and urbanizing on the upper zones of shoreline sandy environments and dunes interrupt their dynamics. This pressure, in addition to the effects of global sea-level rise [72] and other effects of climate change on the coasts [73], makes these features prone to disappearing [74,75].
To summarize, understanding the processes connecting estuarine subsystems is important to enhance current schemes in ecosystem connectivity, as this is key to the long-term functioning of ecosystems [64] and references therein. Complex system transitions, self-organization and connectivity over the years are crucial concepts to understand the dynamics and adaptation of coastal habitats [5,76,77]. This work shows how small-scale features have effects on the functioning of the estuary [78] and mapping is the best tool for its monitoring. Thus, detailed maps, in addition to the adaptation to European and national classifications and regulations, support a better understanding on the functioning of the estuarine systems, including the inner processes occurring therein [16]. Reducing the scale of working features provides stakeholders with the basis for drawing up management plans adapted to their specific systems. This work demonstrates that the reduction of the spatial scale helps to understand the causes of change and facilitates pressures to be identified in the short, medium and long term, and therefore improves the efficiency of estuarine management practices [79].
At present, increasingly available high-quality image datasets combined with automatic classification, i.e., on a pixel-by-pixel basis, are used to analyse changes in different landscapes [18,60,80]. However, for retrospective studies, it is necessary to use aerial orthophotograph libraries. The quality of the aerial orthophotos is an issue to be considered as it may limit accuracy in the identification of feature limits ( Table 2). A comparison between orthophotos and UAS images (Figure 9) illustrates this issue. It is clear that a higher spatial resolution (UAS images) facilitates the visual distinction of the feature limits (see Figure 9, left panes). However, features with microtopograhies, such as saltmarsh plant areas versus seagrass ones (Figure 9, right panes), require the availability of accurate digital elevation models, together with field topographic data, to improve the quality of future analyses.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 19 of 26 Figure 9. Examples of limitations to feature classification arising in the accuracy assessment for the SVB estuary. The images show the differences between acquisition methods. Left images show the limit between beach and embryo dune habitats; Right images show the limit between saltmarsh (darkest green) and mudflat (covered by macro-algae and seagrasses) features. Red dots indicate test points in the field.
The main error sources in previous studies [81] have been related to the quality of the orthophoto and the bias of the analyst. Error sources may differ between cartographies. Therefore, accuracy assessment should be specific to each cartography, and ground truthing campaigns should be performed the same year as the orthophoto to make conditions comparable. With limitations, this study goes beyond previous studies, analysing the behaviour of the estuary as a whole over a period long enough to identify clear trends between geomorphological and ecological subsystems.

Conclusions
The study of estuaries along the Iberian coast under different oceanographic and management conditions has made it possible to establish a common classification for estuaries at the system and subsystem levels, which aims to lay the foundations for the future development of estuarine maps. This classification is the first step for standardizing a methodology to assess the conservation state of Figure 9. Examples of limitations to feature classification arising in the accuracy assessment for the SVB estuary. The images show the differences between acquisition methods. Left images show the limit between beach and embryo dune habitats; Right images show the limit between saltmarsh (darkest green) and mudflat (covered by macro-algae and seagrasses) features. Red dots indicate test points in the field.
The main error sources in previous studies [81] have been related to the quality of the orthophoto and the bias of the analyst. Error sources may differ between cartographies. Therefore, accuracy assessment should be specific to each cartography, and ground truthing campaigns should be performed the same year as the orthophoto to make conditions comparable. With limitations, this study goes beyond previous studies, analysing the behaviour of the estuary as a whole over a period long enough to identify clear trends between geomorphological and ecological subsystems.

Conclusions
The study of estuaries along the Iberian coast under different oceanographic and management conditions has made it possible to establish a common classification for estuaries at the system and subsystem levels, which aims to lay the foundations for the future development of estuarine maps. This classification is the first step for standardizing a methodology to assess the conservation state of estuarine systems in Spain, also applicable in Portugal in the case of the GUA estuary.
This work describes the general relationships between the subsystems composing an estuary, namely, sandy environments (shoreline sandy environments and dune subsystems), muddy environments (saltmarshes/vegetated wetland and mudflats/submerged meadow subsystems), with a distinction between meso-and micro-tidal environments, and the drainage network, comprising the main channel and the associated channel network. The main processes shaping the inner features are also described, highlighting the influence they have on the others. Thus, the drainage network acts as the main connection vector for the entire system, impulse by fluvial and tidal currents and local oceanographic conditions. The hydrophysical conditions control the nature of the features and their associated communities. Sandy environments provide a sediment budget to the estuarine system via the channel network, which promotes the exchange of sediment, organic matter and energy in the entire system. Similarly, in muddy environments, the drainage network subsystem plays a crucial role in developing and maintaining these tide-dependent environments. Additionally, sandy features serve as protection for the muddy ones, where the sandy spits promote the development of saltmarshes and wetlands, respectively. These results suggest that the degradation of one of the subsystems may influence the functioning of the others, harming the ecosystem services provided by the whole system.
Specifically, the San Vicente de la Barquera estuary showed a strong negative trend in saltmarsh surface as a consequence of a possible reduction in sediment input in the estuary and sea level rise. The Guadiana River estuary presented variations in shape over the studied years due to the dynamics of external and intertidal sand bars but with a constant sediment input which allowed for the subsystems progradation, both vertically and horizontally, and also the formation of new features. This system seems to be the more stable one. The Ebro River estuary exhibited a strong positive trend on sandy environments on the left side of the channel but a negative trend on the right side, which will undoubtedly affect the entire system in future years if present erosion rates continue.
Finally, the accuracy assessment for the San Vicente de la Barquera estuary, as an example, has allowed the identification of accuracy weaknesses when mapping estuarine environments, highlighting the need for high-resolution imagery and elevation data series to monitor estuarine changes.
In summary, the development of thematic eco-geomorphological maps improves the capacity to assess the resilience of estuarine systems under present and future pressures by integrating historical and present information on estuaries. In future work, the combination of aerial orthophotos, which provide the historical data for the study of landscapes, with other methods, such as hyperspectral information from airborne and satellite sensors, will facilitate the reconstruction and monitoring of estuarine changes at the landscape level, reducing time and effort. Funding: This research was funded by an EMAS grant (Earth and Marine Sciences joint doctorate programme between the University of Cádiz (Spain) and the University of Ferrara (Italy)). This work is a contribution to the Andalusia Research Groups (P.A.I.) RNM328 and RNM214.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.