Between the tides: Modelling the elevation of Australia's exposed intertidal zone at continental scale

The intertidal zone represents a critical transition between marine and terrestrial ecosystems, supporting a complex mosaic of highly productive and biologically diverse habitats. However, our understanding of these important coastal environments is limited by a lack of spatially consistent topographic data, which can be extremely challenging and costly to obtain at continental-scale. Satellite remote sensing represents an important resource for monitoring extensive coastal zones. Previous approaches to modelling the elevation of the intertidal zone using earth observation (EO) data have been restricted to small study regions or have relied on manual image interpretation, thus limiting their ability to be applied consistently over large geographic extents. In this study, we present an automated open-source approach to generate satellite-derived elevation data for over 15,387km 2 of intertidal terrain across the entire Australian coastline. Our approach combines global tidal modelling with a 30-year time series archive of spatially and spectrally calibrated Landsat satellite data managed within the Digital Earth Australia (DEA) platform. The resulting National Intertidal Digital Elevation Model (NIDEM) dataset provides an unprecedented three-dimensional representation of Australia's vast exposed intertidal zone at 25m spatial resolution. We validate our model against LiDAR, RTK GPS and multibeam bathymetry datasets, finding that modelled elevations are highly accurate across sandy beach (±0.41m RMSE) and tidal flat environments (±0.39m RMSE). Model performance was least accurate (±2.98m RMSE) within rocky shores and reefs and other complex coastal environments with extreme and variable tidal regimes. We discuss key challenges associated with modelling intertidal elevation including tidal model performance and biased observations from sun-synchronous satellites, and suggest future directions to improve the accuracy and utility of continental-scale intertidal elevation modelling. Our model can be applied to tidally-influenced coastal environments globally, addressing a key gap between the availability of sub-tidal bathymetry and terrestrial elevation data.


Introduction
The intertidal zone -the area of coastline periodically exposed and inundated by tides -represents a critical transition between marine and terrestrial ecosystems. Intertidal environments support a complex mosaic of extremely productive and biodiverse habitats ranging from extensive tidal mudflats, sandy beaches, fringing coral reefs and steep rocky cliffs (Banks et al., 2005;Luijendijk et al., 2018). Due to the influence of tidal processes, organisms inhabiting the intertidal zone are typically adapted to extremely dynamic conditions and display strong zonation by elevation along the environmental gradient from permanently to occasionally inundated terrain (Bearup and Blasius, 2017). Because of their high ecological diversity and productivity, intertidal zones serve as key feeding grounds for many endangered shorebird species that use them as critical 'stop-over' points while undertaking cross-continental migrations (Murray et al., 2015;Xia et al., 2017). In addition to their ecological value, intertidal zones also provide many economically significant ecosystem services, including nutrient cycling and carbon storage (Chmura et al., 2003;Billerbeck et al., 2006), storm surge protection (Temmerman et al., 2013), and natural resources for recreational and commercial use (Barbier et al., 1997;Chen et al., 2016). Intertidal zones, however, are also among the world's most vulnerable and threatened ecosystems, with land reclamation, changes in river sediment balances and coastal erosion representing key threatening processes responsible for a global reduction in intertidal extent of up to 16% between 1984 and 2016 (Murray et al., 2018). T low relief, intertidal zones are also likely to be disproportionately affected by global sea-level rise, which is expected to accelerate throughout the 21st century (Galbraith et al., 2002;Kirwan and Megonigal, 2013;Li and Gong, 2016).
Understanding, predicting and managing the impacts of these processes on intertidal ecosystems requires detailed data on the distribution and structure of these habitats across ecologically relevant spatial extents. However, the three-dimensional topography of the intertidal zone remains poorly mapped globally (Eakins and Grothe, 2014;Tseng et al., 2017). Due to the impermeability of water to radar transmission and a lack of repeated observations over tidally-influenced terrain, global digital elevation models (DEMs) produced using Synthetic Aperture Radar (e.g. TerraSAR-X/TanDEM-X WorldDEM and Shuttle Radar Topography Mission DEM or SRTM) or stereo-pair optical imagery (e.g. Aster GDEM) are typically restricted to the terrestrial domain (Eakins and Grothe, 2014;Tseng et al., 2017). Similarly, intertidal zones are regularly omitted from bathymetric models produced using acoustic techniques due to the difficulty of surveying safely from vessels in shallow coastal waters (Hogrefe et al., 2008;Eakins and Grothe, 2014;Weatherall et al., 2015). When combined with hazardous ground survey conditions caused by dynamic tidal processes, this has resulted in a significant gap in the elevation data available across the land-sea interface (Hogrefe et al., 2008;Eakins and Grothe, 2014). Airborne Light Detection and Ranging (LiDAR) bathymetry surveys have shown promise as an approach to address this gap, being capable of rapidly generating accurate, high-resolution bathymetry data at depths of up to 70 m in clear water (Su et al., 2008). However, the reliability of these systems can be strongly influenced by turbidity and breaking white water that frequently affect shallow coastal zone waters, and their high acquisition cost usually limits applications to local-or regional-scales (Su et al., 2008;Gao, 2009;Klemas, 2011).
Temporal 'waterline' methods based on satellite remote sensing represent an alternative approach to modelling the elevation of the intertidal zone (Mason et al., 1997;Chen and Rau, 1998). The rise and fall of the ocean can be used to describe the three-dimensional topography of the coastline by mapping the location of the waterline as a series of topographic contours that cover a range of tidal stages (Zhao et al., 2008;Tseng et al., 2017). Assuming that each waterline represents a constant elevation relative to mean sea level (MSL), these contours can be tagged with tide heights and then interpolated to produce a DEM covering the elevation range between the highest and lowest observed tide (Ryu et al., 2008). Satellites such as the United States Geological Survey (USGS) Landsat mission have continuously monitored coastal zones globally since 1972, providing the temporal depth and resolution required to obtain dense observations across the full tidal range (Boak and Turner, 2005;Gens, 2010). Accordingly, temporal stacks of satellite imagery have been combined with tidal modelling to produce tidally-tagged time series of the coastline for DEM generation (e.g. Mason et al., 1997;Chen and Rau, 1998;Ryu et al., 2008;Zhao et al., 2008;Chen and Chang, 2009;Liu et al., 2013aLiu et al., , 2013bXu et al., 2016;Tseng et al., 2017). However, due to challenge of extracting waterline contours from large numbers of remotely-sensed images, previous approaches have extracted waterlines from a limited selection of images using manual digitisation and visual interpretation (e.g. Chen and Rau, 1998;Zhao et al., 2008;Liu et al., 2013b;Chen et al., 2016). Although these manually digitised contours can produce highly accurate DEMs (e.g. less than 0.4 m RMSE compared to LiDAR; Liu et al., 2013b), this manual process introduces subjectivity, is impractical to apply at a continental-scale, and is restricted by availability of high quality observations covering the entire tidal range.
There is a recognised need for more objective and robust approaches to waterline extraction and intertidal DEM generation that can be applied consistently across space and time (Boak and Turner, 2005). Such approaches are likely to require leveraging the full temporal record of available global-extent earth observation (EO) data. While satellite data with long temporal records such as the Landsat archive were formerly prohibitively expensive to apply at continental-scale, the USGS freedata policy in 2008 has significantly lowered barriers to obtaining EO data (Woodcock et al., 2008;Wulder et al., 2012). More recently, the creation of archives of analysis-ready data or ARD (Dwyer et al., 2018) combined with high-performance computing platforms have provided unprecedented access to petabytes of geometrically and spectrally calibrated satellite imagery . ARD archives and analysis platforms such as Digital Earth Australia (Dhu et al., 2017), Google Earth Engine (Gorelick et al., 2017) and the upcoming Copernicus Data and Information Access Services provide satellite observations that can be analysed and compared consistently across space and time. This has driven a new wave of scientific applications that leverage multiple decades of data to analyse extremely large areas of the Earth's surface (e.g. Mueller et al., 2016;Hermosilla et al., 2017;Bugnot et al., 2018;Egorov et al., 2018;Luijendijk et al., 2018). These developments make it practical for the first time to implement automated waterline extraction and intertidal elevation modelling at the continental-scale.
In this study, we present an open-source workflow for deriving elevation data for the intertidal zone of Australia. We leverage the full 30-year archive of analysis-ready Landsat data managed within the Digital Earth Australia (DEA) platform that provides spatially and spectrally calibrated EO data, enabling time-series analysis on a perpixel basis across the entire Australian continent. We combine this archive with a newly developed multi-resolution tidal modelling framework that accurately associates each satellite observation with modelled tide heights. The resulting National Intertidal Digital Elevation Model (NIDEM; Bishop-Taylor et al., 2018) is a continental-scale dataset providing the first 25 m resolution, three-dimensional representation of Australia's exposed intertidal environments including tidal flats, sandy beaches and shores, and rocky shores and reefs. We anticipate that NIDEM will complement existing intertidal extent products and support a new suite of use cases that require a more detailed understanding of the three-dimensional topography of the intertidal zone, such as hydrodynamic modelling, coastal risk management and ecological habitat mapping projects.

