ABSTRACT

The application of Global Climate Model (GCM) output to a hydrologic model allows for comparisons between simulated recent and future conditions and provides insight into the dynamics of hydrology as it may be affected by climate change. A previously developed numerical model of the Suwannee River Basin, Florida, USA, was modified and calibrated to represent transient conditions. A simulation of recent conditions was developed for the 372-month period 1970-2000 and was compared with a simulation of future conditions for a similar-length period 2039-2069, which uses downscaled GCM data. The MODFLOW groundwater-simulation code was used in both of these simulations, and two different MODFLOW boundary condition “pack- ages” (River and Streamflow-Routing Packages) were used to represent interactions between surface-water and groundwater features. The hydrologic fluxes between the atmosphere and landscape for the simulation of future conditions were developed from dynamically downscaled precipitation and evapotranspiration (ET) data generated by the Community Climate System Model (CCSM). The downscaled precipitation data were interpolated for the Suwannee River model grid, and the downscaled ET data were used to develop potential ET and were interpolated to the grid. The future period has higher simulated rainfall (10.8 percent) and ET (4.5 percent) than the recent period. The higher future rainfall causes simulated groundwater levels to rise in areas where they are deep and have little ET in either the recent or future case. However, in areas where groundwater levels were originally near the surface, the greater future ET causes groundwater levels to become lower despite the higher projected rainfall. The general implication is that unsaturated zone depth could be more spatially uniform in the future and vegetation that requires a range of conditions (substantially wetter or drier than average) could be detrimentally affected. This vegetation would include wetland species, especially in areas inland from the coast.

Keywords:

Groundwater, Climate Model, River System, Downscaling

1. Introduction

Global Climate Models (GCMs) are important tools for simulating historical climate and projecting future climate, including precipitation [1] . The precipitation predicted by some GCMs has been “downscaled”; that is, converted to a finer resolution. The downscaling method can be statistical or dynamic, and both methods have varying de­grees of uncertainty. Comparisons of several models over the conterminous United States indicate that root-mean square errors of precipitation predictions differ by less than 0.1 mm/day between statistical and dynamic methods [2] . Analyses of uncertainty in statistical downscaling methods indicate significant variations between stochastic and regression-based techniques [3] . Dynamic downscaling involves embedding a smaller-scale regional climate model within the GCM [4] . This approach resolves atmospheric processes on a smaller scale and with physically consistent processes, but is computationally intensive and sensitive to uncertainties in the GCM-derived boundary conditions.

The Florida State University Center for Ocean-Atmospheric Prediction Studies (COAPS) has dynamically downscaled GCM-simulated precipitation for the southeastern United States [5] . The downscaled precipitation rates are bias-corrected by using the quantile-matching approach [6] , which assigns corrections by percentile of the precipitation’s cumulative distribution function. COAPS has downscaled other GCM variables including evapotranspiration (ET). Precipitation and ET are both essential to hydrologic simulations, and a predictive model can use the downscaled GCM values to represent future conditions [7] [8] [9] . The downscaled data can be used to compare recent and future simulations to estimate changes in hydrology attributable to climate change.

The Suwannee River Basin, occupying about 9950 mi² in north-central Florida and southern Georgia (Figure 1), is an area where the effects of future precipitation changes are of concern. Forested and agricultural lands account for much of the current land use. Parts of the basin are subject to future population growth and increases in groundwater withdrawals [10] . These factors increase concerns about the effects of precipitation changes and water availability in the future. The Suwannee River, originating in the Okefenokee Swamp in southeastern Georgia (Figure 1) and terminating at the Gulf of Mexico, 12 mi northwest of Cedar Key, is a dominant surface-water feature in the basin (Figure 1). The Alapaha and northern Withlacoochee Rivers are major tributaries to the Suwannee River. The other major tributary is the Santa Fe River, which flows westward from its headwaters to join the Suwannee near Branford, Florida (Figure 1). The lower reaches of these tributaries are incised into the highly transmissive Upper Floridan aquifer, as is most of the Suwannee River (the reach downstream from White Springs), and substantial interactions occur between the aquifer and these river reaches (and the numerous springs occurring along these reaches).

