Hydrometeorological Trends in a Low-Gradient Forested Watershed on the Southeastern Atlantic Coastal Plain in the USA

: Hydrology and meteorological data from relatively undisturbed watersheds aid in identifying effects on ecosystem services, tracking hydroclimatic trends, and reducing model uncertainties. Sustainable forest, water, and infrastructure management depends on assessing the impacts of extreme events and land use change on flooding, droughts, and biogeochemical processes. For example, global climate models predict more frequent high-intensity storms and longer dry periods for the southeastern USA. We summarized 17 years (2005–2021) of hydrometeorological data recorded in the 52 km 2 , third-order Turkey Creek watershed at the Santee Experimental Forest (SEF), Southeastern Coastal Plain, USA. This is a non-tidal headwater system of the Charleston Harbor estuary. The study period included a wide range of weather conditions; annual precipitation (P) and potential evapotranspiration (PET) ranged from 994 mm and 1212 mm in 2007 to 2243 mm and 1063 in 2015, respectively. The annual runoff coefficient (ROC) varied from 0.09 in 2007 (with water table (WT) as deep as 2.4 m below surface) to 0.52 in 2015 (with frequently ponded WT conditions), with an average of 0.22. Although the average P (1470 mm) was 11% higher than the historic 1964–1976 average (1320 mm), no significant ( α = 0.05) trend was found in the annual P ( p = 0.11), ROC ( p = 0.17) or runoff ( p = 0.27). Runoff occurred on 76.4% of all days in the study period, exceeding 20 mm/day for 1.25% of all days, mostly due to intense storms in the summer and lower ET demand in the winter. No-flow conditions were common during most of the summer growing season. WT recharge occurred during water-surplus conditions, and storm-event base flow contributed 23–47% of the total runoff as estimated using a hydrograph separation method. Storm-event peak discharge in the Turkey Creek was dominated by shallow subsurface runoff and was correlated with 48 h precipitation totals. Estimated precipitation intensity–duration–frequency and flood frequency relationships were found to be larger than those found by NOAA for the 1893–2002 period (for durations ≥ 3 h), and by USGS regional frequencies (for ≥ 10-year return intervals), respectively, for the same location. We recommend an integrated analysis of these data together with available water quality data to (1) assess the impacts of rising tides on the hydroperiod and biogeochemical processes in riparian forests of the estuary headwaters, (2) validate rainfall–runoff models including watershed scale models to assess land use and climate change on hydrology and water quality, and (3) inform watershed restoration goals, strategies, and infrastructure design in coastal watersheds.


Introduction
This study analyzed the most recent 17 years' (2005-2021) of hydro-meteorological data in the Turkey Creek watershed, a 52 km 2 forested area in coastal South Carolina that lies adjacent to long-term study areas within the Santee Experimental Forest (SEF) and the Francis Marion National Forest.The purpose was to synthesize water budget data and derive relationships for precipitation intensity-duration-frequencies (PIDF) and flow duration and frequencies.Long-term hydrology, climate, and species-specific land-cover data from headwater forested watersheds are important in the geo-environmental and ecological sciences as well as in modeling for environmental impact studies.In September of 1989, a major hurricane (Hugo) affected this area, such that forest regeneration and the hydrology conditions were reportedly recovered only by 2004 [1].We used data from this watershed-scale monitoring system to answer two key questions.For example, (a) how will a changing climate affect the hydrology (particularly streamflow and ET) of this relatively less-intensively managed forest ecosystem?and (b) how have runoff and peak discharges, which affect stormwater and road infrastructure design, responded, and what could change in the future when high intensity storm events are forecast for the region?

Rationale
The U.S. Environmental Protection Agency (EPA) estimates that 5260 km 2 (2029 mi²) of forest will be cleared for urban development in South Carolina by the year 2050 [2].Emmett [3] noted that harvesting is the predominant agent of forest canopy cover loss in the southeastern US, affecting 42% of the region's forest area from 1986-2010, based on an analysis of the North American Forest Dynamics-Attribution dataset (NAFD-ATT) [4].Urbanization, in some cases after forest clearing, results in more impervious cover which increases surface runoff due to reduced soil-water infiltration and groundwater recharge [2,5,6] as well as evapotranspiration (ET) [7].Schueler [8] noted in their review of urban stream research that as little as 10% of impervious cover in a watershed may alter hydrological processes and impair water quality.Nagy and Lockaby [9] noted that rising land values, increased market concentrations, increased road accessibility, and ultimately higher opportunity costs for maintaining forested land contributed to the conversion of forested to urban land and degradation of aquatic and terrestrial systems in Georgia, USA.These cultural and socioeconomic factors associated with increased urbanization are common globally and are likely to continue in the southeastern United States [10].The Charleston Metropolitan Area in South Carolina is one such urban center, and development within which threatens adjacent, less developed watersheds in coastal South Carolina [11].For example, long-term watersheds at the US Forest Service Santee Experimental Forest within the Francis Marion National Forest are at the wildland-urban interface of the urbanizing Charleston area [12].
Long-term hydrologic, climatic, and species-specific land-cover data have a multitude of uses and applications.These include modeling and quantifying ecohydrological processes controlling water flow paths and residence time [13][14][15], assessment of the impacts of land use change and forest disturbance on flooding and droughts [16] and conservation of ecosystems [17,18], risk assessment tools for extreme events [19,20], and a consideration of societal needs [15,21,22].
Past studies have used long-term data from small, paired, upland forested watersheds to examine ecohydrologic processes including hillslope process-based runoff generation mechanisms (for example see [23,24]).However, some observational studies have been carried out at only limited sites on the low-gradient southeastern Coastal Plain forested watersheds, consisting of wetland forests and pine flatwoods, and where the position of a shallow water table driven by rainfall and ET generally dictates runoff [25][26][27][28][29][30][31][32][33][34][35].Most of those studies were focused on relatively small first-order watersheds on the Atlantic Coastal Plain and analyses using long-term data from the larger watersheds are scarce but are becoming increasingly important as the recent focus on water quantity and quality has rather shifted to a landscape or regional scales [36][37][38].In addition, results reported in literature for effects of afforestation and deforestation on streamflow and peak flow in those small watersheds cannot be simply extrapolated to large watersheds [38,39] because of more complexities within land use (e.g., managed, and natural forests, wetlands, lakes, open lands, urban use) and their interactions [38].

Background and Related Work
Previous studies focused on the water budget [40,41], water budget modeling [40,42], event rainfall-runoff analysis using historical data [43,44], and the water quality [41,45,46] of a third-order, relatively undisturbed, Turkey Creek watershed at the Santee Experimental Forest.The major disturbance to the forest occurred due to Hurricane Hugo that decimated the tree stands [47].Most of these earlier studies presumed the watershed area was much larger than the more recent LiDAR-based estimate (7260 ha vs. 5240 ha) [48].Amatya et al. [49] presented additional ecohydrological and socioeconomic characteristics and discussed key elements for watershed management including water resource assessment, modeling integrated water resources management, environmental and social impact assessment, and land use planning on such a typical coastal forest watershed.Similarly, additional short-term ecohydrological studies on this and adjacent sites were also synthesized by Amatya et al. [50].
Amatya et al. [49] also highlighted several key scientific questions left for future potential studies that would help address multiple hydrology and water resources issues on this rapidly urbanizing forested landscape.One question was "what are the uncertainties associated with current knowledge of surface and groundwater (SGW) interactions and their partitioning that affect the key ecosystem services like water supply, and ultimately water quality, and biogeochemical processes and transport."The SGW part was partially addressed in studies by Callahan et al. [51] and Amatya et al. [52], which used a modified version of a Curve Number-based rainfall-runoff model [53] using storm event data recorded in two first-order sub-catchments nested within Turkey Creek and other similar headwater catchments in the region.The authors concluded that the soil saturation coefficient, as a proxy of antecedent moisture condition, is a key parameter used to describe the partitioning of runoff into overland surface and shallow subsurface runoff.Similarly, Amatya et al. [25] attempted to answer the questions "what are the mechanisms and rates of feedback between topography, hydrology, ecology (vegetation) and how do these vary as a function of local climate using a paired watershed data?"Furthermore, two studies addressed the questions "what is/are the cause(s) of the high methylization efficiency of mercury (Hg) in these streams, and how do forest management practices affect Hg transport and methylization [54]?" Silvicultural practices such as prescribed burning and mechanical thinning which are routinely undertaken in this study watershed interfere with Hg storage in forests and potentially increase Hg export downstream.Those areas can be hotspots for transformation into methylmercury (MeHg) facilitated by anaerobic microbes.Therefore, an accurate quantification of streamflow is needed for estimating impacts of the forest management on Hg export, as a function of streamflow, in these forested watersheds.Lastly, a project focused on the question "how does the interaction of downstream tidal energy with the Turkey Creek watershed flow affect the hydrological and biogeochemical processes including net primary productivity, carbon balance and greenhouse gas emissions?"[55].The motivation was that tidally mediated processes downstream are dependent upon the seasonal flow dynamics of this headwater system.

