Statistical modeling of the long-range-dependent structure of barrier island framework geology and surface geomorphology

Shorelines exhibit long-range dependence (LRD) and have been shown in some environments to be described in the wave number domain by a power-law characteristic of scale independence. Recent evidence suggests that the geomorphology of barrier islands can, however, exhibit scale dependence as a result of systematic variations in the underlying framework geology. The LRD of framework geology, which influences island geomorphology and its response to storms and sea level rise, has not been previously examined. Electromagnetic induction (EMI) surveys conducted along Padre Island National Seashore (PAIS), Texas, United States, reveal that the EMI apparent conductivity (σ_a) signal and, by inference, the framework geology exhibits LRD at scales of up to 10¹ to 10² km. Our study demonstrates the utility of describing EMI σ_a and lidar spatial series by a fractional autoregressive integrated moving average (ARIMA) process that specifically models LRD. This method offers a robust and compact way of quantifying the geological variations along a barrier island shoreline using three statistical parameters (p, d, q). We discuss how ARIMA models that use a single parameter d provide a quantitative measure for determining free and forced barrier island evolutionary behavior across different scales. Statistical analyses at regional, intermediate, and local scales suggest that the geologic framework within an area of paleo-channels exhibits a first-order control on dune height. The exchange of sediment amongst nearshore, beach, and dune in areas outside this region are scale independent, implying that barrier islands like PAIS exhibit a combination of free and forced behaviors that affect the response of the island to sea level rise.

1 Introduction

Barrier island transgression in response to storms and sea level rise depends to varying degrees on preexisting geologic features. The traditional assumption of uniform sand at depth and alongshore cannot explain many observations. Models of barrier island evolution are required to ascertain the degree to which the island is either free (such as a large sand body) or forced (i.e., constrained) by the underlying geology. Despite growing evidence that the underlying geological structure, otherwise termed framework geology, of barrier islands influences nearshore, beach, and dune morphology (e.g., Belknap and Kraft, 1985; Houser, 2012; Lentz and Hapke, 2011; McNinch, 2004; Riggs et al., 1995), this variable remains largely absent from shoreline change models that treat the geology as being uniform alongshore (e.g., Dai et al., 2015; Plant and Stockdon, 2012; Wilson et al., 2015). Spatial variation in the height and position of the dune line impacts the overall transgression of the island with sea level rise (Sallenger, 2000). Transgression is accomplished largely through the transport and deposition of beach and dune sediments to the backbarrier as washover deposits during storms (Houser, 2012; Morton and Sallenger, 2003; Stone et al., 2004).

1.1 Framework geology controls on barrier island evolution

The dynamic geomorphology of a barrier island system is the result of a lengthy, complex, and ongoing history that is characterized by sea level changes and episodes of deposition and erosion (e.g., Anderson et al., 2015; Belknap and Kraft, 1985; Rodriguez et al., 2001). Previous studies demonstrate that the framework geology of barrier islands plays a considerable role in the evolution of these coastal landscapes (Belknap and Kraft, 1985; Evans et al., 1985; Kraft et al., 1982; Riggs et al., 1995). For example, antecedent structures such as paleo-channels, ravinement surfaces, offshore ridge and swale bathymetry, and relict transgressive features (e.g., overwash deposits) have been suggested to influence barrier island geomorphology over a wide range of spatial scales (Hapke et al., 2010, 2016; Houser, 2012; Lentz and Hapke, 2011; McNinch, 2004). In this study, the term framework geology is specifically defined as the topographic surface of incised valleys, paleo-channels, and/or the depth to ravinement surface beneath the modern beach.

As noted by Hapke et al. (2013), the framework geology at the regional scale (> 30 km) influences the geomorphology of an entire island. Of particular importance are the location and size of glacial, fluvial, tidal, and/or inlet paleo-valleys and paleo-channels (Belknap and Kraft, 1985; Colman et al., 1990; Demarest and Leatherman, 1985), and paleo-deltaic systems offshore or beneath the modern barrier system (Coleman and Gagliano, 1964; Frazier, 1967; Miselis et al., 2014; Otvos and Giardino, 2004; Twichell et al., 2013). At the regional scale, nonlinear hydrodynamic interactions between incident wave energy and nearshore ridge and swale bathymetric features can generate periodic alongshore variations in beach–dune morphology (e.g., Houser, 2012; McNinch, 2004) that are superimposed on larger-scale topographic variations as a result of transport gradients (Tebbens, et al., 2002). At the intermediate scale (10–30 km), feedbacks between geologic features and relict sediments of the former littoral system (e.g., Honeycutt and Krantz, 2003; Riggs et al., 1995; Rodriguez et al., 2001; Schwab et al., 2000) act as an important control on dune formation (Houser et al., 2008) and offshore bathymetric features (e.g., Browder and McNinch, 2006; Schwab et al., 2013). Framework geology at the local scale (≤ 10 km) induces mesoscale (∼ 10¹–10² m) to microscale (< 1 m) sedimentological changes (e.g., Murray and Thieler, 2004; Schupp et al., 2006), variations in the thickness of shoreface sediments (Brown and Macon, 1977; Miselis and McNinch, 2006), and spatial variations in sediment transport across the island (Houser and Mathew, 2011; Houser, 2012; Lentz and Hapke, 2011).

To date, most of what is known regarding barrier island framework geology is based on studies performed at either intermediate or local scales (e.g., Hapke et al., 2010; Lentz and Hapke, 2011; McNinch, 2004), whereas few studies exist at the regional scale for United States coastlines (Hapke et al., 2013). The current study focuses on barrier islands in the United States and we do not consider work on barrier islands in other regions. Assessments of framework geology at regional and intermediate spatial scales for natural and anthropogenically modified barrier islands are essential for improved coastal management strategies and risk evaluation since these require a good understanding of the connections between subsurface geology and surface morphology. For example, studies by Lentz and Hapke (2011) and Lentz et al. (2013) at Fire Island, New York, suggest that the short-term effectiveness of engineered structures is likely influenced by the framework geology. Extending their work, Hapke et al. (2016) identified distinct patterns of shoreline change that represent different responses alongshore to oceanographic and geologic forcing. These authors applied empirical orthogonal function (EOF) analysis to a time series of shoreline positions to better understand the complex multiscale relationships between framework geology and contemporary morphodynamics. Gutierrez et al. (2015) used a Bayesian network to predict barrier island geomorphic characteristics and argue that statistical models are useful for refining predictions of locations where particular hazards may exist. These examples demonstrate the benefit of using statistical models as quantitative tools for interpreting coastal processes at multiple spatial and temporal scales (Hapke et al., 2016).

1.2 Statistical measures of coastline geomorphology

It has long been known that many aspects of landscapes exhibit similar statistical properties regardless of the length or timescale over which observations are sampled (Burrough, 1981). An often-cited example is the length L of a rugged coastline (Mandelbrot, 1967), which increases without bound as the length G of the ruler used to measure it decreases, in rough accord with the formula L(G)∽G^1−D, where D≥1 is termed the fractal dimension of the coastline. Andrle (1996), however, has identified limitations of the self-similar coastline concept, suggesting that a coastline may contain irregularities that are concentrated at certain characteristic length scales owing to local processes or structural controls. Recent evidence from South Padre Island, Texas (Houser and Mathew, 2011), Fire Island, New York (Hapke et al., 2010), and Santa Rosa Island, Florida (Houser et al., 2008), suggests that the geomorphology of barrier islands is affected to varying degrees by the underlying framework geology and that this geology varies, often with periodicities, over multiple length scales. The self-similarity of the framework geology and its impact on the geomorphology of these barrier islands was not examined explicitly.