To support analyses of the groundwater system in the Suwannee River Basin, a numerical model was constructed using the MODFLOW code that simulates groundwater flow in the Upper Floridan aquifer and the interaction of the river system with the aquifer [11] . The Suwannee River Basin Model (SRBM) simulates a single aquifer layer under steady-state conditions approximate to September 1990 [11] . Data for calibration came from September 1990 measurements of groundwater levels, discharge at gaging stations for the Suwannee, Alapaha, northern Withlacoochee, Santa Fe, Fenholloway, Aucilla, Econfina, and Steinhatchee Rivers, and spring flows at seven first-magnitude springs [11] (Figure 2). It should be noted that the study area includes two distinct Withlacoochee Rivers that are described herein as “northern Withlacoochee,” referring to the major tributary with headwaters in Georgia, and “southern Withlacoochee,” referring to the river discharging to the Gulf of Mexico approximately 30 mi. southeast of the Suwannee River mouth (Figure 2).

This paper documents the modification and refinement of the SRBM to create a transient version of the model that is suitable for implementing downscaled GCM rainfall and ET data. The calibration of the transient simulations is discussed, along with the adaptation of the downscaled future precipitation to the SRBM, and the repre­sentation of future ET. Because the objective for the modeling was to compare the effects of future climate conditions to relatively recent simulated hydrologic conditions, no attempt was made to predict and simulate future pumping scenarios. Analysis of river and spring flow contributions to the Suwannee Estuary can also be used to predict the effects of changes in precipitation and ET on total freshwater input and salinity.

2. Study Area

The Suwannee, Santa Fe, Alapaha, and northern Withlacoochee Rivers are interconnected and drain the Suwannee River Basin. Outside of the Suwannee River Basin are the Aucilla, Econfina, Fenholloway, and Steinhatchee Rivers to the west and the Waccasassa and southern Withlacoochee Rivers to the south, all of which provide substantial groundwater drainage to the Gulf of Mexico (Figure 2). Large areas of the Suwannee River Basin are devoid of channelized or surface drainage, and most of the drainage directly infiltrates the karst topography of the area, which is generally flat and contains numerous sinkholes. The highly permeable aquifer allows numerous springs, including several major springs and spring groups to augment streamflow (Figure 2). The Upper Floridan aquifer is unconfined in much of the study area, but is generally covered by surficial sediments that confine or partially confine the Upper Floridan aquifer in the northern and eastern parts of the study area.

The climate of the Suwannee River Basin ranges from temperate to humid subtropical. Temperatures typically range from 39˚F to 50˚F in the winter and from 77˚F to 95˚F in the summer. Average annual precipitation in the study area ranges from about 51˚F to 59˚F inches per year (in/yr); about half of the annual precipitation occurs from June through September. Summer precipitation is generally associated with localized thunderstorms that can produce small-scale, intense rainfall. Winter precipitation is generally associated with cold fronts and is more evenly distributed geographically. For the purposes of this study, precipitation is equivalent to rainfall.

3. Methods

The hydrologic modeling for recent and future conditions in the Suwannee River Basin is implemented with MODFLOW-2000 [12] and is based on the original steady-state model by Planert [11] . Modifications to the original model include the addition of parameters for transient simulations and future climatic conditions. Parameter estimation techniques are used to define the spatial and temporal distribution of parameter values required for transient simulation, and potential ET is defined by several methods for the recent and future simulations. The datasets and model files and code are available at Swain [22] .

3.1. Groundwater Model Parameterization

The groundwater model was discretized with 1 layer, 163 rows, and 148 columns and has a uniform grid-spacing of 5000 ft (Figure 3). The Suwannee, Santa Fe, Alapaha, northern Withlacoochee, Fenholloway, Waccasassa, and the southern Withlacoochee Rivers are represented by the MODFLOW River Package RIV1 as in the original steady-state simulation [11] . The upper Fenholloway, Aucilla, Ecofina, and Steinhatchee Rivers and Spring Warrior Creek are represented with the MODFLOW Drain Package (Figure 2 and Figure 3). A total of 696 groundwater withdrawal wells are simulated in the study area. Details of the original steady-state model parameterization can be found in Planert [11] .

3.2. Groundwater and Surface-Water Boundaries

The groundwater head boundaries follow the configuration in Planert [11] . The northern boundary represents the effects of a potentiometric high near Valdosta, Georgia (Figure 4) and is represented by the specified-head package [12] . This representation allows the groundwater head to vary linearly over each stress period of the simulation. A large part of the eastern boundary coincides with the flow lines originating on the Valdosta potentiometric high and is therefore defined as a no-flow boundary. The boundary at the southeast corner has some of the lowest groundwater heads (Figure 2), and is represented as a head-dependent flux boundary by the General-Head Boundary Package [12] . Groundwater head values for the southeastern boundary are developed

