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 101 to 102 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.


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;Mc-Ninch, 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).
Published by Copernicus Publications on behalf of the European Geosciences Union.B. A. Weymer et al.: Long-range-dependent structure of barrier island framework geology

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(Hapke et al., , 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.
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).

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 selfsimilar 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 powerlaw 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).

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 paleochannel (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 (0, d, 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

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 pointwise 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.

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.
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 crosssectional 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.

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.

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 0 -10 5 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 col-lected 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.4m 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)  ison with the σ a measurements.We used a Sensors and Software pulseEKKO Pro ™ system for this purpose.A surveygrade 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 −1 ), 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.

Geomorphometry
Topographic information was extracted from aerial lidar data that were collected by the U.S.  (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.

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(Beran, , 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 σ 2 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, 1992Beran, , 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 2 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 obser-vations 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 σ 2 ).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 (1) It is important to note that all solutions of φ(x 0 ) = 0 and ψ(x) = 0 are assumed to lie outside the unit circle.Additionally, let t (t = 1, 2, . ..) be independent, and identically distributed normal variables with zero variance σ 2 such that an ARMA (p, q) process is defined by the stationary solution of where B is the backward shift operator BX t = X t−1 , B 2 X t = X t−1 , . . .and, specifically, the differences can be expressed in terms of B as Alternatively, an ARIMA (p, d, q) process X t is formally defined as where Eq. ( 3) holds for a dth difference (1 − B) d X 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 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 antipersistent.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 = 0, the model does not exhibit significant short-range effects.If p, d, q = 0, the model contains a combination of both shortand long-memory effects.

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 charac-terization 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).
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 paleochannels, 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.
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 paleochannel that is evident in the satellite imagery in the upperleft 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.

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 www.earth-surf-dynam.net/6/431/2018/Earth Surf.Dynam., 6, 431-450, 2018 the southern zone where σ a increases substantially.However, the lowest σ a values are located within the region of paleochannels 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 paleochannel 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.

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 cyclic-ity 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 2 = 0.98) and 0.95 (r 2 = 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 |σ a (f )| 2 ∼ 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 |β| > 1 indicates a series Earth Surf.Dynam., 6, 431-450, 2018 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 β = −0.97 and β = −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 iswww.earth-surf-dynam.net/6/431/2018/Earth Surf.Dynam., 6, 431-450, 2018 land.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.

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 2 ∼ 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.

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 shortand 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  1).
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.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.

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.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.

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 geophysi-cal 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 shortrange 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 longrange 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 paleochannel 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 largescale 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 andHapke, 2011), andPensacola, 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 andMcNinch (2006), andSchupp 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 paleochannels 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.

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 beachdune 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 www.earth-surf-dynam.net/6/431/2018/Earth Surf.Dynam., 6, 431-450, 2018 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.

Figure 1 .
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.
B. A. Weymer et al.: Long-range-dependent structure of barrier island framework geology

Figure 2 .
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.

Figure 3 .
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.

Figure 4 .
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.

Figure 5 .
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 .
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.

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).
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 σ a d 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

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