Many lines of evidence suggest that geological formations in general are inherently rough (i.e., heterogeneous) and contain multiscale structure (Bailey and Smith, 2005; Everett and Weiss, 2002; Radliński et al., 1999; Schlager, 2004). Some of the underlying geological factors that lead to self-similar terrain variations are reviewed by Xu et al. (1993). In essence, competing and complex morphodynamic processes, influenced by the underlying geological structure, operate over different spatiotemporal scales, such that the actual terrain is the result of a complex superposition of the various effects of these processes (see Lazarus et al., 2011). Although no landscape is strictly self-similar on all scales, Xu et al. (1993) show that the fractal dimension, as a global morphometric measure, captures multiscale aspects of surface roughness that are not evident in conventional local morphometric measures such as slope gradient and profile curvature.

With respect to coastal landscapes, it has been suggested that barrier shorelines are scale independent, such that the wave number spectrum of shoreline variation can be approximated by a power law at alongshore scales from tens of meters to several kilometers (Lazarus et al., 2011; Tebbens et al., 2002). However, recent findings by Houser et al. (2015) suggest that the beach–dune morphology of barrier islands in Florida and Texas is scale dependent and that morphodynamic processes operating at swash (0–50 m) and surf-zone (< 1000 m) scales are different than the processes operating at larger scales. In this context, scale dependence implies that a certain number of different processes are simultaneously operative, each process acting at its own scale of influence, and it is the superposition of the effects of these multiple processes that shapes the overall behavior and shoreline morphology. This means that shorelines may have different patterns of irregularity alongshore with respect to barrier island geomorphology, which has important implications for analyzing long-term shoreline retreat and island transgression. Lazarus et al. (2011) point out that deviations from power-law scaling at larger spatial scales (tens of kilometers) emphasizes the need for more studies that investigate large-scale shoreline change. While coastal terrains might not satisfy the strict definition of self-similarity, it is reasonable to expect them to exhibit long-range dependence (LRD). LRD pertains to signals in which the correlation among observations decays like a power law with separation, i.e., much slower than one would expect from independent observations or those that can be explained by a short-memory process, such as an autoregressive moving average (ARMA) with small (p,q) (Beran, 1994; Doukhan et al., 2003).

1.3 Research objectives

This study performed at Padre Island National Seashore (PAIS), Texas, United States, utilizes electromagnetic induction (EMI) apparent conductivity σ_a responses to provide insight into the relation between spatial variations in framework geology and surface morphology. Two alongshore EMI surveys at different spatial scales (100 and 10 km) were conducted to test the hypothesis that, like barrier island morphology, subsurface framework geology exhibits the LRD characteristic of scale independence. The σ_a responses, which are sensitive to parameters such as porosity and mineral content, are regarded herein as a rough proxy for subsurface framework geology (Weymer et al., 2015a). This assumes, of course, that alongshore variations in salinity and water saturation, and other factors that shape the σ_a response, can be neglected to first order. A corroborating 800 m ground-penetrating radar (GPR) survey, providing an important check on the variability observed within the EMI signal, confirms the location of a previously identified paleo-channel (Fisk, 1959) at ∼ 5–10 m of depth. The overall geophysical survey design allows for a detailed evaluation of the long-range-dependent structure of the framework geology over a range of length scales spanning several orders of magnitude. We explore the applicability of autoregressive integrated moving-average (ARIMA) processes as models that describe the statistical connections between EMI and light detection and ranging (lidar) spatial data series. This paper utilizes a generalized fractional ARIMA ($\mathrm{0},d,\mathrm{0}$) process (Hosking, 1981) that is specifically designed to model LRD for a given data series using a single differencing non-integer parameter d. The parameter d can be used in the present context to discriminate between forced, scale-dependent controls by the framework geology, i.e., stronger LRD (d → 0.5), and free behavior that is scale independent, i.e., weaker LRD (0←d). In other words, it is the particular statistical characteristics of the framework geology LRD at PAIS that we are trying to ascertain from the EMI σ_a signal, with the suggestion that σ_a measurements can be used similarly at other sites to reveal the hidden LRD characteristics of the framework geology.

2 Background and regional setting

2.1 Utility of electromagnetic methods in coastal environments

Methods to ascertain the alongshore variability in framework geology, and to test long-range dependence, are difficult to implement and can be costly. Cores provide detailed point-wise geologic data; however, they do not provide laterally continuous subsurface information (Jol et al., 1996). Alternatively, geophysical techniques including seismic and GPR provide spatially continuous stratigraphic information (e.g., Buynevich et al., 2004; Neal, 2004; Nummedal and Swift, 1987; Tamura, 2012), but they are not ideally suited for LRD testing because the data combine depth and lateral information at a single acquisition point. Moreover, GPR signals attenuate rapidly in saltwater environments whereas seismic methods are labor-intensive and cumbersome. Conversely, terrain conductivity profiling is an easy-to-use alternative that has been used in coastal environments to investigate fundamental questions involving instrument performance characteristics (Delefortrie et al., 2014; Weymer et al., 2016), groundwater dynamics (Stewart, 1982; Fitterman and Stewart, 1986; Nobes, 1996; Swarzenski and Izbicki, 2009), and framework geology (Seijmonsbergen et al., 2004; Weymer et al., 2015a). Previous studies combining EMI with either GPR (Evans and Lizarralde, 2011) or coring (Seijmonsbergen et al., 2004) demonstrate the validity of EMI measurements as a means to quantify alongshore variations in the framework geology of coastlines.

In the alongshore direction, Seijmonsbergen et al. (2004) used a Geonics EM34^™ terrain conductivity meter crossing a former outlet of the Rhine River, Netherlands, to evaluate alongshore variations in subsurface lithology. The survey was conducted in an area that was previously characterized by drilling and these data were used to calibrate the σ_a measurements. The results from the study suggest that coastal sediments can be classified according to σ_a signature and that high σ_a values occur in areas where the underlying conductive layer is thick and close to the surface. Although Seijmonsbergen et al. (2004) propose that EMI surveys are a rapid, inexpensive method to investigate subsurface lithology, they also acknowledge that variations in salinity as a result of changing hydrologic conditions, storm activity, and/or tidal influence confound the geological interpretation and should be investigated in further detail (see Weymer et al., 2016).

The challenge on many barrier islands and protected national seashores is obtaining permission for extracting drill cores to validate geophysical surveys. At PAIS, numerous areas along the island are protected nesting sites for the endangered Kemp's ridley sea turtle and migratory birds while other areas comprise historic archeological sites with restricted access. Thus, coring is not allowed and only noninvasive techniques, such as EMI–GPR, are permitted.

2.2 Regional setting

North Padre Island is part of a large arcuate barrier island system located along the Texas Gulf of Mexico coastline. The island is one of 10 national seashores in the United States and is protected and managed by the National Park Service, a bureau of the Department of the Interior. PAIS is 129 km in length, and is an ideal setting for performing EMI surveys because there is minimal cultural noise to interfere with the σ_a signal, which as stated earlier we regard as a proxy for alongshore variations in framework geology (Fig. 1). Additionally, there are high-resolution elevation data available from a 2009 aerial lidar survey. The island is not dissected by inlets or navigation channels (excluding Mansfield Channel separating the North and South Padre islands) or modified by engineered structures (e.g., groynes, jetties) that often interfere with natural morphodynamic processes (see Talley et al., 2003). The above characteristics make the study area an exceptional location for investigating the relationships between large-scale framework geology and surface morphology.

As described in Weymer et al. (2015a; Fig. 3), locations of several paleo-channels were established by Fisk (1959) based on 3000 cores and seismic surveys. More than 100 boreholes were drilled to the top of the late Pleistocene surface (tens of meters of depth) providing sedimentological data for interpreting the depth and extent of the various paleo-channels. These cores were extracted ∼ 60 years ago, but the remnant Pleistocene and Holocene fluvial–deltaic features described in Fisk's study likely have not changed over decadal timescales.