by interpolating field-measured groundwater heads, and conductance values are based on local hydraulic conductivities. The head in offshore areas is simulated as a head- dependent flux boundary by assigning RIV cells (Figure 4) with the river stage set at sea level and riverbed conductance values computed with a vertical conductivity value of 100 ft/d and a streambed thickness of 50 ft [11] . The bottom of the aquifer is represented as a no-flow boundary at an assumed depth of 1000 ft below the top of the Upper Floridan aquifer [11] . The hydraulic conductivity is spatially varied to account for the variability in transmissivity, which is due to variability in both hydraulic conductivity and aquifer thickness.

The pumping rates for the 696 groundwater withdrawal wells simulated in the model are based on permit compliance files from the Suwannee River Water Management

District (SRWMD) and the Southwest Florida Water Management District (SWFWMD); from the water-use databases at the St. Johns River Water Management District (SJRWMD) and at the U.S. Geological Survey (USGS) South Atlantic Water Science Center; from water-use files at the Florida Department of Environmental Protection and the Georgia Environmental Protection Division of the Department of Natural Resources; and from agricultural water-use estimates from irrigated area and crop net-irrigation requirement data. Mean monthly values from the 1990-2000 period are used as representative monthly well withdrawals for all simulation years. No attempt was made to modify these withdrawals to represent potential changes in the future or other years.

Stage in the Suwannee, Santa Fe, Alapaha, and northern Withlacoochee Rivers was interpolated from eight stream gages (Figure 4). Stage and river-bottom elevations were estimated for the lower Fenholloway, Waccasassa, and the southern Withlacoochee Rivers by using topographic maps. The difference in bottom and shoreline topography was used to estimate river stage. In order to predict river flows for both the recent and future simulations, a method was developed (described below) that used measured rainfall to generate boundary flows at the Withlacoochee River near Pinetta, Georgia; the Alapaha River at Statenville, Georgia; the Suwannee River at White Springs, Florida; and the Santa Fe River near High Springs, Florida (Figure 4). The net recharge or discharge to the aquifer for the river reaches downstream from the gages on these reaches was added or subtracted to these boundary values to determine river discharge from simulation results.

3.3. Transient Simulations

The SRBM was modified from the original steady-state version [11] to become a transient simulation tool, necessary for the simulation of different rainfall time series. For the recent rainfall simulations, the January 1970-December 2000 (372 months) period is used, and for the future rainfall simulation, the January 2039-December 2069 (372 months) period is used. A 1-month time step was chosen for the calibration of the transient model for the recent period. For comparisons between recent and future periods, the time step was reduced to 1 day. As well as providing better temporal resolution, a shorter time step was chosen for stability concerns when ET is transiently computed in the future simulations, as described below.

For the purposes of determining river stages for the future simulation period, the MODFLOW Streamflow-Routing (SFR2) Package [13] was applied to allow discharge in the rivers to be simulated and stages to be calculated (rather than being specified as input to the model). Only the Suwannee, Santa Fe, Alapaha, and northern Withlacoochee Rivers are represented by SFR2 for two 113-month periods, August 1981-Decem- ber 1990 and August 2050-December 2059. The differences in computed stages are used to estimate river stages in the RIV1 Package for the future January 2039-December 2069 period.

3.4. Precipitation and Evapotranspiration in the Recent Period

Rainfall in the recent period was based on data from 15 rainfall stations (Figure 5). For each station, monthly average rainfall and temperature were computed from the field data. These monthly values were used to compute an initial estimate of net recharge (the difference between rainfall and ET) for each station based on the Thornthwaite method [14] as follows. The potential ET is calculated first using the Hamon equation [15] :

(1)

where PET is potential evapotranspiration in millimeters per month (mm/month); d is the number of days in the month, D is the mean monthly hours of daylight in units of 12 hrs, and W_t is a saturated water vapor density term:

(2)

where T is the monthly mean temperature. In order to develop a more realistic actual ET, at times when the precipitation P is less than computed PET, the actual ET is set to the precipitation plus the available soil moisture (STW):

where ST_i−1 is the soil-moisture storage for the previous month and STC is the soil- moisture storage capacity.

Because the Thornthwaite method net-recharge values do not take into account localized topographic dynamics such as losses due to rapid surficial runoff or limitations in ET due to unsaturated zone depth, the net-recharge values were further adjusted through model experimentation. The discharge values at the Ellaville, Branford, Fort White, and Wilcox stream gages (Figure 4) and groundwater levels at monitoring wells (Figure 6) were used as criteria in a trial and error method to estimate multipliers for the net recharge at each rainfall station zone (Figure 5). A multiplier greater than 1.0 could indicate that ET is limited due to unsaturated zone depth, whereas a multiplier less than 1.0 could indicate losses due to runoff. A reasonably good fit was found with the following recharge multipliers:

1.5 for Tallahassee, Monticello, Madison, High Springs, Ocala, Lake City, Quitman, Jasper, Mayo, and Live Oak;

1.0 for Valdosta and Folkston;

0.5 for Steinhatchee;

0.35 for Usher Tower;

and 0.3 for Perry.

The net recharge is calculated for the rainfall zones surrounding each rainfall station (Figure 5) by multiplying the initial estimate of recharge for the zone by its respective net-recharge multiplier. As expected, the lowest net recharge values occur along the coast (Perry, Steinhatchee, and Usher Tower; Figure 5) where coastal runoff losses would predominate and reduce the water available for net recharge.

3.5. Parameter Estimation

Adjusting the model-input hydraulic conductivity values is functionally identical to adjusting hydraulic transmissivity because the model aquifer thickness is defined as constant. Therefore, hydraulic conductivity is the MODFLOW input for calibration, but field values of hydraulic transmissivity are used for initializing and comparison, and the results are discussed in terms of transmissivity. The SRBM transient groundwater simulations were calibrated to reasonable values of hydraulic conductivity (or transmissivity) and storage properties (specific yield and specific storage) based on the limited aquifer testing available and properties derived from the previous steady-state modeling, using the parameter-estimation techniques implemented in the PEST software suite [16] . PEST accesses MODFLOW input and output to evaluate the sensitivity of computed values to perturbations in model-input parameters [17] . Planert [11] calibrated the steady-state version of the SRBM through trial-and-error comparisons of measured and computed groundwater wells, springs, and river discharge. Through experimentation, and guided by the known hydrostratigraphy, Planert [11] divided the aquifer into zones of differing transmissivity and specific storage. For the purpose of further calibrating the SRBM transient groundwater simulation, the parameter-estima­tion application is based on this same zonation, with seven zones chosen for adjustment by PEST (Figure 6). An additional four zones of varying size were not adjusted in the parameter-estimation application and retain the same aquifer parameters as the steady- state model [11] .

The targets for the parameter estimation were measured aquifer heads at the groundwater wells shown in Figure 6. In order to reasonably limit the number of PEST parameters and accounting for previous calibration efforts, river discharges and spring flows were not used as targets for PEST. Upper and lower bounds for estimating transmissivity were defined differently for each zone (Table 1), but the specific storage estimation is given bounds of 0.000005 to 0.000200 for all zones (Table 2). The estimated transmissivity values range from 22,000 to 10,000,000 ft²/d (Table 1), whereas transmissivity values in the nonparameterized zones are as low as 1000 ft²/d (Figure 7). The range of transmissivity field values from pumping tests [11] is the same order of magnitude as the estimated values (Table 1) with the exception of zone 6 where the estimated value is 257,000 ft²/d and the single aquifer test in the zone yielded 25,000 ft²/d. The estimated specific storage values, which were not needed for the earlier steady-state model, range from 0.000030 to 0.000200 (Table 2).

Comparisons of measured and simulated groundwater levels at the wells labeled in Figure 6 indicate that the closest matches are for wells in the central areas near the

rivers (S021516001, S051311001, S061629001) and the poorest match is for a well near the coast (S111117007) where the simulated values are too high (Figure 8). The actual coastal groundwater exchange to offshore may be larger than the simulated value, assuming that the specified values of net recharge and the computed values of direct ET from groundwater are reasonably accurate.

3.6. Incorporation of Downscaled Global Climate Model Data and Simulation of Future Conditions

To gain insight into possible effects of changes in future rainfall, dynamically downscaled GCM rainfall data were scaled appropriately for input to the SRBM. The GCM used for this application is the Community Climate System Model (CCSM) [18] , which was developed by the University Corporation for Atmospheric Research (UCAR). The coupled components of the CCSM include an atmospheric model (Community Atmosphere Model), a land-surface model (Community Land Model), an ocean model

(Parallel Ocean Program), and a sea ice model (Community Sea Ice Model). The A2 increasing-greenhouse scenario for terrestrial carbon emissions was also used [19] .

3.6.1. Future Rainfall and Evapotranspiration

In order to represent ET for future conditions, the preprocessed values representing present conditions were not considered useful. However, the COAPS effort also downscaled latent heat flux values that can be converted directly to potential ET (PET) by dividing by the latent heat of vaporization of water (2450 kJ/kg at 20^°C). This was done for the SRBM area, and the resulting monthly PET values are used in the MODFLOW Evapotranspiration Package. The ET rate is set to the PET value for water levels 2 ft or less below land surface. Typical ET extinction depths, to which the ET rates linearly decrease to zero, for grass and forest on loamy soils can be 20 ft [20] , which was used due to the substantial percentage of forested land. This value would normally not be as deep in the less arid conditions of this study area, but the high connectivity of the karst aquifer warrants a deep extinction depth.