Relative intertidal extents
The foundation of NIDEM is the continental-scale Intertidal Extents Model (ITEM v1.0) developed by Sagar et al. (2017), derived from a 30year time series archive of Landsat observations. The Landsat archive is managed within the Digital Earth Australia (DEA) platform, combining high performance computing with the concept of spatiotemporally consistent ARD (Dhu et al., 2017). These data are processed to surface reflectance utilising a standardised atmospheric and geometric correction workflow to enable operational analysis . The ITEM process is based on sorting all observations in the Landsat archive by tide height, binning observations into ten percent intervals of the observed tidal range, then mapping the typical location of the waterline across a range of tidal stages using Normalised Difference Water Index composite images (NDWI;McFeeters, 1996). NDWI is a remote sensing index designed to detect open water by taking advantage of the high reflectance of visible green light and low reflectance of near-infrared radiation (NIR) by water, and the high reflectance of NIR by dry soil and terrestrial vegetation: Although NDWI has been used extensively to monitor the intertidal zone (e.g. Murray et al., 2012;Dhanjal-Adams et al., 2016;Fan et al., 2018;Luijendijk et al., 2018), NIR bands can be affected by white water near the land-water boundary (Kelly and Gontz, 2018;Pardo-Pascual et al., 2018). Other indices such as the Modified Normalised Difference Water Index (MNDWI; Xu, 2006) have shown promise for monitoring the intertidal zone by utilising short-wave infrared (SWIR) in place of NIR (e.g. Wang et al., 2018;Xu, 2018). However, SWIR-based indices have been found to significantly misrepresent the location of the waterline across exposed tidal flats where water remains during ebb tides (Ryu et al., 2002(Ryu et al., , 2008. To ensure waterlines could be consistently extracted from imagery across tidal stages and intertidal environments that included extensive tidal flats, NDWI was selected and applied to all Landsat observations within each ten percent tidal interval. These individual NDWI layers were then combined into composites for each tidal interval by taking the median NDWI value per pixel, and classified to produce water vs. non-water layers that represented typical waterlines locations from the lowest to the highest observed tides (Fig. 2a).
In this study, we further develop the initial model (ITEM v1.0) described in Sagar et al. (2017) by leveraging improvements in the tidal modelling framework that underpins the ITEM modelling process. The tidal modelling used in ITEM v1.0 was constrained spatially into 1°× 1°r esolution image cells (approximately 110 × 110 km at the equator), with a single tidal height assigned to the centre of each Landsat image based on its time of acquisition. This spatially consistent image cell grid implicitly assumed that tidal heights did not vary across the 1°× 1°cell extent, and that a single modelled tide height per cell could adequately reflect complex tidal dynamics operating at a range of spatial scales. This assumption resulted in sometimes severe model discontinuities at some cell boundaries and poor modelling performance in complex estuaries and other coastal areas characterised by high tidal flux (Sagar et al., , 2017. ITEM v2.0 was developed by replacing the spatially consistent image cell grid with a multi-resolution tidal framework developed by Sagar et al. (2018). The framework uses partitioning methods to allow spatial variability in the tidal model to drive the size and locations of a Voronoi polygon mesh. The 306 resulting tidal modelling polygons ( Fig. 1) are then used as analysis units for the ITEM modelling process, with tide height predictions defined at the nodes of each Voronoi cell. We used the Oregon State University Tidal Prediction Software (OTPS) TPX08 model Erofeeva, 2010, 2002) to predict tide heights. OTPS tidal modelling has been used successfully for continental-scale intertidal modelling applications in both Asia and Australia (Murray et al., 2012;Dhanjal-Adams et al., 2016). The model consists of a multi-resolution bathymetric grid with a 1/6°resolution (∼18 × 18 km at the equator) solution in the global open ocean and a 1/30°local resolution (∼4 × 4 km) solution to improve modelling in complex shallow-water environments, and has been validated to ∼12 cm root mean square error (RMSE) misfit against the Australian Hydrographic Office AusTides tide gauge records (Rogers et al., 2017). The study area for ITEM v2.0 and NIDEM has been extended to cover the Great Barrier Reef and the entire Australian coastline, including Tasmania. The use of the multi-resolution tidal model has produced significant improvements in the coverage and resolution of the relative intertidal extents model, most notably in the offshore regions of northern Australia and the Kimberley coast in north-western Australia.

Unfiltered absolute elevation
ITEM v2.0 details the relative extent of the intertidal zone at intervals of the observed tidal range. As such, it provides topography of the exposed intertidal surface but not an absolute elevation measure. To derive absolute elevations, we used the find_contours function from scikit.measure (Van der Walt et al., 2014) to extract waterline contours ( Fig. 2b) along the boundary of each of the ten percent tidal interval boundaries in ITEM v2.0 (Fig. 2b). This contour extraction method uses the 'marching squares' algorithm (Lorensen and Cline, 1987) to identify precise contour boundaries in a two-dimensional array by linearly interpolating between adjacent pixel values. For each interval boundary,

Fig. 1. Map of the continental-scale
Australian study area, showing the location of the 306 tidal modelling polygons and validation sites for tidal flats (9 sites; circles), sandy beaches and shores (5 sites; triangles), and rocky shores and reefs (7 sites; stars). Labels highlight locations referred to in this work. Image underlay is a composite of all Landsat 8 satellite observations across Australia .
we computed the median tidal height of the ensemble of all Landsat observations originally used to derive the corresponding median composite NDWI layer, and assigned this height as the contour's z-value. To convert contours to continuous elevation data rasters, we first extracted x, y and z point data for each contour vertex, and used these points as inputs for interpolation. We used the griddata interpolation function from scipy.interpolate (Jones et al., 2014), and selected the 'linear' interpolation method to ensure that interpolated elevation values preserved the tidal interval boundaries of ITEM v2.0. This interpolation method computes a Triangulated Irregular Network or Delaunay triangulation of the input data using the Quickhull algorithm (Barber et al., 1996), before performing linear barycentric interpolation on each triangle to estimate new values for each pixel (Fig. 2c). We set the output resolution of the interpolation to 25 × 25 m, producing twodimensional elevation arrays that matched the cell size and extent of the original ITEM v2.0 layers and input Landsat observations. As elevations for the lowest and highest ITEM v2.0 intervals could not be correctly interpolated because they had no lower or upper bounds, we constrained the interpolated elevation arrays to the observed intertidal zone by masking out pixels located within consistently inundated (ITEM v2.0 interval 0) and consistently non-inundated terrain (interval 9). This resulted in a set of 306 'unfiltered' intertidal elevation datasets with elevations in metre units relative to modelled MSL (approximately equivalent to the Australian Height Datum or AHD). To facilitate future re-analysis such as the application of alternative interpolation methods, waterline contours were exported as shapefiles.

Filtered absolute elevation
The accuracy of intertidal elevation modelling approaches based on waterline extraction and tidal modelling is dependent on the ability of a tidal model to correctly assign and sort satellite observations by tidal height (Liu et al., 2013b;Sagar et al., 2017). This assumption may not hold in areas that have undergone significant geomorphological change across the 30-year time series, or where modelled tidal heights differ from actual tidal heights due to inherent model error (Ryu et al., 2008;Sagar et al., 2017). To identify potentially invalid elevation values, we used a 'confidence' layer developed as part of ITEM v2.0. The ITEM confidence layer is based on the per-pixel variance in NDWI values for the ensemble of images used to produce each ten percent tidal interval composite. We used a conservative maximum threshold of 0.25 NDWI standard deviation to mask out over 1980 km 2 of intertidal pixels (5.5% of the total 'unfiltered' dataset) where tidal processes poorly explained inundation patterns across the 30-year time series.
The robust median compositing method used to combine NDWI images in ITEM 2.0 ensured that the input relative intertidal extent layers were relatively free of artefacts including outliers, poorly removed cloud edges and sunglint that commonly affect coastal remote sensing imagery (Sagar et al., 2017;White et al., 2014). However, composite layers produced from multiple satellite observations can still be susceptible to artefacts in regions with high cloud cover or fewer observations in the Landsat archive (Flood, 2013;White et al., 2014;Roberts et al., 2017). For NIDEM, this typically manifested as pixels located either inland of the coastal zone or in areas of deeper water incorrectly mapped as intertidal terrain in data poor areas of the Australian coastline. To remove these artefacts, we restricted NIDEM layers to a 50 m elevation range centred on MSL. This large range relative to the maximum Australian tidal range of over 11 m (Solihuddin et al.,

Fig. 2.
The process followed to extract and interpolate waterline contours for NIDEM for an example area (Forestier Bay in the Pilbara region of Western Australia). Boundaries of each (a) ten percent interval of the observed tidal range from ITEM v2.0 (colours from blue to red indicate pixels inundated at increasingly high tide) were used to (b) extract contours depicting the typical location of the land-water boundary across the 30-year Landsat time series. These contours were assigned the median of all modelled tide heights attributed to the ensemble of Landsat images used to generate each tidal interval. The resulting tidally-tagged contours were used to generate (c) continuous elevation surfaces using triangulated irregular network (TIN) interpolation. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) 2016) was selected to remove obvious false positives without eliminating true intertidal terrain. We used elevation data from Australia's SRTM-derived 1 Second Digital Elevation Model (Gallant et al., 2010) to mask out approximately 402 km 2 of terrain with elevations greater than 25 m above MSL (1.1% of the total 'unfiltered' dataset). Bathymetry data from the national 250 m Australian Bathymetry and Topography Grid (Whiteway, 2009) and the 30 m gbr30 and nthaus30 highresolution depth models for the Great Barrier Reef and Northern Australia (Beaman, 2018a(Beaman, , 2018b were used to mask out pixels located at depths of greater than −25 m in all three datasets. This removed over 18,325 km 2 of misidentified intertidal terrain (50.7% of the total 'unfiltered' dataset) which was largely concentrated (i.e. 9571 and 4574 km 2 , or 26.5 and 12.7%) in two data-poor tidal modelling polygons located off the Archipelago of the Recherche in southern Western Australia. All layers were reprojected to the resolution (25 × 25 m) and projection system (GDA94 Australian Albers, EPSG:3577) of the NIDEM layers using bilinear resampling prior to masking. The resulting masked layers were exported as "filtered" elevation datasets, while the combined mask (i.e. ITEM confidence, bathymetry and elevation limits) was exported as "mask" layers to facilitate re-analysis through the application of customised masking criteria.

Tidal range coverage
To evaluate the representativeness of NIDEM data compared to the full tidal range, we compared the spread of tidal heights coincident with the input Landsat imagery against the full range of modelled tide heights present within each tidal modelling polygon. We calculated three indices: spread (the proportion of the full modelled tidal range observed by Landsat), low tide offset (the proportion of the lowest tidal heights not observed by Landsat) and high tide offset (the proportion of the highest tidal heights not observed by Landsat). The variations in these indices in different regions relate to the sun-synchronous nature of observations from satellites such as the Landsat series, and are discussed in the following section.

Validation
To assess the accuracy of NIDEM, we validated 'filtered' elevations at multiple sites and intertidal environments along the Australian coastline. We used the nationally-consistent coastal 'Smartline' geomorphic and stability map (Sharples et al., 2009) to identify the dominant intertidal landform type for each validation site. These landforms were summarised into three distinct intertidal categories: sandy beaches (including sandy and mixed sand beach and shores), tidal flats (including sand, mud and undifferentiated tidal flats) and rocky reefs and shores (including hard bedrock shores). We obtained validation data from three different elevation and bathymetry data sources: the DEM of Australia derived from LiDAR 5 Metre Grid (Geoscience Australia, 2015), point elevation data collected from Real Time Kinematic (RTK) GPS surveys (Danaher and Collett, 2006;HydroSurvey Australia, 2009), and 1.0 m resolution multibeam bathymetry surveys (Solihuddin et al., 2016). Validation datasets were processed as described below and pooled across sample sites for each intertidal environment type. We compared modelled (i.e. NIDEM) and observed (i.e. validation) data using the pandas Python package (McKinney, 2011) by calculating RMSE accuracy and two correlation coefficients (Pearson's and Spearman's). Pearson's correlation was used to assess whether NIDEM accurately modelled absolute elevation by evaluating the linear relationship between modelled and observed values. Spearman's correlation assessed whether modelled elevations were monotonically related to validation elevations, allowing us to evaluate if at a minimum the relative -not absolute -topography of the intertidal zone was captured by NIDEM. For example, this could assess to what extent actual low tide terrain was correctly identified as low tide terrain in NIDEM even if absolute elevations from the tidal model were incorrect.

LiDAR validation data
The bare-earth LiDAR 5 Metre Grid covers over 245,000 km 2 of predominantly coastal terrain across eastern and northern Australia with accuracy of better than 0.30 m (95% confidence, AHD vertical datum; Geoscience Australia, 2015). However, the dataset has sporadic coverage of the intertidal zone, with some coastal regions affected by elevation discontinuities caused by aerial surveys flown at different tidal stages (e.g. Fig. 3a). To extract usable validation data for the intertidal zone in areas affected by these artefacts, we developed a novel LiDAR point-cloud tidal tagging approach to identify non-inundated intertidal terrain in the DEM. Fourteen validation sites (Fig. 1) were identified within randomly selected tidal modelling polygons based on the availability of extensive intertidal terrain in both NIDEM and LiDAR datasets. LiDAR 5 Metre Grid data was extracted for each of these sites and reprojected using average resampling to match the NIDEM resolution and projection system with the gdalwarp image reprojection and warping utility (GDAL/OGR contributors, 2018). For eight of the fourteen validation sites affected by tidal artefacts, we used the las2txt tool from LAStools software (Isenburg, 2018) to extract a 1% sample of xyz points and associated metadata from the . las format LiDAR pointcloud datasets used to generate the LiDAR 5 Metre Grid (Fig. 3b). For each extracted LiDAR point, we used the OTPS tidal model Erofeeva, 2002, 2010) to compute a tidal height relative to MSL based on the point's time of acquisition timestamp. By rasterising all points with z-values (elevations) greater than the instantaneous tidal height at the time of LiDAR acquisition (plus a 0.15 m buffer to account for point Fig. 3. The process followed to extract usable validation data for the intertidal zone from the (a) bare-earth LiDAR 5 Metre Grid dataset (Geoscience Australia, 2015) for an example area affected by tidal-stage artefacts (Gladstone in central Queensland; dark colours indicate greater depth). For each validation site, (b) xyz LiDAR point clouds were extracted and tagged with tidal heights using the OTPS tidal model for the exact moment each point was acquired. Tidally tagged points located above the water level at the time of LiDAR acquisition were rasterised and used to produce a mask indicating the location of (c) non-inundated intertidal terrain in the LiDAR DEM.
cloud noise), we produced a binary mask representing non-inundated terrain in the LiDAR DEM. This mask was cleaned to remove remaining noise by using the scipy.ndimage mathematical morphology binar-y_opening and binary_closing tools (Jones et al., 2014). Binary opening removed isolated pixels by 'eroding' (shrinking the data area by one pixel) and subsequently 'dilating' (expanding the remaining data area by one pixel) the mask layer, while binary closing filled gaps by first dilating then eroding the array (Serra, 1983). The resulting cleaned mask was then used to extract intertidal validation data from the original LiDAR 5 Metre Grid (Fig. 3c).