Figure 1 Location map and DEM of the study area at Padre Island National Seashore (PAIS), Texas, United States. Elevations for the DEM are reported as meters above sea level (m a.s.l.). Approximate locations of field images (red dots) from the northern (N), central (C), and southern (S) regions of the island showing alongshore differences in beach–dune morphology. Note that views are facing south for the central and southern locations, and the northern location view is to the north. Images taken in October 2014.

Geologic interpretations based on the Fisk (1959) data suggest that the thickness of the modern beach sands is ∼ 2–3 m, and they are underlain by Holocene shoreface sands and muds to a depth of ∼ 10–15 m (Brown and Macon, 1977; Fisk, 1959). The Holocene deposits lie upon a Pleistocene ravinement surface of fluvial–deltaic sands and muds and relict transgressive features. A network of buried valleys and paleo-channels in the central segment of the island, as interpreted by Fisk (1959), exhibits a dendritic, tributary pattern. The depths of the buried valleys inferred from seismic surveys range from ∼ 25 to 40 m (Brown and Macon, 1977). These channels have been suggested to have incised into the Pleistocene paleo-surface and became infilled with sands from relict Pleistocene dunes and fluvial sediments reworked by alongshore currents during the Holocene transgression (Weise and White, 1980). However, the location and cross-sectional area of each valley and paleo-channel alongshore is not well constrained. It is also possible that other channels exist other than those identified by Fisk (1959). As suggested in Weymer et al. (2015a), minima in the alongshore σ_a signal are spatially correlated with the locations of these previously identified geologic features. This observation provides an impetus for using EMI to map the known, and any previously unidentified, geologic features alongshore.

3 Methods

A combination of geophysical, geomorphological, and statistical methods are used in this study to quantify the relationships between framework geology and surface geomorphology at PAIS. A description of the EMI, GPR, geomorphometry and statistical techniques is provided in the following sections.

3.1 Field EMI and GPR surveys

Profiles of EMI σ_a responses are typically irregular and each datum represents a spatial averaging of the bulk subsurface electrical conductivity σ, which in turn is a function of a number of physical properties (e.g., porosity, lithology, water content, salinity). The “sensor footprint”, or subsurface volume over which the spatial averaging is performed, is dependent on the separation between the transmitter–receiver (TX–RX) coils (1.21 m in this study) and the transmitter frequency. The horizontal extent, or radius, of the footprint can be more or less than the step size between subsequent measurements along the profile. The sensor footprint determines the volume of ground that contributes to σ_a at each acquisition point, and as will be discussed later, the radius of the footprint has important implications for analyzing LRD. The footprint radius depends on frequency and ground conductivity, but is likely to be of the same order as, but slightly larger than, the intercoil spacing. Two different station spacings were used to examine the correlation structure of σ_a as a function of spatial scale. An island-scale alongshore survey of ∼ 100 km length was performed using a 10 m station spacing (station spacing ≫ footprint radius) such that each σ_a measurement was recorded over an independently sampled volume of ground. Additionally, a sequence of σ_a readings was collected at 1 m of spacing (station spacing < footprint radius) over a profile length of 10 km within the Fisk (1959) paleo-channel region of the island. This survey design allows for comparison of the long-range-dependent structure of the framework geology over several orders of magnitude (10⁰–10⁵ m).

The 100 km long alongshore EMI survey was performed during a series of three field campaigns, resulting in a total of 21 (each of length ∼ 4.5 km) segments that were collected during 9–12 October 2014, 15–16 November 2014, and 28 March 2015. The EMI σ_a profiles were stitched together by importing GPS coordinates from each measurement into ArcGIS^™ to create a single composite spatial data series. The positional accuracy recorded by a TDS Recon PDA equipped with a Holux^™ WAAS GPS module was found to be accurate within ∼ 1.5 m. To reduce the effect of instrument drift caused by temperature, battery, and other systematic variations through the acquisition interval, a drift correction was applied to each segment and the segments were then stitched together, following which a regional linear trend removal was applied to the composite dataset. An additional 10 km survey was performed along a segment of the same 100 km survey line in one day on 29 March 2015. This second composite data series consists of eight stitched segments.

The same multifrequency GSSI Profiler EMP-400^™ instrument was used for each segment. All transects were located in the back-beach environment ∼ 25 m inland from the mean tide level (MTL). This location was chosen to reduce the effect of changing groundwater conditions in response to nonlinear tidal forcing (see Weymer et al., 2016), which may be significant closer to the shoreline. As will be shown later, there is not a direct correlation between high tide and high σ_a values. Thus, we assume the tidal influence on the EMI signal can be neglected over the spatial scales of interest in the present study. Nevertheless, the duration and approximate tidal states of each survey were documented in order to compare with the EMI signal. Tidal data were accessed from NOAA's Tides and Currents database (NOAA, 2015b). Padre Island is microtidal and the mean tidal range within the study area is 0.38 m (NOAA, 2015a). A tidal signature in EMI signals may become more significant at other barrier islands with larger tidal ranges.

For all surveys, the EMI profiler was used in the same configuration and acquisition settings as described in Weymer et al. (2016). The transect locations were chosen to avoid the large topographic variations (see Santos et al., 2009) fronting the foredune ridge that can reduce the efficiency of data acquisition and influence the EMI signal. Measurements were made at a constant step size to simplify the data analysis; for example, ARIMA models require that data are taken at equal intervals (see Cimino et al., 1999). We choose herein to focus on data collected at 3 kHz, resulting in a depth of investigation (DOI) of ∼ 3.5–6.4 m over the range of conductivities found within the study area (Weymer et al., 2016; Table 1). Because the depth of the modern beach sands is ∼ 2–3 m or greater (see Brown and Macon, 1977; p. 56, Fig. 15), variations in the depth to shoreface sands and muds is assumed to be within the DOI of the profiler, which may not be captured at the higher frequencies also recorded by the sensor (i.e., 10 and 15 kHz) .

An 800 m GPR survey was performed on 12 August 2015 across one of the paleo-channels previously identified by Fisk (1959) located within the 10 km EMI survey for comparison with the σ_a measurements. We used a Sensors and Software pulseEKKO Pro^™ system for this purpose. A survey-grade GPS with a positional accuracy of 10 cm was used to match the locations and measurements between the EMI–GPR surveys. Data were acquired in reflection mode at a nominal frequency of 100 MHz with a standard antenna separation of 1 m and a step size of 0.5 m. The instrument settings resulted in a DOI of up to 15 m. Minimal processing was applied to the data and includes a dewow filter and migration (0.08 m ns⁻¹), followed by automatic gain control (AGC) gain (see Neal, 2004). The theory and operational principles of GPR are discussed in many places (e.g., Everett, 2013; Jol, 2008) and will not be reviewed here.

3.2 Geomorphometry

Topographic information was extracted from aerial lidar data that were collected by the U.S. Army Corps of Engineers (U.S. ACE) in 2009 as part of the West Texas Aerial Survey project to assess post-hurricane conditions of the beaches and barrier islands along the Texas coastline. This dataset is the most recent publicly available lidar survey of PAIS and it provides essentially complete coverage of the island. With the exception of Hurricane Harvey, which made landfall near Rockport, Texas, as a category 4 storm in late August, 2017, Padre Island has not been impacted by a hurricane since July 2008, when Hurricane Dolly struck South Padre Island as a category 1 storm (NOAA, 2015a). The timing of the lidar and EMI surveys in this study precede the impacts of Hurricane Harvey, and it is assumed that the surface morphology across the island at the spatial scales of interest (i.e., 10¹–10² km) did not change appreciably between 2009 and 2015.

