Geostatistic Analysis of Long-Term Subsidence effects on Isostasy of Long Island Applying Continuous GPS Displacement Rates using R Parag, A. S.1, Marsellos, A.E.1, Tsakiri, K.G.2 Dept. Geology, Environment and Sustainability, Hofstra University, NY, USA. . 2 Dept. Information Systems, Analytics, and Supply Chain Management, Rider University, NJ, USA. 1
Abstract:
Method/Materials:
Long Island has experienced subsidence since the last ice age when the large Laurentide ice sheet covered most of New York State and readjusted after it melted. This process of readjustment of the land is known as isostasy, the equilibrium position where a body of land should sit between the surface above and below. However, the region is not uplifting after the weight of the ice melted but rather continues to subside even after this event took place. Modern human environmental interferences, particularly groundwater overpumping, causes uneven land deformation that heightens risks of landslides, sinkholes, and other geohazards (U.S. Geological Survey, n.d.). Central Long Island, with its many residential properties, has a high demand for groundwater, leading to the overpumping of aquifers in certain areas. This excessive withdrawal may contribute to gradual subsidence in parts of the region. In contrast, areas within New York City and Montauk Point exhibit significant subsidence rates despite lower residential groundwater extraction. This is likely due to reliance on upstate water sources, heavy infrastructure, and high population density. While Montauk Point is less populated than central and western Long Island, it still experiences notable subsidence, possibly due to reduced aquifer pumping, which allows underlying geological formations to compact under their own weight. To capture these spatial variations in land motion, we analyzed Figure 1: These geostatistical graphs display the up components of each station and their vertical continuous GPS records from multiple stations movements within their various ranges. Some stations range from 1995-1998 and variations within the 2000s. The x-axis is time in years and the y-axis is the vertical displacement in meters (m). across differing time intervals. Rather than restricting our study to a common observation window, we considered each station’s full operational timeline to assess overall trends. A GIS-interpolated map of average displacement rates reveals localized areasGeostatistical of Analysis: significant subsidence, illuminating the ➢ Displays vertical movements of stations combined effects of natural isostatic ➢ Long Island been subsiding stresses.over time (1995-1998 & 2000s). processes andhashuman-induced since the last ice dueforto strategic Findings support the age need Subsidence Trends: groundwater management and infrastructure glacial isostatic adjustment. ➢ Highest subsidence (~-2 mm/yr) near NYC & planning reduce causes further land land to subsidence ➢ Glacial to isostasy and its associated in this vulnerableeastern Long Island; central Long Island subside under hazards ice weight and coastal setting. shows lower subsidence with occasional uplift when the ice melts. uplift. ➢ Unlike expected uplift, the region Figure 2: These graphs display the up, north, and east components of each station and their oscillations (various colored lines) within their different continues to subside, suggesting Station Observations: ranges. Some stations range from 1995-1998 and variations within the 2000s. additional contributing factors. The x-axis is the period (in days) and the y-axis is the vertical displacement in ➢ NYBR (NYC): Strong long-term subsidence. meters (m). ➢ Groundwater overpumping may ➢ MNP1 (Eastern Long Island): Gradual lead to uneven land deformation, downward trend but higher subsidence increasing geohazard risks like rates than central stations. landslides and sinkholes formed ➢ NYEL, NYPD, ZNY1 (Central Stations): Lower by underground erosion and subsidence rates with slight uplift periods. sediment collapse. Wavelet Spectrum Analysis: ➢ This study uses GPS data to detect ongoing vertical land ➢ Up Component: Multiple peaks, indicating movements beyond 2005–2020 periodic oscillations in vertical research. displacement. Clustering of peaks suggests Figure 3: This GIS map displays the vertical displacement rates of each station and reveals potential areas of subsidence or uplift. The light red and white ➢ GIS-interpolated maps highlight seasonal and sub-annual cycles. colored areas reveal areas of lower subsidence, the dark red colored areas localized subsidence patterns, ➢ North Component: Low but visible longreveal areas of uplift, and the blue colored areas reveal areas of high
Background of Study:
Results:
subsidence. Long Island map is behind the vertical displacement rates.
➢ GPS displacement data from 12 Long Island stations were retrieved using R from the Nevada Geodetic Laboratory. (Fig. 1) ➢ Data processing and cleaning used R packages. ➢ Custom function (GPS_dwd) structured data into northward, eastward, and upward displacement components. ➢ Missing values were interpolated, and stations with excessive data issues (MOR6, QYNS) were excluded. ➢ Outliers were removed using Z-score normalization (values beyond ±3 standard deviations). ➢ Noise reduction applied using Kolmogorov-Zurbenko (KZ) filtering (365day window, 3 replications). ➢ Wavelet Spectrum (WaveletComp) detected periodic signals in Up, North, and East displacement components. (Fig. 2) ➢ Linear regression (lm function) calculated vertical displacement rates from smoothed data. ➢ Geostatistical interpolation (IDW method) mapped subsidence patterns using the gstat and sp packages. ➢ GIS visualization (leaflet package) created an interactive map showing GPS station locations and displacement rates. (Fig. 3)
Acknowledgements:
I would like to thank Dr. Antonios E. Marsellos and the Geology, Environment and Sustainability department at Hofstra University for providing me with the necessary tools to conduct this research. I would also like to thank Dr. Tsakiri for providing me with insights on KZ filtering script.
Credit Authorship Contribution Statement:
Parag, A.: Literature, Editing, RStudio Coding, Writing- Abstract, Introduction, Methods, Results-Figures, Discussion, Conclusion; Marsellos, A.E.: Revision, Editing, RStudio Coding; Tsakiri, K.G.: KZ script insights
Discussion/ Conclusion:
➢ NYC has the highest subsidence due to infrastructure weight and limited groundwater extraction. ➢ Central Long Island shows slight uplift, likely from groundwater overpumping. ➢ Eastern Long Island subsides despite low population, possibly due to retained groundwater weight. ➢ Subsidence may be influenced by groundwater and human activity, not tectonic factors. ➢ GPS data analysis confirms sinking rates, averaging 1.6 mm/year, with some areas exceeding 3 cm/decade. ➢ KZ filtering (365-day) removes noise, preserving long-term displacement trends. ➢ Subsidence varies across regions, with NYC and Montauk Point sinking more than central Long Island. ➢ Wavelet spectrum analysis identifies periodic land displacement possibly influenced by groundwater, temperature, and subsidence. ➢ Vertical (Up) motion shows stronger
References:
Badger, T., Parag, A., Pelletier, A., DeRocchis, S., Marsellos, A. E., & Tsakiri, K. G. (2024). Investigating groundwater and landslides relationship using R: A study on Long Island’s north shore. In 31st Conference on the Geology of Long Island and Metropolitan New York (Abstract). Stony Brook University. Baker, A., et al. (2005). The Impact of Climate Change on Groundwater Resources: A Global Review. Global and Planetary Change, 47(2-4), 221-242. Elsevier. https://doi.org/10.1016/S0921-8181(01)00150-3. Barrett, T., Dowle, M., Srinivasan, A., Gorecki, J., Chirico, M., Hocking, T., & Schwendinger, B. (2024). data.table: Extension of data.frame. R package version 1.16.4. Retrieved from https://CRAN.R-project.org/package=data.table. Bishop, P. (2007). Long-term landscape evolution: linking tectonics and surface processes. Earth Surface Processes and Landforms, 32, 329-365. https://doi.org/10.1002/esp.1493. Bivand, R., Pebesma, E., & Gomez-Rubio, V. (2013). Applied Spatial Data Analysis with R (2nd ed.). Springer, NY. Retrieved from https://asdar-book.org/. Blewitt, G., & Lavallée, D. (2002). Effect of annual signals on geodetic velocity. Journal of Geophysical Research, 107(B7). Bloom, A. L. (1967). Pleistocene Shorelines: A New Test of Isostasy. GSA Bulletin, 78(12), 1477–1494. https://doi.org/10.1130/0016-7606(1967)78[1477:PSANTO]2.0.CO;2. Cheng, J., Schloerke, B., Karambelkar, B., & Xie, Y. (2024). leaflet: Create Interactive Web Maps with the JavaScript 'Leaflet' Library. R package version 2.2.2. Retrieved from https://CRAN.R-project.org/package=leaflet. Close, B., Zurbenko, I., & Sun, M. (2020). kza: Kolmogorov-Zurbenko Adaptive Filters. R package version 4.1.0.1. Retrieved from https://CRAN.R-project.org/package=kza. Davis, G. H. (1987). Land subsidence and sea level rise on the Atlantic Coastal Plain of the United States. Environmental Geology and Water Sciences, 10, 67–80. https://doi.org/10.1007/BF02574663. Dong, D., Fang, P., Bock, Y., Cheng, M. K., & Miyazaki, S. (2002). Anatomy of apparent seasonal variations from GPSderived site position time series. Journal of Geophysical Research, 107(B4). Flavelle, C., & Rojanasakul, M. (2024, February 13). Flooding Isn’t Just About Rising Seas. Groundwater Is a Threat, Too. The New York Times. Retrieved from https://www.nytimes.com/interactive/2024/02/13/climate/floodingsea-levels-groundwater.html?auth=login-google1tap&login=google1tap. GPS data provided by the University of Nevada, Reno Geodetic Laboratory. Retrieved from https://geodesy.unr.edu/gps_timeseries/tenv3/IGS14/. Gräler, B., Pebesma, E., & Heuvelink, G. (2016). Spatio-Temporal Interpolation using gstat. The R Journal, 8(1), 204218. Grolemund, G., & Wickham, H. (2011). Dates and Times Made Easy with lubridate. Journal of Statistical Software, 40(3), 1-25. Retrieved from https://www.jstatsoft.org/v40/i03/. Hijmans, R. J. (2025). raster: Geographic Data Analysis and Modeling. R package version 3.6-31. Retrieved from https://CRAN.R-project.org/package=raster. Hope, J. H., Marsellos, A. E., & Tsakiri, K. G. (2024). High-frequency groundwater level fluctuations, underground erosion, potential sinkhole occurrences across Long Island, NY. In 31st Conference on the Geology of Long Island and Metropolitan New York (Abstract). Stony Brook University. Larson, K. M., & van Dam, T. (2000). Measuring postglacial rebound with GPS and absolute gravity. Geophysical Research Letters, 27(23), 3925-3928. Moore, M., Brown, R., Thomas, M., Santella, K., Ballato, B., Marsellos, A. E., & Tsakiri, K. G. (2023). What can the local GPS network tell us about the Glacial Isostatic Adjustment of Long Island, NY? Retrieved from https://ehs.stonybrook.edu/commcms/geosciences/about/_LIG-Past-Conference-abstract-pdfs/2021Abstracts/Moore.pdf. Moritz, S., & Bartz-Beielstein, T. (2017). imputeTS: Time Series Missing Value Imputation in R. The R Journal, 9(1), 207-218. https://doi.org/10.32614/RJ-2017-009. NASA Jet Propulsion Laboratory. (2023, September 27). NASA-Led Study Pinpoints Areas of New York City Sinking, Rising. Retrieved from https://www.jpl.nasa.gov/news/nasa-led-study-pinpoints-areas-of-new-york-citysinking-rising/. National Oceanic and Atmospheric Administration. (n.d.). Sea level rise viewer. NOAA Office for Coastal Management. Retrieved from https://coast.noaa.gov/digitalcoast/tools/slr.html. Pebesma, E. J. (2004). Multivariable geostatistics in S: the gstat package. Computers & Geosciences, 30, 683-691. Pebesma, E., & Bivand, R. (2005). Classes and methods for spatial data in R. R News, 5(2), 9-13. Retrieved from https://CRAN.R-project.org/doc/Rnews/. R Core Team. (n.d.). The Comprehensive R Archive Network (CRAN). R Foundation for Statistical Computing. Retrieved March 31, 2025, from https://cran.r-project.org/ Roesch, A., & Schmidbauer, H. (2018). WaveletComp: Computational Wavelet Analysis. R package version 1.1. Retrieved from https://CRAN.R-project.org/package=WaveletComp. Shepard, C. C., Agostini, V. N., Gilmer, B., et al. (2012). Assessing future risk: quantifying the effects of sea level rise on storm surge risk for the southern shores of Long Island, New York. Natural Hazards, 60, 727–745. https://doi.org/10.1007/s11069-011-0046-8. Taniguchi, M. (2001). Evaluation of the Groundwater Recharge Rate through the Unsaturated Zone Based on the Temporal Change in Temperature. Ecological Modelling, 138(1-3), 291-301. https://doi.org/10.1016/S09218181(01)00150-3. Tsakiri, K.G., Zurbenko, I.G., 2011a. Effect of noise in principal component analysis. J. Stat. Math. 2, 40. Tsakiri, K.G., Zurbenko, I.G., 2011b. Prediction of ozone concentrations using atmospheric variables. Air Qual. Atmos. Health 4, 111–120. Tsakiri, K., Marsellos, A., Kapetanakis, S., 2018. Artificial Neural Network and Multiple Linear Regression for Flood