Purpose and Scope

This report documents the development of projected DDF curves (2050–89) for the SFWMD. As part of this study, an ensemble method was used to determine median change factors for precipitation depths as well as inter-model variability at locations throughout central and south Florida. Change factors were determined for durations of 1, 3, and 7 days, and return periods of 5, 10, 25, 50, 100, and 200 years. The 1- and 3-day durations are required for the design and permitting of stormwater systems in the State of Florida. The 7-day duration was included to capture potential future changes in precipitation at the multiday timescale, in particular the possibility of high precipitation caused by stalling storms. Change factors were computed from DDF curves fit to precipitation data from downscaled climate datasets for 40-year periods representing future projected climate for the period 2050–89 (centered on 2070) and historical (retrospective) climate for 1966–2005 for the study area. The start of the historical period coincides with the approximate start of the global warming trend from the 1970s onward (Rahmstorf and others, 2017), and the end was selected to coincide with the end of the Coupled Model Intercomparison Project Phase 5 (CMIP5) historical simulations (2005). The downscaled climate datasets used to derive the change factors as part of this study include the Coordinated Regional Downscaling Experiment (CORDEX), Localized Constructed Analogs (LOCA), Multivariate Adaptive Constructed Analogs (MACA), and JupiterWRF.

Terando and others (2020) recommend that, whenever feasible, the entire range of representative concentration pathway (RCP)-based scenarios be considered when assessing potential future impacts. As part of the FPLOS Program, and with the purposes of planning for future stormwater infrastructure projects with a long design lifetime, the SFWMD is interested in the medium-low (RCP4.5) and high (RCP8.5) emission scenarios that are widely available in downscaled climate datasets. Other climate scenarios such as RCP2.6 are generally not available in downscaled climate datasets and are, therefore, not considered as part of this study. These scenarios are described in more detail in the “Datasets" section.

Perica and others (2013) derived historical DDF curves based on historical observations at meteorological stations in central and south Florida, which are published in National Oceanic and Atmospheric Administration (NOAA) Atlas 14. The DDF curves based on partial-duration series (PDS) from NOAA can be multiplied by the change factors derived in this study to determine potential future extreme precipitation depths for events of a given duration and return period.

Background

The climate of Florida and its relation to remote climate oscillations via atmospheric teleconnections have been studied extensively. In the next sections, factors influencing precipitation extremes in the State as well as projected changes in these extremes, will be discussed.

Projected Changes in Precipitation Extremes: Characteristics and Causes

As discussed by Giorgi (2010), precipitation exhibits a much more complex spatiotemporal variability than temperature, and its changes in response to global warming are highly dependent on how regional circulations change under increased GHG emissions. Regional and local forcings such as land use changes can also affect the precipitation change signal, as demonstrated by Bukovsky and others (2021). For this reason, precipitation projections are much more uncertain than temperature projections at regional and local scales, as evidenced by relatively low consensus in the sign of future precipitation changes across climate models. Projected changes in extreme precipitation tend to be even more uncertain. The response of mean sea levels to global warming is slower than the atmospheric response; therefore, sea-level rise projections over the short-term (on the order of a few decades) tend to have a relatively small range. For example, the Southeast Florida Regional Climate Compact’s (2020) Unified Sea Level Rise Projections for Southeast Florida indicates that local sea-level rise will range between 10–17 inches (in.) by 2040 and 21–54 in. by 2100 compared to year 2000 mean sea level at Key West, Florida. However, in the long-term, sea-level rise projections range between 40 and 136 in. by 2120. The much larger projected range of sea-level rise, especially beyond 2070, stems from uncertainties in future GHG emissions and resulting geophysical effects, most notably the response of the Antarctic and Greenland ice sheets to increased levels of warming.

Seneviratne and others (2012) report that global and regional studies indicate it is likely that the frequency of heavy precipitation or the proportion of total precipitation from heavy precipitation events will increase in the 21st century over many global areas as a result of climate change. They also state that changes in extreme precipitation cannot easily be related to changes in total precipitation because changes can be of opposite sign in some cases. This finding was documented in the Intergovernmental Panel on Climate Change’s (IPCC’s) Fifth Assessment Report (Collins and others, 2013) and further validated in their Sixth Assessment Report (IPCC, 2021, p. 15), which states, “There will be an increasing occurrence of some extreme events unprecedented in the observational record with additional global warming, even at 1.5 degrees Celsius [°C] of global warming. Projected percentage changes in frequency are higher for rarer events (high confidence).”

Using data from observing stations worldwide, Pendergrass and Knutti (2018) show that the median number of the wettest days in a year over which half of annual precipitation falls is 12 days, and that this asymmetric distribution of precipitation will become even more asymmetric because of climate change, with one fifth of the projected increase in rainfall occurring during the wettest 2 days of the year. Collins and others (2013) summarize the science of extreme precipitation events and state that changes in extreme precipitation will be driven by two major mechanisms: (1) a thermodynamic mechanism based on the Clausius-Clapeyron (CC) relation (Wallace and Hobbs, 2006), which specifies that as the air temperature increases, the amount of water vapor in the air at saturation also increases (approximately 7-percent increase per degree Celsius of warming); and (2) a dynamic mechanism that states that extreme precipitation events are controlled by variations in horizontal moisture convergence and associated convection, which would change in a more complicated way as a result of global warming. The conjecture under the thermodynamic mechanism is that, as the climate warms and in the absence of moisture limitation, specific humidity increases according to the CC relation, whereas relative humidity remains constant on a global scale (Koutsoyiannis, 2020). Precipitation intensity is believed to be proportional to surface atmospheric moisture content. In particular, the intensity of precipitation extremes, which tend to occur when the atmosphere is close to saturation, is often proportional to the moisture holding capacity of surface air, which increases with temperature according to the CC relation (Wang and others, 2017). However, Wang and others (2008b) found that air temperature and water vapor trends do not support the assumption of unchanging relative humidity over North America. Koutsoyiannis (2020) analyzed monthly data from two global reanalysis datasets and found a decrease in relative humidity with an increase in temperature. On the basis of global observations, reanalysis data, and climate model output, Wang and others (2017) found a peak-shaped relation between daily temperature and daily precipitation extremes, with precipitation extremes increasing at the CC rate for the low-medium range of local temperatures but decreasing at higher local temperatures. After the peak of the temperature-precipitation relation, specific humidity increases more slowly with temperature and as a result, the relative humidity decreases. Wang and others (2017) discuss how this can be the result of moisture limitations (that is, a dynamic mechanism) and temperature responses to precipitation (as opposed to precipitation responses to temperature).

Martinkova and Kysely (2020) reviewed observational studies on the CC relation scaling for middle latitudes. They discuss how several studies have identified super-CC scaling (that is, scaling above the CC relation) for temperatures above 10–15 °C with some studies showing a 2CC relation above this range of temperatures (that is, approximately 14-percent increase per degree Celsius of warming). This range of temperatures is above normal daily high temperatures during the wet season in south Florida and central Florida when most extreme precipitation events tend to occur because of convection and tropical systems. Therefore, a super-CC scaling relation is likely to dominate in the region. Some studies have also identified a higher temperature threshold above which sub-CC scaling (that is, scaling below the CC relation) or even negative scaling (that is, a decrease in precipitation with increasing temperature) can occur; however, this threshold varies with location. According to Martinkova and Kysely (2020), although sub-CC scaling at very high temperatures can result from a lack of moisture, it may instead be caused by smaller sample sizes resulting in underestimation of high quantiles. In their paper, Martinkova and Kysely (2020) summarize the factors for which super-CC scaling regime has been attributed, which include (1) dynamical feedbacks resulting from additional latent heat being released during condensation and increased near-surface humidity resulting in enhanced convection, (2) microphysical processes, (3) the size of cloud cells and their merging, (4) quasi-stationary convective storms, (5) convective events versus large-scale stratiform precipitation events, among others. More recently, Ali and others (2022) found that using quality-controlled hourly precipitation and daily dewpoint temperature data for the continental United States results in higher median scaling that is closer to the CC rate than when using raw data. They also found higher scaling for higher measurement precision of hourly precipitation. A method that removes seasonality in the precipitation and the dewpoint temperature also gives higher scaling. This highlights the importance of using high-quality observations and robust methods in estimating the temperature-humidity scaling relation.

The thermodynamic mechanism is believed to dominate in high latitudes, whereas the dynamic mechanism is believed to dominate at low latitudes such as in the tropics. Collins and others (2013) also state that projections of relative changes in future extreme precipitation may exceed projections of relative changes in future mean precipitation at the regional scales; however, natural variability in extremes is also larger than in the mean, resulting in decreased signal-to-noise ratios, especially for the most extreme events. Large-scale changes in global and regional circulations (for example, storm tracks) of both natural and anthropogenic origin may also dominate over the thermodynamic and dynamic effects for certain regions and events (Collins and others, 2013). For example, Hall and Kossin (2019) have found evidence of increased hurricane stalling along the North American coast since the mid-20th century but do not attribute these changes to anthropogenic climate forcing, stating that these trends could be due to low-frequency natural variability. Bhatia and others (2019) identified recent increases in tropical cyclone intensification rates in the Atlantic basin over the period 1982–2009 with a positive contribution from anthropogenic forcing.

Dougherty and Rasmussen (2019) evaluated the frequency of various flood types by season during the period 2002–13 across the continental United States. They classified flood events as flash floods, slow-rise floods, and hybrid floods that exhibit characteristics of both flash and slow-rise floods. For the State of Florida as a whole, they found that all three types of flood events can occur at various times of the year, predominantly in the summer and fall. However, a larger number of events that could be classified as slow-rise floods were identified as having occurred in the fall (September–November) and a larger number of hybrid floods were identified as having occurred during the summer (June–August) for the State of Florida. Events classified as flash floods in Florida had a mean precipitation duration of 10–20 hours and durations as long as 2 days in south-central Florida, whereas events classified as slow-rise floods had durations of about 7 days in south Florida and durations as long as 10 days in central Florida.

Kharin and others (2013) found that across future emission scenarios, the multimodel median return value for the 20-year daily precipitation event increases by a global mean of 5.3 percent per degree Celsius of warming. Collins and others (2013) show a multimodel ensemble median percent increase in 20-year return values of daily precipitation for the Florida region ranging from about 2.5 to 5 percent per degree Celsius for 2081–2100 with respect to 1986–2005. Over land areas, Kharin and others (2018) found that the probability of an annual maximum 1-day precipitation extreme with a 20-year return period for the current climate is projected to increase by 17 and 36 percent at 1.5 and 2.0 °C of warming, respectively. Kharin and others (2018) also found that the probability of the 50-year return period daily event in the current climate is projected to increase by 20 and 43 percent over land at 1.5 and 2.0 °C of warming, respectively.

The GCMs that the IPCC reports rely on have spatial resolutions that are generally too coarse to be able to provide projections of climate changes at scales that are relevant for impact studies and planning. In particular, regional and local changes in extremes cannot be adequately captured by GCMs, because they are not only affected by teleconnections to global phenomena but are also highly dependent on local microclimatic conditions. To overcome some of these limitations, GCM output can be downscaled to regional and local scales. Trend analysis of historical observations of relevant meteorological and hydrologic variables also provides important information on the local response to climate change. General increases in extreme precipitation are already being observed in the southeastern United States (DeGaetano, 2009), but trends over Florida are mixed in terms of their sign or significance.

Irizarry-Ortiz and others (2013) evaluated historical changes in precipitation at 32 weather stations over the State of Florida with the longest records ranging from 50 to more than 100 years. They only found upward trends in precipitation maxima in the dry season when considering the entire period of record at these stations, although most of these upward trends were found to be insignificant at the 0.05 significance level. For the period after 1950, most stations, except for those in south Florida, exhibited increases in November–January maxima. They highlighted how the multidecadal variability in precipitation related to changes in the AMO complicates the determination of secular trends and their attribution. Mahjabin and Abdul-Aziz (2020) evaluated precipitation data at 24 Florida stations using the Mann-Kendall nonparametric trend test with prewhitening. They found locally significant (p < 0.1) upward trends in the magnitude of 1- to 12-hour extreme precipitation events for the period 1950–2010 and in 6-hour to 7-day extreme precipitation for the period 1980–2010. Trends in precipitation of 2- to 12-hour duration, specifically 0.26 millimeter per year on average for 2-hour duration increasing to 0.53 millimeter per year on average for 12-hour duration, were found to also be globally significant (that is, trends had field significance, where the number of trends found exceeded the number of trends that could occur by chance) for the period 1950–2010 when accounting for cross-correlations across stations using bootstrap resampling. Mahjabin and Abdul-Aziz (2020) also found globally significant downward trends in the annual number of 1- to 3-hour, 1- to 6-hour, and 3- to 6-hour extreme precipitation events for the periods 1950–2010, 1960–2010, and 1980–2010, respectively. Trends in the number of 1- to 7-day extreme precipitation events were found to be mixed and lacked global significance. The Fourth National Climate Assessment (Carter and others, 2018) found an overall increase in the number of heavy precipitation days (defined as precipitation greater than 3 inches per day) interspersed with interdecadal variability for the southeastern United States for the period 1900–2016. For Florida specifically, the study found a mix of upward and downward trends in the number of heavy precipitation days for the period 1950–2016.

The SFWMD (2021) recently evaluated water and climate resiliency metrics districtwide including long-term trends in hydrologic and meteorological observations. A daily peaks-over-threshold (POT) analysis for occurrences above the 1-in-10 and 1-in-25 return frequency event magnitude was performed for six stations with more than 50 years of available precipitation data. The results showed statistically significant upward trends at the 5-percent significance level at two stations but no significant trends at the other four stations. The SFWMD (2021) also performed regional trend analysis using gridded daily precipitation data from the SFWMD’s precipitation “Super-grid” at 2- × 2-mile (mi) resolution (fig. 2B), grouped by SFWMD ArcHydro Enhanced Database (AHED) rain area (as of November 2020; fig. 3). Contact the SFWMD for more information about the November 2020 AHED rain areas. The results from this analysis showed (1) significant upward trends in daily annual maxima for the Eastern Everglades Agricultural Area, Martin-St. Lucie, and Upper Kissimmee AHED rain areas for a 5-year return period; (2) no significant trends in 3-day, 5-year annual maxima; and (3) a significant downward trend in the 5-day, 5-year annual maxima in the Broward AHED rain area. Results of a POT analysis on AHED rain areas is also included in SFWMD (2021).