A 1 m resolution DEM was created from 2009 lidar point clouds available from NOAA's Digital Coast (NOAA, 2017). The raw point cloud tiles were merged to produce a combined point cloud of the island within the park boundaries of PAIS. The point clouds were processed into a continuous DEM using the ordinary kriging algorithm in SAGA GIS, which is a freely available open-source software (http://www.saga-gis.org, last access: 13 January 2018), and subsequent terrain analysis was conducted using an automated approach involving the relative relief (RR) metric (Wernette et al., 2016). Several morphometrics including beach width, dune height, and island width were extracted from the DEM by averaging the RR values across window sizes of 21 m × 21 m, 23 m × 23 m, and 25 m × 25 m. The choice of window size is based on tacit a priori knowledge and observations of the geomorphology in the study area. A detailed description of the procedure for extracting each metric is provided in Wernette et al. (2016).

Each DEM series is paired with the σ_a profile by matching the GPS coordinates (latitude and longitude) recorded in the field by the EMI sensor. Cross-sectional elevation profiles oriented perpendicular to the shoreline were analyzed every 10 m (y coordinate) along the EMI profile to match the same 10 m sampling interval of the σ_a measurements. The terrain variations along each cross-shore profile are summed to calculate beach and island volume based on the elevation thresholds mentioned above. Dune volume is calculated by summing the pixel elevations starting at the dune toe, traversing the dune crest, and ending at the dune heel. In total, six DEM morphometrics were extracted as spatial data series to be paired with the EMI data, each with an identical sample size (n = 9694), which is sufficiently large for statistical ARIMA modeling.

3.3 Statistical methods

Although the procedures for generating the EMI and lidar datasets used in this study are different, the intended goal is the same: to produce spatial data series that contain similar numbers of observations for comparative analysis using a combination of signal processing and statistical modeling techniques. The resulting signals comprising each data series represent the spatial averaging of a geophysical (EMI) or geomorphological elevation variable that contains information about the important processes that form relationships between subsurface geologic features and island geomorphology that can be teased out by means of comparative analysis (Weymer et al., 2015a). Because we are interested in evaluating these connections at both small and large spatial scales, our first approach is to determine the autocorrelation function and Hurst coefficient (self-similarity parameter) H and hence verify whether the data series are characterized by short- and/or long-range memory (Beran, 1992; Taqqu et al., 1995). LRD occurs when the autocorrelation within a series, at large lags, tends to zero like a power function, and so slowly that the sums diverge (Doukhan et al., 2003). LRD is often observed in natural time series and is closely related to self-similarity, which is a special type of LRD.

The degree of LRD is related to the scaling exponent H of a self-similar process, where increasing H in the range 0.5 < H ≤ 1.0 indicates an increasing tendency towards such an effect (Taqqu, 2003). Large correlations at small lags can easily be detected by models with short memory (e.g., ARMA, Markov processes) (Beran, 1994). Conversely, when correlations at large lags slowly tend to zero like a power function, the data contain long-memory effects and either fractional Gaussian noise (fGn) or ARIMA models may be suitable (Taqqu et al., 1995). The R∕S statistic is the quotient of the range of values in a data series and the standard deviation (Beran, 1992, 1994; Hurst, 1951; Mandelbrot and Taqqu, 1979). When plotted on a log∕log plot, the resulting slope of the best-fit line gives an estimate of H, which is useful as a diagnostic tool for estimating the degree of LRD (see Beran, 1994).

It has been suggested that R∕S tends to give biased estimates of H, too low for H > 0.72 and too high for H < 0.72 (Bassingthwaigthe and Raymond, 1994), which was later confirmed by Malamud and Turcotte (1999). Empirical trend corrections to the estimates of H can be made by graphical interpolation, but are not applied here because of how the regression is performed. The R∕S analysis in this study was performed using signal analysis software AutoSignal^™ to identify whether a given signal is distinguishable from a random, white noise process and, if so, whether the given signal contains LRD. The H value is calculated by an inverse variance-weighted linear least-squares curve fit using the logarithms of the R∕S and the number of observations, which provides greater accuracy than other programs that compute the Hurst coefficient.

Two of the simplest statistical time series models that can account for LRD are fGn and ARIMA. In the former case, fGn and its “parent” fractional Brownian motion (fBm) are used to evaluate stationary and nonstationary fractal signals, respectively (see Eke et al., 2000; Everett and Weiss, 2002). Both fGn and fBm are governed by two parameters: variance σ² and the scaling parameter, H (Eke et al., 2000). A more comprehensive class of time series models that has a similar capability to detect long-range structure is ARIMA. Because fGn and fBm models have only two parameters, it is not possible to model the short-range components. Additional parameters in ARIMA models are designed to handle the short-range component of the signal, as discussed by Taqqu et al. (1995) and others. Because the EMI data series presumably contain both short-range and long-range effects, we chose to use ARIMA as the analyzing technique.

ARIMA models are used across a wide range of disciplines in geoscience and have broad applicability for understanding the statistical structure of a given data series as it is related to some physical phenomenon (see Beran, 1992, 1994; Box and Jenkins, 1970; Cimino et al., 1999; Granger and Joyeux, 1980; Hosking, 1981; Taqqu et al., 1995). For example, Cimino et al. (1999) apply R∕S analysis, ARIMA, and neural network analysis to different geological datasets including tree ring data, Sr isotope data of Phanerozoic seawater samples, and El Niño phenomena. The authors show that their statistical approach enables (1) recognition of qualitative changes within a given dataset, (2) evaluation of the scale (in)dependency of increments, (3) characterization of random processes that describe the evolution of the data, and (4) recognition of cycles embedded within the data series. In the soil sciences, Alemi et al. (1988) use ARIMA and Kriging to model the spatial variation in clay-cover thickness of a 78 km² area in northeastern Iran and demonstrate that ARIMA modeling can adequately describe the nature of the spatial variations. ARIMA models have also been used to model periodicity of major extinction events in the geologic past (Kitchell and Pena, 1984).

In all these studies, the statistical ARIMA model of a given data series is defined by three terms ($p,d,q$), where p and q indicate the order of the autoregressive (AR) and moving average (MA) components, respectively, and d represents a differencing or integration term (I) that is related to LRD. The AR element, p, represents the effects of adjacent observations and the MA element, q, represents the effects on the process of nearby random shocks (Cimino et al., 1999; De Jong and Penzer, 1998). However, in the present study our series are reversible spatial series that can be generated, and are identical, with either forward or backward acquisition, unlike a time series. Both p and q parameters are restricted to integer values (e.g., 0, 1, 2), whereas the integration parameter, d, represents potentially long-range structure in the data. The differencing term d is normally evaluated before p and q to identify whether the process is stationary (i.e., constant mean and σ²). If the series is nonstationary, it is differenced to remove either linear (d = 1) or quadratic (d = 2) trends, thereby making the mean of the series stationary and invertible (Cimino et al., 1999), thus allowing determination of the ARMA p and q parameters.

Here, we adopt the definitions of an ARMA (p,q), and ARIMA ($p,d,q$) process following the work of Beran (1994). Let p and q be integers, where the corresponding polynomials are defined as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\varphi }\left(x\right)=\mathrm{1}-\sum _{j=\mathrm{1}}^p{\mathit{\varphi }}_j{x}^j,\\ \text{(1)}& {\displaystyle }& {\displaystyle }\mathit{\psi }\left(x\right)=\mathrm{1}+\sum _{j=\mathrm{1}}^q{\mathit{\psi }}_j{x}^j.\end{array}$$