Study Objectives
The first objective of this study was to analyze the most recent 17 years of (2005-2021) rainfall, stream flow, peak flow rates, potential ET (PET), ET, and groundwater data for the trends and patterns in their responses to rainfall observed on the study watershed (Turkey Creek), and to compare them with those reported in earlier studies for this watershed to assess the potential effects of climate change, such as more frequent high-intensity storm events.The second objective was to derive the daily rainfall and flow frequency duration curves, PIDFs and flood frequencies that are generally used in watershed and stream restoration projects, low flow monitoring to protect stream aquatic habitats, and to inform the design and management of stormwater, roadway, and drainage infrastructure.

Study Site
The study site is the Turkey Creek watershed (WS78) in the Santee Experimental Forest (Figure 1) originally established by the USDA Forest Service in 1964 and monitored through 1984.Recognizing the importance of reference watershed data in a rapidly changing coastal environment, stream gauging was reestablished in Turkey Creek in 2004 to facilitate a large-scale eco-hydrological monitoring and modeling program [56].The present WS-78 gauging station is approximately 800 m upstream of the original location; it is instrumented with a real-time stream gauge sensor and a rain gauge (http://waterdata.usgs.gov/sc/nwis/uv?site_no=02172035, accessed on 20 February 2024).The gauging station is operated in cooperation with the US Geological Survey (USGS) and the College of Charleston.This paper presents the analysis of data collected from 2005 to 2021.
watershed to assess the potential effects of climate change, such as more frequent highintensity storm events.The second objective was to derive the daily rainfall and flow frequency duration curves, PIDFs and flood frequencies that are generally used in watershed and stream restoration projects, low flow monitoring to protect stream aquatic habitats, and to inform the design and management of stormwater, roadway, and drainage infrastructure.

Study Site
The study site is the Turkey Creek watershed (WS78) in the Santee Experimental Forest (Figure 1) originally established by the USDA Forest Service in 1964 and monitored through 1984.Recognizing the importance of reference watershed data in a rapidly changing coastal environment, stream gauging was reestablished in Turkey Creek in 2004 to facilitate a large-scale eco-hydrological monitoring and modeling program [56].The present WS-78 gauging station is approximately 800 m upstream of the original location; it is instrumented with a real-time stream gauge sensor and a rain gauge (http://waterdata.usgs.gov/sc/nwis/uv?site_no=02172035, accessed on 20 February 2024).The gauging station is operated in cooperation with the US Geological Survey (USGS) and the College of Charleston.This paper presents the analysis of data collected from 2005 to 2021.This third-order watershed is located at 33 • 08 ′ N latitude and 79 • 47 ′ W longitude near the small town of Huger, in Berkeley County of South Carolina (Figure 1), and 60 km northeast of the city of Charleston (metropolitan area population of approximately 900,000).It is one of the headwaters of the East Branch of the Cooper River, a major tributary of the Cooper River, which drains to Charleston Harbor.The Turkey Creek watershed is similar to other low-gradient forest watersheds in the southeastern Atlantic coastal plain where rapid urban development has been ongoing since the 1990s.The elevation of the watershed varies from 3.5 m at the outlet to 11.5 m above mean sea level (a.m.s.l.) [48]; Table 1.The sub-tropical climate is characteristic of the coastal plain, with hot and humid summers and moderate winters.The minimum and maximum air temperatures, based on the 1951-2008 period of records at the Santee Experimental Forest, were −8.5 • C and 37.7 • C, respectively, with an average daily temperature of 18.4 • C. Annual precipitation ranged from 830 mm to 1940 mm; and the average precipitation was 1370 mm [27].Winters are generally wet due to low evapotranspiration (ET) demands, and low intensity, longduration rain events.Summers are influenced by high ET demands and short-duration, high-intensity storms; tropical depression storms affect the region during July and October.
Land use within the watershed is comprised of 88% pine forest (mostly regenerated loblolly (Pinus taeda L.) and long leaf pine (Pinus palustris)), 10% wetlands and water, and 2% agricultural lands, roads, and open areas [40].Detailed description of land use and forest cover and types and management are in Morrison [57]; see Appendix A Figure A1 for a land use and land cover map.The watershed was heavily impacted by Hurricane Hugo in September 1989, and the forest overstory trees were almost completely destroyed [47].The current forests are a mixture of remnant large trees, natural regeneration, and about 1000 ha of planted pine [58] (Table 1).The forests are managed using prescribed fire and thinning to reduce fuel hazards and minimize wildfire risks as well as to manage longleaf pine restoration efforts and wildlife habitat (such as red-cockaded woodpecker, Picoides borealis, an endangered species).Treatments to reduce hazardous fuel mainly include broadcast burning, which was as much as 93% of the watershed in 2012, as well as thinning, low intensity under-burn, stand clearcut, and site preparation.The areas and their percentages of total watershed area where such treatments were implemented on the watershed for 2005 to 2021 period are provided in Appendix A Table A1 [59].However, the areas harvested for timber, which was minimal being below 3% (2019) of total area (e.g., well below the harvesting level in surrounding area), and hydrologic soil groups and their distributions, are provided elsewhere [57,60].The watershed is also used for recreational purposes such as hunting, fishing, bird watching, hiking, canoeing, biking, historical tours, horse riding, all-terrain vehicle (ATV) use, and agriculture [49].
The geology of the study site is defined by the Floridan/Santee aquifer which is partially confined beneath the watershed, where the top of the Santee Limestone is approxi-mately 20 m below ground surface (bgs) on the western side of the watershed and about 13 m bgs in the eastern area (Figure 1).More details are given in Amatya et al. [49].The soils in the watershed are mostly poorly drained soils of Wahee (clayey, mixed, thermic Aeric Ochraquults) and Lenoir (clayey, mixed, thermic Aeric Paleaquults) series [61] (Table 1).The watershed also contains small areas of somewhat poorly and moderately well-drained sandy and loamy soils such as Goldsboro and Lynchburg (both fine loamy sand) as shown in the soils map in Appendix A Figure A2.Details of soil types are given elsewhere [49,57].The watershed contains about 56% (2934 ha) wetlands on planted pine (based on the NWI classification https://www.fws.gov/wetlands/,accessed on 17 March 2021), 26% (1392 ha) of other forested wetlands, and <1% (39.7 ha) of isolated wetlands.

Data Collection
A list of recorded hydrometeorological variables used in the study analysis is shown in Table 2.

Rainfall
There are two automatic tipping bucket rain gauges in the study watershed.One gauge is associated with the Campbell Scientific CR10X weather station located near the middle of the Turkey Creek (TC) watershed; another rain gauge is associated with USGS stream gage 02172035; (https://waterdata.usgs.gov/monitoring-location/02172035/,accessed on 20 February 2024) at the watershed outlet (Figure 1).An additional gauge is on Lotti Road just outside the eastern boundary of the watershed.In this study (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), daily rainfall data from only the TC gauge were used for the analysis (Table 2).Since data from the USGS gauge were partially unavailable for February, March, October, and November of 2005, data from the Lotti gauge were used for the year 2005 instead of the TC gauge when it was still not in operation.Event times and tipping percentages from the Lotti gauge were also used to help fill in gaps in the data.Because the manual rain gauges overflowed during the extreme rain event of 3-4 October 2015, the TC manual/tipping bucket ratio from 28 September 2015 to 2 October 2015 was used to estimate rainfall totals from 2-8 October 2015.After the freezing rain and snow event of 3 January 2018, TC meltwater totals and data from a Pluvio weighing bucket gauge at the distant SEF location [62] were used to estimate precipitation during the 3-8 January 2018 period.Data from each of the automatic rain gauges were verified and calibrated using an adjacent manual rain gauge.Breakpoint event rainfall data from the loggers were processed to obtain daily, monthly, and annual totals for each of the three gauges [62,63].

Streamflow
Streamflow records for the 2005 to 2021 period (Table 2) were obtained from a USGS gauging station with real time stage measurements at 15 min intervals using a Sutron Satlink 2 datalogger connected to a bubbler-type water level sensor anchored on the stream bottom at the gauging station, which is located at the outlet of the watershed near the bridge on Hwy 41 (https://waterdata.usgs.gov/monitoring-location/02172035/,accessed on 20 February 2024) (Figure 1).Flow rates were calculated by the USGS using a stagedischarge relationship.The 15 min flow rates were integrated to obtain the daily, monthly, and annual totals [62,63].Missing daily data in early 2005 were extrapolated using a calibration relationship with nearby watershed (WS80) [30].

Weather Data and Potential Evapotranspiration (PET)
Weather data from two automated stations were available to estimate daily PET for the 2005-2021 study period.The first station (TC), located in the middle of the study site (Figure 1) (Table 2), was installed in October 2005 to measure precipitation, air temperature, soil temperature, relative humidity, wind speed and direction, and solar radiation on a half-hourly basis.The second station on a 27 m tall tower also measures air temperature, relative humidity, net radiation, wind speed and direction on a 15 min basis above a forest canopy on nearby WS80 watershed.Details of the stations are provided elsewhere [62][63][64].

Groundwater
Groundwater data were measured as the depth below the ground surface at the well location, with ground elevations rectified using NAVD88 datum.It was measured using water table wells, shallow piezometers (2.5 m depth below ground), and also deeper piezometers (up to 20 m depth) placed into the top of the Gordon Aquifer (Santee Limestone), an important groundwater resource in the region [51]; see Table 2. Time series data for groundwater dynamics at four sites in the watershed were collected on an hourly basis starting in June 2006 on areas representing four major categories of soils (Rains, Goldsboro, Lenoir, and Lynchburg series) in the uplands on the left (south) bank of the watershed and in a fifth well in an area with Wahee soil (on the right bank) starting in January 2011 (Figure 1) [62,63].Automatically recorded data collection at all sites was discontinued in May 2019.
Water elevation data from shallow water table wells for the 2006 to 2021 period were compared to the piezometric data to estimate groundwater recharge in the unconfined aquifer using the water table fluctuation (WTF) method and an estimate of drainable porosity [51,65], and the Darcy flux method [66,67] during static conditions in the absence of storm events.The WTF method is a short-term estimate that assumes that soil water deficits are satisfied by the storm event and the infiltrating water has a direct impact on the water table.The Darcy flux method is a longer-term estimate for a catchment-scale approximation of groundwater recharge potential.Drainable porosity (i.e., specific yield) of the soil was estimated using water table response data for 13 selected rainfall events from July 2006 to December 2007.Water table and piezometer water level data were used from one of the Lenoir soils sites near the TC weather station near the headwaters.Descriptions of well and piezometer installation, the water level monitoring procedure, and recharge estimates are provided in Callahan et al. [51].

Data Analysis Rainfall, Streamflow, Evapotranspiration, and Weather
Measured data on rainfall and runoff coefficient (ROC), computed as a unitless number by dividing annual runoff by rainfall, and estimated PET and ET for the 2005-2021 period were analyzed for their mean annual (calendar year), standard deviation, and coefficient of variation.Similarly, annual ET was approximated by calculating the difference between the annual rainfall and runoff, though this implies the assumption of a negligible change in soil moisture storage.Weather data collected on a half-hourly basis were processed to obtain daily averages.Weather data for 2005 were obtained from the SEF Headquarters [68] and data from the study site were used for 2006 to 2021.Daily average weather data were then used to estimate daily potential evapotranspiration (PET) using the Priestley-Taylor (P-T) method [69] for a standard grass reference [64].Daily PET estimated using maximum and minimum temperature based daily average temperature from 2006-2010 was decreased by 5% to match the daily average temperature obtained by averaging all 30 min readings, consistent with the rest of the 2011-2021 period.However, daily P-T PET for 2015-2021 period estimated using the weather data from the nearby WS80 tower station was used in this study due to an issue with solar radiation used for computing net radiation.The P-T PET-based daily PET were used to obtain total annual PET in this study.All annual analyses were based on a calendar year.Similar analyses were also conducted using water years (April-March), but for 2006-2021 only, when the soils were relatively wetter and near saturation, and potentially with minimal change in storage year-to-year resulting in less variability in annual ET (Appendix A Table A2).Analyses were also carried out by growing (April-October) and dormant (November-March) seasons (Appendix A Table A3).
Daily (24 h) rainfall records were analyzed for their annual distribution including for various storm sizes.Depth-based daily rainfall as well as streamflow were used to analyze rainfall intensity-duration-frequencies as well as flow duration-frequencies, respectively, for the April 2005-December 2021 period using the Weibull probability distribution method [70].Flow duration curves were also analyzed for relatively dry years (2007, 2012), wet years (2015, 2020), and normal years (2008,2013).
Daily streamflow (only from 2006-2021) was also used to quantify partitioned baseflow following the automated recursive digital filtering procedure reported by Arnold et al. [71].The average of Pass1 and Pass2 filtering techniques was used to compute estimated daily base flow.In addition, relationships between daily peak discharges and corresponding rainfall for various durations were explored.Similarly, both the daily water table depths from the five wells and the water levels in the deeper piezometers were also plotted to examine their response to rainfall, general trends, and seasonal variabilities.Additionally, the Pettitt test [72] was employed to detect change, if any, and its significance in trend in piezometric deep water table data.The Pettitt test is based on the U statistic.The null hypothesis assumes a change point detection in water table time series, with an alternative hypothesis lacking the change point.The null hypothesis is rejected if U is higher than U crit .

Trend Analysis
Interannual trends in all the above hydroclimatic variables were estimated using the Sen's slope estimator [73] and tested for significance based on the modified Mann-Kendall (M-K) trend test [73][74][75] under the null hypothesis of no trend (α = 0.05).The M-K test is a non-parametric test that considers the test statistic, S, to have a zero mean and variance estimated using methods by Kendall [74] and Mann [75].Sen [73] developed a non-parametric procedure for estimating the true slope (change per unit time) given a linear trend is present in a time series.Libiseller and Grimvall [76] developed the modified Mann-Kendall trend test for detection of trend in environmental data in the presence of covariates.In this study, we used the modified Mann-Kendall trend test methodology for testing against the null hypothesis of no trends (at 5% significance level) in our hydrometeorological dataset.
In this study, the GEV methodology is preferred over other extreme value-based modeling approaches, such as the Generalized Pareto Distribution for exceedances, to make the analysis comparable with the one adopted in NOAA-Atlas14 [78].The GEV is a threeparameter distribution comprising of the location (µ), scale (σ), and shape (ξ) parameters, and the theoretical cumulative distribution function of a real-valued random variable, X, can be written as [77] The parameter estimation was performed by using the L-moment methodology [79].L-moment-based estimations are more stable and more robust in the presence of extreme values with small sample size [80,81].This is because L-moments do not have sample size-related bounds [79,82].The p-quantile of the GEV distribution was then estimated as, where (1 − p) is the non-exceedance probability.

Log-Pearson Type III (LPIII) Distribution
Design discharge estimation using 17-year streamflow datasets was performed by fitting annual maxima to log-Pearson type III (LPIII) distribution.The LPIII is a statistical method of fitting frequency distribution values for predicting flood at specific sites and is used extensively in hydrologic applications.We used the LPIII distribution because it is highly recommended for use by the U.S. federal agencies as well as by many hydrologists/practitioners for flood frequency analyses [83].This is consistent with the USGS computer program PEAKFQW version 5.2.0, which fits the annual peak-flow data to an LPIII distribution [84].This method has the advantage of extrapolating event data with return periods well beyond the pragmatic flood events.These frequencies were compared with those obtained by the USGS regional regression equations [84] regionalized to a function of drainage area and 24 h, 50 yr PI derived for the same watershed and rain gauge location, respectively.

Annual Rainfall, Streamflow, ROC, ET, and PET
Annual rainfall data are summarized in Table 3 for the 17-year (2005-2021) study period.Rainfall varied from as low as 994 mm in 2007 to as high as 2243 mm in 2015, with an average of 1470 mm and a COV of 0.23.This COV is greater than the 0.18 observed during 13-year historic 1964-1976 period reported by Amatya and Trettin [85], likely due to low and high extremes that were 27% lower and 64% higher, respectively, than the historic average of 1320 mm for the study site as well as the long-term (1946-2008) average of 1370 mm for the nearby the SEF headquarters rain gauge [27].The highest rainfall in 2015 was an indirect result of Hurricane Joaquin (which stayed in the ocean off the coast) leading to extremely high precipitation that resulted in 470 mm in a 2-day (3-4 October) period, which was close to the 440 mm recorded at Charleston airport [86].Eleven out of the seventeen years, particularly from 2013 to 2020, yielded higher than long-term average rainfall.The estimated annual Priestley-Taylor (P-T)-based PET ranged from 1116 mm in 2009 to as high as 1350 mm in 2019, with an average of 1255 mm (Table 3), similar to the water year-based PET of 1254 mm which also had a COV of 0.05 (Appendix A Table A2).This average was about 10% higher than the average of 1137 mm for the 1946-2008 period obtained by Dai et al. [27] and about 5% higher than the 1194 mm average for the 2011-2021 period, both calculated using the Hargreaves-Samani PET method.In addition, the P-T average of 1255 mm was 16% higher than 1079 mm average obtained by Amatya and Trettin [85] using the Thornthwaite PET method for the 1964-1976 period using data from the SEF weather station.One possible reason for the higher P-T-based PET may also be due to omission of daily soil heat flux in the energy component and inherent differences in the PET methods compared here [63,87].It may also be explained by differences in the in recording periods.
Annual streamflow varied from as low as 23 mm in 2012 to as high as 1517 mm in the extremely wet year of 2015 (Table 3; Figure 2).Note that the relatively dry year of 2012, preceded by an even drier year in 2011, resulted in the lowest flow (23 mm) and lowest ROC (0.02) for the entire period of record.Interestingly, 2007 had a similar dry period yet yielded a higher flow and ROC of 94 mm and 0.09, respectively, compared to 2012 despite 2007 having the lowest rainfall on record.We presume this was due to wetter antecedent conditions in 2006 compared to 2011.In addition, the potential energy demands of 2599 mm PET in 2011-2012 were also higher than the 2510 mm recorded for 2006-2007.This potentially resulted in a higher total 2 year ET of 2080 mm with increased storage compared to only 1903 mm for 2006-2007, with a wetter initial year of 2006 compared to the initial year of 2011, yielding a 133 mm lower flow in 2011-2012 (Table 3).Notably, the water year-based ET of 1023 mm in 2007 (Appendix A Table A2) was 123 mm higher than the calendar year ET of 900 mm (Table 3).The same is true for the water year-based ET estimation for 2011 and 2012 both of which have higher ET than the corresponding calendar year values (Table 3).On an annual basis, the water year-basis ET may be more realistic because it minimizes the storage effects due to ET when its demand is the lowest in April.initial year of 2011, yielding a 133 mm lower flow in 2011-2012 (Table 3).Notably, the water year-based ET of 1023 mm in 2007 (Appendix A Table A2) was 123 mm higher than the calendar year ET of 900 mm (Table 3).The same is true for the water year-based ET estimation for 2011 and 2012 both of which have higher ET than the corresponding calendar year values (Table 3).On an annual basis, the water year-basis ET may be more realistic because it minimizes the storage effects due to ET when its demand is the lowest in April.Average annual flow response to the rainfall, as shown by the ROC of 0.24, ranged from 0.02 to 0.68 in 2015 due to the extreme rainfall event.This average is consistent with the average from the water year-based value (Appendix A Table A2) of 0.24 reported for the 1964-1976 period with only 1320 mm rainfall [85] as well as for the small 155 ha WS77 watershed, with similar management practices, for the 2011-2019 period but with slightly higher average (1496 mm) rainfall [25].Hazardous fuel treatment for understory removal on about 93% of the watershed (Appendix A Figure A2) in the dry year of 2012 might not have affected flows due to the very dry conditions.This is consistent with the authors of [60] who noted that these systems (whether treated or untreated) may respond similarly in extreme wet and/or dry conditions.A double mass curve analysis using cumulative annual flow versus rainfall (not shown) indicates some rising flow trend starting in 2015 as also observed in Figure 2.
Average annual ET of 1077 mm for this period (Table 3) (similar to the water-year average ET of 1079 mm; Appendix A Table A2) was 9.6% higher than the 983 mm obtained for the relatively drier historic 1964-1976 period [84], indicating that the system is moisture limited [88].This was mainly attributed to a higher but insignificant rising trend of rainfall (Figure 3a; Table 4 shown below), providing increased recharge (soil moisture), as well as significantly increasing temperature and PET (Figure 3b,c).The annual ET/Rain ratios varied from 0.32 in the wettest year of 2015 to as high as 0.98 in 2012, with an average of 0.76 (same for the water-year; Appendix A Table A2).This value indicates that about 76% of annual rainfall, on average, is lost to ET, assuming no storage and loss to groundwater.This also indicates that, despite abundant moisture from the rainfall in 2015, the available energy of 1216 mm PET might have limited additional ET loss in this year as Average annual flow response to the rainfall, as shown by the ROC of 0.24, ranged from 0.02 to 0.68 in 2015 due to the extreme rainfall event.This average is consistent with the average from the water year-based value (Appendix A Table A2) of 0.24 reported for the 1964-1976 period with only 1320 mm rainfall [85] as well as for the small 155 ha WS77 watershed, with similar management practices, for the 2011-2019 period but with slightly higher average (1496 mm) rainfall [25].Hazardous fuel treatment for understory removal on about 93% of the watershed (Appendix A Figure A2) in the dry year of 2012 might not have affected flows due to the very dry conditions.This is consistent with the authors of [60] who noted that these systems (whether treated or untreated) may respond similarly in extreme wet and/or dry conditions.A double mass curve analysis using cumulative annual flow versus rainfall (not shown) indicates some rising flow trend starting in 2015 as also observed in Figure 2.
Average annual ET of 1077 mm for this period (Table 3) (similar to the water-year average ET of 1079 mm; Appendix A Table A2) was 9.6% higher than the 983 mm obtained for the relatively drier historic 1964-1976 period [84], indicating that the system is moisture limited [88].This was mainly attributed to a higher but insignificant rising trend of rainfall (Figure 3a; Table 4 shown below), providing increased recharge (soil moisture), as well as significantly increasing temperature and PET (Figure 3b,c).The annual ET/Rain ratios varied from 0.32 in the wettest year of 2015 to as high as 0.98 in 2012, with an average of 0.76 (same for the water-year; Appendix A Table A2).This value indicates that about 76% of annual rainfall, on average, is lost to ET, assuming no storage and loss to groundwater.This also indicates that, despite abundant moisture from the rainfall in 2015, the available energy of 1216 mm PET might have limited additional ET loss in this year as opposed to the ET loss of 1094 mm in 2012, close to the rainfall amount, with abundant energy of 1295 mm PET (Table 3).The ET/PET ratio varied from 0.74 in the year 2007 with the lowest rainfall to as high as 1.05 (e.g., ET > PET) in 2009, with an average of 0.86 (Table 3).These data show that, given the moisture and energy available to evaporate it, the percentage increase in annual ET loss can likely be higher than the percentage loss as drainage to streamflow.Annual ET, calculated here as the difference between total rainfall (not throug after canopy interception) and streamflow assuming no other losses, can be expect be higher than the grass-reference PET during wet years due to higher seas  Annual ET, calculated here as the difference between total rainfall (not throughfall after canopy interception) and streamflow assuming no other losses, can be expected to be higher than the grass-reference PET during wet years due to higher seasonal evaporative losses from rainfall intercepted by the forest canopy [89,90].In addition, longterm annual ET may be close to the PET of these forested wetlands with relatively abundant moisture [25,32,89,91].These results are important indicators for understanding pre-development hydrology and stormwater management practices on these forest lands.Annual rainfall including the year of 2015 with or without the extreme wet month of October both showed increasing, but not significant (p = 0.11 and p = 0.12, respectively), trends (Figure 3a; Table 4).This is consistent with Mizzell et al. [92] who observed gauges across the State of South Carolina as well as the authors of [27] who also found a slight increase but not a significant change in annual rainfall observed at the Santee Experimental Forest gauge over the period of 1946-2008 [27].Increasing annual trends (Sens's slope) are also exhibited by minimum, average, and maximum temperatures (Figure 3b), with the minimum and average showing a statistically significant trend (p-value < 0.05) but not the maximum temperature (p = 0.92) (Table 4).The increasing trend in average temperature is consistent with [27] who also reported a steady increase in temperature at the SEF station during the 1981-2008 and Mizzell et al. [92] who also reported a steady increase in temperatures across the state of SC since the 1970s.Both the estimated p-T-based annual PET and total ET (rainfall-streamflow) show significant increasing trends, with p = 0.003 and p = 0.04, respectively (Table 4).Annual streamflow reflected the temporal pattern of the rainfall, also yielding a non-significant (p = 0.46) trend for the same 17-year period (Figure 3a).The annual ROC, response of streamflow to rainfall, exhibited an increasing, but also non-significant (p = 0.43), linear trend (Figure 3d similar to the streamflow (Figure 3a).
The effects of a significant increase in temperature (Figure 3b) and increasing (p = 0.03) trend in PET (Figure 3c) were likely offset by slightly increasing rainfall (Figure 3a), resulting in non-significant trends in streamflow (p = 0.46) and ROC (p = 0.43).This is similar to Dai et al. [27]'s inference that an increase in temperature will further reduce streamflow due to an increase in ET demands if the increased demand is not offset by additional precipitation.Although some of these results are also consistent with those reported for a long-term pine plantation study site in North Carolina [93], a longer period of data may be needed for a definitive conclusion about the trends at this site.Nonetheless, this information on recent trends for these hydrological variables may be of interest to forest managers, engineers, and land use planners.
On the water year basis (analyzed for 2006-2021), the growing season's (April-October) ROC of 0.30 was 50% higher than that (0.20) of the dormant season (November-March) despite nearly 71% (1041 mm) of the average total rainfall (1469 mm) being recorded in the growing season (Appendix A Table A3).This was likely due to the growing season average ET of 795 mm, which was 179% higher than the dormant season average (285 mm) with its reduced ET demands.The growing season ET was 73.6% of the average annual ET of 1079 mm, consistent with similar other regional studies [26,90].
The mean monthly ROC of 0.24 is consistent with that of 0.22 obtained for the annual period (Table 1).On the water year basis (analyzed for 2006-2021), the growing season's (April-October) ROC of 0.30 was 50% higher than that (0.20) of the dormant season (November-March) despite nearly 71% (1041 mm) of the average total rainfall (1469 mm) being recorded in the growing season (Appendix A Table A3).This was likely due to the growing season average ET of 795 mm, which was 179% higher than the dormant season average (285 mm) with its reduced ET demands.The growing season ET was 73.6% of the average annual ET of 1079 mm, consistent with similar other regional studies [26,90].

Daily Rainfall Frequency Duration
Daily rainfall frequency duration for the 2005-2021 period (6209 days) is shown in Figure 5a.Daily rainfall exceeded 60 mm and 12.5 mm for 1% (62 days) and 10% of the time, respectively.The data show that 24 h rainfall exceeded 100 mm or more 0.2% of the time (12 days) (Figure 5b), most of which occurred since 2015.Most of these large events were associated with tropical storm events (Figure 5a), with a 24 h total of 271 mm and 471 mm over two days on 3-4 October of 2015 (Hurricane Joaquin), 219 mm on 8 October 2016 (Hurricane Matthew), and 206 mm on 5 September 2019 (Hurricane Dorian).The highest recorded total rain in October of 2015 was consistent with Mizzell et al. [86] who reported that all time precipitation records were set across the coastal plain, with totals ranging from 254 mm to over 660 mm.
The data in Figure 5b show the number of daily precipitation events of various sizes recorded during the study period.The number of storm events is increasing, but not significant, for all rain sizes; the steepest trend was for 26-50 mm size storms, with the largest number (19) in 2015 followed by (17) in 2017.

Daily Rainfall Frequency Duration
Daily rainfall frequency duration for the 2005-2021 period (6209 days) is shown in Figure 5a.Daily rainfall exceeded 60 mm and 12.5 mm for 1% (62 days) and 10% of the time, respectively.The data show that 24 h rainfall exceeded 100 mm or more 0.2% of the time (12 days) (Figure 5b), most of which occurred since 2015.Most of these large events were associated with tropical storm events (Figure 5a), with a 24 h total of 271 mm and 471 mm over two days on 3-4 October of 2015 (Hurricane Joaquin), 219 mm on 8 October 2016 (Hurricane Matthew), and 206 mm on 5 September 2019 (Hurricane Dorian).The highest recorded total rain in October of 2015 was consistent with Mizzell et al. [86] who reported that all time precipitation records were set across the coastal plain, with totals ranging from 254 mm to over 660 mm.
The data in Figure 5b show the number of daily precipitation events of various sizes recorded during the study period.The number of storm events is increasing, but not significant, for all rain sizes; the steepest trend was for 26-50 mm size storms, with the largest number (19) in 2015 followed by (17) in 2017.
However, how these storm event rainfall sizes are translated into storm runoff event dynamics are yet to be examined.For example, using the rainfall-runoff data from this site for the 1964-1976 historic period, La Torres Torre [43] found that runoff-rainfall ratios were greater for wet (winter-spring) periods compared to dry (summer-fall) periods and runoff response was related to antecedent soil moisture conditions for wet and dry conditions.The second part of the runoff response was, however, partially studied by Epps et al. [28] for streamflow and direct runoff as response to storm events response linking them also to position of the pre-event water table position.In addition, recently Amatya et al. [52] linked the direct runoff response to pre-event water table and a soil profile saturation parameter by using data from first-order low-gradient coastal catchments including two nested ones from this study site.However, how these storm event rainfall sizes are translated into storm runoff event dynamics are yet to be examined.For example, using the rainfall-runoff data from this site for the 1964-1976 historic period, La Torres Torre [43] found that runoff-rainfall ratios were greater for wet (winter-spring) periods compared to dry (summer-fall) periods and runoff response was related to antecedent soil moisture conditions for wet and dry conditions.The second part of the runoff response was, however, partially studied by Epps et al. [28] for streamflow and direct runoff as response to storm events response linking them also to position of the pre-event water table position.In addition, recently Amatya et al. [52] linked the direct runoff response to pre-event water table and a soil profile saturation parameter by using data from first-order low-gradient coastal catchments including two nested ones from this study site.

Daily Flow Duration Frequency Curves
Daily flow duration curves are presented in Figure 6 for the 2005-2021 period since regeneration of forest stands 15 years after Hurricane Hugo impacted the area in 1989 [1].The smaller plot on the right of Figure 6 shows the 1% to 100% exceedance probability behavior for stream flow (discharge).The data showed that the daily flow exceeded zero

Daily Flow Duration Frequency Curves
Daily flow duration curves are presented in Figure 6 for the 2005-2021 period since regeneration of forest stands 15 years after Hurricane Hugo impacted the area in 1989 [1].The smaller plot on the right of Figure 6 shows the 1% to 100% exceedance probability behavior for stream flow (discharge).The data showed that the daily flow exceeded zero discharge 75.7% of the time, compared to the value of 79% of the time observed by Amatya and Radecki-Pawlik [95] using 13 years (1964-1976) of historic data from this watershed.However, no significant trend in annual outflows was observed (Figure 3a).The frequency plot also shows that this third-order coastal watershed had no flow for about a quarter of the year, on average.Flows of 10.0 mm day −1 or higher occurred <1.21% of the time, with the highest daily flow of 256.8 mm on 5 October 2015 by 242 mm on 4 October, 205 mm on 6 October, and 116.7 mm on 7 October in 2015, as a result of nearly 500 mm rain during 3-4 October (Figure 5a).Similarly, Hurricane Matthew produced 96.2 mm of flow on 8 October 2016, as a result of a 200 mm storm total.This trend is consistent with other studies reporting potentially higher and more intensive storms and peak flows [27,93].Flows less than 5.0 mm day −1 occurred 95.7% of time, most of which (59%) occurred during winter (December-February) and spring (March-May) periods.Flows larger than 15.0 mm day −1 were observed mostly during the fall (September-November) season.These results have implications in managing infrastructure for flow for flood resilience as well as in planning for ecological restorations [96,97].
the highest daily flow of 256.8 mm on 5 October 2015 by 242 mm on 4 October, 205 mm on 6 October, and 116.7 mm on 7 October in 2015, as a result of nearly 500 mm rain during 3-4 October (Figure 5a).Similarly, Hurricane Matthew produced 96.2 mm of flow on 8 October 2016, as a result of a 200 mm storm total.This trend is consistent with other studies reporting potentially higher and more intensive storms and peak flows [27,93].Flows less than 5.0 mm day −1 occurred 95.7% of time, most of which (59%) occurred during winter (December-February) and spring (March-May) periods.Flows larger than 15.0 mm day −1 were observed mostly during the fall (September-November) season.These results have implications in managing infrastructure for flow for flood resilience as well as in planning for ecological restorations [96,97].An analysis of daily baseflow computed using an automatic digital filtering technique for the 2006-2021 period provided a mean daily baseflow of 0.22 m 3 s −1 , which is about 35% of the measured mean daily runoff of 0.634 m 3 s −1 .Using a baseline separation method on Florida coastal watersheds, La Torre Torres et al. [44] found a baseflow of 23%, on average, of total runoff for 41 storm events analyzed from this study watershed for the 1964-1976 historic period.This average was lower than the 47% reported by Morrison [57] using a different partitioning method for 15 events of the recent 2011-2015 period.These results show a wide variation in base flow estimates depending upon the analyzed period, temporal scale, and method of partitioning.

Daily Flow Rate Frequencies for Dry, Wet, and Normal Years
Daily flow rate duration curves for two relatively dry (2007,2012), two relatively normal (2008,2013), and two relatively wet (2015, 2020) years in the 2006-2021 period are shown in Figure 7.During dry conditions, both curves for 2007 and 2012 are plotted lower compared to curves for normal and wet years.The flow rates for all ranges of frequency are lowest and have small variability compared to the wet and normal years.The small amount of variability in flow during dry years is caused by higher water storage and a lower water table (as shown below under the Water Table section) in forested watersheds with high ET demand and influence of groundwater sources in total streamflow.Moreover, large areas of poorly drained Wahee soil series on the right streambank (Fig. 1) may play a dominant role in recharging stream during dry season.Flow during dry years was An analysis of daily baseflow computed using an automatic digital filtering technique for the 2006-2021 period provided a mean daily baseflow of 0.22 m 3 s −1 , which is about 35% of the measured mean daily runoff of 0.634 m 3 s −1 .Using a baseline separation method on Florida coastal watersheds, La Torre Torres et al. [44] found a baseflow of 23%, on average, of total runoff for 41 storm events analyzed from this study watershed for the 1964-1976 historic period.This average was lower than the 47% reported by Morrison [57] using a different partitioning method for 15 events of the recent 2011-2015 period.These results show a wide variation in base flow estimates depending upon the analyzed period, temporal scale, and method of partitioning.

Daily Flow Rate Frequencies for Dry, Wet, and Normal Years
Daily flow rate duration curves for two relatively dry (2007,2012), two relatively normal (2008,2013), and two relatively wet (2015, 2020) years in the 2006-2021 period are shown in Figure 7.During dry conditions, both curves for 2007 and 2012 are plotted lower compared to curves for normal and wet years.The flow rates for all ranges of frequency are lowest and have small variability compared to the wet and normal years.The small amount of variability in flow during dry years is caused by higher water storage and a lower water table (as shown below under the Water Table section) in forested watersheds with high ET demand and influence of groundwater sources in total streamflow.Moreover, large areas of poorly drained Wahee soil series on the right streambank (Figure 1) may play a dominant role in recharging stream during dry season.Flow during dry years was observed only for about 20 to 60% of analyzed period, which can have an adverse impact on stream ecosystem health potentially due to stressful conditions caused by low oxygen and concentrated chemicals.The wet years yielded the highest low flow rates, compared to normal and dry years, that are important for ecosystems to ensure optimal biodiversity of macroinventerebrates [98,99].The extremely high flow rates during the wet years were caused by extreme rainfall and very shallow water table elevations causing soil saturation (Figure 5a).Curves for normal years, showing the highest variability in flow, lie between the curves for wet and dry conditions.Compared to dry years, wet years resulted in a lower number of zero flow days (about 10 of 20% of analyzed period), with the rest yielding low flow rates, consistent with Figure 6.
of macroinventerebrates ( [98,99]).The extremely high flow rates during the wet years were caused by extreme rainfall and very shallow water table elevations causing soil saturation (Figure 5a).Curves for normal years, showing the highest variability in flow, lie between the curves for wet and dry conditions.Compared to dry years, wet years resulted in a lower number of zero flow days (about 10 of 20% of analyzed period), with the rest yielding low flow rates, consistent with Figure 6.

Peak Discharge versus Rainfall of Various Duration
The data in Figure 8 show fitted lines using exponential relationships for observed peak discharges with their corresponding rainfall amounts accumulated over 24, 48, 72, and 96 h durations prior to the peak discharge for the study period, except for the extremely large event of 3-5 October 2015.An exponential increase in runoff for storm events with the water table near or above the surface was reported in earlier studies for nearby watersheds [30,100].

Peak Discharge versus Rainfall of Various Duration
The data in Figure 8 show fitted lines using exponential relationships for observed peak discharges with their corresponding rainfall amounts accumulated over 24, 48, 72, and 96 h durations prior to the peak discharge for the study period, except for the extremely large event of 3-5 October 2015.An exponential increase in runoff for storm events with the water table near or above the surface was reported in earlier studies for nearby watersheds [30,100].Apparently, the strongest relationship of the peak discharge on this 52 km 2 watershed was found for the 48 h duration with the highest R 2 of 0.86 followed by the 72 h duration (R 2 = 0.79).Coincidentally, the USGS estimated peak discharge of approximately 158 m 3 sec −1 for that excluded October extreme event (https://nwis.waterdata.usgs.gov/nwis/peak?site_no=02172035&agency_cd=USGS&format=html;accessed on 20 February 2024) was the result of a 48 h rainfall amount exceeding 470 mm.There was no relationship with the 24 h (daily) duration, indicating that the time of con- Apparently, the strongest relationship of the peak discharge on this 52 km 2 watershed was found for the 48 h duration with the highest R 2 of 0.86 followed by the 72 h duration (R 2 = 0.79).Coincidentally, the USGS estimated peak discharge of approximately 158 m 3 s −1 for that excluded October extreme event (https://nwis.waterdata.usgs.gov/nwis/peak?site_no=02172035&agency_cd=USGS&format=html; accessed on 20 February 2024) was the result of a 48 h rainfall amount exceeding 470 mm.There was no relationship with the 24 h (daily) duration, indicating that the time of concentration (t c ) (defined as the time of flow travel from the most remote point on the watershed to its outlet) used for identifying rainfall intensity in event hydrologic models [19,53] for computing flood discharges for water infrastructure design on watersheds of these sizes should exceed 24 h or more.This implies an effect of a larger storm cell size yielding a longer t c as the watershed size increases [19].The t c of a near 3 h storm duration was estimated for computing design discharges for a much smaller 160 ha WS80 watershed [11,19] adjacent to this study site.Most large peak discharges likely resulted for events when most watershed had saturated/ponded soil conditions or near surface water table (WT) depths (25-27 October 2008; 4-6 October 2015; 8 October 2016; see the following section) possibly contributing to overland surface or shallow subsurface runoff on these low-gradient poorly drained landscape [52].

Water Table
Time series data of WT dynamics from wells in four dominant soil types are presented together with daily rainfall in Figure 9. Earlier work by Amatya et al. [60] using data through 2016 showed no significant trends for daily water table position (p = 0.84 for Rains; =0.28 for Lenoir; 1.00 for Goldsboro, and 0.48 for Lynchburg) at well locations on the study site.No data analysis was carried out for the Wahee well due to limited data.
The mean WT depth for that period varied from about 46 cm for the wells in Rains and Lenoir soils to 122 cm for the well in the well-drained Goldsboro soil, with WT depths of 85 cm for the well in moderately drained Lynchburg soil.A similar pattern for mean water table depths was obtained for the data extended through 2021, with a WT as deep as 2.5 m for the Goldsboro location to even surface ponding of as much as 46 cm for the Lenoir soil during the extreme event of 3-4 October 2015, when all the other wells also responded with a WT at the surface [94].However, the duration of saturation and ponding depended on the soil type, as shown by Williams and Amatya [101] who referred to soil drainage class and soil taxonomy in their analysis.Any additional rainfall during pre-saturated conditions may result in localized as well as watershed-wide overland runoff [52] with a potential for large peak discharge flooding the area, as observed during the 3-4 October 2015 extreme event discussed above.
In a recent storm-event modeling study, Amatya et al. [52] used data from two catchments within the Turkey Creek watershed, and a third adjacent small watershed (WS80), each <210 ha in area.Modeling results showed that a threshold rainfall amount of 110 mm generated shallow surface runoff in addition to shallow subsurface flow due to saturationexcess conditions.The pre-event water table elevation with a soil saturation coefficient (α) model parameter described the saturated depth necessary to partition the storm event volume into overland surface and subsurface runoff.Variability in event runoff was attributed to seasonal trends in water table elevation fluctuation as regulated by soil moisture and evapotranspiration [28,44,102].During the extreme rainfall events of 3-4 October 2015, and 8 October 2016, water tables across the watershed responded similarly, rising to the surface with sustained ponding in some well locations, likely contributing to surface runoff and peak discharge.Note that the water table response (Figure 8) shows several saturation-excess periods of varying durations, especially for the more poorly drained sites with Lenoir, Lynchburg, and Rains soil series.The extreme storm event of 3-4 October 2015 created saturated or near-saturated conditions for several months at those sites.The mean WT depth for that period varied from about 46 cm for the wells in R and Lenoir soils to 122 cm for the well in the well-drained Goldsboro soil, with WT de of 85 cm for the well in moderately drained Lynchburg soil.A similar pattern for m water table depths was obtained for the data extended through 2021, with a WT as d as 2.5 m for the Goldsboro location to even surface ponding of as much as 46 cm for A comparative analysis of groundwater position data in the shallow well (screened from surface to 3 m deep) and deeper subsurface piezometer (screened at 14.5 m deep) at the same location as the water table well in the Lenoir soil (upper Turkey Creek) (Figure 10) provided an average recharge estimate to the surficial aquifer of 114 ± 60 mm y −1 [51].The main factor influencing recharge estimates was antecedent water table level, which in turn was influenced by landscape position and soil texture.The shallow water table conditions at this site support a large range of natural wetlands and create management challenges across the region [103].Modest changes in the position of the water table can lead to either groundwater flooding and concomitant management challenges for silvicultural activities, or to ecosystem stresses related to dry conditions in wetlands during times of below-normal precipitation or because of groundwater withdrawal. .The vertical separation of the two data series indicates a consistent albeit small downward hydraulic gradient between the water table and deeper groundwater system.

4.9.Precipitation Intensity-Duration-Frequencies and Design Flood Frequencies
Compared to the mean onsite 24 h and 48 h precipitation depths for the 2005-2021 period, the NOAA-Atlas14 values for the 1893-2002 period considerably underestimate the precipitation for the 10-, 25-, 50-, and 100-year return intervals for the same rain-gauge location (Figure 11a,b).As expected, the single on-site precipitation depths exhibit a wider range of uncertainty depicted by the larger confidence intervals mainly due to their relatively brief data record period.We did not use data from some other on-site gauges within its vicinity but highly recommend that this be explored in future.Similar underestimates, by 10-60% for durations ≥2 h, were observed by [19] with errors increasing for longer return intervals when compared to estimates using data from the adjacent WS80 watershed.Based on these results, it would likely be useful to re-evaluate the hydro-geomorphologic risk analysis of culverts located within this study site, that were noted by [20] as being vulnerable to NOAA-based 100 yr 30 min precipitation intensities.
The data in Figure 12 show that the 2-, 5-, and 10-year design flow rates estimated by the regional USGS regression equations for the coastal plain region [105] seem to match with those obtained using the Log Pearson Type III (LPIII) fitted flow rates for the on-site  5. Analysis of Pettitt's test using a critical value depends on the number of observations (77 and 79 values for water table and piezometer water depth in this case) in the timeseries and the assumed significant level.Based on Jaiswal et al. [104], the critical values analyzed for the 0.05 significance level were determined as 459 and 478 for the water table and piezometer data, respectively.The analysis indicated that the break point could have occurred on 9 December 2012 or 1 April 2010.However, the computed values of U for the Pettitt test were lower than the U crit , indicating a lack of trend and lack of break-change points in both the water table and piezometric data.As expected, the single on-site precipitation depths exhibit a wider range of uncertainty depicted by the larger confidence intervals mainly due to their relatively brief data record period.We did not use data from some other on-site gauges within its vicinity but highly recommend that this be explored in future.Similar underestimates, by 10-60% for durations ≥2 h, were observed by [19] with errors increasing for longer return intervals when compared to estimates using data from the adjacent WS80 watershed.Based on these results, it would likely be useful to re-evaluate the hydro-geomorphologic risk analysis of culverts located within this study site, that were noted by [20] as being vulnerable to NOAA-based 100 yr 30 min precipitation intensities. .The vertical separation of the two data series indicates a consistent albeit small dow hydraulic gradient between the water table and deeper groundwater system.

4.9.Precipitation Intensity-Duration-Frequencies and Design Flood Frequencies
Compared to the mean onsite 24 h and 48 h precipitation depths for the 2005 period, the NOAA-Atlas14 values for the 1893-2002 period considerably underes the precipitation for the 10-, 25-, 50-, and 100-year return intervals for the same rainlocation (Figure 11a,b).As expected, the single on-site precipitation depths exhibit a wider range of tainty depicted by the larger confidence intervals mainly due to their relatively brie record period.We did not use data from some other on-site gauges within its vicin highly recommend that this be explored in future.Similar underestimates, by 10-6 durations ≥2 h, were observed by [19] with errors increasing for longer return int when compared to estimates using data from the adjacent WS80 watershed.Bas these results, it would likely be useful to re-evaluate the hydro-geomorphologic risk ysis of culverts located within this study site, that were noted by [20] as being vuln to NOAA-based 100 yr 30 min precipitation intensities.
The data in Figure 12 show that the 2-, 5-, and 10-year design flow rates estima the regional USGS regression equations for the coastal plain region [105] seem to with those obtained using the Log Pearson Type III (LPIII) fitted flow rates for the o observed annual maximum peak discharge data.The data in Figure 12 show that the 2-, 5-, and 10-year design flow rates estimated by the regional USGS regression equations for the coastal plain region [105] seem to match with those obtained using the Log Pearson Type III (LPIII) fitted flow rates for the on-site observed annual maximum peak discharge data.The USGS regression equation derived flow rates for 25-year and longer return intervals seem to be underestimated, especially for the 100-year return interval compared to the on-site data.Like with the precipitation intensities discussed above, the larger uncertainty bands for the on-site LPIII estimates (Figure 12) are likely due to use of the shorter recording period than the USGS methods.
In an earlier study using 13 years (1964-1976) of historic data with just the Pearson Type III model, Amatya and Radecki-Pawlik [95] reported similar design discharges to those reported by the USGS method [84] for return periods of ≤10-years but lower values, with underpredictions by as much as 29%, than the USGS ones for the higher return periods (Figure 12).Results from the recent 2005-2021 period show much larger mean design flow rates than those obtained from the current USGS method updated by Feaster et al. [105], which suggests increasing frequencies and sizes of larger peak discharges (Figures 6 and 8), consistent with related increases in PIs and sizes (Figures 5 and 11).These results emphasize that planners and designers should take precautionary measures while designing climate resilient stormwater and road drainage infrastructure using the published USGS regional regression equations on such watersheds.The USGS regression equation derived flow rates for 25-year and longer return intervals seem to be underestimated, especially for the 100-year return interval compared to the on-site data.Like with the precipitation intensities discussed above, the larger uncertainty bands for the on-site LPIII estimates (Figure 12) are likely due to use of the shorter recording period than the USGS methods.
In an earlier study using 13 years (1964-1976) of historic data with just the Pearson Type III model, Amatya and Radecki-Pawlik [95] reported similar design discharges to those reported by the USGS method [84] for return periods of ≤10-years but lower values, with underpredictions by as much as 29%, than the USGS ones for the higher return periods (Figure 12).Results from the recent 2005-2021 period show much larger mean design flow rates than those obtained from the current USGS method updated by Feaster et al. [105], which suggests increasing frequencies and sizes of larger peak discharges (Figures 6 and 8), consistent with related increases in PIs and sizes (Figures 5 and 11).These results emphasize that planners and designers should take precautionary measures while designing climate resilient stormwater and road drainage infrastructure using the published USGS regional regression equations on such watersheds.

Summary and Conclusions
This article described 17 years' worth (2005-2021) of hydro-meteorological data in the Turkey Creek watershed, a 52 km 2 area within the Francis Marion National Forest in coastal South Carolina.Conditions during the period of study ranged from very dry to extremely wet.Consistent with other studies, annual mean and annual minimum air temperature, observed annual precipitation, streamflow, runoff coefficient, and ET (as the difference between rainfall and streamflow) all showed increases, but with considerable variability and no statistically significant trends.Both the mean annual rainfall (1470 mm) and streamflow (362 mm) were larger than those during the 1964-1976 period with 1320 mm rainfall and 337 mm of streamflow.These findings on mean annual streamflow are similar (391 mm) to those reported for an adjacent smaller watershed (310-400 mm) with similar management and for relatively less forested coastal plain watersheds in Georgia reported by Bosch et al. [106].Although annual rainfall and their daily intensities and frequencies as well as storm events were shown to be increasing with increased flooding potential, the watershed did not have flow approximately 25% of time, similar to the 22% found during the historic period.Either increased variability in the recent summer-fall precipitation and/or management practices may explain that small difference.Some rainfall differences were observed from June to September, months characterized by high intensity, short duration storms as reported earlier by Haley [40] and La Torre Torres [43].These observations may have implications in managing streamflow for ecological restorations in the face of rising temperatures observed in this study.Daily peak discharge was highly correlated to 48 h duration rainfall, suggesting the effects of storm sizes and possibly the storage on flooding in the watershed.In addition, the on-site rainfall intensity-duration-frequencies as well as the peak discharges for the frequencies in the ≥25-year return intervals were substantially higher than the published NOAA data, the USGS regional equations, and the historic data-based flood frequencies for this location.One reason is both the published NOAA PIDFs [78] and the USGS design discharge data [84], which uses the NOAA PIDFs, extend only up to the year 2013 in contrast with this on-site data from the year for 2005-2021.These frequencies, developed using <30-years of on-site data, will be useful for stormwater infrastructure design, post-event assessments, and mitigation of environmental impacts [107]; however, the results should be cautiously interpreted as minimum conditions.
Annual ET was a significant component of the water budget; however, the assumed negligible change in storage was contradicted by the on-site water table data.For example, in 2011, the water table depth (WTD) was 45 cm at the beginning of January and ended at 166 cm at the end of December in the Rains soil site (Figure 8), indicating a large soil water deficit.Thus, the actual ET in 2011 was likely more than the estimated ET of 986 mm (Table 3) but less than the PET of 1304 mm.In 2012, WTD ended in December at 69 cm (Figure 8), indicating a large gain in soil water storage, in which case the estimated ET of 1094 mm may also be uncertain.In other words, the assumption of negligible change in storage on a year-to-year basis, particularly at the calendar year interval, may result in large uncertainties in ET estimates (Table 1; Appendix A Figure A2).The potential implications of this are understanding vegetation productivity, water availability and demand, local hydrometeorology, closure of the surface energy budget [108], and several other roles that ET plays [109].The recent installations of eddy covariance-based flux towers accompanied by new soil moisture sensors on the adjacent paired watershed at Santee Experimental Forest should provide more accurate ET estimates including its components of transpiration and understory evaporation as well as sub-annual ET estimates for these coastal forests.

Future Directions
Studies on precipitation extremes (very wet and very dry) in natural ecosystems of similar watersheds are limited, but there is a pressing need to identify (1) how these events affect ecosystem processes and functions for sustainable habitat management and restoration, and (2) how scientists can improve predictions.For example, increased precipitation giving rise to elevated water table and increased soil moisture may impact soil water processes while encouraging invasive plant infestations [110], thus reducing biodiversity and limiting recreational access.Extreme dry conditions may also result in plant infestations [111], reduce base flows into groundwater-supported ecosystems such as aquatic habitats [112], and increase the frequency and severity of wildfire [113].Hydrology data (particularly, event runoff ratios, runoff components, and evapotranspiration) can be used as "reference" information for the "pre-development" scenario design/analysis on watersheds with similar characteristics.Ongoing studies on carbon and greenhouse gas emission rates in tidally mediated freshwater forested riparian wetlands [55] are focused on whole-scale ecosystem shifts associated with sea level rise and climate change.
Resource issues are becoming more complex in forested watersheds like those in the coastal plain region of the southeastern U.S that provide a multitude of ecosystem services.However, these forests are increasingly facing various types of disturbances, with the potential to affect hydrology and water quality of streams draining these systems.For example, Francis Marion National Forest acts as a buffer to urbanization from the surrounding Charleston metropolitan area [50].Urban development is intensifying and advancing towards the Francis Marion National Forest boundary as shown in Appendix A Figure A3.However, within the national forest boundary, the land use is predominately forest, with insignificant changes in land use type [114] based on an analysis of the National Land Cover Database from 2001-2016 [115].Additionally, the forest canopy within Francis Marion National Forest, and particularly the Turkey Creek Watershed, is highly stable compared to the surrounding unprotected forested area as shown in Appendix A Figure A4 [116].Hydrological studies of the Turkey Creek watershed, therefore, provide critical "pre-development" reference conditions for the increasingly urbanizing and dynamic surrounding landscape.
Over the past decade, increased rainfall has led to water surplus conditions, reflected in the shallow water table conditions and longer duration of saturation-excess runoff conditions.This leads to longer hydroperiods in depression and floodplain wetlands, supporting those ecosystems which had previously endured multiple droughts during the 2000-2015 period.However, saturated soils create challenges, such as disrupted timber and road management and erosion and increased runoff to receiving waterways [19].Climate change impacts in the region, such as sea level rise and related saltwater intrusion (both in surface water and groundwater) are affecting the overall drainage of watersheds and placing stress on drinking water resources [117,118].The Gordon Aquifer (hydrogeologically equivalent to the Floridan Aquifer) is locally confined by a low permeability unit in the form of the Parkers Ferry Formation in the western portion of the Turkey Creek watershed; this unit is absent beneath the eastern half of the watershed.Previous research showed an average of about 50-175 mm per year of recharge from analysis of groundwater level data in the eastern section [103].Forested watersheds such as Turkey Creek serve as an important water reservoir and provide services to groundwater-dependent ecosystems [119].
Relatively long-term studies inform decisions in rapidly urbanizing areas affected by similar environmental conditions.Multi-site studies may be warranted to accurately understand these complex forest systems in the face of changing land cover/use and climate.Multi-site studies can clarify the linkage between meteorological conditions and localized extreme storm events driven by atmospheric rivers, and extratropical cyclones, which are common phenomena in this region [120].Adequate resources are critical for securing high quality long-term data from successful multi-purpose monitoring like this one for proper management of land and water in an integrated, sustainable way [15,21].Because long-term monitoring is often resource restricted, new data-driven models of varying complexities and scales could be developed and tested [121] and existing ecohydrological models like SWAT [41,122] and others could be improved with watershed data to understand processes and evaluate impacts of various management treatments, including large-scale forest restoration (i.e., longleaf pine [123]).

Figure 1 .
Figure 1.Location of the Turkey Creek watershed in the Francis Marion National Forest in the lower coastal plain of South Carolina.The larger outline is the watershed boundary (yellow) with streams (blue) delineated on true-color satellite image.Two nested catchments (in green and maroon lines) are also shown in Morrison [57].Groundwater wells on four major soil series (Gb-Goldsboro; Le-Lenoir; Ly-Lynchburg; and Ra-Rains) are also shown together with stream monitoring stations.The hatched circle marks the Turkey Creek watershed outlet, USGS stream gaging station 02172035 including a rain gauge.TC is the weather station in the middle of the watershed.

Figure 1 .
Figure 1.Location of the Turkey Creek watershed in the Francis Marion National Forest in the lower coastal plain of South Carolina.The larger outline is the watershed boundary (yellow) with streams (blue) delineated on true-color satellite image.Two nested catchments (in green and maroon lines) are also shown in Morrison [57].Groundwater wells on four major soil series (Gb-Goldsboro; Le-Lenoir; Ly-Lynchburg; and Ra-Rains) are also shown together with stream monitoring stations.The hatched circle marks the Turkey Creek watershed outlet, USGS stream gaging station 02172035 including a rain gauge.TC is the weather station in the middle of the watershed.

Figure 2 .
Figure 2. Annual rainfall, streamflow, and runoff coefficients for the 2005-2021 period.Broken horizontal lines are the average for the rain (1470 mm) and for the streamflow (362 mm).Long-term (1946-2008) average is 1370 mm.Sloped lines are general trends.

Figure 2 .
Figure 2. Annual rainfall, streamflow, and runoff coefficients for the 2005-2021 period.Broken horizontal lines are the average for the rain (1470 mm) and for the streamflow (362 mm).Long-term (1946-2008) average is 1370 mm.Sloped lines are general trends.

Figure 3 .
Figure 3. Annual trends in (a) rainfall, and streamflow, (b) minimum, average, and maximum perature, (c) potential evapotranspiration, and total evapotranspiration, and (d) ratio of flow rain (ROC) for 2006-2021.The p-values estimated for non-parametric M-K test are also pres for each variable.(d) The scatter points indicate the yearly values and the lines indicate linear t in each of the plots.

Figure 3 .
Figure 3. Annual trends in (a) rainfall, and streamflow, (b) minimum, average, and maximum temperature, (c) potential evapotranspiration, and total evapotranspiration, and (d) ratio of flow and rain (ROC) for 2006-2021.The p-values estimated for non-parametric M-K test are also presented for each variable.(d) The scatter points indicate the yearly values and the lines indicate linear trends in each of the plots.

4. 2 .
Trends in Annual Rainfall, Temperature, Streamflow, Runoff Coefficient (ROC), PET, and ET for 2005 to 2021 period.The plot clearly indicates that the study period yielded more wet months, exceeding 200 mm rain, starting in 2014 (21 months) compared to the previous 9-year period (15 months), with extended no flow months, particularly in 2006-2007 and 2011-2012.Those four years yielded very low annual ROC values as discussed above.The wet year of 2015 had 726 mm of flow during October due to an indirect effect of Hurricane Joaquin

Figure 4 .
Figure 4. Monthly rainfall (thin light blue bars) and streamflow (thick dark blue line) for 2006-2021.Horizontal lines are the average for the rain (122.7 mm) (broken brown) and the flow (solid yellow) 31.1 mm).

Figure 4 .
Figure 4. Monthly rainfall (thin light blue bars) and streamflow (thick dark blue line) for 2006-2021.Horizontal lines are the average for the rain (122.7 mm) (broken brown) and the flow (solid yellow) 31.1 mm).

Figure 5 .
Figure 5. (a) Daily rainfall frequency-duration curves and (b) annual trends in the number of storm sizes for three different categories of rainfall: rainfall storm totals > 100 mm, storm totals between 51-100 mm, and storm totals between 26-50 mm.

Figure 5 .
Figure 5. (a) Daily rainfall frequency-duration curves and (b) annual trends in the number of storm sizes for three different categories of rainfall: rainfall storm totals > 100 mm, storm totals between 51-100 mm, and storm totals between 26-50 mm.

Figure 6 .
Figure 6.Daily flow frequency duration curves.A closer view of the lower exceedance values is shown at right.

Figure 6 .
Figure 6.Daily flow frequency duration curves.A closer view of the lower exceedance values is shown at right.

Figure 7 .
Figure 7. Daily flow rates frequency duration curves for dry years of 2007 and 2012, normal years of 2008 and 2013, and wet years of 2015 and 2020.

Figure 7 .
Figure 7. Daily flow rates frequency duration curves for dry years of 2007 and 2012, normal years of 2008 and 2013, and wet years of 2015 and 2020.

Hydrology 2024 , 34 Figure 8 .
Figure 8. Relationship of peak discharge with rainfall amounts for 24, 48, 72-, and 96 h durations observed during the 2005-2021 period without the estimated discharge of extreme event on 5 October 2015.

Figure 8 .
Figure 8. Relationship of peak discharge with rainfall amounts for 24, 48, 72-, and 96 h durations observed during the 2005-2021 period without the estimated discharge of extreme event on 5 October 2015.

Hydrology 2024 , 34 Figure 10 .
Figure 10.Groundwater head data (2006-2020) for the upper Turkey Creek location.Water table and piezometric head for a 14.5-m deep piezometer are shown.Lenoir soil series (clayey Aeric Paleaquults).The vertical separation of the two data series indicates a consistent albeit small downward hydraulic gradient between the water table and deeper groundwater system.

Figure 11 .
Figure 11.NOAA-Atlas14 and on-site-based (a) 24 h precipitation depths and (b) 48 h precipitation depths corresponding to specific return intervals.The onsite precipitation depths are calculated using the GEV analysis.The circles indicate the best estimates, and the error bars indicate the 5-95% confidence intervals of precipitation depth.

Figure 10 .
Figure 10.Groundwater head data (2006-2020) for the upper Turkey Creek location.Water table and piezometric head for a 14.5-m deep piezometer are shown.Lenoir soil series (clayey Aeric Paleaquults).The vertical separation of the two data series indicates a consistent albeit small downward hydraulic gradient between the water table and deeper groundwater system.Plots for both the water table and piezometric data inFigure 10 indicate a possible change in the trend or a break point somewhere between 2010 and 2012.Results of Pettitt's test for change detection for the significance level α of 0.05.are presented in Table5.

4. 9 .
Precipitation Intensity-Duration-Frequencies and Design Flood Frequencies Compared to the mean onsite 24 h and 48 h precipitation depths for the 2005-2021 period, the NOAA-Atlas14 values for the 1893-2002 period considerably underestimate the precipitation for the 10-, 25-, 50-, and 100-year return intervals for the same rain-gauge location (Figure 11a,b).

Figure 10 .
Figure 10.Groundwater head data (2006-2020) for the upper Turkey Creek location.Water ta piezometric head for a 14.5-m deep piezometer are shown.Lenoir soil series (clayey Aeric quults).The vertical separation of the two data series indicates a consistent albeit small dow hydraulic gradient between the water table and deeper groundwater system.

Figure 11 .
Figure 11.NOAA-Atlas14 and on-site-based (a) 24 h precipitation depths and (b) 48 h precip depths corresponding to specific return intervals.The onsite precipitation depths are calcula ing the GEV analysis.The circles indicate the best estimates, and the error bars indicate the confidence intervals of precipitation depth.

Figure 11 .
Figure 11.NOAA-Atlas14 and on-site-based (a) 24 h precipitation depths and (b) 48 h precipitation depths corresponding to specific return intervals.The onsite precipitation depths are calculated using the GEV analysis.The circles indicate the best estimates, and the error bars indicate the 5-95% confidence intervals of precipitation depth.

Figure A3 .
Figure A3.Change in development map for the Charleston Metropolitan Area, Francis Marion National Forest (black boundary), and Turkey Creek watershed (WS-78; yellow boundary) in South Carolina.Based on analysis of the National Land Cover Database 2001 and 2016 land cover layers.[115].

Figure A3 . 34 Figure A4 .
Figure A3.Change in development map for the Charleston Metropolitan Area, Francis Marion National Forest (black boundary), and Turkey Creek watershed (WS-78; yellow boundary) in South Carolina.Based on analysis of the National Land Cover Database 2001 and 2016 land cover layers.[115].Hydrology 2024, 11, x FOR PEER REVIEW 29 of 34

Figure A4 .
Figure A4.Change in forest canopy phenology map for Francis Marion National Forest and Turkey Creek watershed (WS-78; yellow boundary) in South Carolina.Map shows a gross measure of the instability of landscape dynamics from 2008-2017 based on change in phenological composition and characteristics [116].(Source: LanDAT Geospatial Data.2019.USDA Forest Service.Available online: https://landat.org/maps;accessed on 4 October, 2021).

Table 2 .
List of recorded hydro-meteorologic variables and their characteristics.

Table 4 .
Results of Modified MK Test.The results showing statistically significant (p-value < 0.05) trends are indicated in italics.

Table 5 .
Results of Pettitt test for water table and piezometer water depth.