RTK GPS validation data
RTK GPS transect elevation survey data covering tidal flats at East Point in Mindal Bay, Darwin and Moreton Bay, Queensland (Fig. 1) were collected by HydroSurvey Australia (2009) and Danaher and Collett (2006). Point data with elevations in AHD (stated vertical accuracy of ± 0.02 m) from both surveys were rasterised to match the NIDEM resolution and projection system by taking the average elevation value when multiple GPS points fell within a single Landsat pixel.

Multibeam bathymetry validation data
Multibeam bathymetry validation data were obtained for five sites in the Buccaneer Archipelago in the southern Kimberley region (Cockatoo Island, east Tallon Island, west Tallon Island, Waterflow and Irvine/Bathurst Islands; Fig. 1). The data consisted of 1 m resolution gridded AHD reef elevations measured using an Odom ES3 multibeam echo sounder with a Trimble RTX satellite subscription ( ± 0.02 m positioning accuracy; Solihuddin et al., 2016). All multibeam data were resampled to the NIDEM resolution and projection system using average resampling.

Results and discussion
In this study, we present an automated approach to generating satellite-derived elevation data for over 15,387 km 2 of exposed intertidal Fig. 4. The output NIDEM layers at continental-scale (top row; aggregated to 2500 m pixel size for visualisation) and local-scale (bottom row; focused on the Capricornia-Bunker Group of reefs in central Queensland). The NIDEM 'mask' layer (a, d) highlights pixels flagged as exhibiting poor tidal model performance or significant geomorphic change across the 30-year time series based on the ITEM confidence layer (blue), pixels located above 25 m elevation (red) or pixels located below −25 m depth (orange). The effect of applying this mask can be seen by comparing the NIDEM 'unfiltered' layer (b, e) against the NIDEM 'filtered' layer (c, f) which was cleaned by masking by NIDEM 'mask'. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Model performance across sandy beaches and shores
Elevation values from NIDEM agreed well with LiDAR 5 Metre Grid validation datasets across sandy beach and shore intertidal environments. Validation elevations for 121,725 paired pixels pooled across five sites along the Australian coastline were strongly correlated with NIDEM (Pearson's ρ = 0.92, Spearman's ρ = 0.93), resulting in a low RMSE of ± 0.41 m which approached the vertical accuracy of the contributing LiDAR surveys (i.e. < 0.30 m; Fig. 6a, Table 1). Our results compare favourably to more manual, non-automated approaches to waterline delineation, including the construction of intertidal DEMs for a 1,267 km 2 portion of the Dongsha sandbank in China's Jiangsu Province using on-screen digitisation of MODIS and Landsat imagery (Liu et al., 2013a(Liu et al., , 2013b. The Dongsha sandbank DEMs exhibited similarly strong correlations with LiDAR validation data (R 2 = 0.92), with vertical accuracies of between ± 0.45-0.62 m RMSE depending on the number of input satellite images used to generate the DEMs. Tseng et al.