It is important to note that all solutions of ϕ(x₀)=0 and ψ(x)=0 are assumed to lie outside the unit circle. Additionally, let ${\mathit{ϵ}}_t(t=\mathrm{1},\mathrm{2},\mathrm{\dots })$ be independent, and identically distributed normal variables with zero variance ${\mathit{\sigma }}_{\mathit{ϵ}}^{\mathrm{2}}$ such that an ARMA (p,q) process is defined by the stationary solution of

$$\begin{array}{}\text{(2)}& \mathit{\varphi }\left(B\right)X_t=\mathit{\psi }\left(B\right){\mathit{ϵ}}_t,\end{array}$$

where B is the backward shift operator $B{X}_t=X_{t-\mathrm{1}},B^{\mathrm{2}}{X}_t=X_{t-\mathrm{1}},\mathrm{\dots }$ and, specifically, the differences can be expressed in terms of B as $X_t-X_{t-\mathrm{1}}=(\mathrm{1}-B)X_t,(X_t-X_{t-\mathrm{1}})-(X_{t-\mathrm{1}}-X_{t-\mathrm{2}})=(\mathrm{1}-B{)}^{\mathrm{2}}{X}_t\mathrm{\dots }$ Alternatively, an ARIMA ($p,d,q$) process X_t is formally defined as

$$\begin{array}{}\text{(3)}& \mathit{\varphi }\left(B\right){\left(\mathrm{1}-B\right)}^d{X}_t=\mathit{\psi }\left(B\right){\mathit{ϵ}}_t,\end{array}$$

where Eq. (3) holds for a dth difference (1−B)^dX_t.

As mentioned previously, a more general form of ARIMA ($p,d,q$) is the fractional ARIMA process, or FARIMA, where the differencing term d is allowed to take on fractional values. If d is a non-integer value for some −0.5 < d < 0.5 and X_t is a stationary process as indicated by Eq.(3), then the model by definition is called a FARIMA process where d values in the range 0 < d < 0.5 are of particular interest herein because geophysically relevant LRD occurs for 0 < d < 0.5, whereas d > 0.5 means that the process is nonstationary but nonintegrable (Beran, 1994; Hosking, 1981). A special case of a FARIMA process explored in the current study is ARIMA (0d0), also known as fractionally differenced white noise (Hosking, 1981), which is defined by Beran (1994) and others as

$$\begin{array}{}\text{(4)}& X_t=(\mathrm{1}-B{)}^{-d}{\mathit{ϵ}}_t.\end{array}$$

For 0 < d < 0.5, the ARIMA (0d0) process is a stationary process with long-range structure and is useful for modeling LRD. As shown later, different values of the d parameter provide further insight into the type of causative physical processes that generate each data series. When d < 0.5, the series X_t is stationary, which has an infinite MA representation that highlights long-range trends or cycles in the data. Conversely, when d > −0.5, the series X_t is invertible and has an infinite AR representation (see Hosking, 1981). When −0.5 < d < 0, the stationary, and invertible, ARIMA (0d0) process is dominated by short-range effects and is anti-persistent. When d=0, the ARIMA (000) process is white noise, with zero correlations and a constant spectral density. Identification of an appropriate model is accomplished by finding small values of elements $p,d,q$ (usually between 0 and 2) that accurately fit the most significant patterns in the data series. When a value of an element is 0, that element is not needed. For example, if d=0 the series does not contain a significant long-range component, whereas if $p=q=\mathrm{0}$, the model does not exhibit significant short-range effects. If $p,d,q\ne \mathrm{0}$, the model contains a combination of both short- and long-memory effects.

Figure 2 The 100 km (a) and 10 km (b) alongshore EMI surveys showing DEMs of the study area and previously identified paleo-channel region by Fisk (1959). Channels are highlighted in red and green, where the green region indicates the location of the 10 km survey. The 25 ft (7.6 m) contour intervals are highlighted with depths increasing from yellow to red and the center of the channels are represented by the black-dotted lines. For each survey, raw σ_a and zero-mean drift-corrected EMI responses are shown in gray and black, respectively. Tidal conditions during each EMI acquisition segment are shown below each panel. Low (lt) and falling tides (ft) are indicated by blue and light blue shades, respectively. High (ht) and rising tides (rt) are highlighted in red and light red, respectively.

4 Results

4.1 Spatial data series

4.1.1 EMI and GPR surveys

The unprocessed (raw) EMI σ_a responses show a high degree of variability along the island. High-amplitude responses within the EMI signal generally exhibit a higher degree of variability (multiplicative noise) compared to the low-amplitude responses. Higher σ_a readings correspond to a small sensor footprint and have enhanced sensitivity to small-scale near-surface heterogeneities (see Guillemoteau and Tronicke, 2015). Low σ_a readings suggest the sensor is probing greater depths and averaging over a larger footprint. In that case, the effect of fine-scale heterogeneities that contribute to signal variability is suppressed.

The 10 km alongshore survey is located within an inferred paleo-channel region (Fisk, 1959), providing some a priori geologic constraints for understanding the variability within the EMI signal (Fig. 2b). Here, the sample size is n=10 176, permitting a quantitative comparison with the 100 km long data series since they contain a similar number of observations. Unlike the 100 km survey, successive footprints of the sensor at each subsequent measurement point overlap along the 10 km survey. The overlap enables a fine-scale characterization of the underlying geological structure because the separation between the TX and RX coils (1.21 m), a good lower-bound approximation of the footprint, is greater than the step size (1 m).

Figure 3 Comparison of EMI σ_a responses from the 100 km survey with 100 MHz GPR data within one of the Fisk (1959) paleo-channels. The 800 m segment (A–A^′) crosses a smaller stream within the network of paleo-channels in the central zone of PAIS. The DOI of the 3 kHz EMI responses is outlined by the red box on the lower GPR radargram and the interpretation of the channel base (ravinement surface) is highlighted in yellow.

The overall trend in σ_a for the 10 km survey is comparable to that of the 100 km survey, where regions characterized by high- and low-amplitude signals correspond to regions of high and low variability, respectively, implying that multiplicative noise persists independently of station spacing. The decrease in σ_a that persists between ∼ 2.5 and 6 km along the profile (Fig. 2b) coincides in location with two paleo-channels, whereas a sharp reduction in σ_a is observed at ∼ 8.2 km in close proximity to a smaller channel. Most of the known paleo-channels are located within the 10 km transect and likely contain resistive infill sands that should generate lower and relatively consistent σ_a readings (Weymer et al., 2015a). The low σ_a signal caused by the sand indirectly indicates valley incision since it is diagnostic of a thicker sand section, relatively unaffected by the underlying conductive layers. Thus, it is reasonable to assume that reduced variability in the signal is related to the framework geology within the paleo-channels, which we now compare with a GPR profile.

Figure 4 DEM metrics extracted from aerial lidar data. The sampling interval (step size) for each data series is 10 m and the coordinates are matched with each EMI acquisition point. Each panel corresponds to (a) beach width, (b) beach volume, (c) dune height, (d) dune volume, (e) island width, (f) island volume, and (g) EMI σ_a. The island is divided into three zones (red vertical lines) roughly indicating the locations within and outside the known paleo-channel region. A Savitzky–Golay smoothing filter was applied to all data series (lidar and EMI) using a moving window of n=250 to highlight the large-scale patterns in each signal.