The rainfall and PET values were defined at the dynamic downscaled points shown in Figure 9. These monthly values were bilinearly interpolated to the Suwannee River model cell values though the scheme:

(3)

where cellval is the value of rainfall or PET at the model cell, val_i is the value of rainfall or PET at location i out of 4 downscaled points around the model cell, and dist_i is the distance from the downscaled point i to the model cell.

The Suwannee River model incorporating the ET computation previously discussed had substantial stability problems when running on a monthly time step. The ET rate is a strong function of the groundwater depth, and the iterative scheme tends to oscillate: a high water level produces a high ET rate, causing a lower computed water level for the next iteration, resulting in a low ET rate, which causes the cycle to repeat and the model to not converge. Large changes in ET volumes are possible over the monthly time step; therefore the time step was reduced to 1 day (while maintaining a uniform recharge rate for all time steps within a given month), allowing for consistent solution convergence. The forcing functions (stress periods) are defined in monthly increments in the future simulation of conditions for the years 2039-2069, just as in the 1970-2000 simulation of recent historical conditions.

3.6.2. Future Surface-Water Inflows

Measured values for the surface-water boundary flows at the Withlacoochee River near Pinetta, Georgia, the Alapaha River at Statenville, Georgia, the Suwannee River at White Springs, Florida, and the Santa Fe River near High Springs, Florida, stream gages (Figure 4) are available for the recent-period simulation (1970-2000), but can also be used to develop relations between rainfall and the measured flows to generate boundary inflows for future periods. A catchment basin is associated with each rainfall station, with one or more rainfall stations with catchments representing inflows to each of the river boundary locations (Figure 10 and Table 3). The monthly rainfall measured at the station is multiplied by the catchment area and a runoff coefficient to produce the boundary flows. The runoff coefficients are adjusted to calibrate the computed flows to

the measured boundary river flow. Mean errors after calibration are 99 ft³/s at Pinetta, 57 ft³/s at Statenville, 109 ft³/s at White Springs, and 426 ft³/s at High Springs (only 20 months of record available at High Springs).

To predict future (2039-2069) runoff flow in these catchments, the time series of average GCM-generated future rainfall in the SRBM area was used. Attempting to use a more localized predicted rainfall time series to drive these river boundaries was considered unnecessary because of the inherent uncertainty in the GCM predictions. The computed river boundary inflows for recent and future simulation periods are shown in Figure 11.

3.6.3. Future Surface-Water Stages

The primary objective for further developing the SRBM under transient conditions was to quantify the difference in streamflow during relatively recent times with those predicted and resulting from future climatic conditions. To contrast with the RIV1 Package used in the original SRBM, which computes groundwater/surface-water exchanges based on the difference between user-defined surface-water levels and model-simulated aquifer heads, the SFR2 Package was applied to help define differences in river stage between recent and future conditions. The SFR2 Package does not require user-defined surface-water levels (it computes them as part of the model solution), but is more complex and prone to stability problems than RIV1. Accordingly, the SRBM is applied for a shorter simulation period (113 months), and the resulting mean differences in stage between recent and future conditions are used to develop the user-defined stages that are required for RIV1 in the full-length future-conditions simulation. The SFR2 Package computes discharge with the steady uniform equations of flow and mass continuity [13] , which can improve understanding of river discharge. The flow at the end of a stream reach is equal to the inflow plus or minus any groundwater leakage. Depth at the beginning and end of reaches is computed with Manning’s equation in the form:

(4)

where y = depth, Q = flow rate, n = Manning’s friction factor, C = Manning’s constant (1.0 in SI units or 1.49 in English units), w = channel width, and S_o is slope of the stream channel.

The periods of simulation for the SFR2 Package are 113 months each―August 1981 through December 1990 for recent conditions and August 2050 through December 2059 for future conditions. This recent period has sufficient measured river flow data to allow a reasonable simulated-to-measured comparison. The SFR2 Package was applied only to the Suwannee, Santa Fe, Withlacoochee, and Alapaha Rivers (Figure 1), and the RIV1 Package was still used in the Fenholloway River, the Waccasassa River, and the lower Withlacoochee River. River-bottom elevations and widths, and surface-water/ groundwater leakage coefficients are the same as used in the RIV1-only simulations.

The model application with SFR2 was compared with field measurements of stage and discharge at the Ellaville, Branford, and Wilcox stream gages on the Suwannee River and the Fort White stream gage on the Santa Fe River (Figure 4). The simulated variations in stage exceeded the upper range of the defined cross-sectional geometry; therefore, for the Suwanee River south of Branford (Figure 4), the cross-sectional