Terminology Used in This Report

In describing our processing and analysis, the definitions of several terms are important to clarify. In this report, the term “reanalysis” will be used to refer to dynamical climate model simulations that are based on observed boundary conditions and are meant to match observed weather as precisely as possible. The term “historical observations” will be used to refer to weather observations used in statistical downscaling. The term “historical simulations” will be used to refer to GCM and associated downscaled models that are tuned to preindustrial conditions (around 1850 or 1750 depending on reference; see Schurer and others, 2017 for a discussion) and possibly to more recent historical conditions (Mauritsen and others, 2012; and Hourdin and others, 2017), but are otherwise not constrained to precisely match daily weather conditions from 1850 to the present. The term “future projections” refers to GCM and associated downscaled models of future climate that assume a specific GHG emission scenario or concentration trajectory.

Observational Datasets

Three observational precipitation datasets were used in this study to define base historical conditions for comparing DDF curves fitted to precipitation output from downscaled climate datasets for the historical period and for evaluating the performance of climate models in reproducing historical climate extreme indices. The observational datasets include NOAA Atlas 14, volume 9 (Perica and others, 2013), the Parameter-elevation Regressions on Independent Slopes Model (PRISM; Daly and others, 2008, 2021) and the South Florida Water Management District’s Precipitation “Super-grid” (SFWMD, 2005).

Downscaled Climate Datasets