To corroborate the capability of the EMI data to respond to the variable subsurface geology, an 800 m GPR survey confirms the location of a previously identified paleo-channel (Fisk, 1959) at ∼ 5–10 m depth (Fig. 3). A continuous undulating reflector from ∼ 150 to 800 m along the profile is interpreted to be the surface mapped by Fisk (1959), who documented a paleo-channel at this location with a depth of ∼ 8 m. Although the paleo-surface is within the detection limits of the GPR, it is likely that the DOI of the EMI data (∼ 3–6 m) is not large enough to probe continuously along the contact between the more conductive ravinement surface and the more resistive infill sands. Along the transect at shallower depths highlighted by the red box in the lower radargram (Fig. 3), low EMI σ_a values correspond to fine stratifications in the GPR section, which is common for beach sands with little clay content that are not saline-saturated. The EMI highs between ∼ 450 and 530 m coincide with parts of the GPR section that do not have the fine stratification and this may indicate the presence of clay or saline water. Here, the high conductivity zone for both the GPR and EMI is located within a recovering washover channel overlying the paleo-channel that is evident in the satellite imagery in the upper-left panel of Fig. 3. The overwash deposits consisting of a mix of sand and more finely grained backbarrier sediments likely mask the EMI sensors' ability to probe greater depths. Nonetheless, the high-conductivity zone represents a smaller ∼ 100 m segment within the ∼ 500 m wide paleo-channel, suggesting that variations in the EMI responses outside this zone are directly related to variations in the framework geology imaged by GPR.

4.1.2 Lidar-derived DEM morphometrics

The lidar-derived elevation data series along the 100 km transect are presented in Fig. 4. Each data series is shown with respect to the areal DEM of the study area where the approximate locations of each closely spaced paleo-channel are highlighted in gray. This visualization allows a qualitative analysis of the spatial relationships among paleo-channels, subsurface information encoded in the σ_a signal, and surface morphology over the entire length of the barrier island.

The morphology of the beach–dune system, as well as island width, changes substantially from north to south. In the paleo-channel region, beach width decreases in the central channel (∼ 37–42 km) and is more variable outside this region. Beach width generally increases towards the northern section of the island. The volume of the beach tends to be lowest in the northern zone, varies considerably in the central part of the island, then stabilizes and gradually decreases towards the south. These zones correspond to the southern (0–30 km), central (30–60 km), and northern (60–100 km) sections of the island. Alongshore dune heights are generally greater in the south, become slightly more variable in the paleo-channel region, and decrease in the north except for the area adjacent to Baffin Bay. Dune volume is lowest in the northern section, intermittently increases in the central zone, and slightly decreases towards the south. The island is considerably narrower between Mansfield Channel and Baffin Bay (see Fig. 2a), increasing in width in the northern zone; island volume follows a similar trend. Overall, σ_a values are lower northward of the paleo-channel region compared to the southern zone where σ_a increases substantially. However, the lowest σ_a values are located within the region of paleo-channels inferred by Fisk (1959), supporting previous findings in the study area by Weymer et al. (2015a) and Wernette et al. (2018) that suggest a potential geologic control on alongshore geomorphic features.

Each spatial data series (Fig. 4a–g) represents a different superposition of effects caused by physical processes operating across a wide range of temporal and length scales (Weymer et al., 2015a). Short-range fluctuations represent small-scale heterogeneities, whereas long-range components capture variations in each metric at broader length scales. There is a high degree of variability within each signal that is directly related to the geological and geomorphological structure along the island. Within and outside the paleo-channel region, general associations between EMI σ_a responses and DEM metrics are visually subtle, motivating the statistics we now show with ARIMA modeling. To conduct the ARIMA analysis, we chose to divide the island into three zones based on the location of the known paleo-channels. As will be discussed later, the tripartite zonation allows for a quantitative analysis of LRD at three spatial scales (regional, intermediate, local) within and outside the area containing paleo-channels. It is important to note, however, that the framework geology is likely to exhibit LRD regardless of the length-scale over which it is observed.

Figure 5 Autocorrelations of σ_a for the 100 km (a) and 10 km EMI surveys (d). R∕S analysis for the 100 km (b) and 10 km surveys (e). PSD plots for the 100 km (c) and 10 km surveys (f).

Figure 6 Examples of the worst-fit (a, c) and best-fit (b, d) ARIMA models for the 100 and 10 km EMI surveys. Model results are shown for the processed (drift-corrected) σ_a data. Residuals (RMSE) listed for each model give the standard deviation of the model prediction error. For each plot, original data are in red and fitted (model) data are in blue.

4.2 Tests for LRD

4.2.1 Tests for LRD in EMI data series

Both EMI spatial data series appear to be nonstationary since the mean and variance of the data fluctuate along the profile. A closer visual inspection reveals, however, that cyclicity is present at nearly all spatial frequencies (Fig. 6), with the cycles superimposed in random sequence and added to a constant variance and mean (see Beran, 1994). This behavior is typical for stationary processes with LRD, and is often observed in various types of geophysical time series (Beran, 1992), for example records of Nile River stage minima (Hurst, 1951). A common first-order approach for determining whether a data series contains LRD is through inspection of the autocorrelation function, which we have computed in AutoSignal^™ signal analysis software using a fast Fourier transform (FFT) algorithm (Fig. 5a, d). Both EMI signals exhibit large correlations at large lags (at kilometer and higher scales), suggesting the σ_a responses contain LRD, or long-memory effects in time series language. Results from a rescaled range R∕S analysis (Fig. 5b, e) indeed show high H values of 0.85 (r²=0.98) and 0.95 (r²=0.99) for the 100 and 10 km surveys, indicating a strong presence of LRD at both regional and local spatial scales.

The manner in which different spatial frequency (i.e., wave number) components are superposed to constitute an observed EMI σ_a signal has been suggested to reveal information about the causative multiscale geologic structure (Everett and Weiss, 2002; Weymer et al., 2015a). For example, the lowest-wave-number contributions are associated with spatially coherent geologic features that span the longest length scales probed. The relative contributions of the various wave number components can be examined by plotting the σ_a signal power spectral density (PSD). A power law of the form $\left|{\mathit{\sigma }}_{\mathrm{a}}\right(f){|}^{\mathrm{2}}\sim$ f^β over several decades in spatial wave number is evident (Fig. 5c, f). The slope β of a power-law-shaped spectral density provides a quantitative measure of the LRD embedded in a data series and characterizes the heterogeneity, or roughness, of the signal. A value of $\left|\mathit{\beta }\right|$ > 1 indicates a series that is influenced more by long-range correlations and less by small-scale fluctuations (Everett and Weiss, 2002). For comparison, a pure white noise process would have a slope of exactly β=0, whereas a slope of β ∼ 0.5 indicates fractional Gaussian noise, i.e., a stationary signal with no significant long-range correlations (Everett and Weiss, 2002). The β values for the 100 and 10 km surveys are $\mathit{\beta }=-\mathrm{0.97}$ and $\mathit{\beta }=-\mathrm{1.06}$, respectively. These results suggest that both the 100 and 10 km EMI signals contain long-range correlations. However, there is a slightly stronger presence of LRD within the 10 km segment of the paleo-channel region compared to that within the segment that spans the entire length of the island. This indicates that long-range spatial variations in the framework geology are more important, albeit marginally so, at the 10 km scale than at the 100 km scale. It is possible that the variability within the signal and the degree of long-range correlation is also a function of the sensor footprint, relative to station spacing. This is critically examined in Sect. 4.3.

4.2.2 Tests for LRD in surface morphometrics

Following the same procedure as applied to the EMI data, we performed the R∕S analysis for each beach, dune, and island metric. The calculated H values for the DEM morphometrics range between 0.80 and 0.95 with large values of r² ∼ 1, indicating varying but relatively strong tendencies towards LRD. Beach width and beach volume data series have H values of 0.82 and 0.86, respectively. Dune height and dune volume H values are 0.83 and 0.80, whereas island width and island volume have higher H values of 0.95 and 0.92, respectively. Because each data series shows moderate to strong evidence of LRD, the relative contributions of short- and long-range structure contained within each signal can be further investigated by fitting ARIMA models to each dataset.

4.3 ARIMA statistical modeling of EMI