geometry data were extended vertically by 10 ft, which expands the top width in the overbank from 620 ft to 720 ft. However, the cross-sectional geometry in the river remained the same as in RIV1. A better fit was obtained at all stream gages, especially for discharge, although the stage at the Wilcox stream gage on the Suwannee River and the Fort White stream gage on the Santa Fe River remained high (Figure 12). These errors may be due to nonuniform flow dynamics that the SFR2 cannot represent.

The SFR2 Package did not produce particularly accurate fluctuations in simulated stages nor was the full future 2039-2069 simulation period included, so the SFR2 simulated differences in stages between the recent and future simulation were used as a guide to modify the recent-period measured stages to estimate future-period stages for RIV1. The SFR2 simulation was applied to both the recent (1981-1990) and future (2050-2059) periods, and the mean difference in stage between the two simulations was determined at four locations (Figure 4): −0.89 ft at Fort White; 0.64 ft at Ellaville; 0.21 ft at Branford; and −0.29 ft at Wilcox. Sea-level differences are not considered, so zero change in stage at the coastline is defined for the interpolation across the Suwannee, Withlacoochee, and Santa Fe Rivers, producing the change values shown in Figure 13. The surface-water depths in other rivers were not modified from the recent simulation.

4. Results

The recent (1970-2000) simulation period and the future (2039-2069) simulation period model results were compared to evaluate the predicted effects of future rainfall and ET changes. No efforts were made to represent other potential differences between recent and future conditions, such as groundwater withdrawal, sea level, or land cover changes.

4.1. Comparison of Recent and Future Rainfall and Evapotranspiration Rates

Actual ET rates for the future simulation were calculated during the simulation and can be tabulated from the model output. However, the recent simulation used ET computed from the Thornthwaite method summed with rainfall to be input as net recharge. Determining the actual rainfall and ET values in the recent simulation is complicated by the multipliers applied to the net recharge to correct for uncertainties, as discussed previously. When determining the rainfall and actual ET for the recent simulation, it was assumed that the multipliers account for uncertainty in the actual ET and the rainfall does not change, so the actual ET is determined by:

Rearranging:

The measured and downscaled data indicate that the average annual rainfall increases from 55.3 to 61.3 inches (about 11 percent) between the recent and future periods

(Table 4). Simulated average annual ET between these periods increases from 36.2 to 37.9 inches; this 4.5-percent increase is due to higher mean temperatures and latent heat predicted by the CCSM. Although both the actual ET and the rainfall are higher in the future period, the ET as a percentage of rainfall decreased slightly from 65 percent to 62 percent with a 4.4-inch increase in net annual recharge (Table 4).

A higher spatial variation in the average annual net recharge is seen in the future period when compared to that in the recent period (Figure 14(a) and Figure 14(b)), but this can be reasonably attributed to the methods used in each simulation. The net recharge is preprocessed in the recent simulation, but in the future simulation the actual ET is computed during the simulation and must be combined with the preprocessed rainfall to determine net recharge. The recent simulation does not have any negative net recharge values as are seen in the future simulation along the coast and in the far

west (Figure 14(a) and Figure 14(b)). Because the average net recharge is higher in the future simulation (Table 4), the locations with the highest net future recharge areas have values substantially higher than in the recent simulation, mostly in the north and northeast portions of the model domain (Figure 14(c)). This is due to very little ET in these areas even though the average future ET is higher than the recent ET. The effect on groundwater levels is discussed in the next section.

4.2. Comparison of Simulated Recent and Future Groundwater Levels

8

Smaller differences in groundwater levels between recent and future simulations are seen nearer the river system. In both simulations, the groundwater levels generally slope towards the river systems, indicating that the rivers drain the groundwater system. Water levels near most springs (Figure 2) do not differ appreciably between the simulations. The exception is with the westernmost spring near the Aucilla River (Figure 2 and Figure 15), which experiences lower groundwater levels; simulated mean spring flow decreased from 392 ft³/s to 188 ft³/s.

An important trend is observed when the simulated groundwater levels are expressed as a depth below land surface (Figure 16). These groundwater levels are consistently shallower in the future period in areas where they were deepest in the recent period. Conversely, groundwater levels that were shallowest in the recent period become deeper in the future simulation period. Simulation results suggest future groundwater levels will have a greater conformity to land surface (Figure 16); the standard deviation (σ) of water-table depths changed from σ = 36.2 ft in the recent simulation to σ = 28.8 ft