Various efforts to downscale global-climate predictions to local and regional scales have been initiated by climate research groups in the United States and abroad; however, it is notable that most of these downscaled datasets have not been tuned specifically to Florida climatic conditions. In this report, historical and future projections of precipitation based on statistical, dynamical, and hybrid downscaling of GCMs were used to develop future DDF curves for central and south Florida. The GCMs were developed as part of the World Climate Research Programme’s CMIP5 and CMIP6. The CMIP5 model data are a concatenation of historical (retrospective) GCM simulations covering the period 1850–2005 and future projections for the period 2006–2100 (https://pcmdi.llnl.gov/mips/cmip5/data-access-getting-started.html). The CMIP6 historical simulation period is 1950–2014 and future projections are for the period 2015–2100 (https://pcmdi.llnl.gov/CMIP6/Guide/index.html). CMIP6 data are only used by the JupiterWRF downscaled climate dataset that will be discussed in the next sections. See Taylor and others (2012) for an overview of the CMIP5 experimental design and Eyring and others (2016) for CMIP6. It is important to note that the historical CMIP5 and CMIP6 model simulations are not intended to reproduce the precise sequence of historical climate variability. GCMs are often run starting with different initial conditions to create an ensemble of possible trajectories for the historical and future climates that may result because of changes in natural variability (unforced or internal variability, also called “climate noise” in Taylor and others, 2012). This helps separate the climate change “signal” from the “climate noise.” However, often only one ensemble member is used for downscaling, which limits the range of events that are downscaled.

The global performance of the CMIP5 GCMs is summarized in the IPCC’s Fifth Assessment Report (Flato and others, 2013, p. 741–866). The overall model performance for the CMIP5 models in terms of the features of the 20th century climate most relevant to this study is summarized in table 2. For example, the Fifth Assessment Report (Flato and others, 2013, fig. 9.44) states that there is high confidence that the CMIP5 model performance in simulating ENSO is medium. Confidence in the validity of a finding is expressed qualitatively (see Mastrandrea and others, 2010).

The CMIP5 future projections are based on four different RCPs corresponding to low (RCP2.6), medium-low (RCP4.5), medium high (RCP6.0), and high (RCP8.5) year 2100 total radiative forcing values with respect to the preindustrial period (circa 1750; IPCC, 2013). The number next to the RCP label indicates the approximate increase in total radiative forcing, in watts per meter squared, caused by GHG emissions. RCP2.6 represents scenarios in which global mitigation of GHG emissions is significant, and the radiative forcing increase at 2100 from GHG emissions is 2.6 watts per meter squared. RCPs with higher numbers are the result of higher emissions and result in larger changes in global temperatures. RCP8.5 represents future conditions that could result from limited or no climate change mitigation, whereas the two middle scenarios are roughly equally spaced between the low and high scenarios (Terando and others, 2020). It is important to note that RCP8.5 is not necessarily a business-as-usual nor a worst-case scenario, but representative of a plausibly higher level of GHG concentrations. These RCP scenarios are developed using a risk-based framework and are highly uncertain. For this reason, individual RCP scenarios have not been formally assigned a likelihood of occurrence but represent plausible outcomes. Similar to the historical simulations, the CMIP5 future projections are also not intended to simulate the precise sequence of actual future variations in climate and may not capture the timing of future shifts in natural cycles.

Terando and others (2020) recommend that, whenever feasible, the entire range of RCP-based scenarios be considered when assessing potential future impacts. As part of the FPLOS Program, and with the purposes of planning for future stormwater infrastructure projects with long design lifetime, the SFWMD is interested in the evaluation of the middle-low and high emission scenarios RCP4.5 and RCP8.5. The choice of scenarios is also limited by their availability in downscaled climate datasets. The inherent assumption here is that RCP4.5 and RCP8.5 may result in higher extreme precipitation amounts for durations and return periods of interest compared to RCP2.6 because of increased atmospheric warming and higher atmospheric water-holding capacity resulting from these selected emission scenarios. This general assumption does not account for potential modulating factors that could affect precipitation variability and extremes in Florida, such as future changes or shifts in large-scale atmospheric circulations, teleconnection effects, or changes in tropical cyclone intensity. For example, Infanti and others (2020) analyzed precipitation from the Bias-Corrected Spatially Disaggregated (Maurer and others, 2007; Bureau of Reclamation, 2013) dataset for central and south Florida. After subsetting for the best performing models in terms of their ability to capture large-scale sea surface temperature anomalies and regional 2-meter temperature, they found that areas south of Lake Okeechobee may receive less wet-season precipitation in the future under both RCP4.5 and RCP8.5 scenarios. In particular, wet events (those with a Standardized Precipitation Index, SPI, above 0.5) during the wet season would become less wet in the middle-term (2046–72) under RCP4.5 and in the long-term (2073–99) under RCP8.5. The Sea Level Solutions Center at Florida International University (2021) evaluated output from the BCCA downscaled climate dataset (Maurer and others, 2007; Bureau of Reclamation, 2013), for changes in the seasonality of precipitation in Florida and found a decrease in wet-season precipitation and an increase in early dry-season precipitation over the NOAA Climate Division 5 (Everglades; fig. 1). Whether these findings hold for precipitation extremes of interest in this study remains to be seen.

The Coupled Model Intercomparison Project Phase 3 (CMIP3) defined emission scenarios as based on the following notation: economic (A) versus environmental (B) focus, and global (1) versus regional (2) responses (IPCC, 2000). Technological emphasis is added by means of an additional set of letters: fossil-fuel intensive (FI), nonfossil-fuel energy sources (T), or balance across all energy sources (B). The A1B scenario is considered similar to “medium” emissions scenario RCP6.0 from CMIP5 (Walsh and others, 2014). Based on the A1B intermediate emissions scenario simulated with the previous generation of climate models as part of CMIP3, Misra and others (2011) found that a decrease in June–August precipitation and an increase in September–November precipitation is projected in south Florida for the period 2080–2100 with respect to 2000–20. They attribute this decrease in June–August precipitation to the broad summer drying expected for the Caribbean region in CMIP3 models as well as a potential decrease in sea-breeze frontogenesis because of the potential future inundation of the Shark River Slough in Everglades National Park (fig. 1). Potential drying of the Caribbean region is also observed in CMIP5 models with a simulated global warming of 2.0–2.5 °C (Taylor and others, 2018) and is likely by the end of the 21st century under the RCP8.5 scenario (Collins and others, 2013). The simulated drying is associated with a narrowing of the Intertropical Convergence Zone resulting in more intense convection organized over narrower regions and a drier subtropical atmosphere caused by an enhanced and widened subsidence region (Byrne and others, 2018). On the basis of the ensemble mean of six perturbed physics dynamical downscaling simulations of the HadCM3 CMIP3 model under the A1B scenario, Campbell and others (2021) show that the Caribbean drying extends into south Florida during February–April at 1.5 °C global warming, during November–January at 2.0 °C, and during May–July at 2.5 °C. However, these prior studies did not examine whether extreme events might follow the same trends as changes in the mean and in the less extreme quantiles analyzed in these studies.

Most of the downscaled climate datasets used in this study rely on CMIP5 climate data except for the analog resampling and statistical scaling method developed by Jupiter Intelligence (2021), which also uses CMIP6 data. The CMIP6 projection scenarios are described by Eyring and others (2016) and are based on the concept of “shared socioeconomic pathways” (SSPs), which have been developed by the energy modeling community. Updated versions of the four CMIP5 RCP scenarios are available in CMIP6 and called “SSP1-2.6,” “SSP2-4.5,” “SSP4-6.0,” and “SSP5-8.5,” with each of these having a similar change in 2100 radiative forcing levels as their CMIP5 RCP counterpart. However, even though the 2100 radiative forcing is similar, the pathways of emissions, GHG mix, and land uses vary over time between the RCP and corresponding SSP scenario. For this reason, results from CMIP5 and CMIP6 may not be directly comparable. The use of different GCMs or different versions of GCM models in the two projects also precludes a direct comparison. CMIP6 also includes four additional SSPs representing intermediate levels of forcing between the four original scenarios. The CMIP6 model data are a concatenation of historical (retrospective) GCM simulations covering the period 1850–2014 and future projections for the period from 2015 to at least 2100. A description of the GCMs downscaled by the different downscaled climate datasets evaluated in this study can be found in table 3.

Many GCMs have poor skill in simulating extremes because of their coarse spatial resolution and because of the difficulties in capturing subgrid scale physics (Misra and others, 2011). Statistical and dynamical downscaling are two methods used to generate high-resolution climate projections based on large-scale fields simulated by GCMs.

Localized Constructed Analogs (LOCA)

The University of California at San Diego has developed the LOCA statistical-downscaling technique (http://loca.ucsd.edu) to downscale 32 GCMs from the CMIP5 archive at 1/16th degree (approximately 6.6 km) spatial resolution, covering North America from central Mexico through southern Canada. The historical period for this dataset is 1950–2005, and two future projected scenarios are available: RCP4.5 and RCP8.5 over the period 2006–2100, although data are only available through 2099 for some models. Only one ensemble member from each GCM is downscaled by the LOCA method. This means that potential alternative trajectories for the historical and future climates that may result from changes in natural variability are not considered. The interaction of natural variability and secular trends caused by climate change and their effects on precipitation at various scales, including extremes, is therefore incompletely sampled. The Fourth National Climate Assessment (Avery and others, 2018) relied on the LOCA dataset as a source of downscaled climate information. For this study, daily precipitation data for the State of Florida were downloaded from the USGS Center for Integrated Data Analytics THREDDS data server (USGS, 2020a, b). The CMIP5 models downscaled by using the LOCA method and used in this study are listed in table 4. The LOCA grid in central and south Florida is shown in figure 4A.

Table 4.    

General circulation models downscaled by the Localized Constructed Analogs (LOCA) statistical-downscaling method and evaluated in this study.

[HIST, historical; RCP, representative concentration pathway; ensemble member r1i1p1 used for each model for both RCP4.5 and RCP8.5 with the exception of CCSM4_r6i1p1 for RCP4.5, CCSM4_r6i1p1 for RCP8.5, EC-EARTH_r8i1p1for RCP4.5, EC-EARTH_r2i1p1for RCP8.5, GISS-E2-H_r6i1p3, GISS-E2-R_r6i1p1, GISS-E2-H_r2i1p1for RCP8.5, and GISS-E2-R_r2i1p1for RCP8.5]

Table 4.    General circulation models downscaled by the Localized Constructed Analogs (LOCA) statistical-downscaling method and evaluated in this study.

General circulation model	Scenarios available
ACCESS1.0	HIST, RCP4.5, RCP8.5
ACCESS1.3	HIST, RCP4.5, RCP8.5
BCC-CSM1.1	HIST, RCP4.5, RCP8.5
CanESM2	HIST, RCP4.5, RCP8.5
CCSM4	HIST, RCP4.5, RCP8.5
CESM1 (BGC)	HIST, RCP4.5, RCP8.5
CESM1 (CAM5)	HIST, RCP4.5, RCP8.5
CMCC-CM	HIST, RCP4.5, RCP8.5
CMCC-CMS	HIST, RCP4.5, RCP8.5
CNRM-CM5	HIST, RCP4.5, RCP8.5
CSIRO-Mk3.6.0	HIST, RCP4.5, RCP8.5
EC-EARTH	HIST, RCP4.5, RCP8.5
FGOALS-g2	HIST, RCP4.5, RCP8.5
GFDL-CM3	HIST, RCP4.5, RCP8.5
GFDL-ESM2G	HIST, RCP4.5, RCP8.5
GFDL-ESM2M	HIST, RCP4.5, RCP8.5
GISS-E2-H	HIST, RCP4.5, RCP8.5
GISS-E2-R	HIST, RCP4.5, RCP8.5
HADGEM2-AO	HIST, RCP4.5, RCP8.5
HADGEM2-CC	HIST, RCP4.5, RCP8.5
HADGEM2-ES	HIST, RCP4.5, RCP8.5
IPSL-CM5A-LR	HIST, RCP4.5, RCP8.5
IPSL-CM5A-MR	HIST, RCP4.5, RCP8.5
MIROC-ESM	HIST, RCP4.5, RCP8.5
MIROC-ESM-CHEM	HIST, RCP4.5, RCP8.5
MIROC5	HIST, RCP4.5, RCP8.5
MPI-ESM-LR	HIST, RCP4.5, RCP8.5
MPI-ESM-MR	HIST, RCP4.5, RCP8.5
MRI-CGCM3	HIST, RCP4.5, RCP8.5
NorESM1-M	HIST, RCP4.5, RCP8.5

The LOCA method is a statistical downscaling technique that uses past historical observations to add improved fine-scale detail to global climate models (Pierce and others, 2014). The historical observational gridded precipitation dataset used for LOCA is the Livneh and others (2015) dataset over the period 1950–2005. Pierce and others (2021) found that daily precipitation extremes are muted by about 30 percent in the Livneh and others (2015) dataset over most areas of the continental United States, including Florida, and attribute this to the way Livneh and others (2015) split daily precipitation measurements across 2 days, depending on time of observation. A new gridded daily precipitation dataset over the continental United States that preserves extremes has been developed by Pierce and others (2021); however, the LOCA dataset has not been updated on the basis of this new dataset.

Figure 5, which is adapted from Lopez-Cantu and others (2020), shows an overview of the statistical downscaling methodology used for LOCA and described by Pierce and others (2014), which is summarized below. Before downscaling, LOCA performs bias correction on the daily precipitation data. First, LOCA uses a preconditioning approach to correct the annual cycle of daily values. Then, LOCA applies the Preserves Ratio (PresRat) bias-correction scheme to daily precipitation data, as described by Pierce and others (2015). The PresRat scheme consists of applying a multiplicative quantile delta mapping (MQDM) technique to the daily precipitation (described in the “Multiplicative Quantile Delta Mapping” section herein), special treatment of zero-precipitation days, followed by a final correction factor that tries to preserve the mean precipitation change predicted by the GCM as a percentage of the GCM’s historical climatology. As shown in Pierce and others (2015), the final corrections required to maintain the GCM-predicted mean precipitation change tend to be small and between 0.95 and 1.05 for most of Florida during the months of January, April, July, and October. After PresRat bias correction, LOCA also applies a frequency-dependent bias-correction scheme to reduce inaccuracies in the GCMs’ spectra. The spectrum of a time series describes the strength of the variation of a time series over different frequencies, such as daily, annual, and decadal timescales. The frequency-dependent bias-correction method is implemented as a digital filter in the frequency domain to make the model spectrum better match that of the observations.

Implementing the LOCA downscaling method involves the following process. The continental United States is split into regions, and a point near the center of each region is selected and designated as an analog pool point. A spatial weighting mask around each analog pool point is then developed that includes all areas having a positive temporal correlation with the data at the analog pool point. The weighting mask is computed from coarsened observations and extends beyond the region that is represented by the analog pool point. The weighting mask varies by season and is used to limit the region over which the model field for the variable being downscaled is compared to coarsened observations of the variable. Analog days are then selected at the regional scale. For each day that is to be downscaled, a pool of 30 candidate observed analog days is chosen at each analog pool point by matching the model field for the variable being downscaled to days from coarsened observations that have the lowest root-mean-square difference (RMSD). This is done over the masked region corresponding to the analog pool point. This method contrasts with other analog downscaling methods in which the same analog days are chosen for the entire domain that is being downscaled. Only analog days within 45 days of the day of the year being downscaled are allowed. Then the single candidate analog day representing the best match within the local area around the fine-scale grid cell being downscaled is chosen as the single analog day to use for that location. This local matching is done by comparing the modeled field to the coarsened observations for the analog days after interpolating them to the fine-scale grid by using bicubic interpolation. The analog day with the lowest RMSD over a region of 21 fine-scale grid cells around the fine-scale grid cell being downscaled is then selected as the analog day for that point. In the case of precipitation, the selected analog day is scaled by the ratio of the interpolated model field to the interpolated analog day. For most grid cells, only the single locally selected analog day is used in downscaling. However, to reduce discontinuities, locations for which neighboring cells have a different analog day use a weighted combination of the center and adjacent analog days. In contrast to LOCA, other constructed analog methods typically use a weighted mean of the same 30 analog days for the entire domain. LOCA reduces this averaging and is therefore expected to (1) produce better estimates of extreme days, (2) more realistically depict the spatial coherence of the downscaled field (which other analog methods tend to overestimate), and (3) reduce the problem of producing drizzle (light-precipitation) days resulting from other analog downscaling methods. It is worth noting that GCMs themselves have a tendency for drizzle, as documented by Stephens and others (2010) and Pendergrass and Hartmann (2014b), and it is possible that some of this drizzle tendency may be persistent in downscaled datasets.

Deficiencies in the Livneh and others (2015) historical observational dataset used in LOCA have recently been documented by Wootten and others (2021), and Pierce and others (2021). These deficiencies arise because of the way the Livneh and others (2015) split observations of daily precipitation on the basis of their time of observation, where the precipitation total is prorated by the number of hours overlapping the date of the observation. For example, if the observations are made at 8 a.m., one-third (8/24) of the observed precipitation for a given day is assigned to the current day, whereas two-thirds (16/24) of the observed precipitation is assigned to the previous day. This method differs from other gridded datasets such as Daymet and PRISM, which use daily precipitation amounts recorded at any time during a day as the daily total for the 24-hour period from 12 a.m. to 12 a.m. and 7 a.m. to 7 a.m. local time, respectively. The apportioning process in Livneh increases the frequency of wet events, increases the wet spell length, and reduces the intensity of events (Oyler and Nicholas, 2018). Pierce and others (2021) found that daily precipitation extremes are muted by about 30 percent in Livneh over most areas of the continental United States, including Florida. This affects single-day events more than multiday events in both Livneh (Wootten and others, 2021) and the derived LOCA dataset. Previous analyses of DDF curves fit to annual maximum series of precipitation from LOCA’s historical period by Irizarry and others (2016), indicate that return levels for daily maximum precipitation are underestimated by about 40 percent in LOCA and that the underestimation is reduced to about 10 percent for 7-day duration return levels. These findings are generally consistent with those of Pierce and others (2021) and Wang and others (2020), which indicate that a large portion of the bias in daily extremes in LOCA may be due to the deficiencies in the Livneh observational dataset. For Florida in particular, Behnke and others (2016) found that out of the seven gridded datasets evaluated, Livneh performed the worst in capturing extreme precipitation climate indices at meteorological stations that are part of the Florida Automated Weather Network (FAWN, http://fawn.ifas.ufl.edu).

Extreme Value Theory

Extreme value theory is a branch of statistics that deals with the stochastic behavior of extreme and rare events found in the tails of probability distributions and includes events larger or smaller than usual. The probability of occurrence of a climate or weather variable can be described by a probability distribution function. Under global climate change, small changes in the mean of a distribution may result in changes in the frequency of extremes at both ends of a distribution (Cubasch and others, 2013). Furthermore, changes in the variance and skewness or shape of the distribution may also result in substantial changes in the frequency and magnitude of extremes. For stormwater planning in particular, the concern is how the upper tail of the probability distribution, representing the larger values of extreme precipitation, will change under future conditions. Extreme value theory provides a class of models that make it possible to extrapolate from observed or modeled extremes to unobserved or unmodeled extremes (Coles, 2001).

Extreme value theory has developed under two main approaches. In the classical block maxima (or minima) approach, the probability distribution of maxima (minima) over a given block of time, typically a year, is modeled. This is herein referred to as the “annual maxima approach.” The peaks-over-threshold (POT) approach consists of extracting, from a continuous record, the peak values that exceed (or fall below) a particular threshold and modeling their distribution. Both approaches can be implemented under stationary and nonstationary frameworks. Stationarity implies that, given any subset of variables, their joint distribution stays the same through time, which is the approach taken here. First-order stationarity was validated using the Mann-Kendall nonparametric trend test with a significance level (alpha) of 0.05. In contrast, nonstationary processes have characteristics that change systematically through time, for example, as seasonal effects or in the form of trends (Coles, 2001). These changes are reflected as changes in their probability distribution through time. Both extreme value approaches will be discussed hereafter in the context of the stationary framework.

Areal Reduction Factors (ARFs)

ARFs are used in stormwater infrastructure design to convert point precipitation for a given duration and return period into areal precipitation over an area for the same duration and return period. Typically, the area of interest is a watershed. ARFs are also useful when comparing simulated precipitation between datasets with different resolutions. As demonstrated by Chen and Knutson (2008), precipitation output from climate models should be interpreted as areal means rather than as point values, and their use and analysis should remain consistent with that interpretation. Otherwise, large differences in precipitation fields among datasets with different resolutions may be misinterpreted as bias. If the comparison is performed at consistent spatial scales, the true bias may be smaller or even of a different sign. In the computation of DDF change factors, described in the “Multiplicative Quantile Delta Mapping” section, current and future ARFs are assumed to be constant so that they would cancel out. However, ARFs are still useful when assessing model performance by comparing modeled extremes against observed extremes.

The most commonly known ARF curves were developed by the U.S. Weather Bureau (now the National Weather Service) and are documented in Technical Paper No. 29 (U.S. Weather Bureau, 1958). The Weather Bureau’s ARF curves can be used to transform point precipitation to areal precipitation for durations from 1 to 24 hours and areas covering as much as 1,000 km². The curves are based on annual maximum series data at gages in various networks throughout the United States. The gages used in computing ARF had relatively short record lengths (less than or equal to 16 years), which were averaged to obtain areal means of annual maximum series. The Weather Bureau ARFs consist of ratios of mean areal annual maximum series to station-based mean annual maximum series and represent ARFs for an approximately 2-year return period. In some of these earliest studies, ARFs were found to not exhibit systematic regional variations and to be independent of storm magnitude. Since then, ARFs have been studied in more detail and derived on the basis of longer and more spatially complete precipitation datasets for many regions of the United States and the world. ARFs have been found to vary according to the precipitation characteristics of the region under study. Svensson and Jones (2010) document some of these factors. For example, convective storms tend to be more localized than storms associated with frontal systems, resulting in smaller ARFs. In Florida, the seasonal pattern of these precipitation-generating mechanisms results in seasonal variations in ARFs. In some areas, ARFs have been found to decrease with increasing return period (Svensson and Jones, 2010; Kao and others, 2020), because the least frequent (most intense) storms tend to be convective and more localized in these areas. Studies have shown that ARFs decrease with decreasing duration (Olivera and others, 2006; Kao and others, 2020), because short-duration events tend to be more localized. For this reason, ARFs also tend to decrease with area more steeply for shorter durations than for longer durations. ARFs have also been found to vary depending on catchment characteristics, such as shape, topography, and urbanization, although this factor seems to be small (Svensson and Jones, 2010). Proximity to the coast and large water bodies may also affect ARFs.

The characteristics of the data used in developing ARFs, such as period of record, station density, interpolation method, and the methodology used in ARF estimation, will result in varying ARF estimates. Two main approaches for ARF computation are typically used: geographically fixed area ARFs, and storm-centered ARFs. Svensson and Jones (2010) present a review of these different ARF approaches. Storm-centered approaches define the point precipitation as the precipitation at the point having the maximum precipitation in a storm, which varies from storm to storm (Omolayo, 1993). Because of the complexity of this approach in view of storms with multiple centers and the fact that extreme point precipitation and extreme areal precipitation may not be produced by the same storm (Omolayo, 1993), fixed-area methods are most widely utilized and will be used here. In addition, storm-centered ARFs tend to be smaller than fixed-area ARFs and are, therefore, less conservative (Svensson and Jones, 2010).

To evaluate the performance of modeled precipitation extremes from datasets of various resolutions by comparison against observational datasets for a chosen historical period, grid-cell (areal) scale DDF curves were divided by ARFs to approximate station-scale (point-scale) DDF curves. The model-derived station-scale DDF curves for the historical period can then be compared to station-scale DDF curves from NOAA Atlas 14 volume 9 (Perica and others, 2013) or other sources. Because the station data used in developing the NOAA Atlas 14 DDF curves may include data for the period 1840–2008, depending on the station, all the precipitation data available from the historical period (1950–2005) of the downscaled climate datasets were used for this comparison. The two sets of historical DDF curves may not be completely comparable due to the differences inherent in the data and methods used to develop them, which include (1) different periods of record, and (2) different methods for DDF curve fitting (regional frequency analysis on annual maximum series for NOAA Atlas 14 versus CML for the downscaled climate model datasets). Therefore, differences between the two are to be expected.

Pavlovic and others (2016) describe four fixed-area methods that are commonly used in deriving ARFs. For this study, methods M1 and M4 from Pavlovic and others (2016) were tested. Method M1 is similar to that developed by the U.S. Weather Bureau, and results in ARFs that are a function of duration and spatial aggregation scale. Method M4 uses a frequency analysis approach whereby ARFs are estimated as the ratio of quantiles of GEV fits to areal and point annual maximum series for a given return period. Functional forms relating GEV parameters to area and (or) duration can be fitted for smoothing the resulting GEV fits, as done by Overeem and others (2009, 2010) and Pavlovic and others (2016), and ARFs can be derived from these fits. These ARFs will be a function of duration, spatial aggregation scale, and return period. In their review of various ARF methods, Svensson and Jones (2010) describe the advantages of using the Bell (1976) empirical ARF method, which is similar to method M4. The Bell method involves fitting a frequency distribution to areal and point annual maximum series and taking their ratio for the same return period.

The ARF approaches used in this study use the PRISM precipitation dataset discussed in the “Parameter-elevation Regressions on Independent Slopes Model (PRISM)” section. On the basis of the PRISM daily precipitation data from 1981 to 2019, annual maximum series were computed at various durations and spatial scales, which are required for most ARF methods. The procedure for computing the annual maximum series is similar to that of Pavlovic and others (2016); however, the spatial and temporal aggregation steps were reversed, which results in equivalent time series in the end. For PRISM grid cells in central and south Florida, spatial aggregation of the gridded data was performed first. This step consisted in first identifying areas of 1 × 1, 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11, and 13 × 13 grid cells centered at each grid cell of interest and computing the spatial mean of the daily precipitation. These correspond to areas of approximately 18.9; 170.2; 472.8; 926.7; 1,531.9; 2,288.4; and 3,196.1 km², respectively. If any of the n × n cells in an area had no precipitation data (for example, ocean cells in the PRISM dataset), then that area was excluded from subsequent computations. Therefore, areal means were obtained for a smaller and smaller number of cells as the spatial aggregation scale increased. The second step consisted of temporal aggregation of the areal mean precipitation to the durations of interest (1, 3, and 7 days). For durations of 3 and 7 days, rolling sums of the areal mean daily precipitation at each spatial scale were computed. The rolling sum time series can be used in implementing a CML approach to fitting the GP to develop ARF curves by a variation of method M4. Finally, the annual maximum value was extracted for each duration and spatial aggregation scale for each year in the period 1981–2019. Because of the order of operations, the annual maximum series for a given aggregation scale is the mean of concurrent precipitation values; however, annual maximum series at different aggregation scales may not be concurrent. Corrections from constrained to unconstrained annual maximum series were not applied, because these are assumed to cancel out in the computations.

In method M1, ARFs were computed from the temporal mean of the annual maximum series for each duration and spatial aggregation scale at each central grid cell according to

In theory, ARF values should be less than one; however, occasionally, the ARF computed by equation 28 ends up slightly greater than one, in which case, that cell is eliminated for further analysis for the duration and spatial aggregation level where it occurs. This is a result of the mean annual maximum series at a grid cell being smaller than the mean annual maximum series for the n × n area centered at the grid cell because of the fact that the annual maximum series at different spatial scales may not be concurrent. Once the ARF is computed at each central grid cell based on equation 28, the mean of the ARF values for the aggregation scales and durations of interest are computed for climate regions defined in the State of Florida:

Climatologically homogeneous precipitation areas determined by Pathak and others (2009) on the basis of rain-gage adjusted Next-Generation Radar (NEXRAD) data were initially used to define climate regions for ARF averaging. These homogeneous precipitation areas were defined on the basis of their mean precipitation and the spatial correlation of precipitation during the wet season in south Florida (June–October). One limitation of these homogeneous precipitation areas is that their spatial extent does not include areas outside the SFWMD boundary in central Florida. In addition, Welch’s t-tests for the null hypothesis of equal means of ARF across homogeneous precipitation areas showed that the mean ARF for some of these homogenous precipitation areas was very similar for the durations and spatial aggregation scales evaluated. Lumley and others (2002) state that in the case of large samples, the t-test is valid for any distribution, not just for normal distributions. They show that the test is valid for sample sizes as low as 65, even in the case of very skewed distributions. All of the homogeneous precipitation areas evaluated here contain between 97 and 315 PRISM grid cells at which ARF is calculated; therefore, the sample size is considered adequate for use of the Welch’s t-test for equality of means.

According to Pathak and others (2009), convection in the center of the southern Florida peninsula results either from local daytime heating or from the collision of thunderstorm outflow boundaries propagated away from storms associated with the east and west sea breezes and the Lake Okeechobee lake breeze. These mechanisms likely result in more consistent and less variable precipitation in these interior areas, relative to other areas, and similar ARFs. On the basis of the t-test results, some of these homogeneous precipitation areas were merged into larger climate regions. For areas outside the SFWMD, the NOAA National Centers for Environmental Information (NCEI) Florida climate divisions (Guttman and Quayle, 1996) were used to define climate regions for ARF averaging in central Florida (fig. 7). Shapefiles of the boundaries of the NOAA NCEI U.S. Climate Divisions in the State of Florida are available from NOAA (2011).

A modified version of ARF method M4 in Pavlovic and others (2016) was also tested by fitting the GP to rolling sums of durations 1, 3, and 7 days at once using the CML approach for each spatial aggregation level of interest. The ARFs for this method consist of ratios of DDF curves at a spatial aggregation level of interest to the 1 × 1 aggregation level and are a function of duration, return period, and spatial aggregation level.

Multiplicative Quantile Delta Mapping

Quantile mapping (QM), a CDF-matching method (Panofsky and Brier, 1968), is typically applied to bias-correct precipitation time series from climate model simulations, but the method can similarly be used to bias-correct annual maximum series of precipitation or DDF curves. It can also be used on its own as a form of statistical downscaling that attempts to bridge the scale mismatch between the model grid-cell estimates and the point observations. However, this can result in misrepresentation of the spatiotemporal structure of the corrected time series and overestimation of area-mean extremes (Maraun, 2013). The expression for QM is given by

The CDFs are typically developed from data spanning decades and centered around some year of interest. QM only uses information from the historical period to correct for future biases; therefore, it assumes that biases are stationary and that they will persist into the future (Cannon and others, 2015). In other words, QM assumes that F_o,p (the CDF of the observations in the future projection period) is approximated by F_o,h so that as the mean changes, the variance and skew do not, which may not be the case under climate change (Zhang and Zwiers, 2013; Pendergrass and Hartmann, 2014a; Pendergrass and others, 2017). Furthermore, if a future projected value is outside the historical range, then some sort of extrapolation is required and that affects extremes (Wootten, 2018).

To avoid the limitations of traditional QM, other methods have been developed such as quantile delta mapping (QDM). As shown in Cannon and others (2015), QM tends to inflate trends in precipitation extreme indices projected by GCMs, whereas QDM is not as prone to this problem. QDM preserves model-projected changes in quantiles, while simultaneously correcting for systematic biases across quantiles (Cannon and others, 2015). QDM also attempts to bridge the gap between point estimates for the observations versus grid cell estimates in the model. However, it is important to note that changes in the mean may not be adequately preserved by QDM.

QDM can be applied in an additive form, which corresponds to the Equidistant Quantile Mapping algorithm of Li and others (2010), or a multiplicative form, which corresponds to the Maximum Intensity Percentile Based Method of Wang and others (2013). Additive QDM is better suited to correcting biases in air-temperature variables, whereas multiplicative quantile delta mapping (MQDM) is better suited to correcting variables such as precipitation where preserving relative changes is important in order to respect the Clausius-Clapeyron (CC) relation (Wallace and Hobbs, 2006), which relates the amount of atmospheric moisture to temperature changes. Lanzante and others (2021) confirmed the superiority of MQDM over additive QDM by means of a “perfect model” experiment.

In this study, the MQDM version of QM was chosen to bias-correct projected future DDF curves (fig. 8). Multiplicative QDM is given by the following alternative equations and is applied to each duration of interest independently:

The term $F_{o,h}^{-1}\left(F\right)/F_{m,h}^{-1}\left(F\right)$ in equation 32 is the bias-correction factor and the term $F_{m,p}^{-1}\left(F\right)/F_{m,h}^{-1}\left(F\right)$ in equation 33 is the change factor for the particular nonexceedance probability (or return period), and duration. The inherent assumption is that the relative change between historical and future projected conditions can be indicative of future changes in DDF curves even if the models are biased. The ultimate goal of this study is to obtain change factors from the historical period (1966–2005) to the future projection period (2050–89) across the locations of NOAA Atlas 14 stations in the SFWMD. These change factors could then be applied to the PDS-based DDF curves from NOAA Atlas 14 to obtain future projected DDF curves.

NOAA Atlas 14 stations in the Florida Keys (08-2418, 08-4570, 08-5035, and 08-8841 in appendix 1, table 1.1) are inactive in most downscaled model grids (fig. 4A–E). Data from the closest mainland grid cells to these four stations were used in developing DDF curves at these four locations. Therefore, change factors developed for these four locations are highly uncertain and should be used with caution.

Derivation of Change Factors

The following is a description of the analytical steps used to derive change factors based on CORDEX, LOCA, and MACA downscaled climate datasets (fig. 6). For each dataset, the closest grid cells to each of the 174 NOAA Atlas 14 stations in central and south Florida were identified for each of its models. For each of these grid cells, rolling sums of daily precipitation were computed using R package RcppRoll (Ushey, 2018) to obtain the cumulative 3- and 7-day precipitation totals. This process was repeated for each climate model run for historical and future projection periods for each RCP scenario in each dataset. Similar to the annual maxima approach used in developing NOAA Atlas 14 volume 9 by Perica and others (2013), it is expected that all extreme values might be downward-biased because of being constrained to a clock-based interval (that is, a consistent daily timestep). Therefore, rolling sums for the durations of interest were multiplied by correction factors (table 1; Perica and others, 2013) to convert constrained precipitation totals to unconstrained precipitation totals for each duration, that is, to the maximum precipitation for any time period of the given duration. In the annual maxima case, once constrained annual maxima are computed, a consistency check is typically made to ensure that constrained annual maxima for a given duration are larger than for a smaller duration. When that is not the case, an ad hoc small offset (usually on the order of 0.1 in.) is applied to the shorter duration annual maxima to obtain the annual maxima for a longer duration. However, the same process was not followed for the unconstrained rolling precipitation time series from which exceedances were computed because durations of interest (1, 3, and 7 days) are not a multiple of each other and, therefore, longer duration events do not contain a discrete number of shorter duration events. For example, a 7-day event overlaps two 3-day events and a 1-day event (or part of a 3-day event) at the beginning or end of the 7-day event. Despite this, checks were made to ensure that the selected exceedances for each duration are larger than for shorter durations overall by comparing their CDFs.

A POT approach was used to fit the GP probability distribution to exceedances above a high threshold as described in the “Peaks-Over-Threshold” section. This requires the selection of a high threshold of precipitation, which tends to be subjective and is difficult to automate to fit thousands of GP models without user intervention. Therefore, thresholds were defined for every location (grid cell) on the basis of a high percentile of the precipitation totals for each duration that would result in 2–3 exceedances per year (λ) for each duration after declustering. The chosen percentiles were 99th, 98th, and 97th for 1-, 3-, and 7-day durations, respectively. Because of the use of rolling sums of precipitation for various durations (d), there is a high degree of dependence between the d-day precipitation values and their exceedances over a threshold. This is in addition to the inherent temporal correlation in precipitation totals within a region. That is, the exceedances above a given threshold tend to cluster, especially as the threshold is reduced. Therefore, a declustering technique was implemented, and the goal was for the declustered time series to have 2–3 exceedances per year on average. The intervals method of Ferro and Segers (2003) was attempted to determine the ideal run lengths for declustering. However, the run lengths obtained from this method varied substantially across stations, ranging from about 1 to 30 or more days depending on duration. Because a physically based reason for the large variation in run lengths could not be identified, the runs declustering method in the extRemes package from R (Gilleland and Katz, 2016) was used with run lengths equal to the duration in days, plus 1 day. In other words, run lengths of 2, 4, and 8 days were chosen for 1-, 3-, and 7-day durations, respectively. Our choice of run length for the 1-day event is similar to that of NOAA (2022a), which found that dependent events in the time series of exceedances do not affect the precipitation frequency estimates significantly and that declustering the daily exceedances using a 1-day separation period between exceedance events would be adequate. The run length is the number of threshold deficits to be considered as starting a new cluster. These run lengths were chosen to exclude neighboring (clustering) rolling sum precipitation values that include the same high precipitation days. The GP is then fit to cluster maxima, hence the POT approach.

The extRemes package from R (Gilleland and Katz, 2016) was used to compute the interval-based extremal index of the exceedances declustered by the runs declustering method with prespecified run lengths. The objective of this exercise is to confirm that the exceedances declustered on the basis of prespecified run lengths are reasonably independent. A cutoff value of 0.7 for the interval-based extremal index, θ, was assumed to be indicative of adequate tail independence.

Subsequently, the GP was fitted one duration at a time using the traditional maximum likelihood (ML) approach and for all durations at once using the CML approach at each model grid cell. The resulting GP fits based on the CML approach are used to derive DDF curves by using the quantile function of the GP for each duration and converting the return periods of interest into nonexceedance probabilities. The CML approach results in consistent DDF curves across durations at each model grid cell but slightly compromises GOF compared to the traditional ML approach of fitting one duration at a time; however, the CML approach also reduces the number of parameters that are being fit. This tradeoff between model fit and parameter parsimony is quantified by means of the AIC. In addition, a bootstrapping approach was used to determine GOF for the GP.

The resulting DDF curves apply at the grid-cell scale and must be divided by appropriate ARFs to obtain station-scale values for both historical and future projection periods. Assuming the historical ARF applies to the future projection period, the ARF would cancel out in the computation of change factors. However, the ARF allows for comparison of the station-scale DDF curves fitted to simulated historical precipitation data as part of this study against the NOAA Atlas 14 PDS-based DDF curves that were fitted to historical observations of precipitation. Finally, the MQDM method was used to compute change factors as the ratio of the future projected DDF precipitation depth to the historical DDF precipitation depth for each station grid cell in each downscaled climate dataset.

Model Culling

Model culling refers to the process by which some models are eliminated from further analysis, and the remaining subset of climate models (called “best” models, hereafter) are selected to inform future changes in impact-relevant climate variables. The “best” models are defined on the basis of a series of performance measures of relevance to the region and variables of interest. The performance measures are typically based on how well the models reproduce historical observations at weather stations or in comparison to observational gridded products in the region of interest. For this study, precipitation extremes and their return periods in central and south Florida are of the utmost importance. The use of the word “best” within double quotes is to emphasize that these models are the ones that best reproduce historical observations, but we do not imply that these models would necessarily perform best in simulating the future climate or result in the most accurate change factors.

The Expert Team on Climate Change Detection and Indices (ETCCDI; http://etccdi.pacificclimate.org/list_27_indices.shtml) has defined climate extremes indices based on daily precipitation data. The original ETCCDI indices and variations thereof (Zhang and others, 2011) have been used by many researchers worldwide to assess the performance of GCMs and downscaled climate products (see, for example, Donat and others, 2013; Sillmann and others, 2013a, b; and Srivastava and others, 2020). Several ETCCDI precipitation extreme indices were selected to quantify model performance, as documented in table 9. The ETCCDI precipitation extreme indices were supplemented with additional indices of relevance to this study listed in table 9. Preliminary evaluation of precipitation extreme indices for the various downscaled climate datasets used in this study showed correlation between some of the indices and that only a very small number of models remained when model culling was based on all the indices listed in table 9. Based on the desire to include as many models as possible to inform the potential range of future change factors, a decision was made to only use four precipitation extremes indices in selecting the best models. These selected indices consist of annual maxima of consecutive precipitation for durations of 1, 3, 5, and 7 days, and are the most relevant to the estimation of precipitation extremes.

Two different reference observational datasets were used to evaluate historical precipitation extremes indices in this study: the SFWMD “Super-grid” and PRISM. Various studies have found that the model performance in capturing climate extremes indices varies depending on the reference dataset chosen (Gleckler and others, 2008; Sillmann and others, 2013a, b; Donat and others, 2014; Wang and others, 2020). The SFWMD “Super-grid” and PRISM observational datasets were chosen for validation of precipitation extreme indices over the gridded observational precipitation datasets used in bias correction of the downscaled data used in this study (Livneh, Daymet, and gridMET). This choice was made for several reasons. PRISM was chosen as a reference dataset on the basis of Behnke and others (2016), who found that out of the seven gridded datasets evaluated (including Daymet and Livneh), PRISM best captured daily precipitation statistics and climate extremes indices at meteorological stations from the Florida Automated Weather Network (FAWN). The R95p and R99p indices (table 9), which are defined as the annual total precipitation from days with precipitation above the 95th and 99th percentiles of precipitation in the base period, respectively, are exceptions where PRISM performed more poorly than the other gridded datasets. Zhang and others (2018) also evaluated five gridded precipitation products for the State of Florida including PRISM, Real-Time Mesoscale Analysis, and three interpolation methods for FAWN data. They found that PRISM resulted in better fit with the daily FAWN precipitation dataset when evaluated at stations. The SFWMD “Super-grid” dataset is considered by the SFWMD to be the most complete gridded precipitation dataset for south Florida. Here, the common 25-year historical period between PRISM and the downscaled datasets (1981–2005) was chosen as the base period for computation of percentiles in the percentile-based indices and for calculation and performance evaluation of the precipitation extremes indices.

The precipitation extremes indices are calculated using Python code developed by Gibson (2021). The Python code relies on Python bindings (Max Planck Institute for Meteorology, 2022a) to the Climate Data Operators (CDO; Max Planck Institute for Meteorology, 2022b). The Python code was evaluated to make sure it followed the ETCCDI climate extreme index definitions in table 9 and cross-validated against the CDO ECA climate indices package. All precipitation extremes indices were calculated at the native resolution of the downscaled climate model and observational datasets.

The CDO utility remapnn was used to remap the climate indices computed for the CORDEX, LOCA, and MACA downscaled climate datasets to the grids of the PRISM and SFWMD ”Super-grid” datasets using nearest-neighbor interpolation. Nearest-neighbor interpolation was used so as not to add information to the downscaled climate datasets that was not already there to begin with. In essence, the assumption was made that the node values from CORDEX, LOCA, and MACA models represent mean grid cell values, and the node value was assigned to all PRISM and SFWMD “Super-grid” locations having their centroid within a model grid cell. It is important to note that although the resolution of the indices may be the same for the reference and downscaled datasets after remapping, differences in climate indices are to be expected because of differences in the native scale of each dataset and the order of steps followed in remapping (that is, climate index computation first, followed by remapping, or vice versa). For example, climate indices are calculated on the basis of the approximately 4-km spatial mean of daily station precipitation that was used for each PRISM grid cell. In contrast, climate indices are calculated from the 25- to 50-km daily precipitation simulated for each CORDEX grid cell (which is representative of the mean daily precipitation at those scales as described by Chen and Knutson, 2008) and then assigned to all PRISM cells or SFWMD ”Super-grid” cells located within the area of the CORDEX grid cell.

To summarize the spatial pattern match of each climate index in the downscaled datasets to the reference dataset, various performance measures from Gleckler and others (2008), Sillmann and others (2013a), and Srivastava and others (2020) were used. These performance measures were evaluated separately for two distinct climate regions in central and south Florida (fig. 2). The climate regions were defined as the NOAA (2011) U.S. Climate Divisions for the State of Florida with Florida climate divisions 5, 6, and 7 merged into a single region named “south Florida” and called “climate region 5” hereafter.

The root-mean-square error (RMSE) statistic is used to summarize the performance of each model in capturing the climatology (long-term mean for 1981–2005) of each index calculated from the reference dataset at the model resolution:

The best models will have a small RMSE_m,I. The relative performance of a model in simulating the climatological mean of the reference observational dataset can be computing by normalizing the RMSE_m,I as follows:

The median is used for normalization to avoid undue influence by outlier models having unusually large errors. The best models will have a large negative NRMSE_m,I, indicating that the model performs better than most models. A large positive NRMSE_m,I indicates that the model performs worse than most models. To quantify the overall performance of a model in simulating the climatological mean of all indices, the Model Climatology Index (MCI) from Srivastava and others (2020) was used, but it was modified to use the mean of the NRMSE_m,I values over all indices for a particular model as opposed to the median.

In addition to the spatial pattern in the climatological mean, it is important to quantify how different models capture the spatial pattern of the inter-annual variability of each climate index from gridded observations. For this purpose, the inter-annual variability skill score (Chen and others, 2011; Jiang and others, 2015) was used, which is defined as

The best models will have a small IVSS_m,I. The relative performance of a model in simulating the interquartile range of the reference observational dataset can be computed by normalizing the IVSS_m,I as follows:

The best models will have a large negative NIVSS_m,I, indicating that the model performs better than most models. A large positive NIVSS_m,I indicates that the model performs worse than most models. To quantify the overall performance of a model in simulating the inter-annual variability of all indices, the Model Variability Index (MVI) from Srivastava and others (2020) was used but was modified to use the mean of the NIVSS_m,I values over all indices for a particular model as opposed to the median. Finally, the MCI and MVI for each model can be plotted as a scatterplot, and the best performing models will have the largest negative MCI and MVIs.

Constrained Maximum Likelihood

The CML approach was used to fit the GP distribution to declustered precipitation excesses for durations of 1, 3, and 7 days, for the periods 2050–89 and 1966–2005 at every grid cell closest to each of the 174 NOAA Atlas 14 stations in central and south Florida. Overall, after declustering the precipitation excesses using run lengths equal to the duration plus 1 day, it was determined that more than 95 percent of grid cells had an extremal index (as calculated from the intervals method of Ferro and Segers, 2003) of 0.8 or greater in both periods and for all models in all downscaled climate datasets. This increases confidence that the declustering was adequate for the vast majority of cells. This approach was applied to the CORDEX, LOCA, and MACA downscaled climate datasets. Figure 11 is an example of the declustering of precipitation excesses above the 99th percentile threshold (1.49 in.) for a LOCA grid cell centered at lat. 25°5ʹ37.5ʺ N, long. 80°5ʹ37.5ʺ W. The figure shows the threshold corresponding to the 99th percentile of the daily precipitation at this grid cell (1.49 in.) as well as the daily precipitation data. The original data above the threshold line were eliminated from the GP fitting by the declustering method because they are considered to be part of the same event. The declustered exceedances correspond to the maximum values for each cluster and were used for GP fitting.

Figure 12 shows the mean annual cycle of the number of 1-day declustered threshold exceedances for the historical period (1966–2005) and the future period (2050–89) as simulated in each downscaled climate dataset for the two main climate regions of interest in the SFWMD (fig. 2), south-central Florida (climate region 4) and south Florida (climate region 5). The 1-day threshold exceedance events were selected on the basis of the 99th percentile of the daily precipitation time series for each period and location independently; therefore, the actual threshold value used in defining the events varies between periods. The range of the mean annual cycle of the number of 1-day threshold exceedances in the historical simulations from different climate models in each downscaled climate dataset tends to match the annual cycle from PRISM, the SFWMD “Super-grid,” and the datasets used in bias correction reasonably well, with the exception of CORDEX. Each downscaled climate dataset was bias-corrected to different gridded observational datasets; therefore, as expected, they match the bias-correction datasets better than they match PRISM and the SFWMD “Super-grid” in the historical period. The exceptions are some CORDEX models that show a more uniform annual distribution of the number of 1-day threshold exceedance events than the observational and the other downscaled climate datasets. An investigation into CORDEX RCM runs driven by ERA-Interim reanalysis (Simmons and others, 2007) boundary conditions showed that the annual distribution of the number of 1-day threshold exceedances was also quite uniform in those runs, especially in south-central Florida. This indicates a problem with the way some RCMs simulate the timing of extremes, and this issue does not appear to be completely corrected by the bias-correction procedure. This may be a consequence of the inability of many raw CORDEX models to capture the annual cycle of precipitation discussed in the “Results” section and the bias-correction algorithm amplifying large simulated events during months when the observational datasets tend to show fewer extreme events. Of the RCMs driven by ERA-Interim reanalysis that were evaluated, the WRF model at 25-km and 50-km resolutions best captures the seasonality of 1-day threshold exceedances. However, other models such as CanRCM4 are better at capturing the mean annual cycle of precipitation than WRF.

Comparison of figures 9 and 12 shows a similar seasonality in the mean annual cycle of precipitation and the mean annual cycle of 1-day threshold exceedance events. However, the amplitude of the annual cycle of 1-day threshold exceedance events is smaller than that of precipitation. A decline in the mean precipitation and an even more pronounced decline in the number of 1-day threshold exceedance events is observed for July compared to June and August. This pattern is likely due to precipitation during July being more driven by local convection as opposed to tropical storms. This explanation is consistent with Winsberg (2020), who found that for most of the State, the monthly frequency distribution of torrential precipitation is bimodal, peaking in June and September. Winsberg (2020) explained the June peak as resulting from the Intertropical Convergence Zone moving north and having more of an influence on the local weather and the September peak as being caused by the greater frequency of tropical storms reaching the State during that month, which corresponds to the peak of the Atlantic hurricane season.

Figure 13 shows the median change in the mean annual cycle of the number of declustered threshold exceedance events from the historical period (1966–2005) to the future period (2050–89) for durations of 1, 3, and 7 days, as simulated in each downscaled climate dataset for south-central Florida (climate region 4) and south Florida (climate region 5). In the case where the percentile-based threshold is allowed to vary between the historical and future periods, the median projected change in the mean annual cycle of the number of declustered threshold exceedance events in all datasets (fig. 13, solid lines) indicates a future decrease in the number of 1-day threshold exceedances from April to September in CORDEX, and from June to September in LOCA and MACA in both climate regions. MACA also shows a decrease in May for climate region 5 in south Florida (fig. 13B). A future increase in the number of 1-day threshold exceedances is predicted by all datasets in October and to a lesser extent, in some dry-season months. These future changes are consistent with those of the mean annual cycle of precipitation, although more months are affected. However, these conclusions are based on the threshold changing from the historical to the future period based on predefined percentiles, consistent with the way the threshold is defined in fitting the GP distribution in each individual period. When the percentile-based historical threshold is used consistently to define exceedance events for the historical and future periods, the median change in the number of threshold exceedance events (fig. 13, dashed lines) is generally higher than in the case where the threshold varies between periods (fig. 13, solid lines). This means that the percentile-based historical threshold will be exceeded more often in the future period than in the historical period according to the models; viewed another way, it means that the threshold defined on the basis of percentiles within a period increases from the historical to the future period. The main exception to this is LOCA in climate region 5 (south Florida) for all three durations, where on an annual basis, the median change in the number of threshold exceedance events is slightly higher in the case when the threshold varies between periods (fig. 13B, D, F, “LOCA varying threshold”) than in the case where the historical threshold is used in both periods (fig. 13B, D, F, “LOCA historical threshold”). This indicates that the percentile-based thresholds decrease slightly from the historical to future periods for all three durations in LOCA. Another exception is the 7-day duration from MACA in climate region 5 (south Florida), where on an annual basis, the median change in the number of threshold exceedance events is slightly higher in the case when the threshold varies between periods (fig. 13F, “MACA varying threshold”) than in the case where the historical threshold is used in both periods (fig. 13F, “MACA historical threshold”). Overall, the median change in the number of future exceedance events for all durations (fig. 13) follows a similar seasonal pattern as the projected median changes to the mean annual cycle of precipitation (fig. 10), regardless of whether the future events are defined on the basis of the historical threshold or using the future data to define the percentile-based future threshold. In both cases, the number of future exceedance events is generally projected to increase during the dry season (particularly in October) and decrease during the wet season.

Lastly, figure 14 shows the GP fits obtained from the CML method at the LOCA grid cell centered at lat. 25°5ʹ37.5ʺ N, long. 80°5ʹ37.5ʺ W. The corresponding DDF and intensity-duration-frequency curves are shown in figure 15.

Figure 16 shows boxplots of the distribution of the delta AIC (ΔAIC) for all grid cells in all models for each downscaled climate dataset in the historical (1966–2005) and future projection (2050–89) periods. The ΔAIC is being calculated with respect to the traditional ML model (that is, $AI{C}_{min}$ in equation 22 is the AIC of the traditional ML model). In this case, a ΔAIC value of −4 (−2 × 2, where 2 is the number of fewer parameters in the CML model) indicates that the negative log-likelihood of the CML and traditional ML models are the same and that GP fits have no degradation when using the CML model as opposed to the ML model. Values of ΔAIC greater than −4 indicate some degradation in the CML model compared to the traditional ML model, as reflected in larger negative log-likelihood values. Overall, median values of ΔAIC range from −3 to −2 for all downscaled climate datasets and periods. It is evident that more than 95 percent of grid cells have values below the ΔAIC thresholds of 2 and 7 for which the (simpler) CML model has substantial evidence and support compared to the traditional ML model, according to Burnham and Anderson (2002) and Burnham and others (2011), respectively. Very few cells exceed the $\text{Δ}AIC$ threshold of 7 units. The tradeoff between a slightly higher negative log-likelihood and a simpler model that results in consistent DDF curves automatically seems justifiable on the basis of these results.

Goodness of Fit

Figure 17 shows bar graphs indicating the percentage of model grid cells with p-values less than 0.05 for the lag-1 KACF and the MK trend test on excess values in the historical (1966–2005) and future projection (2050–89) periods. The null hypothesis that precipitation excess data are not autocorrelated is measured by the KACF, whereas the null hypothesis that it does not exhibit monotonic trends is quantified by the MK test. The nominal rate line in figure 17 represents the expected percentage of rejections of the null hypothesis when it is true, which is the significance level or expected nominal value of 5 percent. Also shown in this figure is the 95-percent prediction interval for the significance level under multiple comparison testing and under the assumptions that the data are uncorrelated and that the number of false rejections of the null hypotheses follows a binomial distribution as in Serinaldi and Kilsby (2014). The prediction interval is slightly wider for CORDEX than for LOCA and MACA, reflecting the smaller number of grid cells with NOAA Atlas 14 stations in them because of its lower resolution compared to LOCA and MACA. The percentage of rejections of the null hypothesis of no autocorrelation and no trend are close to the nominal value for most periods and durations. The main exception is CORDEX for which the null hypothesis is rejected more often than the nominal rate for the 1- and 7-day MK tests in the historical period. The 1- and 3-day MK tests for CORDEX in the future period also show slightly larger rejection rates than the nominal rate. Another exception is MACA for which the null hypothesis rejection rate is slightly larger than the nominal rate for the 3-day MK test in the historical period, for the 1- and 3-day MK test in the future projection period. The 1- and 7-day MK tests for LOCA in the future projection period also have rejection rates slightly larger than the nominal rate. The 7-day KACF test for MACA in the historical period also shows a rejection rate that is slightly larger than expected. It should be noted that the maximum rejection rate obtained for the MK tests is 8.2 percent, which is comparable to the 7–8-percent rejection rates corresponding to the 98th percentile-based threshold chosen by Serinaldi and Kilsby (2014) in fitting the GP to rainfall observations for the period 1970–2011. On the basis of these results, we conclude that it is acceptable to use a stationary approach in fitting the GP for these two 40-year periods.

Individual MACA models show more variability, with higher percentages of grid cells (as much as about 50 percent) for which the null hypothesis is rejected in a sizeable number of models, especially for the MK test. This result implies that the excesses in MACA are showing significant monotonic trends, most notably in the future projection period. This result is consistent with findings that MACA amplifies the change signal in the native GCMs (Lopez-Cantu and others, 2020; Wang and others, 2020), which would invalidate the assumption of stationarity within historical and future projection periods in a subset of MACA models, in which case a nonstationary model may be more appropriate. It is also possible that the presence of positive autocorrelation in the excess time series may be causing an overestimation of the significance of the computed trend, resulting in rejecting the null hypothesis of no trend more often than according to the chosen significance level of 0.05 (von Storch and Navarra, 1995). As discussed by Serinaldi and Kilsby (2015), when using nonstationary models, there is a high degree of subjectivity in defining how the parameters of the extreme value distribution might vary in time which increases uncertainty in return levels. There is also an inherent assumption that these parameter variations will hold for the entire future design life period. It is also difficult to parse out multidecadal natural variability influences from trends induced by climate change in nonstationary models. Therefore, for simplicity, the stationarity assumption is assumed to be valid.

Figure 18 shows bar graphs indicating the percentage of model grid cells for which the null hypothesis that the historical (1966–2005) and the future projected (2050–89) precipitation excess data follow a GP distribution is rejected at the 0.05 significance level for each duration in each downscaled climate dataset. The p-values are bootstrapped on the basis of the CML assumption, as discussed in the “Peaks-Over-Threshold” section and appendix 3. Results are only shown for the model grid cells closest to the 174 NOAA Atlas 14 stations used in this study. The figure shows a nominal rate of 5 percent and a 95-percent prediction interval for the significance level. Overall, CORDEX tends to perform worse for the 1-day duration compared to the other datasets, especially in the future period, with rejection rates that are much higher than the nominal rate of 5 percent that is expected by chance if, in fact, all data follow a GP distribution. For the 1-day duration, only MACA in the historical period frequently shows rejection rates below 5 percent. LOCA tends to perform worse than the other datasets for 3- and 7-day durations. Limited GOF testing based on the traditional ML approach with bootstrapping shows a similar percentage of rejections indicating that the performance deterioration when using CML does not have a large influence on the large test rejection rates for the GOF statistics and that, in fact, the raw excess data do not always seem to follow a GP distribution.

Figure 19 shows L-moment ratio diagrams for precipitation excesses in the historical (1966–2005) and future projection (2050–89) periods for each downscaled climate dataset. JupiterWRF is only included in the 1-day panel because 3- and 7-day durations were not evaluated for this dataset. The empirical L-moment ratios for all datasets shown in figure 19 tend to follow the generalized Pareto (GP) line shown in the figure. Although not shown, the centroid of the empirical L-moments for each individual model in a downscaled climate dataset and the overall centroid for each dataset and time period are located along the GP line. Also, as the duration increases the cloud of points shifts from the portion of the GP curve to the right of the exponential (“EXP”) point shown in the figure, which corresponds to positive shape parameters, to the portion of the GP curve on the left of the exponential point, which corresponds to negative shape parameters. This shift reflects a general decline in shape parameters with duration. Some LOCA points show lower L-skew and L-kurtosis especially for the 1-day duration, indicating a tendency for lower GP shape parameters and lower extremes in both periods. There is also a tendency for higher L-kurtosis in the future projection period (fig. 19B) than in the historical period (fig. 19A), indicating a tendency for higher GP shape parameters and higher extremes in the future. A portion of the cloud of points to the left of the exponential point tends to be located somewhat below the GP curve. The fact that the cloud of data seems to define an area of the L-moment ratio diagram indicates that three- or four-parameter distributions may be more flexible in fitting the data than the GP distribution. The generalized Beta Type 2 (4 parameters; Papalexiou and Koutsoyiannis, 2012; Chen and Singh, 2017) is an alternative distribution for modeling extremes, of which Burr Type III, generalized Gamma, Log-logistic, Burr Type XII, GP, and GEV are special cases. The log-Pearson Type 3 distribution (three parameters) is another commonly used distribution in modeling extremes (Griffis and Stedinger, 2007). These more complex distributions, which have more than one shape parameter and cover an area of the L-moment ratio diagram as opposed to a curve, may be useful in modeling the points in the L-moment ratio diagram that diverge most from the GP curve. However, a more complex distribution will not necessarily give more accurate results because the same amount of data will be used in fitting more parameters, which may result in overfitting. It is also possible that the scatter in the L-moment ratio diagram may be due to the use of a constant percentile to define the threshold for each duration at all locations as opposed to determining an optimal threshold for each specific location.

Figure 19 also shows more points farther to the right of the GP curve in the future projection period compared to the historical period, which reflects an overall tendency for larger shape parameters and larger extremes in the future. Figure 20 shows L-moment ratio diagrams for annual maxima data obtained from each downscaled climate dataset. The results for annual maxima show substantially more scatter of the points around the GEV curve than the peaks-over-threshold (POT) data around the GP curve (fig. 19). This difference may be at least partly due to the smaller sample sizes used to determine L-moments in the annual maxima case compared to the POT case. This interpretation is supported by the fact that the L-moment ratios for POT based on the longer historical period 1950–2005 (not shown here) were located even more compactly around the GP curve than in the shorter historical period 1966–2005. Therefore, it appears that the POT approach is an overall improvement over the annual maxima approach, even when a single percentile value per duration is used in defining the threshold for all grid cells.

Figure 21 bar charts show the percentage of model grid cells for which one of the five best fitting probability distributions is one of the available functions on the positive real line in the gamlss R package. The selection of the best fitting probability distributions was performed using the fitDist function in gamlss, which uses the AIC to balance GOF with distribution parsimony. In this figure, the first four bars (EXP, PARETO2, PARETO2o, and GP) correspond to different parameterizations of the GP and the exponential distribution, which is the special case of the GP having a shape parameter of zero. The bars for the remaining distributions are ordered from higher to lower percentage. It is evident that the GP is an adequate distribution for most locations; however, various parameterizations of the Weibull distribution also show up as some of the best. Using global precipitation data, Serinaldi and Kilsby (2014) show that positive precipitation (that is, equivalent to a POT approach with a threshold of zero) tends to follow a Weibull distribution and reference studies that show that as the threshold is decreased there is a progressive shift from GP to Weibull. This indicates that it is possible that the selected constant percentile-based threshold for each duration may be too low at some locations (that is, that the asymptotic conditions for the GP are not being met) resulting in more rejections of the null hypothesis of a GP distribution than the expected nominal value of 5 percent (significance level). However, figure 19 shows that the empirical L-moment ratios are in better alignment with the GP curve than with the Weibull curve. In fact, the deviation of the cloud of points from the GP curve in the region below the exponential point is away from the Weibull distribution rather than toward it. Therefore, the GP distribution has more support than the Weibull distribution. Note also that although the generalized Gamma and generalized Beta Type 2 distributions have more parameters, they do not feature among the top best distributions according to the AIC. In other words, the improvement in fit is not sufficient to justify the higher number of parameters in these distributions.

Areal Reduction Factors

Figure 22 shows the mean ARFs from method M1 of Pavlovic and others (2016) for 1-, 3-, and 7-day durations by ARF region (fig. 7). Also shown are error bars representing the mean ARF for each ARF region plus and minus one standard deviation of the grid cell values within each ARF region. As expected, the mean ARF increases with increasing duration and its variability decreases. Overall, ARF values decrease from north to south Florida, and ARF values tend to be lower on the west coast of south Florida than the east coast. According to Winsberg (2020), mid-latitude low-pressure systems can bring heavy rain over north Florida in the winter and, as a result, that part of the State has no seasonal concentration of heavy precipitation. The precipitation associated with these low-pressure frontal systems tends to be more spatially uniform than in convective systems. Baigorria and others (2007) found a more widely spread spatial correlation pattern in the northeast-southwest direction around weather stations during the frontal rainy season in the southeastern United States (perpendicular to the path of weather fronts), whereas the convective rainy season is characterized by correlations that decrease rapidly over short distances from each weather station. The influence of winter frontal systems in characterizing some of the annual maximum precipitation events likely explains the higher ARF values obtained for northern Florida. Winsberg (2020) found that, for most of the State, the monthly frequency distribution of torrential precipitation is bimodal, with a peak in June and another peak in September. He explained the June peak as being the result of the Intertropical Convergence Zone moving north and having more of an influence on the local weather, and the September peak as being due to the greater frequency of tropical storms reaching the State during that month.

Marshall and others (2004) explain how the west coast sea-breeze front typically does not penetrate as far inland as the one on the east coast because of its interaction with the large-scale easterly flow. However, this results in stronger convergence and greater precipitation along the west coast sea-breeze front. As shown by numerical simulations by Pielke (1974) and Baker and others (2001), the largest maximum vertical motion forms over the southwestern corner of the Florida peninsula in areas where the convex curvature of the coastline accentuates the horizontal convergence in the thermally driven sea breeze. Paxton and others (2009) have shown that the convergent pattern leads to a more localized circulation that appears to be associated with tornado development in southwest Florida. This more localized convection likely results in higher spatial variability in precipitation and may explain the lower ARFs obtained for the west coast of south Florida than in other regions of the State. Although ARFs vary across ARF regions, the error bars have a large degree of overlap and the differences among ARFs tend to be small. T-tests for the null hypothesis of equal means across ARF regions were also performed. In most cases, the null hypothesis of equal means of ARFs across ARF regions was rejected with the exception of the 1-day mean ARFs across the Central Everglades and Lower East Coast ARF regions.

For method M4 from Pavlovic and others (2016), where ARFs are a function of duration, return period, and spatial aggregation scale, ARFs for the longer return periods were found to be quite uncertain because of the extrapolation of the right tail of the GP distribution. This results in a large number of cases where the ARF, computed as the ratio of DDF precipitation depths at a given spatial aggregation level (1 × 1, 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11, or 13 × 13 grid cells) to the 1 × 1 grid-cell aggregation level, was larger than one. The frequency of local ARF values exceeding one increased with return level and occurs more than 50 percent of the time for the 200-year return period. ARF values greater than one were set to one before evaluating trends in ARFs across duration and return period. For central and south Florida as a whole, ARFs were found to increase with increasing duration for return periods of 5 and 10 years. For return periods longer than 10 years, ARFs for 7-day duration move progressively closer to the ARFs for the 1-day duration, both of which are smaller than the ARFs for the 3-day duration. The ARFs do not vary much with return period for return periods of 5, 10, and 25 years, and decrease with increasing return period for return periods of 50, 100, and 200 years. Possible physical explanations for the behavior of the ARF curves derived from method M4 are beyond the scope of this study and are complicated by the larger number of ARF values that were set to one for the longest return periods. Because of the uncertainties associated with extrapolation of the more extreme quantiles, method M1 was chosen as the preferred method for obtaining ARFs. The reciprocal of the M1 ARFs were used to convert grid-cell (areal-) scale DDF curves to station- (point-) scale DDF curves for comparison against the NOAA Atlas 14 PDS-based historical DDF curves. As shown in figure 22, the ARF correction is close to 1 for most datasets, except for the CORDEX models and in particular for the coarser NAM-44i models.

Historical Bias and Spatial Pattern

Figure 23 shows the overall percentage difference of the model-derived station-scale DDF depths for the entire historical period available in the downscaled climate datasets (1950–2005) compared to the PDS-based DDF depths from NOAA Atlas 14 volume 9 (Perica and others, 2013). The percentage differences are calculated from the mean DDF depths at all 174 NOAA Atlas 14 stations in central and south Florida. The DDF depths from NOAA Atlas 14 volume 9 are based on statistical fits to precipitation observations within the period 1840–2008 depending on the station. Although the two datasets use different methods for DDF fitting and different periods of record, it is still informative to evaluate their overall differences. From this figure, it is evident that the JupiterWRF dataset shows the lowest absolute difference from observed DDF depths for the 1-day duration, on the order of less than 5 percent, followed by CORDEX-Daymet with median differences of about −15 to −23 percent, and CORDEX-gridMET with median differences of about −22 and −27 percent. LOCA seems to perform more poorly than the other datasets for the 1-day duration, with median differences in the range of −36 to −40 percent. This is likely a result (at least in part) of deficiencies in its training dataset (Livneh and others, 2015) mentioned earlier in this report. However, the MACA-Livneh dataset did not show differences as large as those for LOCA, with differences in the range of −28 to −29 percent for the 1-day duration, indicating that the differences in LOCA may also be due to the downscaling method. Two different versions of the Livneh dataset (Livneh and others, 2013, 2015) and two different periods are used in bias correcting LOCA and MACA, which may also explain some of the differences. MACA-gridMET shows larger negative differences than those for MACA-Livneh. MACA-gridMET models always show negative differences, whereas some MACA-Livneh models show positive differences for the 7-day duration. Overall, the longer return periods show more negative differences (with the exception of CORDEX for the 1-day duration) and the differences are more variable across models than for the shorter return periods. The negative differences generally decrease with duration for most datasets. For the 7-day duration, some MACA-Livneh and CORDEX-Daymet models even exhibit positive differences. Overall, the range of negative differences increases with duration. Although CORDEX starts with the smallest 1-day negative difference with respect to observations, its percent difference does not improve with duration quite as much as those for the other two datasets, indicating that the RCMs and (or) the bias-correction scheme may be inadequate in simulating multiday extreme events. Although not shown in this figure, lower-resolution CORDEX models (NAM-44i) tend to have slightly smaller negative differences than the higher-resolution CORDEX models (NAM-22i); however, the difference in median percentage difference is less than 2.3 percent. Overall, the CORDEX-Daymet negative differences are smaller than in CORDEX-gridMET by 3.9–9.2 percent, depending on duration and return period, and the improvement in CORDEX-Daymet compared to CORDEX-gridMET is more pronounced for longer durations and longer return periods.

The differences between extremes fitted to observations and those derived from the downscaled climate datasets likely result from a combination of factors, some of which have been mentioned previously. These factors include the following, among others:

To evaluate the effect of the different periods of record used in comparing the official NOAA Atlas 14 PDS-based DDF curves and model-derived station-scale DDF depths for the historical period 1950–2005, DDF curves were developed by fitting the GEV distribution to annual maxima data provided by NOAA Atlas 14 at 97 (out of 174) stations (fig. 1) having at least 40 years of data within the period 1950–2005. For simplicity, the method of L-moments was used to fit the GEV distribution to annual maxima one duration at a time without corrections for CDF crossing. It was verified that, on average, data from the 97 stations are representative of data from the larger set of 174 stations. This was accomplished by comparing the median of the official NOAA Atlas 14 PDS-based DDF curves for all 174 stations against the median of the official NOAA Atlas 14 PDS-based DDF curves for the subset of 97 stations for each duration and return period. Then the difference between the model-derived station-scale DDF depths for the 1950–2005 historical period and the DDF curves that were fit to observed annual maxima for 1950–2005 was computed for the 97 stations. It was found that the median percent differences between these two sets of DDF depths were reduced by a mean absolute value of 3.9 percent, with a maximum absolute value of about 7.5 percent, from those shown in figure 23 indicating potential nonstationarity in the observed extremes, some of which may date back to the year 1840 at some stations. However, differences between the method used to develop DDF curves for this exercise and the method used to develop the official NOAA Atlas 14 PDS-based DDF curves may also explain the differences, at least in part. To address this question, we fitted a GEV distribution using the method of L-moments to all the available annual maxima data at each of the 174 stations in the region for the period 1840–2008; data availability within the period 1840–2008 varies among these stations. The DDF depths fitted using the GEV distribution and the method of L-moments for 1840–2008 were slightly smaller than the official NOAA Atlas 14 PDS-based DDF depths for the same period. As a result, the absolute median difference between the model-derived 1950–2005 DDF curves and these newly-fitted 1840–2008 DDF curves had a slight decline of 1.1 percent, on average, and a maximum decline of 3.3 percent compared to the difference between the model-derived 1950–2005 DDF curves and the official NOAA Atlas 14 PDS-based DDF curves. Therefore, it appears that the difference in period of record explains a larger portion of the difference between the official NOAA Atlas 14 PDS-based DDF curves and model-derived station-scale DDF depths for the historical period 1950–2005 than the difference in DDF-fitting methodology. Regardless of these findings, the model-derived station-scale DDF depths still appear to be highly negatively biased when compared to DDF curves developed from observations.

Figure 24 shows the overall percentage difference in station-scale precipitation depths obtained by fitting the GP distribution using CML to precipitation data from the datasets used to bias-correct each downscaled climate dataset compared to observations from NOAA Atlas 14 (1840–2008). The periods of record for each bias-correction dataset vary by downscaled climate dataset as described under the applicable subsection in the “Downscaled Climate Datasets” section. The precipitation data were obtained for the closest grid cell to each of the 174 NOAA Atlas 14 stations at the native resolution of each bias-correction dataset. The native resolution was calculated to be approximately 1 km for Daymet, 6.6 km for Livneh and others (2013 and 2015), and 4.4 km for gridMET in central and south Florida. The ARF values corresponding to these resolutions are very close to one (fig. 22). Comparison of figure 24 with the median percentage difference in precipitation depths across models for each downscaled climate dataset in figure 23 shows similar magnitudes and patterns of change with respect to changes in duration and often with respect to changes in return period. These similarities indicate that the overall percent difference of the model-derived station-scale DDF depths for the entire historical period available in the downscaled climate datasets (1950–2005) compared to the PDS-based DDF depths from NOAA Atlas 14 volume 9 is largely a result of the bias-correction datasets not being able to capture the DDF depths from NOAA Atlas 14. A similar exercise was performed for figure 24 as for figure 23. When comparing to the DDF curves developed by fitting the GEV distribution to annual maxima data provided by NOAA Atlas 14 at the 97 stations having at least 40 years of data within the period 1950–2005, the median percent differences shown in figure 24 were reduced by a mean absolute value of 4.5 percent, with a maximum reduction of about 7.4 percent. When comparing to the DDF curves developed by fitting the GEV distribution to annual maxima at all 174 NOAA Atlas 14 stations, the median percent differences shown in figure 24 declined by an absolute value of 1.2 percent, on average, with a maximum decline of 3.2 percent. Again, it appears that the difference in period of record explains a larger portion of the difference between the official NOAA Atlas 14 PDS-based DDF curves and station-scale DDF depths fitted to the bias-correction datasets than the difference in DDF-fitting methodology. Large negative biases remain, however, indicating that the gridded observational datasets used in the training and bias-correction steps of downscaling do not adequately capture the localized station-based extremes used in developing the NOAA Atlas 14 DDF curves. The biases were found to be the lowest for Daymet, which has the highest resolution of the four bias-correction datasets.

Figure 25 shows the overall percentage change in DDF depths for the future projection period (2050–89) compared to the historical period (1966–2005) for all models in each downscaled climate dataset. Comparing this figure with figure 23 indicates that the absolute value of the percent differences from observations for LOCA is much larger than the predicted percentage change from the historical to future period. The negative percent differences from observations for LOCA correspond to bias-correction factors for multiplicative quantile delta mapping (MQDM) that are greater than one, whereas the positive percent changes from present to future in LOCA corresponds to change factors that are greater than one. In the case of LOCA, the bias-correction factor would be much larger than the predicted change factors, especially for the 1-day duration, reducing the confidence in those change factors. The median percentage change for CORDEX is more comparable in magnitude to the median percentage difference than for LOCA. The median percentage change for CORDEX is smaller than the absolute value of the percentage difference with respect to observations for the shorter return periods, and equal to or larger than the absolute value of the percentage difference with respect to observations for the longer return periods. In the case of MACA, the median percentage change is smaller than the absolute value of the percentage difference from observations only for the shorter return periods for the 1-day duration and larger otherwise. Overall, the 5th, 16th, 50th, 84th percentiles and low outliers of percent change from the historical period to the future period are similar for MACA-Livneh and MACA-gridMET however, MACA-Livneh has much higher 95th percentiles and high outliers of percent change, especially for the longer return periods. In general, the median percentage change is similar between the two CORDEX bias-corrected datasets with CORDEX-gridMET only as much as 3.7 percent higher than CORDEX-Daymet. The median percentage change for JupiterWRF is larger than the absolute of the percentage difference from observations for all return periods for the 1-day duration. Figure 26 shows that, in general, the higher-resolution CORDEX models (NAM-22i) predict smaller median percentage changes in DDF depths compared to lower-resolution CORDEX models (NAM-44i). For the 5- and 10-year, 1-day event, CORDEX NAM-22i models predict slightly higher median percentage changes in DDF depths compared to CORDEX NAM-44i models.

Figure 27 shows Taylor diagrams comparing the model-derived station-scale DDF curves to the NOAA Atlas 14 PDS-based DDF curves for 1-day events. According to Taylor (2001), the Taylor diagram characterizes the statistical relation between two fields, a “test” field (often representing a field simulated by a model) and a “reference” field (usually representing “truth,” on the basis of observations). Note that the means of the fields are subtracted out before computing their second-order statistics, so the diagram does not provide information about overall differences, but solely characterizes the “centered pattern error.” For the 5-year event, the LOCA dataset has a much lower pattern correlation (0.3–0.5) than the other datasets. The pattern correlation is the Pearson product-moment coefficient of the linear correlation between two variables that are respectively the values of the same variables at corresponding locations on two different maps. For the 5-year event, the MACA models cluster in two regions, with the points close to the red line having a standard deviation ratio close to one corresponding to MACA-Livneh and the points having a lower standard deviation ratio corresponding to MACA-gridMET. Also, for the 5-year event, most of the CORDEX models cluster close to JupiterWRF and in between the two MACA clusters. The points having the lowest centered RMSD in the diagrams indicate the best performance, which in the case of the 1-day, 5-year event corresponds to some CORDEX models and MACA-gridMET models, even when the standard deviation of MACA-Livneh is closer to the observed standard deviation. As the return period increases for events of 1-day duration, CORDEX and JupiterWRF seem to deteriorate the most of all four datasets. The spatial standard deviation increases beyond that of the observations in JupiterWRF as the return period increases for events of 1-day duration, and its pattern correlation decreases. The spatial standard deviation for CORDEX gets closer to that of the observations as the return period increases; however, the pattern correlation decreases and even becomes negative for some CORDEX models for the 50- to 200-year events. This result is likely due a combination of the coarse model resolution in CORDEX not being able to capture local variations in extremes (despite the use of mean ARFs for particular regions) and the increased uncertainty associated with the estimation of the rarer extremes. In contrast, LOCA and MACA-Livneh’s performance does not deteriorate as much as the return period increases. Overall, MACA-Livneh performs best among the datasets when all return periods are considered. It is worth noting how the standard deviations tend to increase with increasing return period for all downscaled climate datasets. This is likely because the NOAA Atlas 14 DDF curves are derived by means of regional frequency analysis, which reduces the spatial variation of the shape parameter and the return levels as opposed to the CML method, which is applied more locally. Similar conclusions can be made from the Taylor diagrams for 3- and 7-day durations (figs. 28 and 29); however, for those longer durations, the standard deviation ratios are even higher than for the 1-day duration, especially for the long return periods.

Model Culling

Initially, the best models were selected by evaluating the MCI and MVI criteria from Srivastava and others (2020) and considering all the precipitation extremes indices listed in table 9. This process resulted in a very small subset of best models. On the basis of the desire to include as many models as possible to inform the potential range of future change factors, a decision was made to only use four precipitation extremes indices in selecting best models, namely the annual maxima of precipitation for durations of 1, 3, 5, and 7 consecutive days. Tables 10–12 list the best models for each climate region according to the MCI and MVI when each downscaled climate dataset (CORDEX, LOCA, and MACA) is evaluated on its own for the four precipitation extremes indices chosen. Table 13 lists the best models when all datasets are evaluated together. In this case, the median RMSE and IVSS are defined on the basis of all the models in all the downscaled climate datasets. A model was considered to be among the best if it had a negative MCI and a negative MVI when compared to either the PRISM or SFWMD “Super-grid” observational datasets. Because of the small sample size, all models were considered the best models for the JupiterWRF dataset.

In table 10, it is notable that most of the best CORDEX models for climate region 5 (south Florida) are bias-corrected with Daymet. The best models for both climate region 4 (south-central Florida) and climate region 5 include twice as many high-resolution models (NAM-22i) as low-resolution models (NAM-44i). For LOCA (table 11), about half of the best models in climate region 4 are also in climate region 5. For MACA (table 12), in climate region 4 the best models are 12 from MACA-gridMET and 4 from MACA-Livneh, whereas in climate region 5, the best models are 7 from MACA-Livneh and 5 from MACA-gridMET. Only three models are considered best in both climate regions 4 and 5 for MACA: HADGEM2-CC365 (gridMET), IPSL-CM5A-MR (gridMET), and bcc-csm1-1-m (Livneh). When all datasets are evaluated together (table 13), only 4 LOCA models are considered best in climate region 4 and 0 are considered best in climate region 5. For climate region 4, the best models are 16 from CORDEX and 19 from MACA. For climate region 5, the best models are 12 from CORDEX and 31 from MACA. It is important to note that when MCI and MVI are computed using all of the precipitation extremes indices in table 9 considering all datasets together as a single group, no LOCA models had negative values for both indices, indicating that all LOCA models performed worse than the median model across all datasets and regions.

Change Factors

Change factors in DDF precipitation depths for the period 2050–89 with respect to the period 1966–2005 were computed for all model grid cells associated with NOAA Atlas 14 stations in central and south Florida from all four downscaled climate datasets evaluated as part of this study. The change factors generated as part of this study are available in Irizarry-Ortiz and Stamm (2022). Figures 30A, 31A, 32A, 33A show median change factors for climate regions 4 (south-central Florida) and 5 (south Florida) considering all RCPs and all models or only the best models. Median change factors and their range generally increase with increasing return period and are similar across duration for climate region 4. For climate region 5, median change factors and their range increase with increasing return period and mostly decrease with increasing duration. The increase in change factors with return period is consistent with Lopez-Cantu and others (2020), who found that at the continental scale, high-end daily extremes generally increase more than low-end extremes, as determined by evaluating various downscaled climate datasets for the period 2044–99 with respect to 1951–2005. That is, heavy precipitation events get even heavier under future conditions. The exception is JupiterWRF (which is primarily based on CMIP6 GCMs) for which 1-day change factors do not vary much with return period in both climate regions. The Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC, 2021, p. 15) also states that “There will be an increasing occurrence of some extreme events unprecedented in the observational record with additional global warming, even at 1.5 °C of global warming. Projected percentage changes in frequency are higher for rarer events (high confidence).”

Median change factors from LOCA are generally the lowest, followed by those from JupiterWRF (1-day duration only), CORDEX, and MACA, which has the highest change factors. One exception is the median change factors from CORDEX being slightly larger than those from MACA for the 5- and 10-year, 7-day event in climate region 5 (fig. 32A). When only the historically best-performing models are considered in climate region 4 (fig. 31A), the median increases across most durations and return periods for MACA but decreases slightly for CORDEX and LOCA compared to the case for all models (fig. 30A), resulting in a larger inter-model range after culling. If all datasets are considered together, the median increases and gets closer to that of MACA after culling. When only the best models are considered in climate region 5 (fig. 33A), the medians for all datasets are reduced compared to the case for all models (fig. 32A), with the median for LOCA falling slightly below one for the more frequent extreme events. If all datasets are considered together, the median increases to almost that of MACA after culling. This happens because, for climate region 5, when all downscaled datasets are considered together in evaluating precipitation index performance, the best models are 12 from CORDEX and 31 from MACA. Overall, climate region 4 has higher median change factors (1.05–1.55 depending on dataset, duration, and return period) than climate region 5 (1.0–1.4 depending on dataset, duration, and return period) when considering the best models and all RCPs for both areas. The extension of the projected Caribbean drying described in the “Downscaled Climate Datasets” section may be a possible explanation for the lower change factors obtained for south Florida.

Figures 34–38 show maps of median change factors interpolated from the median at-station change factors based on best models from CORDEX, LOCA, MACA, all models from JupiterWRF, and best models selected from the evaluation of precipitation extreme indices when considering all datasets together, respectively. These maps were developed using multiquadric radial basis function interpolation from the geospt R package (Melo and others, 2012). Median change factors at the 10 closest NOAA Atlas 14 stations to each pixel location on an equally spaced 200- × 200-pixel grid were used in the interpolation. The parameter values of the multiquadric radial basis function were chosen so that the interpolation was exact at the locations of the NOAA Atlas 14 stations and so that values were not extrapolated beyond the range of the station values. The choice of interpolation method was subjective because the intent was to be able to show general patterns in median change factors rather than precise values at unsampled locations. These maps show that median change factors based on best models are slightly lower for south Florida than for south-central Florida, as shown in figures 31A and 33A. This is most evident for LOCA as shown in figure 35, which shows some areas with change factors less than one in large portions of south Florida. The exception is JupiterWRF for which a north-south gradient in median change factors is not evident (fig. 37). The fact that MACA contributes the largest proportion of best models when all downscaled datasets are considered together in evaluating precipitation index performance in climate region 5 is evident when comparing the maps of median change factors in south Florida (figs. 36 and 38). The increase in median change factors with return period shown in figures 31A and 33A is most evident in the maps for MACA (fig. 36), followed by the maps for all datasets (fig. 38) and CORDEX (fig. 34).

For the entire continental United States, Lopez-Cantu and others (2020) found that, despite being the least-biased overall compared to observed extremes, the MACA-gridMET multimember mean change factors (for the period 2044–99 with respect to 1951–2005) are also considerably larger than the other datasets evaluated, which include BCCA version 2, LOCA, and CORDEX without bias correction (called “CORDEXnoBC” hereafter). They found that BCCA version 2 results in the lowest change factors for the continental United States, followed by LOCA, low-resolution CORDEXnoBC models, high-resolution CORDEXnoBC models, and MACA-gridMET models. Furthermore, Lopez-Cantu and others (2020) found that the downscaled climate datasets somewhat preserve the pattern of the change signal in the native GCM change factors; however, the magnitude of the signal is reduced in BCCA version 2 and increased in MACA-gridMET. For the southeastern United States, using the CanESM2 GCM as an example, they found for RCP4.5 that LOCA and low-resolution CORDEXnoBC models tend to preserve the GCM change signal the best, whereas MACA-gridMET amplifies the change signal, especially for longer return periods. High-resolution CORDEXnoBC models were found to amplify the change signal more evenly across return periods. Their findings were similar for RCP8.5, but they found that LOCA tends to mute the change signal whereas its amplification in MACA-gridMET and the high-resolution CORDEXnoBC models is slightly lower than in RCP4.5.

As shown in figures 30B–33B, change factors are generally somewhat higher in RCP8.5 than in RCP4.5 for all downscaled climate datasets in both climate regions when all or only the best models are considered, with some exceptions. For climate region 5, figure 32B shows higher change factors for CORDEX models in RCP4.5 than in RCP8.5 for 1- and 3-day events with return periods shorter than 25 years, and for all durations in the 7-day event when all models are considered. However, the number of CORDEX models available for each RCP and summarized in this figure is imbalanced, with only 14 models available for RCP4.5 out of a total of 68 when all models are considered. When the 14 CORDEX models available for both RCPs are compared, RCP8.5 generally shows larger median change factors compared to RCP4.5 for climate regions 4 and 5. For climate region 5, median change factors are as much as 0.11 higher in RCP8.5 than in RCP4.5 for the 1-day duration, between 0.02 lower and 0.09 higher for the 3-day duration, and between 0.03 lower and 0.07 higher for the 7-day duration. For climate region 4, median change factors from these same 14 models are as much as 0.11 higher in RCP8.5 than in RCP4.5 for the 1-day duration, as much as 0.07 higher for the 3-day duration, and between 0.02 lower and 0.04 higher for the 7-day duration. No best CORDEX models are available for RCP4.5; therefore, the RCP4.5 markers for CORDEX are missing from figures 31B and 33B. It is also notable that JupiterWRF results in slightly larger median change factors in RCP4.5 compared to RCP8.5 in both climate regions (figs. 30B and 32B), with differences ranging from 0.02 to 0.03.

NOAA’s Physical Sciences Laboratory provides a Climate Change Web Portal (NOAA, 2022b) that summarizes climate change projections for key variables from the RCP4.5 and RCP8.5 CMIP5 scenarios. The output is provided as time series for predefined regions or as zoomable global maps. Projections of key variables and their anomalies are provided for predefined periods. Using this tool, the median temperature increase from 1956–2005 to 2070–99 expected for south and central Florida, based on the ensemble of CMIP5 models, ranges between 1.4 and 1.6 °C for RCP4.5, and between 2 and 2.25 °C for RCP8.5, with the largest increases in the central Florida region. Temperature increases between 1956–2005 and 2040–69 are very similar. Assuming that the Clausius-Clapeyron (CC) relation holds for precipitation extremes, an approximately 7-percent increase in precipitation would be expected per degree Celsius of warming. This would correspond to change factors of 1.10–1.11 for RCP4.5 and of 1.14–1.16 RCP8.5 based on the median temperature increase projections from the CMIP5 ensemble. The median LOCA change factors shown in figures 30–33 are generally lower than the expected range if the extremes followed the CC relation, indicating a sub-CC relation. Conversely, the median change factors from MACA and CORDEX are generally higher than expected based on the median temperature increase projections from the CMIP5 ensemble, especially for the longest return periods, indicating a super-CC relation. These deviations from CC cannot be attributed to physical factors such as moisture convergence and enhanced convection because they are based on statistically and dynamically downscaled projections, which have been bias-corrected, and therefore any links to physics have been severed. The median change factors from JupiterWRF are more in line with expectations from the CC relation alone; however, this does not imply that they are more accurate than change factors obtained from other downscaled climate datasets. The generally higher median change factors obtained for RCP8.5 than for RCP4.5 are in line with expectations from sole consideration of the CC relation.

Figure 39 shows boxplots of change factors for all models and stations in climate region 5 (south Florida). The variability in change factors increases with increasing return period when all datasets are considered together because of the influence of MACA models. Change factors as high as 8 in south Florida for the 1-day, 200-year event (fig. 39) and as high as 10 in south-central Florida for the 7-day, 200-year event are estimated for individual stations by MACA. CORDEX also shows change factors as high as 6.5 for the 7-day, 200-year event in south Florida (fig. 39) and as high as 6 for the 7-day, 200-year event in south-central Florida. The highest outlier change factors from CORDEX are slightly lower and occur less frequently than in MACA. When these very high change factors are multiplied by NOAA Atlas 14 DDF curves, the resulting future projected DDF values are beyond record high precipitation accumulations being recorded around the globe in the last decade. For example, record precipitation accumulations were observed during Hurricane Harvey in Texas and Louisiana in 2017 with the highest storm-total rainfall reported as 60.6 in. at two locations in Texas exceeding the previously accepted U.S. tropical cyclone storm total rainfall record (NOAA, 2018). The estimated return period for rainfall of this magnitude in southeast Texas is 1,000–10,000 years (Emanuel, 2017). Rainfall during Hurricane Florence in North Carolina in 2018 exceeded 30 in. near landfall with a localized maximum of 35.9 in. that set a State record for tropical cyclone rainfall (NOAA, 2019). The 3- and 4-day rainfall totals for this event surpassed the 1,000-year event at a particular site (Griffin and others, 2019). Hurricane Dorian in 2019 produced storm-total rainfall of 22.8 in. at a location in the Bahamas (NOAA, 2020). For an August 2016 storm in Louisiana, radar rainfall indicated more than 34 in. of rainfall at a location during the storm duration, well over the estimated 1,000-year return interval (Brown and others, 2020). A shift toward more frequent extremes, especially for rarer events, was documented in the Intergovernmental Panel on Climate Change’s (IPCC’s) Fifth Assessment Report (Collins and others, 2013) and further validated in their Sixth Assessment Report (IPCC, 2021). Attribution studies have determined that climate change has increased the intensity of these historical events and made them more likely. (See van der Wiel and others [2017] for the August 2016 Louisiana event; van Oldenborgh and others, [2017] for Hurricane Harvey; Reed and others [2020] for Hurricane Florence; and Reed and others [2021] for Hurricane Dorian.)

Evidence of record-breaking precipitation extremes and their links to climate change has been accumulating over the last few decades. However, given that the very high change factors computed here were generated by statistical downscaling or bias-corrected dynamical downscaling, the confidence in them is not as high as it would be if they had been generated by a purely physically based model. In addition, the highest change factors are associated with very long return periods, which are much longer than the number of years of data used in DDF fitting and are highly uncertain despite the use of a POT approach in this study. The boxplots of median station change factors across models for climate region 5 (south Florida) (fig. 40) show the central tendency of the station change factors across models in the region. Although the influence of individual station outliers is removed when looking at medians across stations as in this figure, it is evident that MACA still results in higher change factors than CORDEX and LOCA, especially for 1- and 3-day durations. From these two figures, it is evident that change factors below one (that is, drier extremes) are possible as indicated by CORDEX, LOCA, and MACA, with the median for LOCA being very close to a change factor of one. Although a similar boxplot is not shown here for climate region 4 (south-central Florida), it can be inferred from figures 30 and 32, that overall, climate region 4 has higher change factors than climate region 5. In particular, the median for LOCA is above one for all durations and return periods in climate region 4. In fact, it is the 16th percentile for LOCA that is very close to one for climate region 4.

Spot checks of the data associated with some of the extremely high change factors in MACA show that they are the result of extremely high outlier events in the future projection period 2050–89. The specific case of a change factor of 10 in MACA is due to 35 in. of precipitation during a 3- to 7-day event in the future projection period at a grid cell in southwest Florida. The DDF fit at this station results in an extrapolation to 164 in. of precipitation for the 200-year event, whereas 35 in. of precipitation are fitted for an event with a return period of between 25 and 50 years. This is a limitation arising from fitting DDF curves for a limited future 40-year projection period based on local data, in that it is possible that this very high event only happened once in a much longer period than the 40 years analyzed. However, these very high precipitation events appear to be much more common in MACA than in the other downscaled climate datasets and are likely a result of the downscaling method.

Wootten (2018) shows how subtle decisions performed in the process of statistical downscaling such as tail adjustment, trace adjustment, and interpolation can affect the skill of the prediction and increase the uncertainty in future projections, especially for extremes and event occurrence. Wang and others (2020) found that MACA projects significant increases in the 20-year daily event in New England compared to LOCA for which the projected changes are inconclusive over most of the area. In LOCA, model consensus is lacking over most of the region, but where models do agree, the projected changes have a smaller magnitude than in MACA. Their finding is generally consistent with this study’s findings for Florida. Wang and others (2020) note that projected changes in the mean annual maximum daily event are quite close between MACA-Livneh and MACA-gridMET but differ significantly between LOCA and MACA-Livneh. As a result, they conclude that differences in the projected changes are primarily caused by differences in the downscaling method, rather than the observational data used in training. However, LOCA and MACA used two different versions of the Livneh dataset for bias correction, which complicates the comparison.

The factors identified as potentially contributing to biases in DDF curves for the historical period, in addition to other factors such as scenario uncertainty, contribute to the total uncertainty in derived change factors. Wootten and others (2017) attempt to characterize the sources of uncertainty from global climate models and downscaling techniques. In particular, they evaluate the sources of uncertainty in the decadal mean of the annual number of days with extreme precipitation (more than 1 in.) for central Florida across various downscaled climate datasets, including some analyzed in this study and others that were not. Based on data from an idealized scenario that only compares 16 common GCMs and two common emission scenarios from MACA and BCCA, they find that, for the 2070s, GCM model uncertainty dominates, followed by natural variability and downscaling uncertainty. The emission scenario contributes the least to the total uncertainty, although its contribution increases with time whereas the contribution of natural variability decreases with time. Although a quantitative uncertainty analysis has not been performed as part of this study, Wootten and others (2017) findings seem qualitatively consistent with this study’s findings for change factors on the basis of an examination of figures 39 and 40. Once the station median is computed, the range of variation across models is reduced. The difference in median between datasets is smaller than between models within a dataset, and the difference between RCPs is even smaller. The DDF fitting methodology introduces additional uncertainty in the estimated change factors.

Because of the large uncertainties in the derived change factors and associated future DDF curves, methods for decision making under deep uncertainty could be used to increase the flexibility in the planning process. These methods are documented by Marchau and others (2019) and include decision scaling, robust decision making, and dynamic adaptive policy pathways, among others. Flexibility in the phasing and design of various flood control adaptation measures could reduce long-term costs while considering the evolution in climate science as well as local and global changes in climate and sea-level rise with consideration for potential tipping points. Research gaps and recommendations have also been outlined by Florida International University’s Sea Level Solutions Center (Florida International University, 2021).