Table 1
Sites used for validating NIDEM 'filtered' intertidal elevation data, including LiDAR data from the DEM of Australia derived from LiDAR 5 Metre Grid (Geoscience Australia, 2015) for 14 sites, elevation data collected from Real Time Kinematic (RTK) GPS surveys in Darwin, Northern Territory (HydroSurvey Australia, 2009) and Moreton Bay, Queensland (Danaher and Collett, 2006), and 1.0 m resolution multibeam bathymetry surveys for five sites in the southern Kimberley region, Western Australia (Solihuddin et al., 2016 Bishop-Taylor, et al. Estuarine, Coastal and Shelf Science 223 (2019) 115-128 22 years of Landsat imagery, although this DEM covered a relatively small extent of 16 km 2 . Poor modelling results occurred at a single site near Robbins Island in north-western Tasmania where accuracies were significantly lower (RMSE of ± 0.57 m) than the four other sites assessed (≤0.34 m; Table 1). This result was driven by a set of outlying values where NIDEM predicted low elevations between 0.0 and 0.5 m for an area of intertidal terrain with actual elevations ranging from 1.0 to 2.0 m (Fig. 6a). The Robbins Island site occurs on the boundary of two large shallow tidal basins (Boullanger Bay and Big Bay) and the intersection of three estuaries (the Duck, Montagu and Welcome Rivers), creating highly variable and unpredictable tidal conditions (Spruzen et al., 2008;Donaldson et al., 2012). This complex tidal regime challenges the application of any continental-scale tidal model, especially given that the individual instantaneous modelled tide heights used to assign elevations to waterlines are likely to rapidly decrease in accuracy with increasing distance from a given modelling point (Donaldson et al., 2012;Bell et al., 2016). In this study, we attempted to moderate the influence of spatially variable tidal conditions by adopting a multi-resolution tidal modelling framework that allowed local complexity in tides to drive the scale and boundaries of the model . This minimised the distance between each satellite image and the locations used to calculate tides, resulting in a closer match between the tidal stage observed from space and the resulting attributed tide height. At the Robbins Island site, however, the resulting automatically derived tidal modelling polygons dissected Boullanger Bay in the west of the study site, causing an area of intertidal terrain within the bay to be assigned tidal heights that were more appropriate for the western Tasmanian open coast (see Appendix A1). Although this result highlights areas where the multi-resolution tidal modelling framework used in this study could be improved, results for other sites indicate that the approach provides a significant improvement over previous artefactprone continental-scale modelling frameworks which used a regular 1 × 1°grid to assign tide heights to imagery .

Model performance across tidal flats
We compared the performance of NIDEM across 99,506 pixels of tidal flat terrain at nine LiDAR 5 Metre Grid and RTK GPS survey sites across northern and southern Australia (Fig. 6b, Table 1). Although overall accuracies were high and consistent with sandy beach and shore sites (RMSE of ± 0.39 m), tidal flats exhibited slightly lower correlations overall (Pearson's ρ = 0.78, Spearman's ρ = 0.81). This result was expected given two interacting characteristics of tidal flat environments: their extensive, low-slope morphology and dynamism. Tidal flats are known to be highly variable, and can undergo significant geomorphic change in response to tidal processes and changes in sediment flux (Mason et al., 2010;Chen et al., 2016). Even small vertical changes in elevation in these gently sloping environments can produce large horizontal shifts in waterline locations (e.g. a 0.1 m change in elevation on a 1:500 slope tidal flat would result in a 50 m horizontal discrepancy in waterline position, equivalent to two Landsat pixels; Mason et al., 1995;Zhao et al., 2008). This presented a particular challenge to our composite-based approach, as coastal change at any point across our 30-year time-series would invalidate or reduce the accuracy of our waterline-derived modelled elevations. Previous studies have generated intertidal DEMs based on shorter temporal extents to reduce the influence of coastal change (e.g. Zhao et al., 2008;Liu et al., 2013b), however these approaches risk reducing the overall accuracy of modelled elevations by limiting available satellite observations and increasing the relative influence of natural variability and noise. To minimise the impact of geomorphic change while maximising our use of the Landsat archive, our approach used NDWI variance to mask individual pixels where tidal modelling poorly explained patterns of inundation (i.e. ITEM 2.0 confidence; Sagar et al., 2017). Although this approach effectively removed areas of significant coastal change, the choice of a single universal masking value (i.e. 0.25 NDWI standard deviation) unavoidably preserved areas of intertidal terrain subject to subtler, long-term geomorphic change. In these areas, NIDEM represents median or 'typical' elevation conditions across the 30-year time series.
Validation results for tidal flat environments revealed a distinct lack of modelled elevation data for the upper portion of the tidal range (e.g. between 0.5 and 1.3 m AHD; Fig. 6b). In tidal flats, the elevation zone above mean sea level is typically occupied by coastal wetlands including mangroves (Bunt et al., 1985; see Appendix A2). Although these vegetated intertidal communities are regularly inundated by tidal flows, patterns of inundation are difficult to observe from satellites due to the presence of dense canopy cover. This was particularly true for the NDWI index used to detect water in our model, which intentionally differentiates between open water and vegetated features based on their low and high reflectance of NIR radiation respectively (McFeeters, 1996). Although areas of exposed tidal flat terrain located behind mangroves were accurately modelled by NIDEM across our nine validation sites (i.e. areas above 1.3 m AHD in Fig. 6b), vegetative resistance to flowing water in mangrove, saltmarsh and other coastal wetland communities can cause significant hydrodynamic attenuation in tidal flow (Rodríguez et al., 2017;Montgomery et al., 2018). This would affect the spatial distribution of extracted waterlines by creating lags between tidal conditions and patterns of observed water. Due to these issues, caution should be applied when interpreting elevation values for intertidal regions on the landward side of coastal wetlands (e.g. Appendix A2), or other environments such as estuaries or areas influenced by artificial hydraulic structures where tidal flows may be similarly modified or restricted (Williams and Watford, 1997;Rodríguez et al., 2017).

Model performance across rocky shores and reefs
Although NIDEM produced accurate modelled elevations in sandy beaches and tidal flat environments, results were less accurate within rocky shore and reef environments. Modelled elevations (n = 2299) had a low Pearson's correlation of 0.46 and a high RMSE of ± 2.98 m vertical accuracy compared to LiDAR 5 Metre Grid and multibeam bathymetry data for seven validation sites (Fig. 6c, Table 1). Spearman's correlation was higher at 0.79, indicating that NIDEM largely captured the relative -but not absolute -topography of the surveyed sites. This poor result was driven by validations across five Buccaneer Archipelago reef sites within the southern Kimberley region, including one site (Cockatoo Island) where NIDEM elevations were negatively correlated with multibeam bathymetry (Pearson's and Spearman's ρ of −0.26 and −0.19; Table 1). The fringing reefs of the Buccaneer Archipelago are exposed to an extreme and dynamic tidal regime, including the largest tidal range of any coral reef system (approaching 11 m), with strong tidal currents that can exceed 18 km/h (Purcell, 2002;Solihuddin et al., 2016). Of particular relevance to the NIDEM approach, Lowe et al. (2015) observed significant asymmetry in tidal patterns on Tallon Island (one of the sites included in the NIDEM validation), with the shallow, elevated intertidal reef platform rapidly inundating over a short ∼2 h period during flood tides, before draining slowly over ∼10 h during ebb tides. The presence of 'trapped' water on the reef platform for extended periods caused areas of shallow exposed reef to appear as permanently inundated terrain in the ensemble of input Landsat observations, and violated the assumptions of the input ITEM model by causing modelled tidal heights to become unsynchronised from actual water levels. Our results indicate that satellite-derived intertidal DEMs such as NIDEM are unlikely to produce accurate elevations in regions that exhibit significant unaccounted-for tidal asymmetry, or where the resolution or local accuracy of the tidal models used to assign tide heights is low. Increasing the accuracy of estimated tide heights in these locations will likely require accounting for finerscale bathymetry and substrate conditions using locally specified hydraulic modelling (e.g. Lowe et al., 2015) that would be challenging to apply at the continental-scale. Future work should focus on developing additional methods for quantifying local tidal asymmetry and poor tidal modelling performance so that the accuracy and validity of satellitederived intertidal DEMs can be quantitatively assessed for each pixel of intertidal terrain.

Potential of continental-scale intertidal elevation modelling
Our results demonstrate that satellite-derived elevation models based on waterline extraction can approach the accuracy of LiDAR for modelling the topography of tidal flats and sandy beach and shore environments. The NIDEM approach is particularly suitable for remote areas of inaccessible coastline (e.g. northern Australia) where higherresolution approaches, such as shallow-water multibeam bathymetry or airborne bathymetric LiDAR, would be cost prohibitive or impractical. By providing an important 'missing link' between existing medium resolution terrestrial topographic and marine bathymetric datasets, we anticipate elevation data from NIDEM will contribute to unified coastal terrain models that combine best-available topography and bathymetry datasets into single continuous 'topobathy' datasets (Hogrefe et al., 2008). Examples are the nthaus30 depth model for Northern Australia (Beaman, 2018b) and the gbr30 depth model for the Great Barrier Reef (Beaman, 2018a) used to test the inclusion of NIDEM as source elevation data (Fig. 7). NIDEM data fills in the shallow reef flats around islands and sand bars within river channels, while subtly reducing the stepping effect across the land-sea interface. NIDEM could thus provide coastal managers with near-seamless elevation data extending from inland to the ocean floor, and support a holistic approach to understanding the physical and ecological processes influencing the coastal zone. This may include providing valuable baseline elevation data for predicting the impact of coastal hazards such as storm surges or tsunami inundation (e.g. Skinner et al., 2015;Smolders et al., 2015), investigating mechanisms of coastal erosion and sediment transport (Ryu et al., 2001(Ryu et al., , 2008Hsu et al., 2013;Gharibreza et al., 2014), or improving modelling of sea-level rise under future climate change scenarios (Galbraith et al., 2002;Thorner et al., 2014). Access to accurate intertidal elevation data is also critical for studying and conserving coastal ecosystems, particularly given the important role intertidal topography and tidal dynamics play in regulating environmental stress and spatial patterns of species richness, abundance and productivity Valdivia et al., 2011). The outputs from NIDEM facilitate whole-of-landscape approaches to ecological modelling and three-dimensional habitat mapping in the intertidal zone that are not artificially restricted to either the marine or terrestrial domain.
A key advantage of our continental-scale approach to intertidal elevation modelling was using composite layers to identify the typical location of the waterline at various tidal heights. Intertidal zones exhibit considerable natural variability in the location of the waterline between identical tides or during ebb and flow stages (Ryu et al., 2001;Boak and Turner, 2005;Li and Gong, 2016). Wave run-up or the effect of wind can affect the location of waterlines by tens of meters on low sloping beaches or flats (Thieler and Danforth, 1994), while long-term sea-level variation, seasonal influences or storm surge events can drive Fig. 7. Hillshaded nthaus30 depth model for Northern Australia (Beaman, 2018b) showing the Joseph Bonaparte Gulf and Victoria River area (a) before and (b) after the inclusion of NIDEM; hillshaded gbr30 depth model for the Great Barrier Reef (Beaman, 2018a) showing the Fitzroy River area of central Queensland (c) before and (d) after NIDEM. Note the NIDEM data subtly improves the stepping effect across the land/ocean interface and helps fill in the shallow reef flats around the islands and sand bars, resulting in better-defined river channels. additional unpredictable variability (Boak and Turner, 2005;García-Rubio et al., 2015). Recent long-term shoreline trend studies have shown that combining multiple observations into composite layers can effectively isolate these factors, and allow waterline contours to be extracted that can be up to twice as spatially accurate as the resolution of the input satellite imagery (Almonacid-Caballer et al., 2016;Hagenaars et al., 2018). Generating median composites at ten percent increments of the tidal range allowed us to both reduce sources of natural variability and extract waterline contours that were most representative of the shoreline position across the full range of observed tidal conditions (Almonacid-Caballer et al., 2016). Using a robust central tendency median compositing method based on good quality Landsat pixel observations combined with an elevation and bathymetry filtering step also greatly reduced sensor artefacts and noise, including false positives in data-poor deeper ocean areas affected by sunglint or clouds. Issues with noise persist in some areas close to the coast, including the Archipelago of the Recherche in southern Western Australia, Port Phillip Bay in Victoria, and the south-eastern coast of Tasmania and King Island (Fig. 1). However, across most of the Australian coastline these approaches effectively eliminated the need for manually-derived masks (e.g. Chen and Chang, 2009;Murray et al., 2012) or the subjective manual selection of clear scenes (Boak and Turner, 2005). This allowed us to leverage the full Landsat archive, and facilitated the automated extraction of waterline contours at a scale that would be impractical based on manual waterline digitisation (e.g. Zhao et al., 2008). Importantly, this automated approach improves reproducibility, allowing NIDEM to be improved over time as new Landsat observations are added to the DEA platform. To facilitate reanalysis, NIDEM is released in both a 'filtered' and 'unfiltered' version, allowing users to modify or apply custom filtering to the raw data depending on their application.

Limitations and future work
The complex behaviour of tides mean that a sun synchronous sensor like Landsat does not observe the full range of the tidal cycle at all locations (Eleveld et al., 2014;Parke et al., 1987). To date, however, few studies modelling the extent and topography of the intertidal zone have addressed issues of tidal bias (e.g. Murray et al., 2012;Chen et al., 2016;Tseng et al., 2017), with the remainder implicitly assuming that satellite observations of the coastline were representative of the full local tidal range. While biases in the proportion of the tidal range observed do not affect the accuracy of absolute elevation models like NIDEM, they can prevent models from providing elevation data for areas of the intertidal zone exposed or inundated at the extremes of the tidal range. This risks giving misleading insights into the true extent of were biased away from low or high tides respectively. For example, a polygon with a spread of 0.7, a high tide offset of 0.05 and a low tide offset of 0.25 indicates that Landsat observed 70% of the tidal range, but did not image the highest 5% or lowest 25% or of tide heights. the intertidal zone, and reduces the comparability of upper and lower elevations at different locations. The portion of the tidal cycle observed by a sun synchronous sensor can be estimated using a tidal model shown in Fig. 8. Across the Australian continent, both the overall spread of the observed tidal range compared to the full tidal range (Fig. 8a) and offsets in this proportion relative to the lowest and highest tidal extremes ( Fig. 8b and c) varied greatly, even across relatively small distances in areas of rapid tidal flux. This evaluation of the observed tidal range at a particular location provides valuable information to users about the 'fitness for purpose' of NIDEM at a given location for their specific application. We strongly recommend that future regional-, continental-or global-scale intertidal zone analyses using EO data consider potential biases associated with the representativeness of their input data relative to the full tidal range. This would serve as an important first step to obtaining a better understanding of how processes driving intertidal variability affect remote sensing analyses in the intertidal zone. Liu et al. (2013b) performed a quantitative assessment of the factors influencing the accuracy of intertidal DEMs, finding that overall DEM accuracy was strongly correlated with the number of available input satellite observations. This increase in accuracy was driven largely by the spatial coverage of waterline contours relative to total intertidal extent, with higher densities of waterline contours leading to reduced RMSE. Interpolation error represents a key source of uncertainty in waterline modelling approaches, with low contour densities increasing the proportion of the study area requiring interpolation (Mason et al., 2001;Liu et al., 2013b). The NIDEM approach based on tenth percentile tidal composites can produce relatively large spacing between waterline contours compared to waterlines digitised from individual observations, particularly in extensive shallow intertidal areas with low topographic relief (Mason et al., 2001). Rather than combining tidally tagged imagery at constant ten percent intervals of the tidal range, future work could develop an adaptive approach to composite generation tailored to the local availability of high quality satellite observations. This may result in a smaller number of higher quality tidal composites in areas with few clear observations, and a denser collection of tidal composites in areas with abundant clear observations throughout the 30-year time-series. By increasing the density of extracted waterline contours, this adaptive approach would improve the ability of the DEM to resolve finer-scale topographic features, and potentially allow high quality elevation estimates to be produced using shorter temporal extents. This may facilitate the comparison of DEMs from different temporal epochs to track and volumetrically estimate rates of coastal change (e.g. Ryu et al., 2008;Mason et al., 2010). This could be a particularly powerful approach for modelling the coastal zone if combined with additional sources of earth observation data available from satellites, such as the European Space Agency's Sentinel-2A and 2B, with a high revisit frequency compared to Landsat (approximately 5 days compared to 16 days; Drusch et al., 2012).

Conclusion
In this study, we have presented an automated approach to modelling the elevation of the intertidal zone based on a 30-year time series archive of Landsat remote sensing data and a multi-resolution tidal modelling framework. The resulting NIDEM dataset is to our knowledge the first continental-scale elevation model of the exposed intertidal zone, and provides an unprecedented three-dimensional representation of Australia's tidal flats, sandy beaches and shores, and rocky shores and reef environments at 25 m spatial resolution. NIDEM is based on freely available EO data and open-source software, making it directly applicable to any tidally influenced coastal environment globally. Future work will focus on integrating our approach with higher spatial and temporal resolution sources of EO data (e.g. Sentinel 2), and developing adaptive approaches to waterline contour extraction that will maximise the tidal resolution of the model based on the local availability of high quality satellite observations. This could enable modelling and comparison of finer-scale intertidal topographic features across time and at the continental-scale.