in the future simulation. This trend toward conformity is from a combined effect of higher rates of rainfall and ET in the future simulation (Table 4). Areas where the groundwater levels are relatively lower in the recent simulation received more rainfall recharge in the future simulation, and ET rates are always lower relative to rainfall at

these greater depths to water. Areas where groundwater levels are relatively higher also get more rainfall recharge, but the higher ET rate near the surface offsets the higher rainfall recharge. Although the average increase in ET for the entire model domain is less than the average increase in rainfall, ET is relatively greater in wetland areas near the coast where water levels are much closer to surface (Figure 14(b)). The resulting future net recharge rates are more negative in coastal areas, which are consequently noticeably drier in the future simulation (Figure 16). Conversely, the easternmost areas of the model domain become substantially wetter in the future compared to those in the recent simulations.

4.3. Comparison of Discharge in Recent and Future Simulation Periods

River levels for the entire 372-month future simulation period were estimated based on testing of the SFR2 Package with the 113-month future period. These future changes in river levels were applied, and the total leakage exchanges for all the river reaches were summed to compute net flows at points along the rivers. A comparison of monthly average flow rates for the recent and future periods at the Ellaville, Branford, and Wilcox stream gages on the Suwannee River and the Fort White stream gage on the Santa Fe River (Figure 4) indicates differing patterns of discharge with similar peak magnitudes (Figure 17). The ratio of the standard deviations (σ) between the future and recent periods (σ_future/σ_recent) was determined at each station to quantify discharge variability between the recent and future simulations: Ellaville = 0.924, Branford = 0.934, Wilcox = 0.981, and Fort White = 1.040. Note that the three sites on the Suwannee River all show a reduction in the simulated discharge variability with more reduction in variability

further upstream, indicating upstream boundary effects (Figure 4). Conversely, Fort White showed a slight increase in variability so there is no single model-wide trend. It is likely that the differences between river discharge fluctuations in the recent and future simulations are primarily due to differences in the rainfall-runoff time series used for the upstream river boundaries.

Cumulative simulated discharge for the recent and future periods also was compared (Figure 18). At all sites, the total flow over the entire 31-year period was greater for the future simulation than for the recent simulation, although for the first 18 years of each simulation there is little difference between the periods. This pattern is consistent with the 4.36-inch increase in simulated net recharge (Table 4) and the higher groundwater levels in the western inland areas (Figure 15) near the northern parts of the Suwannee River (Figure 2).

5. Model Limitations

The SRBM provides a useful tool for estimating the effects of climatic changes as represented by GCM output. As is the case with all model applications, parameters and processes are simplified and estimated, resulting in uncertainties. The major limitations of the model application can be summarized as follows:

1) Future rainfall was dynamically downscaled from the CCSM GCM and is subject to all the uncertainties and assumptions in the CCSM and the downscaling process. Predictive modeling relies on extrapolating known processes into the future by using assumptions for unknown quantities. Model results must therefore be viewed as potential, but not necessarily inevitable, outcomes.

2) Future-condition ET estimates rely on downscaled CCSM GCM data and are subject to inherent limitations in accuracy. The rate at which direct ET from groundwater is reduced at increasing depths below land surface is not based on data from the study area, but rather on a simple standard linear function. Simulated groundwater levels are therefore dependent on this ET rate function.

3) Net recharge rates (precipitation minus actual ET) are calculated with different methods in the recent and future simulations. A value of net recharge is precalculated and input to the recent simulation, but ET is computed during the future simulation and varies with groundwater head. This appears to be a factor in the substantially higher spatial variability in simulated future netrecharge compared to the recent simulation (Figure 14) and may exaggerate the changes in groundwater levels between recent and future simulations (Figure 15).

4) Surface-water levels must be specified by the user for the 372-month simulations, causing prediction uncertainty in the future simulation. The SFR2 Package was applied for recent and future 113-month test periods to develop insight into potential differences in river stage. Although the routing of flow is accounted for when computing stage, these stages are subject to uncertainties. The assumption then is made that average differences between the recent and future river stage indicate a shift that can be applied to the user-defined stage without regard to interseasonal variations. The boundary

surface-water inflows are based on rainfall-runoff relations and do not take into account the complex interactions that were simulated within the study area.

5) Simulation of groundwater levels was simplified as one layer with no vertical gradients in hydraulic head. The change in aquifer flow area with changing head over time is neglected because it is assumed that the saturated thickness of the aquifer is constant. Boundary conditions assume that the potentiometric high near Valdosta, Georgia, behaves the same in the recent and future simulation periods, and the eastern boundary continues to coincide with the flow lines and acts as a no-flow boundary. This can be a factor in the large increase in simulated future groundwater heads in the northeastern model area. Similarly, the boundary at the southeast corner of the study area repre- senting drainage to the Ocklawaha River is the same in the recent and future simulation periods. Sea-level change in offshore areas is not accounted for in the future simulation period.