The results of the tests described in Sect. 4.2.1 for estimating the self-similarity parameter H and the slope of the PSD function suggest that both EMI data series, and by inference the underlying framework geology, exhibit LRD. The goal of our analysis using ARIMA is to estimate the p, d, and q terms representing the order, respectively, of AR, integrated (I), and MA contributions to the signal (Box and Jenkins, 1970) to quantify free vs. forced behavior along the island. For the analysis, the “arfima” and “forecast” statistical packages in R were used to fit a family of ARIMA ($p,d,q$) models to the EMI σ_a data and island morphometrics (Hyndman, 2015; Hyndman and Khandakar, 2008; Veenstra, 2013). Results of 10 realizations drawn from a family of ARIMA ($p,d,q$) models and their residuals (RMSE) are presented in Table 1. The worst fit (ARIMA 001) models are shown for the 100 and 10 km (Fig. 6a, c) surveys. The best fit (ARIMA 0d0) models for both the 100 and 10 km surveys are shown in Fig. 6b and d, respectively. For this analysis, the tests include different combinations of $p,d,$ and q that model either short-range: ARIMA (100; 001; 101; 202; 303; 404; 505), long-range: ARIMA (010; 0d0), or composite short- and long-range processes: ARIMA (111). It is important to note that AR and MA are only appropriate for short-memory processes since they involve only near-neighbor values to explain the current value, whereas the integration (the I term in ARIMA) models long-memory effects because it involves distant values. Note that ARIMA was developed for one-way time series, in which the arrow of time advances in only one direction, but in the current study we are using it for spatial series that are reversible. Different realizations of each ARIMA ($p,d,q$) data series were evaluated, enabling physical interpretations of LRD at regional, intermediate, and local spatial scales. Determining the best-fitting model is achieved by comparing the residual score, or RMSE, of each predicted data series relative to the observed data series, where lower RMSE values indicate a better fit (Table 1).

Table 1 Comparison of residuals (RMSE) of each ARIMA model for the 100 and 10 km EMI surveys.

Table 2 Comparison of residuals (RMSE) of each ARIMA model for all spatial data series. Note that the residuals for each DEM metric correspond to the analysis performed at the regional scale (i.e., 100 km).

Based on the residuals and visual inspection of each realization (Fig. 6), two observations are apparent: (1) both EMI data series are most accurately modeled by an ARIMA (0d0) process with non-integer d, and (2) the mismatch between the data and their model fit is considerably lower for the 10 km survey compared to the 100 km survey. The first observation suggests that the data are most appropriately modeled by a FARIMA process, i.e., a fractional integration that is stationary (0 < d < 0.5) and has long-range dependence (see Hosking, 1981). This implies that spatial variations in framework geology at the broadest scales dominate the EMI signal and that small-scale fluctuations in σ_a caused, for example, by changing hydrological conditions over brief time intervals less than the overall data acquisition interval, or fine-scale lithological variations less than a few station spacings, are not as statistically significant. Regarding the second observation, the results suggest that a small station spacing (i.e., 1 m) is preferred to accurately model both short- and long-range contributions within the signal because large station spacings cannot capture short-range information. The model for the 10 km survey fits better because both p (AR) and q (MA) components increase with a smaller step size since successive volumes of sampled subsurface overlap. On the contrary, the sensor footprint is considerably smaller than the station spacing (10 m) for the 100 km survey. Each σ_a measurement in that case records an independent volume of ground, yet the dataset still exhibits LRD, albeit not to the same degree as in the 10 km survey.

Figure 7 Example of the best fit ARIMA (0d0) models for each lidar-derived DEM metric: (a) beach width, (b) beach volume, (c) dune height, (d) dune volume, (e) island width, and (f) island volume.

Table 3 Summary table showing the computed d parameters that most appropriately model each ARIMA (0d0) iteration (i.e., lowest RMSE).

4.4 ARIMA statistical modeling of island metrics compared with EMI

A sequence of ARIMA ($p,d,q$) models was also evaluated for the elevation morphometrics series to find best fits to the data. The analysis comprised a total of 36 model tests (Table 2). The RMSE values reveal that (1) all data series are best fit by an ARIMA (0d0) process with fractional d, i.e., a FARIMA process; (2) the ARIMA models, in general, more accurately fit the EMI data than the DEM morphometric data likely because the morphology is controlled by more than the framework geology alone; and (3) in all cases, the poorest fit to each series is the ARIMA (001) or MA process. This, in turn, means that the differencing parameter d is the most significant parameter amongst p, d, and q. It is important to note that different values of d were computed based on the best fit of each FARIMA model to the real data. A graphical representation of the FARIMA-modeled data series for each DEM metric is shown in Fig. 7, allowing a visual inspection of how well the models fit the observed data. Because each data series has its own characteristic amplitude and variability, it is not possible to compare RMSE among tests without normalization. The variance within each data series can differ by several orders of magnitude.

Instead of normalizing the data, a fundamentally different approach is to compare the EMI σ_ad values with respect to each metric at regional, intermediate, and local scales (Table 3). Higher positive d values indicate a stronger tendency towards LRD. According to Hosking (1981), {x_t} is called an ARIMA (0d0) process and is of particular interest in modeling LRD as d approaches 0.5 because in such cases the correlations and partial correlations of {x_t} are all positive and decay slowly towards zero as the lag increases, while the spectral density of {x_t} is concentrated at low frequencies. It is reasonable to assume that the degree of LRD may change over smaller intermediate and/or local scales, which implies a breakdown of self-similarity. For a self-similar signal, d is a global parameter that does not depend on which segment of the series is analyzed. In other words, the d values should be the same at all scales for a self-similar structure.

The results of the FARIMA analysis at the intermediate scale vary considerably within each zone of the barrier island (north, central, south) and for each spatial data series (Table 3). In the southern zone (0–30 km), EMI σ_a and beach volume have the strongest LRD (d=0.44), whereas the other metrics exhibit weak LRD (ranging from d ∼ 0 to 0.2), which may be characterized approximately as a white noise process. Within the paleo-channel region (30–60 km), all of the island metrics show a moderate to strong tendency towards LRD (0.3 ≤ d ≤ 4.2); however, the EMI signal does not (d=0.11). In the northern zone (60–100 km) all data series contain moderate to strong LRD with the exception of beach and island width.

A FARIMA analysis was also conducted at the local scale by dividing the island into 10 km segments, starting at the southern zone (0–10 km) and ending at the northern zone of the island (90–100 km). A total of 70 FARIMA model realizations were evaluated and the resulting d values demonstrate that the EMI data segments show a stronger presence of LRD (d > 0.4) within the paleo-channels (30–60 km) and further to the north (60–80 km) in close proximity to the ancestral outlet of Baffin Bay. These findings indicate that there may be local and/or intermediate geologic controls along different parts of the island, but that the framework geology dominates island metrics at the regional scale.

5 Discussion

Although it has long been known that processes acting across multiple temporal and length scales permit the shape of coastlines to be described by mathematical constructs such as power-law spectra and fractal dimension (Lazarus et al., 2011; Mandelbrot, 1967; Tebbens et al., 2002), analogous studies of the subsurface framework geology of a barrier island have not been carried out. This research supports previous studies demonstrating that near-surface EMI geophysical methods are useful for mapping barrier island framework geology and that FARIMA data series analysis is a compact statistical tool for illuminating the long and/or short-range spatial correlations between subsurface geology and geomorphology. The results of the FARIMA analysis and comparisons of the best-fitting d parameters show that beach and dune metrics closely match EMI σ_a responses regionally along the entire length of PAIS, suggesting that the long-range dependent structure of these data series is similar at large spatial scales. However, further evaluation of the d parameters over smaller data segments reveals that there are additional localized framework geology controls on island geomorphology that are not present at the regional scale.