6) Changes in land use and water use due to human activities are not accounted for in the future-conditions simulation. Increased groundwater pumping, for example, can be expected as population increases, but no changes to groundwater withdrawals were simulated in the future-conditions simulation. Changing land use would induce changes in surface runoff and net-recharge rates, which have been shown to have an important effect on groundwater levels.

7) The scope of the present study included a single GCM, the CCSM, to predict future rainfall and ET. Distinct differences in predicted rainfall and ET exist among GCM results because the assumptions needed to predict future climate have inherent uncertainty. Nevertheless, the SRBM simulation procedure described herein can be repeated with downscaled data from other GCMs. Our study was not designed to evaluate the variations between climate models but rather to demonstrate the hydrologic implications of one possible future climate scenario. This study demonstrates how downscaled GCM data can be combined with a groundwater flow model to gain insight into potential effects of climate change.

8) The accuracy of the model results depends on the sensitivity of the model output to uncertainty in the model input parameters. The results of the parameter estimation and the scenario simulations suggest that groundwater heads are notably sensitive to net recharge values; at least for the range of values experienced in the scenario simulations. This indicates that uncertainty in aquifer transmissivity is a smaller source of error compared to factors affecting net recharge, such as the different methods of computing net recharge mentioned in limitation 3 above.

The limitations listed above help guide the interpretation of model results and the relative confidence of the model results. The recharge and boundary variability issues listed in limitations 3 and 5, respectively, indicate that the high simulated groundwater heads in the northeastern model domain are quite uncertain, although the trends in net recharge assure some confidence that a groundwater level increase would be expected. The predicted groundwater differences near the coast have fewer limitations that induce uncertainty, so ultimately the future trend of higher rainfall and higher ET causing the groundwater table depth to become more uniform seems reasonable. Even with the limitations described herein, the model simulations presented in this study can provide useful insight to Suwannee River Basin hydrology under potential future conditions of climate change.

6. Summary

A numerical model of the Suwannee River Basin area was adapted to compare recent hydrology to future hydrologic conditions resulting from predicted changes in rainfall and evapotranspiration (ET). The steady-state model of Planert [11] was modified to represent transient conditions and was calibrated with the PEST parameter-estimation code [16] . The surface-water system was represented using the MODFLOW River Package, which required user-defined river levels.

Some variables for simulation of the future period were developed from dynamically downscaled precipitation and potential ET produced by the Community Climate System Model (CCSM). The precipitation and potential ET data were interpolated for the Suwannee model grid from this downscaled dataset. The conditions simulated in the recent period (1970-2000) represent ET estimated with the Thornthwaite method based on measured data. Downscaled potential ET was applied with a standard function to decrease actual ET at lower groundwater levels under future conditions. This method produced undesirable oscillations at the month-long simulation time step originally used, and consequently the time step was reduced to 1 day.

To help guide the estimation of future river levels, the MODFLOW Streamflow- Routing (SFR2) Package was applied in a shorter recent simulation period (1981-1990) and future simulation period (2050-2059), and average differences in river stage were determined. These differences in mean river stage were then used to create a time series of river stages for the River Package in the 2039-2069 simulation. The recent and future flows were calculated at river locations by adding or subtracting the simulated leakages for all upstream river segments to or from the initial inflows.

Comparisons of simulation results for the recent and future periods provide useful insight into the potential effects of future climate changes. Limitations from the representation of no-flow model boundaries and the spatial distribution of net recharge must be considered when interpreting model results. Even with these limitations, the simulations indicate that increasing rainfall and ET together in the future period caused groundwater levels to rise in areas where they were low in the recent period, and groundwater levels declined in areas where they were higher in the recent period. The general implication is that unsaturated zone depth would be more spatially uniform in the future, and vegetation that requires a range of conditions (substantially wetter or drier than average) would be detrimentally affected. This vegetation would include wetland species, especially in areas near the coast.

Acknowledgements

The authors appreciate the input of Trey Grubbs at the Suwannee River Water Management District in the development of this paper.

Cite this paper

Swain, E. and Davis, J.H. (2016) Applying Downscaled Global Climate Model Data to a Groundwater Model of the Suwannee River Basin, Florida, USA. American Journal of Climate Change, 5, 526-557. http://dx.doi.org/10.4236/ajcc.2016.54037

References