At the intermediate scale, a low EMI d value (d=0.11) suggests there is only a weak framework geologic control on barrier island morphometrics. A possible explanation is that the paleo-channels, located within a ∼ 30 km segment of the island, are not regularly spaced and on average are less than a few kilometers wide. This implies that the framework geology controls are localized (i.e., effective in shaping island geomorphology only at smaller spatial scales). At the local scale, relationships between the long-range dependence of EMI and each metric vary considerably, but there is a significant geologic control on dune height within the paleo-channel region (d > 0.4). It is hypothesized that the alongshore projection of the geometry of each channel is directly related to a corresponding variation in the EMI signal, such that large, gradual minima in σ_a are indicative of large, deep channel cross sections and small, abrupt minima in σ_a represent smaller, shallow channel cross sections. At shallower depths within the DOI probed by the EMI sensor, variability in the σ_a signal may correspond to changes in sediment characteristics as imaged by GPR (Fig. 3). Located beneath a washover channel, a zone of high-conductivity EMI σ_a responses between ∼ 450 and 530 m coincides with a segment of the GPR section where the signal is more attenuated and lacks the fine stratification that correlates much better with the lower σ_a zones. The contrasts in lithology between the overwash deposits and stratified infilled sands were detected by both EMI and GPR measurements.

It is argued herein that differences in the d parameter between EMI σ_a readings (our assumed proxy for framework geology) and lidar-derived surface morphometrics provide a new metric that is useful for quantifying the causative physical processes that govern island transgression across multiple spatial scales. All of the calculated d values in this study are derived from ARIMA (0d0) models that fit the observations, and lie within the range of 0 < d < 0.5, suggesting that each data series is stationary but does contain long-range structure that represents randomly placed cyclicities in the data. For all models in our study, the d values range between ∼ 0 and 0.50, which enables a geomorphological interpretation of the degree of LRD and self-similarity at different spatial scales. In other words, the d parameter not only provides an indication of the scale dependencies within the data, but also offers a compact way of analyzing the statistical connections between forced (stronger d ∼ 0.5) and free (weaker d ∼ 0) behavior that may be more influenced by morphodynamic processes operating at smaller spatial scales.

Alongshore variations in beach width and dune height are not uniform at PAIS (Weymer et al., 2015b) and exhibit different spatial structure within and outside the paleo-channel region (Fig. 5). These dissimilarities may be forced by the framework geology within the central zone of the island but are influenced more by contemporary morphodynamic processes outside the paleo-channel region. This effect could be represented by higher-wave-number components embedded within the spatial data series. Beach and dune morphology in areas that are not controlled by framework geology (e.g., the northern and southern zones) exhibit more small-scale fluctuations representing a free system primarily controlled by contemporary morphodynamics (e.g., wave action, storm surge, wind).

Because variations in dune height exert an important control on storm impacts (Sallenger, 2000) and ultimately large-scale island transgression (Houser, 2012), it is argued here that the framework geology (or lack thereof) of PAIS acts as an important control on island response to storms and sea level rise. This study supports recent work by Wernette et al. (2018) suggesting that framework geology can influence barrier island geomorphology by creating alongshore variations in either oceanographic forcing and/or sediment supply and texture that controls smaller-scale processes responsible for beach–dune interaction at the local scale. The forced behavior within the paleo-channel region challenges shoreline change studies that consider only small-scale undulations in the dune line that are caused by natural randomness within the system. Rather, we propose that dune growth is forced by the framework geology, whose depth is related to the thickness of the modern shoreface sands beneath the beach. This depth is the primary quantity that is detected by the EMI sensor. With respect to shoreline change investigations, improving model performance requires further study of how the framework geology influences beach–dune morphology through variations in wave energy, texture, and sediment supply (e.g., Houser, 2012; McNinch, 2004; Schwab et al., 2013).

Our findings extend previous framework geology studies from Outer Banks, NC (e.g., Browder and McNinch, 2006; McNinch, 2004; Riggs et al., 1995; Schupp et al., 2006), Fire Island, NY (e.g., Hapke et al., 2010; Lentz and Hapke, 2011), and Pensacola, FL (e.g., Houser, 2012), where feedbacks between geologic features and relict sediments within the littoral system have been shown to act as an important control on dune growth and evolution. Nonetheless, most of these studies focus on offshore controls on shoreface and/or beach–dune dynamics at either local or intermediate scales because few islands worldwide exist that are as long and/or continuous as North Padre Island. To our knowledge, few framework geology studies have specifically used statistical testing to analyze correlations between subsurface geologic features and surface morphology. Two notable exceptions include Browder and McNinch (2006), and Schupp et al. (2006), both of which used chi-squared testing and cross-correlation analysis to quantify the spatial relationships among offshore bars, gravel beds, and/or paleo-channels at the Outer Banks, NC. Although these techniques are useful for determining spatial correlations among different datasets, they do not provide information about the scale (in)dependencies between the framework geology and surface geomorphology that FARIMA models are better designed to handle. The current study augments the existing literature in that (1) it outlines a quantitative method for determining free and forced evolution of barrier island geomorphology at multiple length scales, and (2) it demonstrates that there is a first-order control on dune height at the local scale within an area of known paleo-channels, suggesting that framework geology controls are localized within certain zones of PAIS.

Further study is required to determine how this combination of free and forced behavior resulting from the variable and localized framework geology affects island transgression. Methods of data analysis that would complement the techniques presented in this paper might include power spectral analysis, wavelet decomposition, and shoreline change analysis that implicitly includes variable framework geology. These approaches would provide important information regarding (1) coherence and phase relationships between subsurface structure and island geomorphology and (2) nonlinear interactions of coastal processes across large and small spatiotemporal scales. Quantifying and interpreting the significance of framework geology as a driver of barrier island formation and evolution and its interaction with contemporary morphodynamic processes is essential for designing and sustainably managing resilient coastal communities and habitats.

6 Conclusions

This study demonstrates the utility of EMI geophysical profiling as a new tool for mapping the length-scale dependence of barrier island framework geology and introduces the potential of FARIMA analysis to better understand the geologic controls on large-scale barrier island transgression. The EMI and morphometric data series exhibit LRD to varying degrees, and each can be accurately modeled using a non-integral parameter d. The value of this parameter diagnoses the spatial relationship between the framework geology and surface geomorphology. At the regional scale (∼ 100 km), small differences in d between the EMI and morphometrics series suggest that the long-range-dependent structure of each data series with respect to EMI σ_a is statistically similar. At the intermediate scale (∼ 30 km), there is a greater difference among the d values of the EMI and island metrics within the known paleo-channel region, suggesting a more localized geologic control with less contributions from broader-scale geological structures. At the local scale (10 km), there is a considerable degree of variability among the d values of the EMI and each metric. These results all point toward a forced barrier island evolutionary behavior within the paleo-channel region transitioning into a free, or scale-independent, behavior dominated by contemporary morphodynamics outside the paleo-channel region. In a free system, small-scale undulations in the dune line reinforce natural random processes that occur within the beach–dune system and are not influenced by the underlying geologic structure. In a forced system, the underlying geologic structure establishes boundary constraints that control how the island evolves over time. This means that barrier island geomorphology at PAIS is forced and scale dependent, unlike shorelines which have been shown at other barrier islands to be scale independent (Tebbens et al., 2002; Lazarus et al., 2011). The exchange of sediment amongst nearshore, beach, and dune in areas outside the paleo-channel region is scale independent, meaning that barrier islands like PAIS exhibit a combination of free and forced behaviors that will affect the response of the island to sea level rise and storms. We propose that our analysis is not limited to PAIS but can be applied to other barrier islands and potentially in different geomorphic environments, both coastal and inland.